数值代数课设.docx
《数值代数课设.docx》由会员分享，可在线阅读，更多相关《数值代数课设.docx（23页珍藏版）》请在冰点文库上搜索。
数值代数课设
1.用matlab编写一个Gauss消去法的程序求解方程组
=
程序：
>>clear
>>A=[210;34-1;526]
A=
210
34-1
526
>>b=[1;0;1]
b=
1
0
1
>>[L,U]=lu（A）
L=
0.40000.07141.0000
0.60001.00000
1.000000
U=
5.00002.00006.0000
02.8000-4.6000
00-2.0714
>>c=L\b
c=
1.0000
-0.6000
0.6429
>>x=U\c
x=
0.8621
-0.7241
-0.3103
2、列主元Gauss消去法求解方程组
=
程序：
A=input（'请输入线性方程组的增广矩阵A='）;
n=length（A）-1;
x=zeros（n,1）;
aa=zeros（n,1）;
forj=1:
n
fori=1:
（n+1）
AA（j,i）=abs（A（j,i））;
end
end
fork=1:
（n-1）
fori=k:
n
aa（i-（k-1））=AA（i,k）;
end
fori=k:
n
ifAA（i,k）==max（aa）
break
end
end
ifAA（i,k）==0
break
fprintf（'方程组系数矩阵奇异\n'）;
else
forj=k:
（n+1）
jh=A（i,j）;
A（i,j）=A（k,j）;
A（k,j）=jh;
end
end
fenzi=A（k,k）;
forj=k:
（n+1）
A（k,j）=A（k,j）/fenzi;
end
forp=（k+1）:
n
jj=A（p,k）;
forj=k:
（n+1）
A（p,j）=A（p,j）-jj*A（k,j）;
end
end
end
ifk==（n-1）
x（n）=A（n,（n+1））/A（n,n）;
fori=（n-1）:
（-1）:
1
he=0;
forj=（i+1）:
n
he=he+A（i,j）*x（j）;
end
x（i）=A（i,（n+1））-he;
end
end
x
请输入线性方程组的增广矩阵A=[2401;3-511;1781]
x=
0.4032
0.0484
0.0323
3、用LU分解求解线形方程组
=
程序：
>>A=[21-15;164-1;3275;0-182]
A=
21-15
164-1
3275
0-182
>>b=[8;4;2;1]
b=
8
4
2
1
>>[L,U]=lu（A）
L=
0.6667-0.0625-0.66921.0000
0.33331.000000
1.0000000
0-0.18751.00000
U=
3.00002.00007.00005.0000
05.33331.6667-2.6667
008.31251.5000
0002.5038
>>[L,U,P]=lu（A）
L=
1.0000000
0.33331.000000
0-0.18751.00000
0.6667-0.0625-0.66921.0000
U=
3.00002.00007.00005.0000
05.33331.6667-2.6667
008.31251.5000
0002.5038
P=
0010
0100
0001
1000
>>y=L\（P*b）
y=
2.0000
3.3333
1.6250
7.9624
>>x=U\y
x=
-5.3063
2.3333
-0.3784
3.1802
4、用Cholesky分解求解线形方程组
=
程序：
>>clear
>>A=[63-2;351;-217]
A=
63-2
351
-217
>>b=[2;-1;1]
b=
2
-1
1
>>L=chol（A）
L=
2.44951.2247-0.8165
01.87081.0690
002.2783
>>Y=L'\b
Y=
0.8165
-1.0690
1.2332
>>X=L\Y
X=
0.9541
-0.8807
0.5413
5、利用Jacobic迭代法求解方程组
=
，eps=1.0
10
，
=
程序：
function[x,k,index]=Jacobi（A,b,eps,it_max）
ifnargin<4it_max=100;end
ifnargin<3eps;end
n=length（A）;k=0;
x=zeros（n,1）;y=zeros（n,1）;index=1;
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
index=0;return;
end
y（i）=y（i）/A（i,i）;
end
ifnorm（y-x,inf）break;
end
x=y;k=k+1;
end
运行：
>>A=[0.98-0.05-0.02;-0.04-0.90.07;-0.020.090.94]；
>>b=[1;1;1]；
>>[x,k,index]=Jacobi（A,b,1e-6,100）；
x=
0.9904
-1.0628
1.1867
k=
6
index=
1
6、利用G-S迭代法求解方程组
=
，esp=1.0×10
=
程序：
functiony=gaussseidel（A,b,x0,esp）
D=diag（diag（A））;
U=-triu（A,1）;
L=-tril（A,-1）;
G=（D-L）\U;
f=（D-L）\b;
y=G*x0+f;
n=1;
whilenorm（y-x0）>=esp&n<=1000
x0=y;
y=G*x0+f;
n=n+1;
end
y
n答案：
A=[1031;2-103;1310];
b=[14;-5;14];
x0=[0;0;0];
esp=1.0e-3;
gaussseidel（A,b,x0,esp）
y=
1.00003896860354
1.00002773079501
0.99998778390114
n=
6
ans=
1.00003896860354
1.00002773079501
0.99998778390114
7、利用SOR迭代法求解方程组
=
Eps=1.0*10^-6,
取
程序：
function[n,x]=sor（A,b,X,nm,w,ww）
%用超松弛迭代法求解方程组Ax=b
%输入：
A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：
x为求得的方程组的解构成的列向量，n为迭代次数
n=1;
m=length（A）;
D=diag（diag（A））;%令A=D-L-U,计算矩阵D
L=tril（-A）+D;%令A=D-L-U,计算矩阵L
U=triu（-A）+D;%令A=D-L-U,计算矩阵U
M=inv（D-ww*L）*（（1-ww）*D+ww*U）;%计算迭代矩阵
g=ww*inv（D-ww*L）*b;%计算迭代格式中的常数项
%下面是迭代过程
whilen<=nm
x=M*X+g;%用迭代格式进行迭代
ifnorm（x-X,'inf'）disp（'迭代次数为'）;n
disp（'方程组的解为'）;x
return;
%上面：
达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=x;n=n+1;
end
%下面：
如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp（'在最大迭代次数内不收敛!
'）;
disp（'最大迭代次数后的结果为'）;x
运行：
>>A=[0.78-0.02-0.12-0.14;-0.020.86-0.040.06;-0.12-0.040.72-0.08;-0.140.06-0.080.74];
b=[0.76;0.08;1.12;0.68];
c=1000;
d=1e-6;
f=1.03;
k=[0;0;0;0];
g=sor（A,b,k,c,d,f）
迭代次数为
n=
7
方程组的解为
x=
1.5350
0.1220
1.9752
1.4130
g=
7
8、用共轭梯度法求解上题的线性方程组，
Eps=1.0*10^-6,
9、利用幂法求：
的按模最大的特征值和对应的特征向量，限定最大的迭代步骤n=500，Eps=1.0×
，
=
程序：
function[k,lambda,Vk,Wc]=mifa（A,x0,eps,n）
lambda=0;k=1;Wc=1;state=1;V=x0;
while（（k<=n）&（state==1））
Vk=A*V;[mj]=max（abs（Vk））;mk=m;
tzw=abs（lambda-mk）;Vk=（1/mk）*Vk;
Txw=norm（V-Vk）;Wc=max（Txw,tzw）;V=Vk;lambda=mk;state=0;
if（Wc>eps）
state=1;
end
k=k+1;Wc=Wc;
end
if（Wc<=eps）
disp（'迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：
'）
else
disp（'迭代次数k已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：
'）
end
Vk=V;k=k-1;Wc;
主程序：
>>A=[123;234;345];
x0=[1;1;1];
esp=1.0e-6;
n=500;
[k,lambda,xk,Wc]=mifa（A,x0,eps,n）,
[k,lambda,xk,Wc]=mifa（A,x0,0.000001,500）,
[x,D]=eig（A）,Dzd=max（diag（D））,wuD=abs（Dzd-lambda）,wux=x（:
2）./xk,
答案：
迭代次数k已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：
k=
500
lambda=
9.6235
xk=
0.5247
0.7623
1.0000
Wc=
1.7764e-015
迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：
k=
8
lambda=
9.6235
xk=
0.5247
0.7623
1.0000
Wc=
1.5190e-007
x=
0.82770.40820.3851
0.1424-0.81650.5595
-0.54280.40820.7339
D=
-0.623500
0-0.00000
009.6235
Dzd=
9.6235
wuD=
9.2421e-009
wux=
0.7781
-1.0710
0.4082
10、利用反幂法求A=
的按模最小的特征值和对应的特征向量，限定最大的迭代步骤n=500，Eps=1.0×
，
=
程序：
function[k,lambdan,Vk,Wc]=ydwyfmf（A,V0,jlamb,jd,max1）
[n,n]=size（A）;A1=A-jlamb*eye（n）;jd=jd*0.1;RA1=det（A1）;
ifRA1==0
disp（'请注意：
因为A-aE的n阶行列式hl等于零，所以A-aE不能进行LU分解.'）
return
end
lambda=0;
ifRA1~=0
forp=1:
n
h（p）=det（A1（1:
p,1:
p））;
end
hl=h（1:
n）;
fori=1:
n
ifh（1,i）==0
disp（'请注意：
因为A-aE的r阶主子式等于零，所以A-aE不能进行LU分解.'）
return
end
end
ifh（1,i）~=0
disp（'请注意：
因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.'）
k=1;Wc=1;state=1;Vk=V0;
while（（k<=max1）&（state==1））
[LU]=lu（A1）;Yk=L\Vk;Vk=U\Yk;[mj]=max（abs（Vk））;
mk=m;Vk1=Vk/mk;Yk1=L\Vk1;Vk1=U\Yk1;
[mj]=max（abs（Vk1））;mk1=m;Vk2=（1/mk1）*Vk1;tzw1=abs（（mk-mk1）/mk1）;
tzw2=abs（mk1-mk）;Txw1=norm（Vk）-norm（Vk1）;
Txw2=（norm（Vk）-norm（Vk1））/norm（Vk1）;
Txw=min（Txw1,Txw2）;tzw=min（tzw1,tzw2）;Vk=Vk2;
mk=mk1;Wc=max（Txw,tzw）;Vk=Vk2;mk=mk1;state=0;
if（Wc>jd）
state=1;
end
k=k+1;%Vk=Vk2,mk=mk1,
end
if（Wc<=jd）
disp（'A-aE的秩R（A-aE）和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：
'）
else
disp（'A-aE的秩R（A-aE）和各阶顺序主子式值hl、迭代次数k已经达到最大迭代次数max1,按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：
'）
end
hl,RA1
end
end
[V,D]=eig（A,'nobalance'）,Vk;k=k-1;Wc;lambdan=jlamb+1/mk1;
主程序：
A=[120;4-32;514];
x0=[1;1;1];
[k,lambda,Vk,Wc]=ydwyfmf（A,x0,0.2,0.000001,500）
请注意：
因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.
A-aE的秩R（A-aE）和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：
hl=
0.8000-10.5600-21.7280
RA1=
-21.7280
V=
-0.1577-0.55310.3776
-0.3243-0.0508-1.0000
-1.00001.0000-0.1070
D=
5.112700
01.18360
00-4.2964
k=
7
lambda=
1.1836
Vk=
0.5531
0.0508
-1.0000
Wc=
5.3678e-009
升级会员 