病态线性方程组的求解

数值分析（Hilbert矩阵）病态线性⽅程组的求解Matlab程序（Hilbert 矩阵）病态线性⽅程组的求解理论分析表明，数值求解病态线性⽅程组很困难。

考虑求解如下的线性⽅程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n H h ?=，11ij h i j = +-，i ，j = 1,2，…，n 1. 估计矩阵的2条件数和阶数的关系2. 对不同的n ，取(1,1,,1)n x =∈，分别⽤Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，⽐较结果。

3. 结合计算结果，试讨论病态线性⽅程组的求解。

第1⼩题：condition.m %第1⼩题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出，1）在n⼩于等于13之前，图像近似直线log（cond（H））~1.519n-1.8332）在n⼤于13之后，图像趋于平缓，并在⼀定范围内上下波动，同时随着n的增加稍有上升的趋势第2⼩题：solve.m%m第2⼩题主程序N=4000;。

1 / 7 数值分析课程实验报告题目：病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。

考虑求解如下的线性方程组的求解Hx = b ，期中，期中H 是Hilbert 矩阵，()ij n n H h ´=，11ij h i j =+-，i ，j = 1,2，…，n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ，取(1,1,,1,1))nx =Î，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，比较结果。

3.结合计算结果，试讨论病态线性方程组的求解。

解答过程1.估计矩阵的2-条件数和阶数的关系矩阵的2-条件数定义为：1222()Cond A A A-=´，将Hilbert 矩阵带入有：1222()Cond H H H -=´调用自编的Hilbert_Cond 函数对其进行计算，取阶数n = 50，可得从，可得从1阶到50阶的2-条件数，以五位有效数字输出，其中前10项见表1。

表1.前十阶Hilbert 矩阵的2-条件数阶数1 2 3 4 5 2-条件数1 19.281 524.06 1.5514e+004 4.7661e+005 阶数6 7 8 9 10 2-条件数1.4951e+007 4.7537e+008 1.5258e+010 4.9315e+011 1.6025e+013 从表1可以看出，随着阶数每递增1，Hilbert 矩阵的2-条件数都至少增加一个数量级，但难以观察出明显的相依规律。

故考虑将这些数据点绘制在以n 为横轴、Cond （H ）2为纵轴的对数坐标系中（编程用Hilbert_Cond 函数同时完成了这个功能），生成结果如图1。

图1.不同阶数下Hilbert 矩阵的2-条件数分布条件数分布由图可见，当维数较小时，在y-对数坐标系中Cond （H ）2与n 有良好的线性关系；但n 超过10后，线性趋势开始波动，n 超过14后更是几乎一直趋于平稳。

病态方程组的求解理论分析表明，数值求解病态线性方程组很困难，考虑求解如下的线性方程组的求解：Hx=b ，其中H 是Hilbert 矩阵，H=(hij)n×n ，h ij =1i+j−1,ii ,jj =1,2…..nn1.估计矩阵的2-条件数和矩阵阶数的关系；2.对不同的n ，取x =(1,1,1…,1)∈R n ，利用Hx=b 得出b.然后分别用Guass 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代求解方程： Ax=b ；把近似解x (k )与准确解x 对比，计算结果；3.综合计算结果，讨论病态方程组的求解。

解：1.用MATLAB 计算1-9阶Hilbert 矩阵的2-条件数得到的结果如下表所示（cond2（Hnn ）表示n 阶Hilbert 矩阵的2-条件数）：对表中数据进行分析可发现，随矩阵的阶数的增加，2-条件数越来越大，后项与前项之比约为30，2-条件数呈近似指数式增长。

对cond2（H nn ）取以10为底的对数，可发现在前几项log10(cond2(Hn))随n 的增加约呈线性增加关系。

n 为20时，log10(cond2(Hn))与n 的关系与下图所示：此时可发现当n ≤13时，log10(cond2(Hn))与n 接近于线性关系，n >13后，log10(cond2(Hn))变化十分缓慢。

为了对图像的线性进行比较，对n 为1-13部分1-9阶hilbert 矩阵的2-条件数n 12 3 4 5 6 7 8 9 cond2(H nn ) 119.2815 524.0568 1.55E+04 4.77E+05 1.49E+07 4.75E+08 1.53E+10 4.93E+11进行线性拟合，结果如下图所示：从上图可看出，n较小时，log10(cond2(Hn))与n成近似线性关系。

继续增大n，n=200时，log10(cond2(Hn))与n如下图所示：由图可看出，n为200时，log10(cond2(Hn))与n的关系与n=20时相近，在n较小时log10(cond2(Hn))与n的关系接近线性关系，当n大于一定的值后log10(cond2(Hn))变化十分缓慢，但有上升的趋势。

第2次作业病态线性方程组的计算在科学技术、工程和经济等领域中都会遇到求解线性方程组问题。

