数值方法计算实习题信计101班和数应1.docx
《数值方法计算实习题信计101班和数应1.docx》由会员分享，可在线阅读，更多相关《数值方法计算实习题信计101班和数应1.docx（26页珍藏版）》请在冰豆网上搜索。
数值方法计算实习题信计101班和数应1
信计091
要求：
1、用Matlab语言或你熟悉的其他算法语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：
课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：
严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。
一、下表给出了飞行中鸭子的上部形状的节点数据，试用三次样条插值函数（自然边界条件）和20次Lagrange插值多项式对数据进行插值。
用图示出给定的数据，以及
和
。
0.9
1.3
1.9
2.1
2.6
3.0
3.9
4.4
4.7
5.0
6.0
1.3
1.5
1.85
2.1
2.6
2.7
2.4
2.15
2.05
2.1
2.25
7.0
8.0
9.2
10.5
11.3
11.6
12
12.6
13.0
13.3
2.3
2.25
1.95
1.4
0.9
0.7
0.6
0.5
0.4
0.25
解：
>>x=[0.91.31.92.12.63.03.94.44.75.06.07.08.09.210.511.311.61212.613.013.3];
>>y=[1.31.51.852.12.62.72.42.152.052.12.252.32.251.951.40.90.70.60.50.40.25];
（1）三次样条插值法
在MATLAB中编写m文件
function[f,f0]=scyt（x,y,y2_1,y2_N,x0）
%y2_1和y2_N分别为自然边界条件
%插值点x的坐标：
x0
symst;
f=0.0;f0=0.0;
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等'）;
return;
end
fori=1:
n
if（x（i）<=x0）&&（x（i+1）>=x0）
index=i;
break;
end
end
A=diag（2*ones（1,n））;
A（1,2）=1;
A（n,n-1）=1;
u=zeros（n-2,1）;
lamda=zeros（n-1,1）;
c=zeros（n,1）;
fori=2:
n-1
u（i-1）=（x（i）-x（i-1））/（x（i+1）-x（i-1））;
lamda（i）=（x（i+1）-x（i））/（x（i+1）-x（i-1））;
c（i）=3*lamda（i）*（y（i）-y（i-1））/（x（i）-x（i-1））+3*u（i-1）*（y（i+1）-y（i））/（x（i+1）-x（i））;
A（i,i+1）=u（i-1）;
A（i,i-1）=lamda（i）;
end
c
（1）=3*（y
（2）-y
（1））/（x
（2）-x
（1））-（x
（2）-x
（1））*y2_1/2;
c（n）=3*（y（n）-y（n-1））/（x（n）-x（n-1））-（x（n）-x（n-1））*y2_N/2;
m=zgf（A,c）;
h=x（index+1）-x（index）;
f=y（index）*（2*（t-x（index））+h）*（t-x（index+1））^2/h/h/h+y（index+1）*（2*（x（index+1）-t）+h）*（t-x（index））^2/h/h/h+m（index）*（t-x（index））*（x（index+1）-t）^2/h/h-m（index+1）*（x（index+1）-t）*（t-x（index））^2/h/h;
f0=subs（f,'t',x0）;
其中的zgf函数为追赶法，其程序为
functionx=zgf（A,b）
n=rank（A）;
for（i=1:
n）
if（A（i,i）==0）
disp（'Error:
对角有元素为0！
'）;
return;
end
end;
d=ones（n,1）;
a=ones（n-1,1）;
c=ones（n-1）;
for（i=1:
n-1）
a（i,1）=A（i+1,i）;
c（i,1）=A（i,i+1）;
d（i,1）=A（i,i）;
end
d（n,1）=A（n,n）;
for（i=2:
n）
d（i,1）=d（i,1）-（a（i-1,1）/d（i-1,1））*c（i-1,1）;
b（i,1）=b（i,1）-（a（i-1,1）/d（i-1,1））*b（i-1,1）;
end
x（n,1）=b（n,1）/d（n,1）;
for（i=（n-1）:
-1:
1）
x（i,1）=（b（i,1）-c（i,1）*x（i+1,1））/d（i,1）;
end
在MATLAB中输入指令
>>[f,f0]=scyt（x,y,0,0）
f=
1000/729*（27/5*t-1377/100）*（t-39/10）^2+1000/729*（522/25-24/5*t）*（t-3）^2+100/81*（-6396162352027119/288230376151711744*t+19188487056081357/288230376151711744）*（39/10-t）^2-100/81*（-176836856862157557/90071992547409920+4534278381080963/9007199254740992*t）*（t-3）^2
f0=
2.5851得三次样条插值函数
S（x）=1000/729*（27/5*x-1377/100）*（x-39/10）^2+1000/729*（522/25-24/5*x）*（x-3）^2+100/81*（-6396162352027119/288230376151711744*x+19188487056081357/288230376151711744）*（39/10-x）^2-100/81*（-176836856862157557/90071992547409920+4534278381080963/9007199254740992*x）*（x-3）^2
>>xi=0.9:
0.01:
13.3;yi=interp1（x,y,xi,'spline'）;
>>title（'试验一--三次样条插值图示'）
（2）用拉格朗日法插值
%定义Lagrange程序
functionf=Language（x,y,x0）
symst;
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等'）;
return;
end
f=0.0;
for（i=1:
n）
l=y（i）;
for（j=1:
i-1）
l=l*（t-x（j））/（x（i）-x（j））;
end;
for（j=i+1:
n）
l=l*（t-x（j））/（x（i）-x（j））;
end;
f=f+l;
simplify（f）;
if（i==n）
if（nargin==3）
f=subs（f,'t',x0）;
else
f=collect（f）;
f=vpa（f,6）;
end
end
end
>>Language（x,y）
ans=
52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2
二、已知Wilson矩阵
，且向量
，则方程组
有准确解
。
⑴用Matlab内部函数求
，
的所有特征值和
；
⑵令
，解方程组
，并求出向量
和
，从理论结果和实际计算结果两方面分析方程组
解的相对误差
与
的相对误差
的关系；
⑶再改变扰动矩阵
（其元素的绝对值不超过0.005），重复第2问。
解：
（1）A=[10787;7565;86109;75910];b=[32;23;33;31];
M=det（A）
M=
1
A的所以特征值：
>>D=eig（A）
D=
0.0102
0.8431
3.8581
30.2887
条件数>>n=30.2887/0.0102
n=
2.9695e+003
所以A的行列式为1，cond（A）2=2.9695e+003
（2）>>B=[1078.17.2;7.085.0465;85.989.899;6.994.9999.98];
>>b=[32;23;33;31];
>>[rank（B）,rank（[B,b]）]
ans=
44
>>x=B\b
x=-5.4761
11.4934
-1.4292
2.4838
求
假设X=
>>x=B\b;x1=[1;1;1;1];X=x-x1
X=
-6.4761
10.4934
-2.4292
1.4838
求
>>norm（X）
ans=
12.6552
12.655就是
。
%求
>>norm（X）/norm（x1）
ans=
6.3276
6.3276即为
>>norm（a）
ans=
0.2244
>>norm（A）
ans=
30.2887
>>norm（a）/norm（A）%求
ans=
0.0074
所以
=0.0074
inv（A）%求A的逆矩阵，下求
>>d=inv（A）;
>>norm（d）*norm（a）*norm（x）
ans=
288.4858
>>1-norm（d*（a））
ans=
-12.3693
>>288.4858/-12.3693
ans=
-23.3227
所以
=-23.322
>>norm（X）
ans=
0.0074
所以
>
（3）改变
，取a2=[00.0020.0010.003;0.0010.00400.001;0.003-0.004-0.0010;-0.001-0.0020-0.003]
B2=a2+A;C2=[Bb]
b=[32;23;33;31]
>>rank（B2）
ans=
4
>>rank（C2）
ans=
4
rank（B2）=rank（C2）
所以扰动矩阵有唯一解
>>x2=B2\b
x2=
1.0649
0.8893
1.0272
0.9859
>>x=B\b;x1=[1;1;1;1];X=x-x1%求
（设X=
）
X2=
-6.4761
10.4934
-2.4292
1.4838
norm（X2）%求
>>norm（X2）
ans=
12.6552
12.6552就是
>>norm（X）/norm（x1）%求
>>norm（X2）/norm（x1）
ans=
6.3276
所以
=6.3276
>>norm（a2）
ans=
0.0071
>>norm（a）/norm（A）%求
>>norm（a2）/norm（A）
ans=
2.3336e-004
所以
=2.3336e-004
norm（d）*norm（a2）*norm（x2）
ans=
9.0875
>1-norm（d*（a2））
ans=
0.6943
norm（X2）
ans=
12.6552
9.0875/0.6943
ans=
13.0887
所以
=13.0887
<
三、解三对角线性方程组的追赶法及其应用
⑴编写解三对角线性方程组的追赶法的通用程序，并应用于方程组
，检验程序的正确性；（解为
）
⑵求微分方程边值问题
的数值解（取步长
），并与精确解比较（精确解为
）。
说明：
离散化微分方程时，
解：
clearall;
a=[2,2,2,2,2];
b=[-1,-1,-1,-1];
c=[-1,-1,-1,-1];
r=[1,0,0,0,0];
n=length（a）;
b=[0,b];
u
（1）=r
（1）/a
（1）;
v
（1）=c
（1）/a
（1）;
fork=2:
n-1
u（k）=（r（k）-u（k-1）*b（k））/（a（k）-v（k-1）*b（k））;
v（k）=c（k）/（a（k）-v（k-1）*b（k））;
end
u（n）=（r（n）-u（n-1）*b（n））/（a（n）-v（n-1）*b（n））;
x（n）=u（n）;
fork=n-1:
-1:
1
x（k）=u（k）-v（k）*x（k+1）;
end
fprintf（'Èý¶Ô½Ç·½³Ì×éµÄ½âΪ\n'）
fork=1:
nfprintf（'x（%1d）=%10.8f\n',k,x（k））
end
>>zhuiganfa%调用追赶法
三对角方程组的解为
x
（1）=0.83333333
x
（2）=0.66666667
x（3）=0.50000000
x（4）=0.33333333
x（5）=0.16666667
和结果
很好地吻合。
4、公元1225年，比萨的数学家Leonardo研究了方程
，得到一个根
，没有人知道他用什么方法得到这个值。
对于这个方程，分别用下列方法：
⑴迭代法
；⑵迭代法
；⑶对⑴的Steffensen加速方法；⑷对⑵的Steffensen加速方法；⑸Newton法。
求方程的根（可取
），计算到Leonardo所得到的准确度。
解：
由题意编写m文件如下
function[x0,k,er,x]=diedai（g,x0,wucha,max）
%g是给定的迭代函数
%x0是给定的初始值
%wucha是规定的误差范围
%max是所应许的最大迭代次数
%k是迭代次数加1
%x是不动点近似值
%x（x1，x2...，xn）
X
（1）=x0;
fork=2:
max
X（k）=feval（'g',X（k-1））;
k,er=abs（X（k）-X（k-1））
x=X（k）;
if（erbreak;end
ifk==max
disp（'超出迭代次数'）;
end
end
其中定义的g函数是
functiony=g（x）;y=20/（x^2+2*x+10）;
在命令窗口中输入>>diedai（'g',1,10^（-9）,15）
k=
15
er=
1.410125245193683e-005
超出迭代次数
X=
Columns1through3
1.000000000000001.538461538461541.29501915708812
Columns4through6
1.401825309448601.354209390404291.37529809248738
Columns7through9
1.365929788170651.370086003401821.36824102361284
Columns10through12
1.369059812007481.368696397555521.36885768862873
Columns13through15
1.368786102577991.368817874396091.36880377314363
ans=1
取解为1.36880377314363
对于
，只需修改diedai.m文件中的g，把其改为g1，编写m文件functiony=g1（x）;
y=（20-2*x^2-x^3）/10;
在命令窗口中输入diedai（'g1',1,10^（-9）,30）
…..
k=
24
er=
1.36601568861328
k=
25
er=
1.36860974051282
k=
26
er=
1.36942327571766
取解为1.36860974051282
牛顿法：
编写m文件
function[p1,er,k,y]=ndf（f,df,p0,tol,max）
%f是非线性函数
%df是f的微商
%p0是初始值
%tol是给定的允许误差
%max是迭代的最大次数
%p1是牛顿法求得的近似解
%er是p0的误差估计
%k是迭代次数
%y=f（p1）
p0,feval（'f',p0）
fork=1:
max
p1=p0-feval（'f',p0）/feval（'df',p0）;
er=abs（p1-p0）;
p0=p1;
p1,er,k,y=feval（'f',p1）
if（（erbreak,end
end
定义函数m文件
functiony=f（x）
y=x^3+2*x^2+10*x-20;
定义微商函数m文件
functiony=df（x）
y=3*x^2+4*x+10;
在命令窗口输入
>>ndf（'f','df',1,10^（-9）,10）
p0=k=
11
ans=y=0.91756564217382
-7
p1=p1=1.411764705882351.36933647058824
er=er=
0.411764705882350.04242823529412
k=k=
24
y=y=
0.011148119412453.907985046680551e-014
p1=p1=
1.368808188617531.36880810782137
K=er=
31.776356839400251e-015
y=k=
1.704487072373695e-0065
p1=y=
1.368808107821370
er=ans=
.0796********
由结果知道牛顿法迭代到第三次已经达到所要求的精度
故方程的根为1.36880810782137
3.
编写Steffensen.m文件
i=2;x0=1;%设初始值
f=inline（'20/（x^2+2*x+10）'）;%迭代格式
y=f（x0）;
z=f（y）;
x1=x0-（y-x0）.^2/（z-2*y+x0）;
S.result=[x0;x1];
whileabs（x1-x0）>=1e-9%迭代精度
x0=x1;
y=f（x0）;
z=f（y）;
x1=x0-（y-x0）.^2/（z-2*y+x0）;
i=i+1;
S.result（i）=x1;
end
S.step=[（0:
i-1）]';
fprintf（'迭代步数为:
\t%d\n',i-1）;
forj=1:
i
fprintf（'%10d',S.step（j））;fprintf（'\t'）;
f
在命令窗口输入Steffensen
迭代步数为:
4
01.0000000
11.3708139
21.3688082
31.3688081
41.3688081
分析结果知，Steffensen迭代加速步数减少了，到第四步已经达到了精度要求。
4.把迭代格式改为（20-2*x^2-x^3）/10，保存，在命令窗口输入
Steffensen
迭代步数为:
5
01.0000000
11.3334921
21.3684154
31.3688081
41.3688081
51.3688081
分析结果知，Steffensen迭代加速步数减少了，到第五步已经达到了精度要求。
五、用不同的数值方法计算积分
的近似值，其中
⑴取不同的步长
，分别用复合梯形公式和复合辛普森公式计算积分，比较两个公式的计算效果，是否存在一个最小的
，使得精度不能再被改善？
⑵用龙贝格求积公式，取
，并打印出T-表。
解：
（1）用复合梯形公式，编写fhtx.m文件
functions=fhtx（f,a,b,n）
%f是被积函数
%ab分别为积分的上下限
%n是子区间的个数
%s是梯形总面积
h=（b-a）/n;
s=0;
fork=1:
（n-1）
x=a+k*h;
s=s+feval（'f',x）;
end
s=h*（feval（'f',a）+feval（'f',b））/2+h*s;
编写被积函数文件tf.m
functiony=f（x）;
y=（10/x）^2*sin（10/x）;
在命令窗口输入
>>fhtx（'tf',1,3,10）
ans=
-4.7789
>>fhtx（'tf',1,3,15）
ans=
-2.8604
>>fhtx（'tf',1,3,20）
ans=
-2.2208
用复化辛普森公式，编写fhxps.m文件
functions=fhxps（f,a,b,n）
%f是被积函数
%ab分别为积分的上下限
%n是子区间的个数
%s是梯形总面积
h=（b-a）/（2*n）;
s1=0;s2=0;
fork=1:
n
x=a+（2*k-1）*h;
s1=s1+feval（'f',x）;
end
fork=1:
n-1
x=a+2*k*h;
s2=s2+feval（'f',x）;
end
s=h*（feval（'f',a）+feval（'f',b）+4*s1+2*s2）/3;
在命令窗口输入
>>fhxps（'f',1,3,15）
ans=
-1.4136
>>fhxps（'f',1,3,20）
ans=
-1.4220
取n较大时
>>fhxps（'f',1,3,1000）
ans=
-1.42602475569236
>>fhtx（'tf',1,3,1000）
ans=
-1.42633620719670
>>fhtx（'tf',1,3,2000）
ans=
-1.
升级会员 