基于matlab平台的三种迭代法求解矩阵方程

 学院：电子工程学院

基于matlab平台的三种迭代法求解矩阵方程组

求解系数矩阵由16阶Hilbert方程组构成的线性方程组的解，其中右端项为[2877/851,3491/1431,816/409,2035/1187,2155/1423,538/395,1587/1279,573/502,947/5,1669/1691,15/1717,414/475,337/409,905/1158,1272/1711,173/244].

要求：1）Gauss_Sedel迭代法；

      2）最速下降法；

      3）共轭梯度法；

      4）将结果进行分析对比。

解：根据题目要求，编写了对应算法的matlab程序，求解结果如下：（求解精度为10e-4，最大迭代次数1000）

1、方程的解：如下图1所示

图1 三种方法求解的结果对比

图2 Gause_Sedel算法收敛特性

图3 最速下降法收敛特性

图3 共轭梯度法收敛特性

从图中可以看到，在相同的最大迭代次数和预设求解精度条件下，共轭梯度算法仅需要4次迭代便可求出方程组的解，耗时0.000454秒，而且求出解的精度最高；Gauss_Sedel方法需要465次迭代，耗时0.006779秒，求解精度最差；最速下降法需要398次迭代，耗时0.007595秒，求解精度与共轭梯度算法差不多，因此两者求出的解也几乎相同。从中可以得出结论，共轭梯度算法无论从求解精度还是求解速度上都优于其他两种，最速下降法在求解精度上几乎与共轭梯度算法持平，但求解速度更慢。Gauss_Sedel方法在求解精度和速度两方面都最差。

具体的解为:

Gauss_Sedel迭代法:   (共需465次迭代，求解精度达到9.97e-5)

X=[0.995328360833192    1.01431732497804    1.05286123930011    0.934006974137998    0.931493373808838    0.966508138403066    1.00661848511341.0379********.0518********    1.06215849948572    1.048576731223.028*********    1.01999170162638    0.971831831519515    0.952526166634813    0.916996019179182].

最速下降法:    (共需398次迭代，求解精度达到9.94e-5)

X=[0.998835379744322    1.01507463472900    0.9825093720185     0.980191460759243    0.991245169713628    1.00378022225329    1.01350884374478    1.01928337905816.020*********    1.01930314197028    1.01444777381651    1.0070405297    0.998384452250809    0.9873994044377    0.975767814970912    0.963209150871750].

共轭梯度法:    (共需4次迭代，求解精度达到3.98e-5)

X=[ 0.9972751179456.027*********    0.977623373409853     0.973206695321590    0.986133032967607    1.0012025234    1.01322158496914.020*********.023*********.021*********    1.016780454399    1.00920310863874    0.999772406055155    0.988443827498859    0.976094192496949    0.962844741655005].

Matlab程序

主程序:

clc;clear;

%%  本程序用于计算第二次数值分析作业，关于希尔伯特矩阵方程的解，用三种方法，分析并比较，也可推广至任意n维的矩阵方程 %%

A=hilb(16);   %生成希尔伯特系数矩阵

b=[2877/851;3491/1431;816/409;2035/1187;2155/1423;538/395;1587/1279;573/502;947/5;1669/1691;15/1717;414/475;337/409;905/1158;1272/1711;173/244]; %右端向量

M=1000;   %最大迭代次数

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'];

子函数1(Gause_Sedel算法):

function [x,n,xx,cc,jingdu] = gauss_seidel(A,b,err,M)

%   利用迭代方法求解矩阵方程   这里是高斯赛尔得迭代方法

%   A 为系数矩阵 b 为右端向量  err为精度大小   返回求解所得向量x及迭代次数  

%   M 为最大迭代次数  cc 迭代矩阵普半径  jingdu 求解过程的精度 n 所需迭代次数  xx 存储求解过程中每次迭代产生的解

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('迭代次数太多可能不收敛！');

        break;

    end

  n=k; 

  jingdu(k)=norm(A*x-b);

end

end

子函数2(最速下降算法):

function [x,n,xx,jingdu]=zuisuxiajiangfa(A,b,eps,M)

%   利用迭代方法求解矩阵方程   这里是最速下降迭代方法

%   A 为系数矩阵 b 为右端向量  err为精度大小   返回求解所得向量x及迭代次数  

%   %   M 为最大迭代次数   jingdu 求解过程的精度  n 所需迭代次数  xx 存储求解过程中每次迭代产生的解 

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('迭代次数太多可能不收敛！');

        break;

    end

n=k;

jingdu(k)=norm(r);

end

end

子函31(共轭梯度法):

function [x,xx,n,jingdu]=con_grad(A,b,eps,M)

%   利用迭代方法求解矩阵方程   这里是共轭梯度迭代方法

%   A 为系数矩阵 b 为右端向量  err为精度大小   返回求解所得向量x及迭代次数  

%   M 为最大迭代次数    jingdu 求解过程的精度 n 所需迭代次数  xx 存储求解过程中每次迭代产生的解

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('迭代次数太多可能不收敛！');

        break;

end

   n=k;

jingdu(k)=norm(r);

end   

end