在很多有广泛背景的数学问题中，如样条逼近、微分方程边值问题的差分法和有限元方法都要求求解线性方程组，而且通常会遇到求解大型方程组的问题。

在数值求解大型方程的过程中，一般会遇到很到计算上的问题，矩阵的病态性是其中之一。

矩阵的病态性质一般来自两个方面，一是由于矩阵本身的性质所引起，如希尔伯特矩阵。

这种矩阵，虽然阶数可能不大，但其条件数却非常大；另一类病态性质是由于矩阵的阶数所引起的。

矩阵在阶数比较小时，性质很好；但当阶数升高后，条件数则变得很大。

如用差分法或有限元法离散微分方程后得到的线性方程组。

它们的理论性质很好，甚至在理论上可证明它们是对称正定的。

然而，由于矩阵得阶数很高，可达到上百万阶，因而矩阵得条件数是相当大的。

方程的病态性质给方程组的数值求解带来巨大麻烦，甚至使某些看似理论性质很好的方程如对称正定矩阵，也不太可能得到通常意义下唯一的数值解。

病态矩阵是指矩阵有很小的扰动时，所引起解的变化很大的矩阵。

矩阵的病态性质一般可用矩阵的条件数来指示。

矩阵的条件数记为：-(cond)=AAA虽然矩阵的病态性，如上面提到的，可能来自两个方面,但在与具体问题相联系的求解过程中遇到的问题是基本一致。

下面这个方程组来自对偏微分方程（Stokes方程）的某种有限元方法离散后得到的。

