数值分析上机实习报告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
升级会员 