数值分析上机作业总.docx
《数值分析上机作业总.docx》由会员分享，可在线阅读，更多相关《数值分析上机作业总.docx（22页珍藏版）》请在冰豆网上搜索。
数值分析上机作业总
数值分析上机实验
一、解线性方程组直接法（教材49页14题）
追赶法程序如下：
functionx=followup（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
主程序如下：
functionzhunganfa
A=[2-2000000;-25-200000;0-25-20000;00-25-2000;000-25-200;0000-25-20;00000-25-2;000000-25];
b=[220/27;0;0;0;0;0;0;0];
x=followup（A,b）
计算结果：
x=
8.1478
4.0737
2.0365
1.0175
0.5073
0.2506
0.1194
0.0477
二、解线性方程组直接法（教材49页15题）
程序如下：
functiontiaojianshu（n）
A=zeros（n）;
forj=1:
1:
n
fori=1:
1:
n
A（i,j）=（1+0.1*i）^（j-1）;
end
end
c=cond（A）
d=rcond（A）
当n=5时
c=
5.3615e+005
d=
9.4327e-007
当n=10时
c=
8.6823e+011
d=
5.0894e-013
当n=20时
c=
3.4205e+022
d=
8.1226e-024
备注：
对于病态矩阵A来说，d为接近0的数；对于非病态矩阵A来说，d为接近1的数。
三、解线性方程组的迭代法（教材74页14题）
（1）用Jacobi迭代法求：
Jacobi迭代法程序如下：
function[x,n]=jacobi（A,b,x0,eps,varargin）
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin<3
error
return
elseifnargin==5
M=varargin{1};
end
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
B=D\（L+U）;
f=D\b;
x=B*x0+f;
n=1;
whilenorm（x-x0）>=eps
x0=x;
x=B*x0+f;
n=n+1;
if（n>=M）
disp（'Warning:
迭代次数太多，可能不收敛'）;
return;
end
end
本题主程序如下：
functionyakebidiedai
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];
[x,n]=jacobi（A,b,x0）
计算结果：
x=
1.0000
-2.0000
3.0000
-2.0000
1.0000
n=
67
经过67次迭代，得到最终结果
（2）用Gauss-Seidel迭代法求：
Gauss-Seidel迭代法程序如下：
function[x,n]=gauseidel（A,b,x0,eps,M）
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin==4
M=200;
elseifnargin<3
error
return;
end
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
G=（D-L）\U;
f=（D-L）\b;
x=G*x0+f;
n=1;
whilenorm（x-x0）>=eps
x0=x;
x=G*x0+f;
n=n+1;
if（n>=M）
disp（'Warning:
迭代次数太多，可能不收敛'）;
return;
end
end
本题主程序如下：
functiongaosidiedai
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
b=[12;-27;14;-17;12];
x0=[0;0;0;0;0];
[x,n]=gauseidel（A,b,x0）
计算结果：
x=
1.0000
-2.0000
3.0000
-2.0000
1.0000
n=
38
经过38次迭代，得到最终结果。
四、矩阵特征值与特征向量的计算（教材100页13题）
幂法求最大特征值的程序：
function[l,v,s]=pmethod（A,x0,eps）
ifnargin==2
eps=1.0e-6;
end
v=x0;
M=5000;
m=0;
l=0;
for（k=1:
M）
y=A*v;
m=max（y）;
v=y/m;
if（abs（m-l）l=m;
s=k;
return;
else
if（k==M）
disp（'收敛速度过慢'）;
l=m;
s=M;
else
l=m;
end
end
end
求解本题程序如下：
functionmifa
A=[19066-8430;6630342-36;336-168147-112;30-3628291];
x0=[0001]';
[l,v,s]=pmethod（A,x0）
求解结果：
l=
343.0000
v=
-0.6667
-2.0000
0
1.0000
s=
114
结论：
经过114次迭代，求得此矩阵的最大特征值为343.0000，及其对应特征向量为[-0.6667;-2.0000;0;1.0000]
五、函数逼近（教材164页16题）
本题采用最小二乘法进行拟合：
线性最小二乘法程序如下：
function[a,b]=LZXEC（x,y）
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等'）;
return;
end
A=zeros（2,2）;
A（2,2）=n;
B=zeros（2,1）;
fori=1:
n
A（1,1）=A（1,1）+x（i）*x（i）;
A（1,2）=A（1,2）+x（i）;
B（1,1）=B（1,1）+x（i）*y（i）;
B（2,1）=B（2,1）+y（i）;
end
A（2,1）=A（1,2）;
s=A\B;
a=s
（1）;
b=s
（2）;
首先利用1/y代替y，1/x代替x并采用线性最小二乘法求出a与b：
functionzuixiaoercheng
x1=[2356791011121416171920];
y1=[106.42108.26109.58109.50109.86110.00109.93110.59110.60110.72110.90110.76111.10111.30];
x2=1./x1;
y2=1./y1;
[b,a]=LZXEC（x2,y2）
计算结果：
b=
8.4169e-004
a=
0.0090
绘制图形：
fplot（'x/（0.0090*x+8.4169e-004）',[2,20]）
gridon
title（'最小二乘拟合'）
六、数值微分与数值积分（教材207页26题）
本题采用高斯—勒让德求积公式求解：
高斯—勒让德求积公式程序如下：
functionI=IntGauss（f,a,b,n,AK,XK）
if（n<5&&nargin==4）
AK=0;
XK=0;
else
XK1=（（b-a）/2）*XK+（（a+b）/2）;
I=（（b-a）/2）*sum（AK.*subs（sym（f）,findsym（f）,XK1））;
End
ta=（b-a）/2;
tb=（a+b）/2;
switchn
case0,
I=2*ta*subs（sym（f）,findsym（sym（f））,tb）;
case1,
I=ta*（subs（sym（f）,findsym（sym（f））,ta*0.5773503+tb）+...
subs（sym（f）,findsym（sym（f））,-ta*0.5773503+tb））;
case2,
I=ta*（0.55555556*subs（sym（f）,findsym（sym（f））,ta*0.7745967+tb）+...
0.55555556*subs（sym（f）,findsym（sym（f））,-ta*0.7745967+tb）+...
0.88888889*subs（sym（f）,findsym（sym（f））,tb））;
case3,
I=ta*（0.3478548*subs（sym（f）,findsym（sym（f））,ta*0.8611363+tb）+...
0.3478548*subs（sym（f）,findsym（sym（f））,-ta*0.8611363+tb）+...
0.6521452*subs（sym（f）,findsym（sym（f））,ta*0.3398810+tb）...
+0.6521452*subs（sym（f）,findsym（sym（f））,-ta*0.3398810+tb））;
case4,
I=ta*（0.2369269*subs（sym（f）,findsym（sym（f））,ta*0.9061793+tb）+...
0.2369269*subs（sym（f）,findsym（sym（f））,-ta*0.9061793+tb）+...
0.4786287*subs（sym（f）,findsym（sym（f））,ta*0.5384693+tb）...
+0.4786287*subs（sym（f）,findsym（sym（f））,-ta*0.5384693+tb）+...
0.5688889*subs（sym（f）,findsym（sym（f））,tb））;
end
本题计算程序如下（采用4个节点）：
functionshuzhijifen
a=1;
b=1;
fors=-5:
0.05:
5
q1（1,a）=IntGauss（'cos（1/2*x^2）',0,s,4）;
q2（1,b）=IntGauss（'sin（1/2*x^2）',0,s,4）;
a=a+1;
b=b+1;
end
plot（q1,q2）;
绘制图形：
七、非线性方程及非线性方程组的解法
本题采用牛顿法进行求解：
牛顿法解非线性方程程序如下：
function[r,n]=mulNewton（F,x0,eps）
ifnargin==2
eps=1.0e-4;
end
x0=transpose（x0）;
Fx=subs（F,findsym（F）,x0）;
var=findsym（F）;
dF=Jacobian（F,var）;
dFx=subs（dF,findsym（dF）,x0）;
r=x0-inv（dFx）*Fx;
n=1;
tol=1;
whiletol>eps
x0=r;
Fx=subs（F,findsym（F）,x0）;
dFx=subs（dF,findsym（dF）,x0）;
r=x0-inv（dFx）*Fx;
tol=norm（r-x0）;
n=n+1;
if（n>100000）
disp（'迭代步数太多，可能不收敛'）;
return;
end
end
本题解决方案如下：
首先，绘制此方程的图形，大概确定其与X轴的交点位置。
由于
，可以得出
因此绘制程序如下：
ezplot（'log（（513+0.6651*x）/（513-0.6651*x））-x/（1400*0.0918）',[-772,772,-10,10]）;
gridon
得到图形如下图所示：
经过放大后，发现图形与x轴的交点接近
处。
计算非零根：
令
为牛顿法接非线性方程的初值。
程序如下：
symsx
f=log（（513+0.6651*x）/（513-0.6651*x））-x/（1400*0.0918）;
x0=-765
[r,n]=mulNewton（f,x0）
解得：
r=-767.3861
x0=765
[r,n]=mulNewton（f,x0）
解得：
r=767.3861
结论：
此方程的两个非零根分别为：
八、常微分方程数值解法（教材266页19题）
本题分别采用四阶ADAMS预测—校正算法和经典RK法进行求解：
四阶ADAMS预测—校正算法如下：
functiony=DEYCJZ_yds（f,h,a,b,y0,varvec）
formatlong;
N=（b-a）/h;
y=zeros（N+1,1）;
x=a:
h:
b;
y
（1）=y0;
y
（2）=y0+h*Funval（f,varvec,[x
（1）y
（1）]）;
y（3）=y
（2）+h*Funval（f,varvec,[x
（2）y
（2）]）;
y（4）=y（3）+h*Funval（f,varvec,[x（3）y（3）]）;
fori=5:
N+1
v1=Funval（f,varvec,[x（i-4）y（i-4）]）;
v2=Funval（f,varvec,[x（i-3）y（i-3）]）;
v3=Funval（f,varvec,[x（i-2）y（i-2）]）;
v4=Funval（f,varvec,[x（i-1）y（i-1）]）;
t=y（i-1）+h*（55*v4-59*v3+37*v2-9*v1）/24;
ft=Funval（f,varvec,[x（i）t]）;
y（i）=y（i-1）+h*（9*ft+19*v4-5*v3+v2）/24;
end
经典RK算法程序如下：
functiony=DELGKT4_lungkuta（f,h,a,b,y0,varvec）
formatlong;
N=（b-a）/h;
y=zeros（N+1,1）;
y
（1）=y0;
x=a:
h:
b;
var=findsym（f）;
fori=2:
N+1
K1=Funval（f,varvec,[x（i-1）y（i-1）]）;
K2=Funval（f,varvec,[x（i-1）+h/2y（i-1）+K1*h/2]）;
K3=Funval（f,varvec,[x（i-1）+h/2y（i-1）+K2*h/2]）;
K4=Funval（f,varvec,[x（i-1）+hy（i-1）+h*K3]）;
y（i）=y（i-1）+h*（K1+2*K2+2*K3+K4）/6;
end
其中FUNVAL函数程序如下：
functionfv=Funval（f,varvec,varval）
var=findsym（f）;
iflength（var）<4
ifvar
（1）==varvec
（1）
fv=subs（f,varvec
（1）,varval
（1））;
else
fv=subs（f,varvec
（2）,varval
（2））;
end
else
fv=subs（f,varvec,varval）;
end
本题解决方案如下（步长h=0.1）：
程序：
functionchangweifen
symsxy;
z=-y+2*cos（x）;
yy1=DELGKT4_lungkuta（z,0.1,0,pi,1,[xy]）;
yy2=DEYCJZ_yds（z,0.1,0,pi,1,[xy]）;
a=0:
0.1:
pi;
yy3=cos（a）+sin（a）;
plot（a,yy1,'r*'）;
gridon;
holdon;
plot（a,yy2,'b+'）;
plot（a,yy3,'-go'）;
legend（'RK','Adams','准确解'）
整体图形：
局部放大图形：
继续放大：
数据结果：
yy1（RK）
yy2（ADAMS）
yy3（准确解）
yy1-yy3
yy2-yy3
1
1
1
0
0
1.094837464
1.1
1.094837582
-1.18389E-07
0.005162418
1.178735678
1.189000833
1.178735909
-2.30293E-07
0.010264924
1.250856361
1.266114065
1.250856696
-3.351E-07
0.01525737
1.310478904
1.324215642
1.310479336
-4.32217E-07
0.013736305
1.357007579
1.369446642
1.3570081
-5.21089E-07
0.012438542
1.389977487
1.401242287
1.389978088
-6.012E-07
0.011264199
1.409059202
1.419252042
1.409059875
-6.72089E-07
0.010192167
1.414062067
1.423285143
1.4140628
-7.33353E-07
0.009222343
1.404936093
1.41328163
1.404936878
-7.84658E-07
0.008344753
1.381772465
1.389323897
1.381773291
-8.25739E-07
0.007550607
1.344802625
1.351635461
1.344803481
-8.56415E-07
0.00683198
1.294395964
1.300578527
1.29439684
-8.76584E-07
0.006181686
1.231056128
1.236650238
1.231057014
-8.86229E-07
0.005593224
1.155415987
1.160477584
1.155416873
-8.85423E-07
0.005060711
1.068231314
1.072811013
1.068232188
-8.74325E-07
0.004578825
0.970373228
0.974516833
0.970374081
-8.53183E-07
0.004142752
0.862819494
0.866568452
0.862820316
-8.22334E-07
0.003748136
0.746644754
0.75003657
0.746645536
-7.82199E-07
0.003391034
0.623009788
0.626078402
0.623010521
-7.33279E-07
0.003067882
0.493149914
0.495926042
0.49315059
-6.76156E-07
0.002775452
0.358362651
0.360874088
0.358363262
-6.11484E-07
0.002510826
0.219994747
0.222266651
0.219995287
-5.39985E-07
0.002271364
0.079428728
0.081483866
0.079429191
-4.62442E-07
0.002054676
-0.061930915
-0.060071936
-0.061930535
-3.7969E-07
0.001858599
-0.202671764
-0.200990293
-0.202671471
-2.92614E-07
0.001681178
-0.341387584
-0.339866738
-0.341387382
-2.02132E-07
0.001520644
-0.476692371
-0.475316868
-0.476692262
-1.09196E-07
0.001375394
-0.607234205
-0.605990211
-0.607234191
-1.47751E-08
0.00124398
-0.731708756
-0.730583746
-0.731708836
8.01498E-08
0.00112509
-0.848872314
-0.847854952
-0.848872489
1.74596E-07
0.001017537
-0.95755422
-0.95663424
-0.957554488
2.6759E-07
0.000920247
结论：
通过图形和数据结果可以看出，在本题中利用经典RK方法获得解的精确度比4阶ADAMS方法要高很多。
《数值计算方法》
上机实验
姓名：
刘天博
学号：
升级会员 