數(shù)值分析第二次作業(yè)
學院：電子工程學院
基于matlab平臺的三種迭代法求解矩陣方程組
求解系數(shù)矩陣由16階Hilbert方程組構成的線性方程組的解，其中右端項為[2877/851,3491/1431,816/409,2035/1187,2155/1423,538/395,1587/1279,573/502,947/895,1669/1691,1589/1717,414/475,337/409,905/1158,1272/1711,173/244].要求：1）Gauss_Sedel迭代法；
2）最速下降法；
3）共軛梯度法；
4）將結果進行分析對比。
解：根據(jù)題目要求，編寫了對應算法的matlab程序，求解結果如下：（求解精度為10e-4，最大迭代次數(shù)1000）
1、方程的解：如下圖1所示
圖1
三種方法求解的結果對比
圖2
Gause_Sedel算法收斂特性
圖3
最速下降法收斂特性
圖3
共軛梯度法收斂特性
從圖中可以看到，在相同的最大迭代次數(shù)和預設求解精度條件下，共軛梯度算法僅需要4次迭代便可求出方程組的解，耗時0.000454秒，而且求出解的精度最高；Gauss_Sedel方法需要465次迭代，耗時0.006779秒，求解精度最差；最速下降法需要398次迭代，耗時0.007595秒，求解精度與共軛梯度算法差不多，因此兩者求出的解也幾乎相同。從中可以得出結論，共軛梯度算法無論從求解精度還是求解速度上都優(yōu)于其他兩種，最速下降法在求解精度上幾乎與共軛梯度算法持平，但求解速度更慢。Gauss_Sedel方法在求解精度和速度兩方面都最差。
具體的解為:
Gauss_Sedel迭代法:
(共需465次迭代，求解精度達到9.97e-5)
X=[0.***
1.01431732497804
1.05286123930011
0.***
0.93***38
0.***
1.00661848511341
1.03799789809258
1.05***4
1.062***
1.04857676431223
1.02856199041113
1.01999170162638
0.97***15
0.***
0.***].最速下降法:
(共需398次迭代，求解精度達到9.94e-5)
X=[0.***
1.0***00
0.***
0.***
0.***
1.00378022225329
1.0***78
1.01928337905816
1.02085909665194
1.01930314197028
1.01444777381651
1.00704058989297
0.***
0.***
0.***
0.***].共軛梯度法:
(共需4次迭代，求解精度達到3.98e-5)
X=[
0.***
1.02707840189049
0.***
0.***
0.986***7
1.00128902564234
1.0***14
1.02047386502293
1.02300905060565
1.02163015083975
1.01678089454399
1.00920310863874
0.***
0.***
0.***
0.***].Matlab程序
主程序:
clc;clear;
%%
本程序用于計算第二次數(shù)值分析作業(yè)，關于希爾伯特矩陣方程的解，用三種方法，分析并比較，也可推廣至任意n維的矩陣方程
%%
A=hilb(16);
%生成希爾伯特系數(shù)矩陣
b=[2877/851;3491/1431;816/409;2035/1187;2155/1423;538/395;1587/1279;573/502;947/895;1669/1691;1589/1717;414/475;337/409;905/1158;1272/1711;173/244];
%右端向量
M=1000;
%最大迭代次數(shù)
err=1.0e-4;
%求解精度
[x,n,xx,cc,jingdu]=yakebi_diedai(A,b,err,M);
%
雅克比算法求解
tic;
[x1,n1,xx1,cc1,jingdu1]=gauss_seidel(A,b,err,M);
%
gauss_seidel算法求解
toc;
tic;
[x2,n2,xx2,jingdu2]=zuisuxiajiangfa(A,b,err,M);
%
最速下降法求解
toc;
tic;
[x3,flag,jingdu3,n3]=bicg(A,b,err);
%
matlab內置雙共軛梯度算法求解
toc;
tic;
[x4,xx4,n4,jingdu4]=con_grad(A,b,err,M);
%
教材共軛梯度算法求解
toc;
%%
計算相應結果，用于作圖
%%
num=[1:16]＇;
jie=[num,x1,x2,x4];
%
三者的解對比
%
三者的收斂情況對比
num1=[1:n1]＇;
fit1=[num1,jingdu1＇];
num2=[1:n2]＇;
fit2=[num2,jingdu2＇];
num4=[1:n4]＇;
fit4=[num4,jingdu4＇];
子函數(shù)1(Gause_Sedel算法):
function
[x,n,xx,cc,jingdu]
=
gauss_seidel(A,b,err,M)
%
利用迭代方法求解矩陣方程
這里是高斯賽爾得迭代方法
%
A
為系數(shù)矩陣
b
為右端向量
err為精度大小
返回求解所得向量x及迭代次數(shù)
%
M
為最大迭代次數(shù)
cc
迭代矩陣普半徑
jingdu
求解過程的精度
n
所需迭代次數(shù)
xx
存儲求解過程中每次迭代產(chǎn)生的解
for
ii=1:length(b)
if
A(ii,ii)==0
x=＇error＇;
break;
end
end
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=(D-L)＼U;
cc=vrho(B);
%迭代矩陣普半徑
FG=(D-L)＼b;
x0=zeros(length(b),1);
x=B*x0+FG;
k=0;
xx(:,1)=x;
while
norm(A*x-b)>err
x0=x;
x=B*x0+FG;
k=k+1;
xx(:,k+1)=x;
if
k>=M
disp(＇迭代次數(shù)太多可能不收斂?。?;
break;
end
n=k;
jingdu(k)=norm(A*x-b);
end
end
子函數(shù)2(最速下降算法):
function
[x,n,xx,jingdu]=zuisuxiajiangfa(A,b,eps,M)
%
利用迭代方法求解矩陣方程
這里是最速下降迭代方法
%
A
為系數(shù)矩陣
b
為右端向量
err為精度大小
返回求解所得向量x及迭代次數(shù)
%
%
M
為最大迭代次數(shù)
jingdu
求解過程的精度
n
所需迭代次數(shù)
xx
存儲求解過程中每次迭代產(chǎn)生的解
x0=zeros(length(b),1);
r0=b-A*x0;
t0=r0＇*r0/(r0＇*A*r0);
x=x0+t0*r0;
r=b-A*x;
xx(:,1)=x;
k=0;
while
norm(r)>eps
r=r;
x=x;
t=r＇*r/(r＇*A*r);
x=x+t*r;
r=b-A*x;
k=k+1;
xx(:,k+1)=x;
if
k>=M
disp(＇迭代次數(shù)太多可能不收斂！＇);
break;
end
n=k;
jingdu(k)=norm(r);
end
end
子函31(共軛梯度法):
function
[x,xx,n,jingdu]=con_grad(A,b,eps,M)
%
利用迭代方法求解矩陣方程
這里是共軛梯度迭代方法
%
A
為系數(shù)矩陣
b
為右端向量
err為精度大小
返回求解所得向量x及迭代次數(shù)
%
M
為最大迭代次數(shù)
jingdu
求解過程的精度
n
所需迭代次數(shù)
xx
存儲求解過程中每次迭代產(chǎn)生的解
x0=zeros(length(b),1);
r0=b-A*x0;
p0=r0;
%
t0=r0＇*r0/(r0＇*A*r0);
%
x=x0+t0*r0;
%
r=b-A*x;
%
xx(:,1)=x;
k=0;
x=x0;
r=r0;
p=p0;
while
norm(r)>eps
x=x;
r=r;
p=p;
afa=r＇*r/(p＇*A*p);
x1=x+afa*p;
r1=r-afa*A*p;
beta=r1＇*r1/(r＇*r);
p1=r1+beta*p;
x=x1;
r=r1;
p=p1;
k=k+1;
xx(:,k)=x;
if
k>=M
disp(＇迭代次數(shù)太多可能不收斂?。?;
break;
end
n=k;
jingdu(k)=norm(r);
end
end
點擊咨詢