数值分析实验.docx
《数值分析实验.docx》由会员分享，可在线阅读，更多相关《数值分析实验.docx（29页珍藏版）》请在冰豆网上搜索。
数值分析实验
数值实验
目录
数值实验1
第一章实验1
1.根号5的极限1
2.求pi的近似值2
3．画图y=sinx与其泰勒展开3
4．三维图mesh与surf实现。
4
第二章实验6
1.列主元三角分解A.6
2.追赶法求方程，电路的电流8
3.方程组的性态与矩阵条件数的实验9
4.Wilson矩阵。
求特征值，条件数。
误差12
第三章实验15
1.分别用Jacobi,Gauss-Seild,共轭梯度法解方程15
2.利用共轭梯度法计算矩阵--10^5阶19
3.利用cgs,bicg,bicgstab,等计算矩阵解20
第五章实验22
1.Newton插值f（x）=1/（1+4x^2）.图形，误差22
2.f（x）=1/（1+4x^2）插值24
3.飞机的外形轮廓24
第九章四阶龙格库塔方法（自选）25
1.解算微分方程组25
2.Matlab的四阶Rk方法25
3.作图比较，一定误差27
第一章实验
1.根号5的极限
代码如下:
x0=5;
fori=0:
1:
100
xi=sqrt（x0）;
ei=x0-xi;
fprintf（'NO:
%d,xi=%0.8f,Theerroris:
%0.8f\n',i,xi,ei）;
x0=xi;
ifei<10^（-8）
break;
end
end
可以看出极限是1.
2.求pi的近似值
代码如下：
sum0=0;
forn=1:
1:
50000
y=（-1）^（n+1）/（2*n-1）;
sum1=y+sum1;
pi_1=4*sum0;
pi_2=4*sum1;
error=pi_2-pi_1;
sum0=sum1;
fprintf（'NO:
%d,pi_2=%0.8f,Theerror:
%0.8f\n',n,pi_2,error）;
ifabs（error）<10^（-4）
break;
end
end
求得结果是：
3．画图y=sinx与其泰勒展开
x=0:
pi/100:
2*pi;
y=sin（x）;
y1=0;
y2=0;
y3=0;
fori=0:
2
y1=y1+（-1）^（i）*x.^（2*i+1）/factorial（2*i+1）;
end
fori=0:
5
y2=y2+（-1）^（i）*x.^（2*i+1）/factorial（2*i+1）;
end
fori=0:
10
y3=y3+（-1）^（i）*x.^（2*i+1）/factorial（2*i+1）;
end
plot（x,y,'*r',x,y1,'b',x,y2,'-g',x,y3,'k'）
axis（[02*pi-1.51.5]）
看出n=2发散，n=10几乎与y=sinx重合。
4．三维图mesh与surf实现。
首先用mesh函数.
fori=0:
1:
2;
k=[0.2,0.1,0.05];
[xi,yi]=meshgrid（-10:
k（1,i+1）:
10）;
z=exp（-abs（xi））+cos（xi+yi）+1./（xi.^2+yi.^2+1）;
subplot（2,2,i+1）
title（'mesh'）
mesh（xi,yi,z）
end
再用surf函数.
fori=0:
1:
2;
k=[0.2,0.1,0.05];
[xi,yi]=meshgrid（-10:
k（1,i+1）:
10）;
z=exp（-abs（xi））+cos（xi+yi）+1./（xi.^2+yi.^2+1）;
subplot（2,2,i+1）
title（'mesh'）
mesh（xi,yi,z）
end
第二章实验
1.列主元三角分解A.
clc
A=[11111;12345;1361015;14102035;15153570]
E=eye（5）
[m,n]=size（A）;
ifm~=n
error（'不是方阵'）
return
end
ifdet（A）==0
error（'不能分解'）
end
u=A;p=eye（m）;l=eye（m）;
fori=1:
m
forj=i:
m
t（j）=u（j,i）;
fork=1:
i-1
t（j）=t（j）-u（j,k）*u（k,i）;
end
end
a=i;b=abs（t（i））;
forj=i+1:
m
ifbb=abs（t（j））;
a=j;
end
end
ifa~=i
forj=1:
m
c=u（i,j）;
u（i,j）=u（a,j）;
u（a,j）=c;
end
forj=1:
m
c=p（i,j）;
p（i,j）=p（a,j）;
p（a,j）=c;
end
c=t（a）;
t（a）=t（i）;
t（i）=c;
end
u（i,i）=t（i）;
forj=i+1:
m
u（j,i）=t（j）/t（i）;
end
forj=i+1:
m
fork=1:
i-1
u（i,j）=u（i,j）-u（i,k）*u（k,j）;
end
end
end
l=tril（u,-1）+eye（m）
u=triu（u,0）
B=A/E
2.追赶法求方程，电路的电流
A=[2-2000000;-25-200000;0-25-20000;00-25-2000;
000-25200;0000-25-20;00000-25-2;000000-25]
bb=[220/270000000]
n=size（A,1）;
s=zeros（n,1）;
%-----È¡³öÈý¶Ô½Ç--------
b=diag（A）;
a=diag（A,-1）;
c=diag（A,1）;
d=zeros（n,1）;
u=zeros（n-1,1）;
fori=1:
n-1
d
（1）=b
（1）;
u（i）=c（i）/d（i）;
d（i+1）=b（i+1）-a（i）*u（i）;
end
%-----×·µÄ¹ý³Ì------------
y=zeros（n,1）;
y
（1）=bb
（1）/d
（1）;
fori=2:
n
y（i）=（bb（i）-a（i-1）*y（i-1））/d（i）;
end
%-----¸ÏµÄ¹ý³Ì---------------
s（n）=y（n）;
fori=n-1:
-1:
1
s（i）=y（i）-u（i）*s（i+1）;
end
s
3.方程组的性态与矩阵条件数的实验
clc
n=5;
A=zeros（n）;
C=zeros（n）;
b=zeros（1,n）;
fori=1:
n
x（i）=1+0.1*i;
forj=1:
n
A（i,j）=x（i）^（j-1）;
C（i,j）=1/（i+j-1）;
b（i）=b（i）+A（i,j）;
end
end
A
C
b
d=cond（A,2）
DD=A\b'
（1）当n=5,10,20时，cond（A）=5.3615e+05,8.6823e+11,1.5428e+22。
看出越来越病态了。
（2）当n=5,10,20时，解分别如下：
说明条件数越大，越病态，发散了。
当n=10时就一个值有一定误差。
（3）n=10,解方程（A+delta）x=b.
clc
n=10;
A=zeros（n）;
C=zeros（n）;
b=zeros（1,n）;
fori=1:
n
x（i）=1+0.1*i;
forj=1:
n
A（i,j）=x（i）^（j-1）;
%C（i,j）=1/（i+j-1）;
b（i）=b（i）+A（i,j）;
end
end
A（2,2）=A（2,2）+10^（-8）
A（10,10）=A（10,10）+10^（-8）
A
d=cond（A,2）
x=A\b'
那么结果如下：
可以看出很小的扰动导致解误差变得较大了。
4.Wilson矩阵。
求特征值，条件数。
误差
（1）行列式，条件数，特征值
（2）（A+A0）（x+x0）=b,求此时的x0与||x0||2.
代码如下：
clc
A=[10787;7565;86109;75910]
A1=[107.28.16.9;7.085.076.025;8.25.899.969.01;6.985.048.979.98]
b=[32233331]
A0=A1-A
x=A\b'
x1=A1\b'
x0=x1-x
fanshux0=sqrt（x0（1,1）^2+x0（2,1）^2+x0（3,1）^2+x0（4,1）^2）
fanshux=sqrt（1+1+1+1）
eigAA=eig（A'*A）;
fanshuA=sqrt（max（eigAA））
eigA0A0=eig（A0'*A0）;
fanshuA0=sqrt（max（eigA0A0））
x0_x=fanshux0/fanshux
A0_A=fanshuA0/fanshuA
condA=cond（A）
x0_xx=（condA*A0_A）/（1-condA*A0_A）
结果如下：
然后比较：
然后利用公式计算:
左边=0.8225，右边=1.0358
则有0.8225<1.0358。
第三章实验
1.分别用Jacobi,Gauss-Seild,共轭梯度法解方程
（1）Jacobi迭代法
clc
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]
b=[12-2714-1712]
n=length（A）;k=0;
ep=10^（-7）
it_max=100;
x=zeros（n,1）;y=zeros（n,1）;
fprintf（'NO:
x1x2x3x4x5\n'）
while1
fori=1:
n
y（i）=b（i）;
forj=1:
n
ifj~=i
y（i）=y（i）-A（i,j）*x（j）;
end
end
ifabs（A（i,i））<1e-10|k==it_max
return;
end
y（i）=y（i）/A（i,i）;
end
ifnorm（y-x,inf）break;
end
x=y;
k=k+1;
fprintf（'NO:
%d,%f%f%f%f%f\n',k,x（1,1）,x（2,1）,x（3,1）,x（4,1）,x（5,1））
end
得的结果为
（2）Gauss-Seild迭代法
clc
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]
b=[12-2714-1712]
x0=[00000]';
n=length（A）;
ep=10^（-7）;
N=100;
k=1;
x=zeros（n,1）;
fprintf（'\n**********start**************\n',k）;
fork=1:
N
fprintf（'%d\t',k）;
%calculatex
（1）
x
（1）=b
（1）/A（1,1）;
forj=2:
n
x
（1）=x
（1）-A（1,j）*x0（j）/A（1,1）;
end
fprintf（'%f\t',x
（1））;
%calculatex（i）
fori=2:
n-1
x（i）=b（i）/A（i,i）;
forj=1:
i-1
x（i）=x（i）-A（i,j）*x（j）/A（i,i）;
end
forj=i+1:
n
x（i）=x（i）-A（i,j）*x0（j）/A（i,i）;
end
fprintf（'%f\t',x（i））;
end
%calculatex（n）
x（n）=b（n）/A（n,n）;
forj=1:
n-1
x（n）=x（n）-A（n,j）*x（j）/A（n,n）;
end
fprintf（'%f\t',x（n））;
ifnorm（x-x0,inf）break;
else
k=k+1;
x0=x;
end
fprintf（'\n'）;
end
fprintf（'\n***********end*************\n',k）;
x0
运行结果如下：
可以看出Gauss-Seild只要41次就可以计算结果，而Jacobi需要76次。
（3）共轭梯度迭代法
clc
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]
b=[12-2714-1712]'
x0=[00000]';
N=length（A）;
eps=10^（-7）;
x=zeros（N,1）;%µü´ú½üËÆÏòÁ¿.normÊǾØÕó£¬ÏòÁ¿·¶Êý
r=b-A*x;%Æäʵ±íʾµÄÊǸºËÑË÷·½Ïò¡£fai=1/2x'AX-bx
d=r;
fork=0:
N-1
fprintf（'µÚ%d´Îµü´ú£º',k+1）;
a=（norm（r）^2）/（d'*A*d）
x=x+a*d
rr=b-A*x;%rr=r（k+1）
if（norm（rr）<=eps）||（k==N-1）
break;
end
B=（norm（rr）^2）/（norm（r）^2）;
d=rr+B*d;
r=rr;
end
计算结果如下：
理论上只需要N次得到精确值
2.利用共轭梯度法计算矩阵--10^5阶
clc
n=10^5;
fori=1:
n
forj=1:
n
ifi==j
A（i,j）=3;
elseifabs（i-j）==1
A（i,j）=-1;
elseif（i+j==n+1）&&（i~=n/2）&&（i~=n/2+1）
A（i,j）=1/2;
else
A（i,j）=0;
end
end
end
b（1,1）=2.5;
b（1,2:
10^5/2-1）=1.5;
b（1,10^5/2）=1.0;
b（1,10^5/2+1）=1.0;
b（1,10^5/2+2:
10^5-1）=1.5;
b（1,10^5）=2.5;
b=b';
N=length（A）;
eps=10^（-7）;
x=zeros（N,1）;%µü´ú½üËÆÏòÁ¿.normÊǾØÕó£¬ÏòÁ¿·¶Êý
r=b-A*x;%Æäʵ±íʾµÄÊǸºËÑË÷·½Ïò¡£fai=1/2*x'AX-bx
d=r;
fork=0:
N-1
fprintf（'µÚ%d´Îµü´ú£º',k+1）;
A=（norm（r）^2）/（d'*A*d）;
x=x+A*d;
x'
rr=b-A*x;%rr=r（k+1）
if（norm（rr）<=eps）||（k==N-1）
break;
end
B=（norm（rr）^2）/（norm（r）^2）;
d=rr+B*d;
r=rr;
end
3.利用cgs,bicg,bicgstab,等计算矩阵解
利用cgs（）函数
A=gallery（'wilk',21）
b=sum（A,2）;
tol=1e-12;maxit=15;
M1=diag（[10:
-1:
111:
10]）;
x=cgs（A,b,tol,maxit,M1）
利用bicg（）函数
A=gallery（'wilk',21）
b=sum（A,2）;
tol=1e-12;maxit=15;
M1=diag（[10:
-1:
111:
10]）;
x=bicg（A,b,tol,maxit,M1）
利用bicgstab（）函数
A=gallery（'wilk',21）
b=sum（A,2）;
tol=1e-12;maxit=10;
M1=diag（[10:
-1:
111:
10]）;
x=bicgstab（A,b,tol,maxit,M1）
第五章实验
1.Newton插值f（x）=1/（1+4x^2）.图形，误差
（1）输出插值多项式。
程序：
function[p]=NewtonChazhi（x,y,n）
%UNTITLED3此处显示有关此函数的摘要
%此处显示详细说明
f=zeros（n,n）;
f（:
1）=y';
forj=2:
n
fori=j:
n
f（i,j）=（f（i,j-1）-f（i-1,j-1））/（x（i）-x（i-j+1））;
end
end
p=f（n,n）;
fork=n:
-1:
2
p=conv（p,poly（x（k-1）））;
d=length（p）;
p（d）=p（d）+f（k-1,k-1）;
end
p
disp（'牛顿插值多项式：
'）;
polynomial=poly2sym（p,'x'）
end
x=-5:
1:
5;
y=1./（1+4*x.^2）;
n=length（x）;
[p]=NewtonChazhi（x,y,n）
结果：
牛顿插值多项式：
polynomial=
-（7319042784910035*x^10）/147573952589676412928+（256*x^8）/93425+（5*x^7）/1152921504606846976-（7410620163505401*x^6）/144115188075855872+x^5/9007199254740992+（36624*x^4）/93425+x^3/9007199254740992-（5148893614132309*x^2）/4503599627370496+（19*x）/36028797018963968+1
p=
Columns1through6
-0.000000.00270.0000-0.05140.0000
Columns7through11
1.39200.0000-1.14330.00001.0000
（2）在[-5,5]均匀插入99个节点，计算这些节点上函数f（x）的近似值，在同一张图上画出原函数与插值函数的图形。
程序：
x=-5:
1:
5;
y=1./（1+4*x.^2）;
n=length（x）;
[p]=NewtonChazhi（x,y,n）
X=linspace（-5,5,101）;
y1=1./（1+4*X.^2）;
fori=1:
101
N（i）=polyval（p,X（i））;
E（i）=N（i）-y1（i）;
end
plot（X,y1,X,N）
（4）误差图形
2.f（x）=1/（1+4x^2）插值
3.飞机的外形轮廓
y=[0,-30,-36,-35,-28.44,-9.4,0];
x2=520:
-8.67:
0;
x=[520,280,156.6,78,39.62,3.1,0];y2=spline（x,y,x2）;%用样条插值计算近似值
x1=[0,3.1,39.62,78,156.6,280,520];y1=[0,9.4,28.44,35,36,30,0];
x3=0:
8.67:
520;
y3=spline（x1,y1,x3）;
plot（x2,y2,'r-',x3,y3,'g-'）;title（'机翼外形曲线'）;
第九章四阶龙格库塔方法（自选）
四阶龙格-库塔（Runge-Kutta）方法，因为在本人研究方向上，要计算姿态结算方程组，而经典的方法就是四阶RK方法，故这章自己选题做这个方面的实验。
1.解算微分方程组
导航姿态方程如下：
角度与角速度关系
2.Matlab的四阶Rk方法
3.作图比较，一定误差
黑线表示实验测试计算结果，红线绿线表示matlab计算结果。
两者误差较大。
由于计算中全是三角函数计算，误差大。
升级会员 