数值分析上机实习报告汇编Word文档格式.docx
《数值分析上机实习报告汇编Word文档格式.docx》由会员分享，可在线阅读，更多相关《数值分析上机实习报告汇编Word文档格式.docx（28页珍藏版）》请在冰豆网上搜索。
math.h>
usingnamespacestd;
double*newton（doublea,doubleb,doubleeps）;
//牛顿迭代函数
doublenewtonz（doublex）;
//牛顿迭代子函数
voidmain（）
{
doublea=0.1,b=1.9,eps=0.00001,*result;
//初始数据
cout<
<
"
\n牛顿法解方程：
x^7-28x^4+14=0,在（0.1,1.9）中求近似根，初始值为区间端点，\n误差为0.00001。
\n"
endl;
学号：
2014021966姓名：
徐林\n"
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）//代入a迭代计算
{
k++;
x2=x1;
x1=newtonz（x1）;
//调用牛顿迭代子函数
x0=fabs（x1-x2）;
}re[0]=x1;
x0=b,k=0;
eps）//代入b迭代计算
}re[1]=x1;
returnre;
1.3计算结果打印
1.4MATLAB上机程序
functiony=Newton（f,df,x0,eps,M）
d=0;
fork=1:
M
iffeval（df,x0）==0
d=2;
break
else
x1=x0-feval（f,x0）/feval（df,x0）;
end
e=abs（x1-x0）;
x0=x1;
ife<
=eps&
abs（feval（f,x1））<
=eps
d=1;
end
ifd==1
y=x1;
elseifd==0
y='
迭代M次失败'
;
else
y='
奇异'
functiony=df（x）
y=7*x^6-28*4*x^3;
End
functiony=f（x）
y=x^7-28*x^4+14;
x0=1.9;
eps=0.00001;
M=100;
x=Newton（'
f'
'
df'
x0,eps,M）;
vpa（x,7）
1.5问题讨论
1.需注意的是，要使用Newton迭代法须
满足定理中的条件Ⅰ,Ⅱ,Ⅲ，以及
0。
要用误差范围来控制循环的次数，保证循环的次数和质量，编写程序过程中要注意标点符号的使用，正确运用适当的标点符号，Newton迭代法是局部收敛的，在使用时应先确定初始值，否则所得的解可能不在所要求的范围内。
（3）因为newton法求方程是平方收敛的，所以较为精确,但是要求出函数的导数，且必须有二阶导数。
第二题
2.已知函数值如下表:
1
2
3
4
5
0.69314718
1.0986123
1.3862944
1.6094378
6
7
8
9
10
1.7917595
1.9459101
2.079445
2.1972246
2.3025851
=1
=0.1
试用三次样条插值求
及
的近似值。
2.1理论依据及方法应用条件
三次样条插值函数可定义为：
对于[a,b]上的一个划分∏
a<
x0<
x1<
x2<
…<
xn-1<
xn=b.（n>
=2）
如果定义在[a,b]上的函数S（x）,满足
（1）.在[xi,xi+1]上为3次多项式；
（2）.S（x）,S＇（x），S＂（x）在[a,b]上连续，则称S（x）在[a,b]上划分
的3次样条函数，如果对于
，
还满足
，则称
为
的三次样条插值函数。
其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中1,x,,x2,x3,（x-xj）+3为基函数,而取B样条函数Ω3（（x-xj）/h）为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X-1,XN+1,则任意三次样条函数可用Ω3（（x-xj）/h）线性组合来表示S（x）=
cjΩ3（（x-xj）/h）这样对不同插值问题,若能确定cj由解的唯一性就能求得S（x）。
由s（xi）=yi,I=1,2,…Ns’（x0）=y0’,s’（xN）=yN’可得
S（xi）=
cjΩ3（（xi-xj）/h）=yi
S’（x0）=1/h
cjΩ3’（（x0-xj）/h）=y’0
S’（xN）=1/h
cjΩ3’（（xN-xj）/h）=y’N
2.2计算程序
#include<
stdio.h>
#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;
j<
=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]）;
=8;
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;
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;
j++）/*追赶法求解三弯矩方程*/
B[j]=b[j]-a[j]/B[j-1]*c[j-1];
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;
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]）;
/*解f（x）的值*/
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];
/*解f’（x）的值*/
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）;
2.5问题讨论
1.若要用追赶法求解三对角方程组，三对角阵需要满足：
（i=1,2,…,n）均非奇异，保证
有唯一的Doolittle分解;
2.样条插值效果比Lagrange插值好,三次样条插值的解存在且唯一，近似误差较小.并且没有Runge现象。
第三题
3.用Romberg算法求
（允许误差ε=0.00001）。
3.1理论依据及方法应用条件
数值积分的Romberg算法计算步骤如下:
当
时，就停机。
3.2计算程序
#defineN9
floatf（floatx）/*定义函数f（x）*/
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）
ifnargin<
4,eps=1e-6;
h=b-a;
R（1,1）=（h/2）*（feval（f,a）+feval（f,b））;
n=1;
J=0;
err=1;
while（err>
eps）
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;
fork=1:
J
R（J+1,k+1）=（4^k*R（J+1,k）-R（J,k））/（4^k-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.^x）*（x.^1.4）*（5*x+7）*sin（x*x）;
[T,n]=mromb（f,1,3,1.e-5）
3.5问题讨论
1、Romberge算法的优点是:
a、把积分化为代数运算,而实际上只需求T1（i）,以后用递推可得。
b、算法简单且收敛速度快,一般4或5次即能达到要求。
c、节省存储量，算出的可存入。
2、Romberge算法的缺点是:
a、对函数的光滑性要求较高。
b、计算新分点的值时,这些数值的个数成倍增加。
第四题
4.用定步长四阶Runge-Kutta法求解
打印
4.1理论依据及方法应用条件
Runge-Kutta法的基本思想:
不是按Taylor公式展开，而是先写成
处附近的值的线性组合（有待定系数）,再按Taylor公式展开，然后确定待定常数。
四阶古典Runge-Kutta公式:
4.2计算程序
#include<
intmain（）
inti;
doubleh=0.0005;
doublek1,k2,k3,k4;
doubley1=0.0,y2=0.0,y3=0.0;
=200;
i++）
k1=k2=k3=k4=h*1.0;
y1+=（k1+2*k2+2*k3+k4）/6;
k1=k2=k3=k4=h*y3;
y2+=（k1+2*k2+2*k3+k4）/6;
k1=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）=%fy2（0.025）=%fy3（0.025）=%f"
y1,y2,y3）;
continue;
}
if（i==90）
\ny1（0.045）=%fy2（0.045）=%fy3（0.045）=%f"
if（i==170）
\ny1（0.085）=%fy2（0.085）=%fy3（0.085）=%f"
if（i==200）
\ny1（0.100）=%fy2（0.100）=%fy3（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）;
S（:
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法稳定，精度高，可根据有
变化的情况与需要的精度自动修改步长，误差小且程序简单，存储量少。
2.但是Runge-Kutta法需要每步都计算函数值
四次,在函数较复杂时,工作量就会变得较大,可靠性有待核查。
第五题
5.已知A与b
A=
12.38412
2.115237
-1.061074
1.112336
-0.113584
0.718719
1.742382
3.067813
-2.031743
19.141823
-3.125432
-1.012345
2.189736
1.563849
-0.784156
1.112348
3.123124
15.567914
3.123848
2.031454
1.836742
-1.056781
0.336993
-1.010103
27.108437
4.101011
-3.741856
2.101023
-0.71828
-0.037585
19.897918
0.431637
-3.111223
2.121314
1.784317
9.789365
-0.103458
-1.103456
0.238417
-0.784165
14.713846
3.123789
-2.213474
3.
升级会员 