数值分析上机实习报告Word格式文档下载.docx
- 文档编号:22923647
- 上传时间:2023-02-06
- 格式:DOCX
- 页数:27
- 大小:74.40KB
数值分析上机实习报告Word格式文档下载.docx
《数值分析上机实习报告Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《数值分析上机实习报告Word格式文档下载.docx(27页珍藏版)》请在冰豆网上搜索。
〃初始数据
usingnamespaeestd;
double*newton(doublea,doubleb,doubleeps);
doublenewtonz(doublex);
voidmain()
{
doublea=0.1,b=1.9,eps=0.00001,*result;
eout<
<
"
\n牛顿法解方程:
xA7-28xA4+14=0,在(0.1,1.9)中求近似根,初始值为区间端点,\n误差为0.00001。
\n"
endl;
学号:
2014021966姓名:
徐林\n"
endl;
result=newton(a,b,eps);
if(a<
=result[0]&
&
result[0]<
=b)
cout<
近似根为:
"
=result[1]&
result[1]<
//
结束,按任意键关闭"
getchar();
}//主函数结束
//*
******************************************************************
doublenewtonz(doublex){
doublex1=0.0,t;
//牛顿迭代子函数
t=(7*pow(x,6)-4*28*pow(x,3));
if(t==0)
exit(0);
x1=(x-((pow(x,7)-28*pow(x,4)+14)/t));
returnx1;
}
double*newton(doublea,doubleb,doubleeps){
doublex0=0.0,x1=1.0,x2=0.0,re[2];
intk=0;
x0=a;
while(x0>
eps)
k++;
x2=x1;
x1=newtonz(x1);
x0=fabs(x1-x2);
}re[0]=x1;
//代入a迭代计算
//调用牛顿迭代子函数
x0=b,k=0;
x0=fabs(x1-x2);
}re[1]=x1;
//代入b迭代计算
returnre;
1.3计算结果打印
■1'
C:
\U5er5\AdminiStratoDesl'
lop、.D亡bug\胡亚梅eye"
li.845497
Pi'
essanj/ke^tocontinue
1.4MATLAB上机程序
functiony=Newton(f,df,xO,eps,M)d=0;
fork=1:
M
iffeval(df,xO)==O
d=2;
break
else
x1=xO-feval(f,xO)/feval(df,xO);
end
e=abs(x1-x0);
x0=x1;
ife<
=eps&
abs(feval(f,x1))v=epsd=1;
ifd==1
y=x1;
elseifd==0
y='
迭代M次失败:
奇异’
functiony=df(x)
y=7*xA6-28*4*xA3;
End
functiony=f(x)
y=xA7-28*xA4+14;
>
x0=1.9;
eps=0.00001;
M=100;
x=Newton('
f,'
df,x0,eps,M);
vpa(x,7)
1.5问题讨论
1.需注意的是,要使用Newton迭代法须f(x)x728x414满足定理中的条件I,n,川,以及
f(xc)•f'
(x°
)>
0。
要用误差范围来控制循环的次数,保证循环的次数和质量,编写程序过程中要
注意标点符号的使用,正确运用适当的标点符号,Newton迭代法是局部收敛的,在使用时应先确定
初始值,否则所得的解可能不在所要求的范围内
(3)因为newton法求方程是平方收敛的,所以较为精确,但是要求出函数的导数,且必须有二阶导数。
第二题
2.已知函数值如下表:
X
1
2
3
4
5
f(x)
0.69314718:
1.0986123
1.3862944
1.6094378
6
7
8
9
10
1.7917595
1.9459101
2.079445
2.1972246
2.3025851
(x)
(1)=1
(10)=0.1
试用三次样条插值求f(4.563)及f'
(4.563)的近似值
2.1理论依据及方法应用条件三次样条插值函数可定义为:
对于[a,b]上的一个划分n
a<
xo<
xi<
X2<
・・.<
Xn-i<
Xn=b.(n>
=2)
如果定义在[a,b]上的函数S(x),满足
(1).在[xi,xi+i]上为3次多项式;
(2).S(x),S'
(x),
S"
(x)在[a,b]上连续,则称S(x)在[a,b]上划分的3次样条函数,如果对于f(x)[a,b],
s(x)f(x)还满足s(xjf(x),i0~n,则称s(x)为f(x)的三次样条插值函数。
其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中1,x,,x2,x3,(x-xj)+3为基函数,而取B样条函数Q3((x-Xj)/h)为基函数.由于三次样条函数空间是N+3隹的,故我们把分点扩大到X1,Xn+1则任意三次样条函数可用Q3((X-Xj)/h)线性组合来表示S(x)=N1CjQ3((x-xj)/h)这样对不同插值
j1
问题,若能确定Cj由解的唯一性就能求得S(x)。
由s(xi)=yi,I=1,2,…Ns'
(x0)=y0,s(xN)=yn可得
N1
S(Xi)=CjQ3((Xi-Xj)/h)=yi
)=1/h
CjQ3
(X0-Xj)
/h)=y
(xN)=1/h
(XN-Xj)
d0
d1
dN
A—y'
o)
00
yy.
(j1jh•4h.(
j1j
h.
j
(yN
hN1
yN
yj
hj1
yj1
)6f(Xj1,Xj,Xj1)
/*宏定义*/
/*追赶法求解三弯矩方程*/
hj
2.2计算程序
#include<
stdio.h>
#include<
#defineN10
main()
floats,ds,t;
floatdy0=1,dy9=0.1;
intj;
intx[N]={1,2,3,4,5,6,7,8,9,10};
floaty[N]={0,0.69314718,1.0986123,1.3862944,1.6094378,
1.7917595,1.9459101,2.079445,2.1972246,2.3025851};
intb[N]={2,2,2,2,2,2,2,2,2,2},h[N-1];
floatd[N],u[N-1],v[N-1],a[N-1],c[N-1],B[N],l[N-1],p[N],X[N];
for(j=1;
jv=9;
j++)
h[j-1]=x[j]-x[j-1];
d[0]=6/h[0]*(y[1]/h[0]-y[0]/h[0]-dy0);
d[9]=6/h[8]*(dy9-y[9]/h[8]+y[8]/h[8]);
j<
=8;
j++)d[j]=6/(h[j-1]+h[j])*(y[j+1]/h[j]-y[j]/h[j]-y[j]/h[j-1]+y[j-1]/h[j-1]);
for(u[8]=1,j=0;
=7;
u[j]=h[j-1]/(h[j-1]+h[j]);
for(v[0]=1,j=1;
jv=8;
v[j]=h[j]/(h[j-1]+h[j]);
for(j=0;
a[j]=u[j];
c[j]=v[j];
for(B[0]=b[0],j=1;
B[j]=b[j]-a[j]/B[j-1]*c[j-1];
=9;
l[j]=a[j]/B[j-1];
p[j]=d[j]-l[j]*p[j-1];
X[9]=p[9]/B[9];
for(j=8;
j>
=0;
j--)
X[j]=p[j]/B[j]-c[j]*X[j+1]/B[j];
t=4.563;
/*解f(x)的值
/*解f'
(x)
/*打印结果*/
s=X[3]*pow((x[4]-t),3)/6/h[3]+X[4]*pow((t-x[3]),3)/6/h[3]+(y[3]-X[3]*h[3]*h[3]/6)*(x[4]/h[3]-t/h[3])+(y[4]-X[4]*h[3]*h[3]/6)*(t/h[3]-x[3]/h[3]);
*/
ds=-X[3]*pow((x[4]-t),2)/2/h[3]+X[4]*pow((t-x[3]),2)/2/h[3]-(y[3]-X[3]*h[3]*h[3]/6)/h[3]+(y[4]-X[4]*h[3]*h[3]/6)/h[3];
的值*/
printf("
s=%f\nds=%f\n"
s,ds);
}2.3计算结果打印
2.4MATLAB上机程序
functionQ=san(ssss,p)
Q=zeros(2,1);
x=[1;
2;
3;
4;
5;
6;
7;
8;
9;
10];
y=[0;
0.69314718;
1.0986123;
1.3862944;
1.6094378;
1.7917595;
1.9459101;
2.079445;
2.1972246;
2.3025851];
h=zeros(10,1);
d=zeros(10,1);
u=zeros(10,1);
v=zeros(10,1);
r=zeros(10,1);
l=zeros(10,1);
z=zeros(10,1);
m=zeros(10,1);
fort=1:
1:
h(t)=x(t+1)-x(t);
d
(1)=6/h
(1)*((y
(2)-y
(1))/h
(1)-1);
d(10)=6/h(9)*(0.1-(y(10)-y(9))/h(9));
u(t+1)=h(t)/(h(t)+h(t+1));
v(t+1)=1-u(t+1);
d(t+1)=6/(h(t)+h(t+1))*((y(t+2)-y(t+1))/(x(t+2)-x(t+1))-(y(t+1)-y(t))/(x(t+1)-x(t)));
u(10)=1;
v
(1)=1;
r
(1)=d
(1);
fort=2:
l(t)=u(t)/r(t-1);
r(t)=d(t)-l(t)*v(t-1);
z
(1)=d
(1);
z(t)=d(t)-l(t)*z(t-1);
m(10)=z(10)/r(10);
fort=9:
-1:
m(t)=(z(t)-v(t)*m(t+1))/r(t);
ifp>
=t&
p<
(t+1)
Q(:
1)=feval(ssss,p,t,x,m,h,y);
functionQ=ssss(p,t,x,m,h,y)
Q(1,1)=((power((x(t+1)-p),3)*m(t)+power((p-x(t)),3)*m(t+1))/6+(y(t)-m(t)*h(t)*h(t)/6)*(x(t+1)-p)+(y(t+1)-m(t+1)*h(t)*h(t)/6)*(p-x(t)))/h(t);
Q(2,1)=(-(power((x(t+1)-p),2)*m(t)+power((p-x(t)),2)*m(t+1))/2+(y(t)-m(t)*h(t)*h(t)/6)+(y(t+1)-m(t+1)*h(t)*h(t)/6))/h(t);
end
2.5问题讨论
1.若要用追赶法求解三对角方程组,三对角阵需要满足:
A(i=1,2,…,n)均非奇异,保证A有唯一的
Doolittle分解;
aG丰0;
2.样条插值效果比Lagrange插值好,三次样条插值的解存在且唯一,近似误差较小.并且没有Runge现象。
3.用Romberg法求:
xx1.4(5x7)sinx2dx(允许误差&
3.1理论依据及方法应用条件弹b?
f(a)f(b)数值积分的Romberg算法计算步骤如下:
⑴1(I1)
T1T1
4mTmk)Tm
=0.00001)。
(0)(0)
111
时,
(k1)
Tm1
就停机
ba
211
fa
i1
(2i
m1,2,,lk1,2,,lm1
3.2计算程序
/*定义函数f(x)*/
#include<
#defineN9floatf(floatx)
floaty;
y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);
return(y);
floatT[N+1][N+1],h[N+1],a=1,b=3,m[N+1];
inti,l;
T[1][0]=(b-a)*(f(a)+f(b))/2;
l=1;
while(l<
=N)
m[l]=0;
for(i=1;
i<
=(pow(2,l-1));
i++)
m[l]+=f(a+(2*i-1)*(b-a)/pow(2,l));
T[1][l]=(T[1][l-1]+(b-a)*m[l]/pow(2,l-1))/2;
l++;
i=1;
while(i<
for(l=1;
l<
=N-i+1;
l++)
T[i+1][l-1]=(pow(4,i)*T[i][l]-T[i][l-1])/(pow(4,i)-1);
h[i]=T[i][0]-T[i+1][0];
if(fabs(h[i])<
=1e-5)break;
i++;
Theansweris:
%f"
T[i+1][0]);
3.3计算结果打印
3.4MATLAB上机程序
function[T,n]=mromb(f,a,b,eps)
ifnarginv4,eps=1e-6;
end
h=b-a;
R(1,1)=(h/2)*(feval(f,a)+feval(f,b));
n=1;
J=0;
err=1;
while(err>
J=J+1;
h=h/2;
S=0;
fori=1:
n
x=a+h*(2*i-1);
S=S+feval(f,x);
R(J+1,1)=R(J,1)/2+h*S;
J
R(J+1,k+1)=(4Ak*R(J+1,k)-R(J,k))/(4Ak-1);
err=abs(R(J+1,J+1)-R(J+1,J));
n=2*n;
R;
T=R(J+1,J+1);
formatlong
f=@(x)(3.Ax)*(x.A1.4)*(5*x+7)*sin(x*x);
[T,n]=mromb(f,1,3,1.e-5)
3.5问题讨论
1、Romberge算法的优点是:
a、把积分化为代数运算,而实际上只需求「⑴,以后用递推可得
b、算法简单且收敛速度快,一般4或5次即能达到要求。
c、节省存储量,算出的可存入。
2、Romberge算法的缺点是:
a、对函数的光滑性要求较高。
b、计算新分点的值时,这些数值的个数成倍增加。
第四题
4.用定步长四阶Runge-Kutta法求解
厂dy,/dt1
dyz/dty3
dy3/dt10001000y2100y3
yi(0)0
y2(0)0
Iy3(0)0
h0.0005,打印yi(0.025),y(0.045),%(0.085),%(0.1),(i1,2,3)
4.1理论依据及方法应用条件
Runge-Kutta法的基本思想:
旳1不是按Taylor公式展开,而是先写成tn处附近的值的线性组合(有待定系数),再按Taylor公式展开,然后确定待定常数。
四阶古典Runge-Kutta公式:
Yn1
Yn(K12K
22K3K4)
K1
hF(Xn,Yn)
K2
hF(Xn
£
h,Yn
2K1)
K3
^h,Yn
丄心)
K4
-K3)
4.2计算程序
#include<
intmain()
inti;
doubleh=0.0005;
doublek1,k2,k3,k4;
doubley1=0.0,y2=0.0,y3=0.0;
=200;
k1=k2=k3=k4=h*1.0;
y1+=(k1+2*k2+2*k3+k4)/6;
k仁k2=k3=k4=h*y3;
y2+=(k1+2*k2+2*k3+k4)/6;
k仁h*(1000-1000*y2-100*y3);
k2=h*(1000-1000*y2-100*(y3+0.5*k1));
k3=h*(1000-1000*y2-100*(y3+0.5*k2));
k4=h*(1000-1000*y2-100*(y3+k3));
y3+=(k1+2*k2+2*k3+k4)/6;
if(i==50)
printf("
\ny1(0.025)=%fcontinue;
if(i==90)
\ny1(0.045)=%fcontinue;
if(i==170)
\ny1(0.085)=%fcontinue;
if(i==200)
\ny1(0.100)=%f
y2(0.025)=%f
y2(0.045)=%f
y2(0.085)=%f
y2(0.100)=%f
y3(0.025)=%f"
y1,y2,y3);
y3(0.045)=%f"
y3(0.085)=%f"
y3(0.100)=%f\n\n"
4.3计算结果打印
4.4MATLAB上机程序
functionY=R_K(df1,a,b,h)
m=(b-a)/h;
Y=zeros(3,1);
S=zeros(3,1);
K=zeros(3,4);
x=a;
y1=a;
y2=a;
y3=a;
forn=1:
m
K(:
1)=feval(df1,x,y1,y2,y3);
x=x+0.5*h;
S(:
1)=Y+0.5*h.*K(:
1);
y1=S(1,1);
y2=S(2,1);
y3=S(3,1);
2)=feval(df1,x,y1,y2,y3);
2);
3)=feval(df1,x,y1,y2,y3);
1)=Y+h.*K(:
3);
4)=feval(df1,x,y1,y2,y3);
Y=Y+h.*(K(:
1)+2.*K(:
2)+2.*K(:
3)+K(:
4))/6;
functionZ=df1(x,y1,y2,y3)
Z=zeros(3,1);
Z
(1)=0*x+0*y1+0*y2+0*y3+1;
Z
(2)=0*x+0*y1+0*y2+1*y3;
Z(3)=0*x+0*y1-1000*y2-100*y3+1000;
4.5问题讨论
1.定步长四阶runge-kutta法稳定,精度高,可根据有y'
f(t,y)变化的情况与需要的精度自动修
改步长,误差小且程序简单,存储量少
2.但是Runge-Kutta法需要每步都计算函数值f(t,y)四次,在函数较复杂时,工作量就会变得较大
可靠性有待核查。
第五题
5.已知A与b
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 上机 实习 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)