该微分方程为⎪⎩⎪⎨⎧=Ω∈=Ω∈=∇+∆Ω∂0),( , 0 div ),( , ),(u y x u y x y x f p u 其中)),(),,((),(21y x u y x u y x u =为二维空间的向量函数，为标量函数),(y x p 。

为方便起见，现已经有方程组的系数矩阵和右端向量：这些数据存于数据文件‘stokes_coff_....mat ’中，可用Matlab 的load 命令调入即可。

关于数据文件中数据的说明：CoffA:方程组的系数矩阵；具有如下形式：⎪⎪⎭⎫ ⎝⎛=C B B ACoffA T变量BT 即为T B ；CoffA 为对称正定矩阵。

实验一病态线性代数方程组的求解1.估计Hilbert矩阵2-条件数与阶数的关系运行tiaojianshu.m输入m=10 可以得到如下表的结果2.选择不同维数，分别用Guass消去（LU分解），Jacobi迭代，GS 迭代，SOR迭代求解，比较结果。

说明：Hx=b，H矩阵可以由matlab直接给出，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。

对于Jacobi迭代，GS迭代，SOR迭代，取迭代初值x0为0向量，迭代精度eps=1.0e-6，迭代次数<100000, SOR迭代中w=1.2和0.8分别计算。

取不同的n值，得到如下结果:对于Guass法，可以看出来，随着n的增大，求解结果误差变大，这是因为随着n增大，系数矩阵的条件数变大，微小的扰动就容易造成很大的误差。

最后得不到精确解。

对于Jacobi迭代，计算结果为Inf，说明是发散的。

对于GS迭代和SOR迭代，结果是收敛的，但是可以看出迭代次数比较多，并且对于不同维数GS和SOR收敛速度不一样，有时候GS快，有时SOR快。

对SOR取不同的w迭代速度也不一样，存在一个最优的松弛因子w。

并且可以知道，迭代次数多少跟初值x0也有关系。

3.讨论病态问题求解的算法。

通过上面的实验分析，可以看出，求解病态矩阵的时候要小心，否则可能得不到所要求的精确度。

可以采用高精度运算，用双倍多倍字长，使得由于误差放大而损失若干有效数字位之后，还能保留一些有效位。

另外可以通过对原方程作某些预处理，降低系数矩阵的条件数，因为cond(aA)=cond(A),所以不能通过将每一个方程乘上相同的常数来达到这个目标，可考虑将矩阵的每一行和每一列分别乘上不同的常数，亦即找到可逆的对角阵D1和D2将方程组化为D1AD2y=D1b，x=D2y这称为矩阵的平衡问题，但是这样计算量比原问题本身要多。

或者通过变分原理将求解线性方程组的问题转化为等价的求解无约束函数最优化问题的极小值等等，可以参考[1]郑洲顺，黄光辉，杨晓辉求解病态线性方程组的混合算法实验一所编程序如下：1.求条件数tiaojianshu.mm=input ('input m:=') ;N=[1:m];for i=1:mn=N(i);H=hilb(n);k=cond(H);disp('矩阵的阶数')disp(n)disp('矩阵')disp(H)disp('矩阵的条件数')disp(k)end2.Guass法①Guass.mfunction guass(n)n=input('请输入系数矩阵的维数n=');H=hilb(n);x_exact=ones(n,1);b=H*x_exact;x=Doolittle(H,b)a=input('是否继续,继续请按1,结束请按0:');if a==1guass(n);else end②Doolittle.mfunction x=Doolittle(A,b)% LU分解求Ax=b 的解N=size(A);n=N(1);L=eye(n,n);%L的对角元素为1U=zeros(n,n); %U为零矩阵U(1,1:n)=A(1,1:n)%U第一行L(1:n,1)=A(1:n,1)/U(1,1)%L第一列for k=2:nfor i=k:nU(k,i)=A(k,i)-L(k,1:(k-1))*U(1:(k-1),i)endfor j=(k+1):nL(j,k)=(A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k) endendy=solveDownTriangle(L,b);%调用下三角系数矩阵求解函数x=solveUpTriangle(U,y);%调用上三角系数矩阵求解函数③solveDownTriangle.mfunction x=solveDownTriangle(A,b)%求下三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=1:nif (i>1)s=A(i,1:(i-1))*x(1:(i-1),1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end④solveUpTriangle.mfunction x=solveUpTriangle(A,b)% 求上三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=n:-1:1if (i<n)s=A(i,(i+1):n)*x((i+1):n,1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end3.Jacobi法function [x,n]=jacobi(a,x0)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0');eps=1.0e-6;%解的精度控制m=10000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;D=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);B=D\(L+U);f=D\b;x=x0;n=0;tol=1;while tol>=epsx=B*x0+f;n=n+1;tol=norm(x-x0);x0=x;if (n>=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)4.GS法function [x,n]=gauseidel(a,x0)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0');eps=1.0e-6;%解的精度控制m=100000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;D=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);G=(D-L)\U;f=(D-L)\b;x=x0;n=0;tol=1;while tol>=epsx=G*x0+f;n=n+1;tol=norm(x-x0);x0=x;if (n>=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)5.SOR法function [x,n]=SOR(a,x0,w)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0=');w=input('请输入w=');eps=1.0e-6;%解的精度控制m=100000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;if (w<=0||w>=2)error;return;endD=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=x0;n=0;tol=1;while tol>=epsx=B*x0+f;n=n+1;tol=norm(x-x0);x0=x;if (n>=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)。

数值实验3_3 病态的线性方程组求解
自动化系李琳琳 2004211068
1．选择问题的维数为6，分别用Gauss消去法、J迭代方法、GS迭代方法和SOR 迭代方法求解并与问题的解比较，所得结果见下表
2．逐步增大问题的维数，仍用上述方法求解，结果如下：
由上述计算结果可以看出：
病态方程组的数值求解必须小心进行，否则得不到所要求的准确度或不稳定。

由本题中所做的数值实验可以看出，对于系数矩阵为Hilbert矩阵的病态方程组，Guass(即LU分解)方法和J迭代法都是无效的，结果发散。

而GS迭代和SOR迭代方法则较有效的解决了这个问题，得到较为精确的结果。

另一方面，也可以说明GS方法和SOR方法在收敛性方面更有优势。

其中，当系数矩阵A对称正定时，GS法一定收敛，而J法却不一定；且采用最优松弛因子的SOR方法要比GS法和J法收敛快得多。

1 / 8数值分析实验六：解线性方程组的迭代法2016113 张威震1 病态线性方程组的求解1.1 问题描述理论的分析表明，求解病态的线性方程组是困难的。

实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b 的求解，其中系数矩阵H 为Hilbert 矩阵，,,1(),,,1,2,,1i j n n i j H h h i j n i j ⨯===+-这是一个著名的病态问题。

通过首先给定解（例如取为各个分量均为1）再计算出右端b 的办法给出确定的问题。

实验要求：（1）选择问题的维数为6，分别用Gauss 消去法、列主元Gauss 消去法、J 迭代法、GS 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数（至少到100），仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（3）讨论病态问题求解的算法1.2 算法设计首先编写各种求解方法的函数，Gauss 消去法和列主元高斯消去法使用实验5中编写的函数myGauss.m 即可，Jacobi 迭代法函数文件为myJacobi.m ，GS 迭代法函数文件为myGS.m ，SOR 方法的函数文件为mySOR.m 。

1.3 实验结果1.3.1 不同迭代法球求解方程组的结果比较选择H 为6*6方阵，方程组的精确解为x* = (1, 1, 1, 1, 1, 1)T ，然后用矩阵乘法计算得到b ，再使用Gauss 顺序消去法、Gauss 列主元消去法、Jacobi 迭代法、G-S 迭代法和SOR 方法分别计算得到数值解x1、x2、x3、x4，并计算出各数值解与精确解之间的无穷范数。

Matlab 脚本文件为Experiment6_1.m 。

迭代法的初始解x 0 = (0, 0, 0, 0, 0, 0)T ，收敛准则为||x(k+1)-x(k)||∞<eps=1e-6，SOR方法的松弛因子选择为w=1.3，计算结果如表1。

《科学与工程计算》实验报告学号：姓名：一、实验内容：考虑方程组Hx=b 的求解，其中系数矩阵H 为Hilbert 矩阵，nj i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯这是一个著名的病态问题。

通过首先给定解（例如取为各个分量均为1）再计算出右端b 的办法给出确定的问题。

实验要求：（1）选择问题的维数为6，分别用Jacobi 迭代法、GS 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？ （3）讨论病态问题求解的算法。

二、程序设计的基本思想、原理和算法描述：1、 算法 Jacobi 迭代法若A 为稀疏矩阵，只需遍历非零元素GS 迭代法若A 为稀疏矩阵，只需遍历非零元素每步迭代计算量相当于一次矩阵与向量的乘法；不需要保留上一步的迭代解，与Jacobi 迭代法计算量一样。

SOR 迭代法（稠密矩阵）2、 函数组成double max(double array[100]) 求数组中的最大值函数3、 输入/输出设计对于方程组Hx=b 的求解，系数矩阵H 为Hilbert 矩阵，矩阵中的数由下列函数生成。

nj i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯X*取一个特解[1,1,1, (1)b 数组由矩阵的每行元素相加得来。

4、符号名说明double c[100] 用来存储第k+1次和第k 次迭代结果的差值的绝对值 double x[100] 第k+1次迭代结果，储存解数组 double x0[100] 初始向量double r 第k+1次和第k 次迭代结果的差值的绝对值的最大值 double sum 矩阵方程变换后右侧值的和 int k 迭代次数double a[100][100] 存储Hilbert 矩阵 double b[100] 存储b 向量三、源程序及注释：Jacobi 迭代法#include<iostream> #include<math.h> #include <iomanip> #include <stdio.h> using namespace std; int n;double max(double array[100])//求最大值函数 {double a=array[0]; int i;for(i=1; i<n; i++) {if(a<array[i]) a=array[i]; return a; } }double c[100]= {0.0};double x[100]= {0.0}; //第k+1次迭代结果，储存解数组 double x0[100]= {0.0}; //初始向量double r,sum=0; int main(){double s=0;int i,k,j;//k为迭代次数double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cin>>n;cout<<"输出a数组："<<endl;double max(double array[100]);for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//矩阵中的数精确到六位}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];//矩阵每行的和}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;for(k=1;; k++){for(i=0; i<n; i++){for(j=0; j<n; j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+((b[i]-sum)/a[i][i]);//雅克比迭代方法计算方式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;r=max(c);if(r<s)//输出迭代次数{for(i=0; i<n; i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数："<<k<<endl;return 0;}}return 0;}GS迭代法#include<iostream>#include<math.h>#include <iomanip>#include <stdio.h>using namespace std;int n;double max(double array[100])//求最大值函数{double a=array[0];int i;for(i=1;i<n;i++){if(a<array[i])a=array[i];}return a;}main() {double s=0;double max(double array[100]);double c[100]={0.0};double x[100]={0.0};//第k+1次迭代结果，储存解数组double x0[100]={0.0};//初始向量int i,k,j;double r,sum=0;double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cout<<"输出a数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//矩阵中的数精确到六位}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];//矩阵每行的和}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;for(k=1;;k++){for(i=0;i<n;i++){for(j=0;j<n;j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+((b[i]-sum)/a[i][i]);//gs迭代方法计算方式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0;i<n;i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数："<<k<<endl;return 0;}}}SOR迭代法#include<iostream>#include<math.h>#include <iomanip>#include <stdio.h>using namespace std;int n;double max(double array[100]){double a=array[0];int i;for(i=1;i<n;i++){if(a<array[i])a=array[i];}return a;}int main() {double s=0,w=0;int i,k,j;//k为迭代次数double max(double array[100]);double c[100]= {0.0}; //double x[100]= {0.0}; //第k+1次迭代结果，储存解数组double x0[100]= {0.0}; //初始向量int i,k,j;double r,sum=0;double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cin>>n;cout<<"输出a数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//cout<<a[i][j]<<" ";}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;cout<<"输入松弛因子："<<endl;cin>>w;for(k=1;;k++){for(i=0;i<n;i++){for(j=0;j<n;j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+(w*(b[i]-sum)/a[i][i]);//sor迭代方法的计算公式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0;i<n;i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl; cout<<"迭代次数："<<k<<endl;return 0;}}}四、运行输出结果：JacobiGSSOR五、结果比较分析：说明：Hx=b，H矩阵可以由函数hi,j=1/(i+j-1)直接由程序生成，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。