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

-.z.数值分析第二次作业学院：电子工程学院基于matlab平台的三种迭代法求解矩阵方程组求解系数矩阵由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〕将结果进展分析比照。解：根据题目要求，编写了对应算法的matlab程序，求解结果如下：〔求解精度为10e-4，最大迭代次数1000〕方程的解：如下列图1所示图1三种方法求解的结果比照图2Gause_Sedel算法收敛特性图3最速下降法收敛特性图3共轭梯度法收敛特性从图中可以看到，在一样的最大迭代次数和预设求解精度条件下，共轭梯度算法仅需要4次迭代便可求出方程组的解，耗时0.000454秒，而且求出解的精度最高；Gauss_Sedel方法需要465次迭代，耗时0.006779秒，求解精度最差；最速下降法需要398次迭代，耗时0.007595秒，求解精度与共轭梯度算法差不多，因此两者求出的解也几乎一样。从中可以得出结论，共轭梯度算法无论从求解精度还是求解速度上都优于其他两种，最速下降法在求解精度上几乎与共轭梯度算法持平，但求解速度更慢。Gauss_Sedel方法在求解精度和速度两方面都最差。具体的解为:Gauss_Sedel迭代法:(共需465次迭代，求解精度到达9.97e-5)*=[0.3192 1.804 1.011 0.7998 0.8838 0.3066 1.341 1.258 1.654 1.572 1.223 1.113 1.638 0.9515 0.4813 0.9182].最速下降法:(共需398次迭代，求解精度到达9.94e-5)*=[0.4322 1.900 0.0185 0.9243 0.3628 1.329 1.478 1.816 1.194 1.028 1.651 1.297 0.0809 0.4377 0.0912 0.1750].共轭梯度法:(共需4次迭代，求解精度到达3.98e-5)*=[0.9456 1.049 0.9853 0.1590 0.7607 1.234 1.914 1.293 1.565 1.975 1.399 1.874 0.5155 0.8859 0.6949 0.5005].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/895;1669/1691;1589/1717;414/475;337/409;905/1158;1272/1711;173/244];%右端向量M=1000;%最大迭代次数err=1.0e-4;%求解精度[*,n,**,cc,jingdu]=yakebi_diedai(A,b,err,M);%雅克比算法求解tic;[*1,n1,**1,cc1,jingdu1]=gauss_seidel(A,b,err,M);%gauss_seidel算法求解toc;tic;[*2,n2,**2,jingdu2]=zuisu*iajiangfa(A,b,err,M);%最速下降法求解toc;tic;[*3,flag,jingdu3,n3]=bicg(A,b,err);%matlab置双共轭梯度算法求解toc;tic;[*4,**4,n4,jingdu4]=con_grad(A,b,err,M);%教材共轭梯度算法求解toc;%%计算相应结果，用于作图%%num=[1:16]';jie=[num,*1,*2,*4];%三者的解比照%三者的收敛情况比照num1=[1:n1]';fit1=[num1,jingdu1'];num2=[1:n2]';fit2=[num2,jingdu2'];num4=[1:n4]';fit4=[num4,jingdu4'];子函数1(Gause_Sedel算法):function[*,n,**,cc,jingdu]=gauss_seidel(A,b,err,M)%利用迭代方法求解矩阵方程这里是高斯赛尔得迭代方法%A为系数矩阵b为右端向量err为精度大小返回求解所得向量*及迭代次数%M为最大迭代次数cc迭代矩阵普半径jingdu求解过程的精度n所需迭代次数**存储求解过程中每次迭代产生的解forii=1:length(b)ifA(ii,ii)==0*='error';break;endendD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=(D-L)\U;cc=vrho(B);%迭代矩阵普半径FG=(D-L)\b;*0=zeros(length(b),1);*=B**0+FG;k=0;**(:,1)=*;whilenorm(A**-b)>err*0=*;*=B**0+FG;k=k+1;**(:,k+1)=*;ifk>=Mdisp('迭代次数太多可能不收敛！');break;endn=k;jingdu(k)=norm(A**-b);endend子函数2(最速下降算法):function[*,n,**,jingdu]=zuisu*iajiangfa(A,b,eps,M)%利用迭代方法求解矩阵方程这里是最速下降迭代方法%A为系数矩阵b为右端向量err为精度大小返回求解所得向量*及迭代次数%%M为最大迭代次数jingdu求解过程的精度n所需迭代次数**存储求解过程中每次迭代产生的解*0=zeros(length(b),1);r0=b-A**0;t0=r0'*r0/(r0'*A*r0);*=*0+t0*r0;r=b-A**;**(:,1)=*;k=0;whilenorm(r)>epsr=r;*=*;t=r'*r/(r'*A*r);*=*+t*r;r=b-A**;k=k+1;**(:,k+1)=*;ifk>=Mdisp('迭代次数太多可能不收敛！');break;endn=k;jingdu(k)=norm(r);endend子函31(共轭梯度法):function[*,**,n,jingdu]=con_grad(A,b,eps,M)%利用迭代方法求解矩阵方程这里是共轭梯度迭代方法%A为系数矩阵b为右端向量err为精度大小返回求解所得向量*及迭代次数%M为最大迭代次数jingdu求解过程的精度n所需迭代次数**存储求解过程中每次迭代产生的解*0=zeros(length(b),1);r0=b-A**0;p0=r0;%t0=r0'*r0/(r0'*A*r0);%*=*0+t0*r0;%r=b-A**;%**(:,1)=*;k=0;*=*0;r=r0;p=p0;whilenorm(r)>eps*=*;r=r;p=p;afa=r'*r/(p'*A*p)