数值分析(Hilbert矩阵)病态线性方程组的求解Matlab程序
病态线性方程组的求解
病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n H h ⨯=，11ij h i j =+-，i ，j = 1,2，…，n1. 估计矩阵的2条件数和阶数的关系2. 对不同的n ，取(1,1,,1)nx =∈，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，比较结果。
3. 结合计算结果，试讨论病态线性方程组的求解。
1）估计矩阵的2-条件数和阶数的关系矩阵的2-条件数定义为：1222()Cond A A A-=⨯，将Hilbert 矩阵带入有：1222()Cond H H H -=⨯。
使用MA TLAB 自带的cond2函数进行计算并画出log10（cond2）和阶数n 的关系曲线如下：可见当n 小于13的时候，条件数的对数与阶数有较好的线性关系，但是随着阶数的提高，条件数趋于“稳定”地振荡。
但是事实上，n较大时，H矩阵已经奇异，直接使用cond函数计算结果可能存在不准确性。
原因是对于条件数过大的矩阵使用inv函数求逆矩阵是不可靠的，应该使用invhilb函数进行求逆，并代入定义式中求解，生成的结果如下所示。
线性度较好，可知，Hilbert矩阵的2-条件数会随其阶数n的增加呈指数增大趋势，因此当n较大时Hilbert矩阵将是严重病态的。
对不同的n，采用各种方法求解方程编写程序，选取n=2,3，5,10,15,20，迭代条件为迭代100000次或者是计算精度达到1e-6，若迭代次数少于设定最大值，表示相邻两次迭代达到精度要求，或者是计算结果超出范围。
X0取零向量，w取1.2，计算结果如下所示：由上可见，当n大于2时，Jacobi法已经不收敛，故其迭代次数已经没有意义。
当ｎ＝15时，Gauss消去法已经不收敛。
并且随着阶数的上升，gauss消去法的误差也随之上升。
数值分析中求解线性方程组的MATLAB程序(6种)
数值分析中求解线性方程组的MATLAB程序（6种）1.回溯法（系数矩阵为上三角）function X=uptrbk(A,B)%求解方程组，首先化为上三角，再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));2.系数矩阵为下三角function x=matrix_down(A,b)%求解系数矩阵是下三角的方程组n=length(b);x=zeros(n,1);x(1)=b(1)/A(1,1);for k=2:1:nx(k)=(b(k)-A(k,1:k-1)*x(1:k-1))/A(k,k);end3.普通系数矩阵（先化为上三角，在用回溯法）function X=uptrbk(A,B)%求解方程组，首先化为上三角，再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));4.三角分解法function [X,L,U]=LU_matrix(A,B)%A是非奇异矩阵%AX=B化为LUX=B，L为下三角，U为上三角%程序中并没有真正解出L和U，全部存放在A中[N,N]=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;for p=1:N-1[max1,j]=max(abs(A(p:N,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=R(p);R(p)=R(j+p-1);R(j+p-1)=d;if A(p,p)==0'A is singular.No unique solution'break;endfor k=p+1:Nmult=A(k,p)/A(p,p);A(k,p)=mult;A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endendY(1)=B(R(1));for k=2:NY(k)=B(R(k))-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);for k=N-1:-1:1X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);endL=tril(A,-1)+eye(N)U=triu(A)5.雅克比迭代法function X=jacobi(A,B,P,delta,max1);%雅克比迭代求解方程组N=length(B);for k=1:max1for j=1:NX(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)breakendendX=X';k6.盖斯迭代法function X=gseid(A,B,P,delta,max1);%盖斯算法，求解赋初值的微分方程N=length(B);for k=1:max1for j=1:Nif j==1X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);elseif j==NX(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);elseX(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);endenderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)break;endendX=X';k。
数值分析 第三章 基于MATLAB的科学计算—线性方程组1概述
科学计算—理论、方法及其基于ＭＡＴＬＡＢ的程序实现与分析 三、 解线性方程组（线性矩阵方程）解线性方程组是科学计算中最常见的问题。
所说的“最常见”有两方面的含义：１） 问题的本身是求解线性方程组；２） 许多问题的求解需要或归结为线性方程组的求解。
关于线性方程组B A x B Ax 1-=⇒=（１）其求解方法有两类：１） 直接法：高斯消去法（Gaussian Elimination ）； ２） 间接法：各种迭代法（Iteration ）。
１、高斯消去法１） 引例考虑如下（梯形）线性方程组：()⎪⎩⎪⎨⎧==+==+-=⇒⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--⇔⎪⎩⎪⎨⎧==-=+-5.0141315.3221122004301211214322332321321332321x x x x x x x x x x x x x x x 高斯消去法的求解思路：把一般的线性方程组（１）化成（上或下）梯形的形式。
２）高斯消去法——示例考虑如下线性方程组：⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---⇔⎪⎩⎪⎨⎧=++-=-+-=+-306015129101.2001.221113060129501.2001.221321321321321x x x x x x x x x x x x１） 第一个方程的两端乘12加到第二个方程的两端，第一个方程的两端乘－１加到第三个方程的两端，得⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--3060031110001.0001.00111321x x x２） 第二个方程的两端乘001.010-加到第三个方程的两端，得 ⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--60600311010001.0001.00111321x x x３） 从上述方程组的第三个方程依此求解，得()⎪⎩⎪⎨⎧==+-==+-=600300001.03100024011332321x x x x x x ３）高斯消去法的不足及其改进——高斯（全、列）主元素消去法在上例中，由于建模、计算等原因，系数2.001而产生0.0005的误差，实际求解的方程组为⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---306015129101.20005.22111321x x x ⎪⎩⎪⎨⎧===⇒70.4509.30142.2565321x x x 注：数值稳定的算法高斯列主元素消去法就是在消元的每一步选取（列）主元素—一列中绝对值最大的元取做主元素，高斯列主元素消去法是数值稳定的方法。
matlab演示病态方程组
matlab演示病态方程组
病态方程组是指具有高度敏感性和不稳定性的方程组，即使在输入数据上稍微的变化也会导致输出结果的巨大变化。
在 MATLAB 中演示病态方程组可以通过以下步骤进行：
步骤1，定义病态方程组。
首先，我们需要定义一个病态方程组。
例如，我们可以选择一个已知的病态方程组，如希尔伯特矩阵方程组。
希尔伯特矩阵是一个特定形式的方阵，其元素由公式H(i,j) = 1/(i+j-1)计算得出。
步骤2，求解方程组。
利用 MATLAB 的求解器，如 "\" 运算符或 `linsolve` 函数，对病态方程组进行求解。
在求解过程中， MATLAB 会自动选择适当的数值方法来解决病态方程组。
步骤3，观察数值稳定性。
在求解完成后，可以通过观察结果的数值稳定性来判断方程组
的病态程度。
可以计算残差向量的范数，或者利用条件数等指标来
评估方程组的数值稳定性。
步骤4，参数变化的影响。
进一步，可以尝试对输入数据进行微小的变化，然后观察方程
组解的变化情况。
这可以帮助理解病态方程组的敏感性和不稳定性。
步骤5，使用数值稳定的方法。
最后，可以尝试使用一些数值稳定的方法来解决病态方程组，
如正则化方法、迭代方法等。
这些方法可以在一定程度上减轻病态
方程组带来的数值计算困难。
通过以上步骤，可以在 MATLAB 中演示病态方程组的求解过程，并观察其数值稳定性和参数变化的影响。
这有助于加深对病态方程
组特性的理解，以及选择合适的数值方法来解决这类问题。
实验一用matlab求解线性方程组
实验1.1 用matlab 求解线性方程组第一节 线性方程组的求解 一、齐次方程组的求解rref （A ） %将矩阵A 化为阶梯形的最简式null （A ） %求满足AX ＝0的解空间的一组基，即齐次线性方程组的基础解系【例１】 求下列齐次线性方程组的一个基础解系，并写出通解：我们可以通过两种方法来解： 解法1：>> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果： ans=1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵，得化简后的方程⎪⎩⎪⎨⎧=+--=+--=-+-02200432143214321x x x x x x x x x x x x取x2，x4为自由未知量，扩充方程组为即提取自由未知量系数形成的列向量为基础解系，记所以齐次方程组的通解为解法2： clearA=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2];B=null(A, 'r') % help null 看看加个‘r ’是什么作用,若去掉r ，是什么结果？执行后可得结果： B=1 0 1 0 0 1 0 1⎩⎨⎧=-=-004321x x x x ⎪⎪⎩⎪⎪⎨⎧====44432221x x x x x x x x ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡11000011424321x x x x x x ,00111⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε,11002⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε2211εεk k x +=易见，可直接得基础解系所以齐次方程组的通解为二、非齐次线性方程组的求解 Matlab 命令的基本格式：X ＝A\b %系数阵A 满秩时，用左除法求线性方程组AX ＝b 的解注意：A/B 即为AB -1, 而A\B 即为A -1B.C =[A,b];D ＝rref(C) % 求线性方程组AX ＝b 的特解，即D 的最后一列元素【例2】 求下列非齐次线性方程组的解：,00111⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε,11002⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε⎪⎪⎪⎩⎪⎪⎪⎨⎧=+=++=++=++=+150650650651655454354332121x x x x x x x x x x x x x 2211εεk k x +=解： clearA=[5 6 0 0 0;1 5 6 0 0;0 1 5 6 0;0 0 1 5 6;0 0 0 1 5]; b=[1;0;0;0;1];format rational %采用有理数近似输出格式,比较format short 看看x=A\b执行后可得所求方程组的解. 作业：【第一题】 求下列非齐次线性方程组的通解.A=[1 2 3 1;1 4 6 2;2 9 8 3;3 7 7 2] B=[3;2;7;12] format rational x=A\B x =⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++1227737389222643324321432143214321x x x x x x x x x x x x x x x x42/31/2684838239393950-7/3【第二题】计算工资问题一个木工，一个电工，一个油漆工，三个人相互同意彼此装修他们自己的房子。
数值分析（Hilbert矩阵）病态线性方程组的求解Matlab程序
数值分析（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)nx =∈K ?，分别用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;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max); [GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x); nS_R;xSOR(n)=normS_R; xnSOR(n)=nS_R; normC_G=norm(C_G-x); nC_G;xCG(n)=normC_G; xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n); %消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k); endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多，迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多，迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多，迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多，迭代可能不收敛');break;endend。
数值分析实验报告——Hilbert矩阵的求解
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后更是几乎一直趋于平稳。
MATLAB线性方程组和代数方程数值求解
第10讲 线性方程组和代数方程数值求解 （第8章 MATLAB 方程数值求解）目的：一、 掌握矩阵的分解与线性方程组的解 二、 代数方程数值求解的方法。
—————————————————————————————————————— 一、掌握线性方程组和代数方程数值求解的方法。
（一）矩阵的分解 1、化简矩阵的计算工作量当方程组AX b =的系数矩阵为方阵且可逆时，方程组有唯一解。
1X A b −=，求1A −或者1A b −的方法是将()()1,~,rA E E A−或者()()1,~,rA b E A b −这两个方法都需要经过将矩阵A 化为单位矩阵E 这一过程，而将矩阵A 化为单位矩阵E 通常采用高斯消元法。
假设110a ≠利用11a 将红框内的元素化为0，在计算化简的过程中，蓝框内的元素也要一起同步计算。
重复这一过程，可以将A 从上往下化简为上三角矩阵：可以继续将这个矩阵从下往上化简为对角矩阵，在往上化简时，比如用nn b 化简红框元素，由于主对角线下元素是零，此时蓝框内的元素不会一起同步计算。
112233000000000000nn a b b b，（2）对比从上往下，和从下往上化简的过程可知，化简一般矩阵为单位矩阵的计算工作量远远大于将三角矩阵化简为单位矩阵的计算工作量。
2、方阵的LU 分解与方程组n A X b =的求解矩阵的LU 分解是指将矩阵A 分解为下三角矩阵L 与上三角矩阵U 的乘积，即A L U =⋅，方程组AX b =等同于LUX b =，解得11X U L b −−=，由于,L U 是三角矩阵，求逆时运算量比较小。
所以如果,L U 已知，用11X U L b −−=求解比用1X A b −=更加有效。
但如果,L U 未知，用11X U L b −−=求解方程组，还需要将A 分解为,L U ，如果求出的,L U 只用于求解一个方程组AX b =，显然并不比用1A −求解更有效。
矩阵病态线性代数方程组的求解
实验一病态线性代数方程组的求解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)。
数值分析希尔伯特病态线性方程组
病态线性代数方程组的求解理论的分析表明，求解病态的线性代数方程组是困难的。
考虑方程组Hx = b的求解，其中H为Hilbert矩阵，H=(hij)n⨯n，hij=1，i,j=1,2,...,n i+j-11. 估计Hilbert矩阵2-条件数与阶数的关系；2. 选择问题的不同维数，分别用Gauss消去法，Jacobi迭代，GS迭代和SOR迭代求解，比较结果；3. 讨论病态问题求解的算法。
解：1、取Hilbert矩阵阶数最高分别为n=20和n=100。
采用Hilbert矩阵的2-条件数作为实验的比较对象，画出的曲线如下图所示：lg(cond(Hn))nlg(cond(Hn))~n关系图lg(cond(Hn))n从图中可以看出，在n≤13之前，图像接近直线，在n＞13之后，图像趋于平缓，在一定的范围内上下波动。
为了比较图像的线性部分，作出一条直线与已知曲线进行比较。
比较直线的关系式为：lg(cond(Hn))=1.519n-1.833，结果下图所示。
lg(cond(Hn))~n关系图lg(cond(Hn))n从图2中可以看出，当n较小时，lg(cond(Hn))~n之间近似满足线性关系。
当n 继续增大到100时，lg(cond(Hn))~n关系图下图所示：lg(cond(Hn))~n关系图lg(cond(Hn))n从图中可以看出，图像的走势符合在n=20时的猜想，在n大于一定的值之后，图像趋于平缓，且在一定范围内震荡，同时又有一定上升趋势，但上升速度很慢。
2、选择不同的阶数n，设方程组的精确解为xz=(1,1,…,1)T进行计算，用四种方法解x_Guass1、x_Jacobi1、x_GS1、x_SOR1对比表如下表所示。
Gauss消去法求解：选择问题的阶数为3~8时，用Gauss消去法求得的解与精确解一致，当阶数为9~14时，解开始出现偏差，而且n越大，偏差越大。
用迭代法求解：取初始向量为1.2(1,1,…,1)T.无论n为多少阶，用Jacobi迭代方法迭代出现发散的不稳定现象，无法求解；用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，当迭代次数取得相当大（20000次）时解仍在精确解附近波动；取w=1.5，用SOR迭代方法迭代不发散，能求得解，收敛速度较GS迭代快一些，但仍非常缓慢。
MATLAB与数值分析第二部分—5线性方程组的数值方法(2013) (1)
b2 2 m21 1 109
10 9 0
（大数吃小数）
1 109
1 109
原因：小主元 导致计 算失败。
x2 1,
x1 0
（结果差异很大）
全主元消去法
1 每一步选绝对值最大的元素为主元素，保证 | mik | 。
Step k: ① 选取 | ai
a11 ... a1i det( Ai ) ... ... ... ai 1 ... aii
注：事实上，只要 A 非奇异，即 A1 存在，则可通过逐次 消元及行交换，将方程组化为三角形方程组，求出唯 一解。
问题?
• 消元法是按照系数矩阵的主对角线上的元素（主元） 进行消元。从而可能出现： （1）某个主元为零，导致消元过程无法进行。 （2）当某个主元的绝对值很小时，计算结果误 差很大。
这种方法需要计算n +1个n阶行列式并作 n次除法，而每个 n阶行列式计算需作(n - 1) ×n!次乘法。
计算量十分惊人!!
如n = 30,需2.38×1035 次乘法。理论上是绝对正确， 但n较大时，实际计算不可行。
在实际中如何利用计算机这一强有力工具求解线性方程？
§2.1
消元法
——计算机解线性方程组的一种直接法
9 1 10 例： 1 1
109 2
1 0
109 109
109 x2 1 , x1 0 9 10
在第一行 是个小数， 仍是小主 元 导致计 算失败
注意：这两个方程组 在数学上严格等价。 取a12为主元， 避免大吃小
矩阵三角分解法
用矩阵理论分析消元法：
m 1
数值分析——线性方程组直接解法Hilbert矩阵
数值分析第一次上机实习报告——线性方程组直接解法一、问题描述设 H n = [h ij ] ∈ R n ×n 是 Hilbert 矩阵, 即11ij h i j =+- 对n = 2,3,4，…13,(a) 取11n n x R ⨯⎛⎫ ⎪=∈ ⎪ ⎪⎝⎭，及n n b H x =，用Gauss 消去法和Cholesky 分解方法来求解n n H y b =，看看误差有多大.(b) 计算条件数：2()n cond H(c) 使用某种正则化方法改善(a)中的结果.二、方法描述1. Gauss 消去法Gauss 消去法一般用于系数矩阵稠密且没有任何特殊结构的线性方程组。
设H =[h ij ]，y = (y 1,y 2,…，y n )T . 首先对系数矩阵H n 进行LU 分解，对于k=1,2，…n,交替进行计算：1111),,1,,1(),1,2,,k kj kj kr rj r k ik ik ir rk r kk u h l u j k k n l a l u i k k n u -=-=⎧=-=+⎪⎪⎨⎪=-=++⎪⎩∑∑…… 由此可以得到下三角矩阵L=[l ij ]和上三角矩阵U=[u ij ]. 依次求解方程组Ly=b 和Ux=y ，111,1,2,,1(),,1,,1i i i ir r r n i i ir r r i ii y b l y i n x y u x i n n u -==+⎧=-=⎪⎪⎨⎪=-=-⎪⎩∑∑…… 即可确定最终解。
2. Cholesky 分解法对于系数矩阵对称正定的方程组n n H y b =，存在唯一的对角元素为正数的下三角矩阵L ，使得H=LL T 。
因此，首先对矩阵H n 进行Cholesky 分解，即1122111()1()j jj jj jk k j ij ij ik jk k jj l h l l h l l l -=-=⎧=-⎪⎪⎨⎪=-⎪⎩∑∑ 1,i j n =+… L 的元素求出之后，依次求解方程组Ly=b 和L T x=y ，即1111111(),2,3,i i i ik k k ii b y l y b l y i n l -=⎧=⎪⎪⎨⎪=-=⎪⎩∑… 11(),1,2,n n nn n i i ki k k i nn y x l x y l x i n n l =+⎧=⎪⎪⎨⎪=-=--⎪⎩∑…1 由此求得方程组n n H y b =的解。
数值分析实验-病态线性方程组的算法设计
数值分析课程实验报告 实验名称 病态线性方程组的算法设计
班级
学号 姓名 序号
任课教师 评分 一、 实验目的
1、 初步病态线性方程组的判定。
2、 初步了解常规方法在求解病态线性方程组时遇到的困难。
3、 针对病态问题设计求解算法并验证算法的有效性。
二、用文字或图表记录实验过程和结果
1、Hilbert 矩阵如下：
11/21/1/21/31/(1)()1/1/(1)1/(21)n ij n
n H h n n n ⎡⎤⎢⎥+⎢⎥==⎢⎥⎢⎥+-⎣⎦L L M M M L 其中1(1)ij h i j =+-，它是一个对称正定矩阵，并且()n cond H 随着n 的增加迅速增加，利用Matlab 分析如下：
可以发现在阶数不断增大
Hilbert 矩阵的条件数不断增大，
这样使得求解Hilbert 病态方程
变得非常困难，即使A 或b 有微小
扰动时，即使求解过程是精确进
行的（即没有舍入误差），所得的
解相对于原方程的解也会有很大
的相对误差。
这就需要提出病态
线性方程组的求解方法，对于一
般的方程求解常用的有高斯（选
主元）消去算法、高斯—赛德尔迭代。
本试验先使用用列主元高斯消去算法和高斯-赛德尔迭代算法求解线性方程组：
n
H x b = 其中11(,,),(1,2,,)n T n i ij j b b b b h i n ====∑L L 。
2、高斯列主元消去算法
(1)设计流程图：。
数值分析（宋）第1次大作业Hilbert矩阵病态问题研究
Hilbert 矩阵病态问题研究 （数值分析第一次大作业） 姓名：** 学号：** 班级：**1）Hilbert 矩阵的阶数n 与ln(())n cond H 的关系猜想：ln(())n cond H 与n 呈线性关系，其中()n cond H 按2范数计算。
绘制ln(())n cond H n 曲线。
分别取11050500n ≤≤、、，得到ln(())n cond H n 曲线如图1-1、图1-2及图1-3所示。
程序详见附录1。
图1-1. 110n ≤≤由图1-1可知，110n ≤≤，ln(())n cond H 是n 的线性函数，猜想正确。
图1-2. 150n ≤≤由图1-2知，当15n >时，ln(())n cond H 与n 之间的线性关系已经不存在，而且ln(())n cond H 的值大致在(40,50)内间波动，猜想与实际不完全相符。
图1-3. 1500n ≤≤图1-3进一步说明了ln(())n cond H 与n 之间的变化关系：当n 小于某一值（设该值为k ）时，ln(())n cond H 是n 的线性函数，而当n 大于k 时，随着n 的增大，ln(())n cond H 与n 间的线性关系不再成立，且其值在某一区间内波动。
为进一步确定k 的大小，绘制114n ≤≤时的曲线，如图1-4所示，可知k 的取值应为13。
图1-4. 114n ≤≤2）由n H 至ˆnH 的预处理 绘制ˆln(()/())n n cond H cond H n 曲线。
其中11ˆn nH D H D --=，D 为由n H 的对角元素开方构成的对角矩阵。
条件数按2范数计算。
程序详见附录2。
分别取11350500n ≤≤、、，得到如图2-1、图2-2和图2-3所示曲线。
由曲线图像可知：当Hilbert 矩阵的阶数12n ≤时，ˆln(()/())n ncond H cond H 随n 增大而逐渐减小，而n 继续增大时，ˆln(()/())n n cond H cond H 的取值将在区间(-7,4)内波动，且主要集中在(0,-3)区间内。
数值分析MATLAB科学计算—线性方程组
科学计算—理论、方法及其基于ＭＡＴＬＡＢ的程序实现与分析 三、 解线性方程组（线性矩阵方程）解线性方程组是科学计算中最常见的问题。
所说的“最常见”有两方面的含义：１） 问题的本身是求解线性方程组；２） 许多问题的求解需要或归结为线性方程组的求解。
关于线性方程组B A x B Ax 1-=⇒=（１）其求解方法有两类：１） 直接法：高斯消去法（Gaussian Elimination ）； ２） 间接法：各种迭代法（Iteration ）。
１、高斯消去法１） 引例考虑如下（梯形）线性方程组：()⎪⎩⎪⎨⎧==+==+-=⇒⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--⇔⎪⎩⎪⎨⎧==-=+-5.0141315.3221122004301211214322332321321332321x x x x x x x x x x x x x x x 高斯消去法的求解思路：把一般的线性方程组（１）化成（上或下）梯形的形式。
２）高斯消去法——示例考虑如下线性方程组：⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---⇔⎪⎩⎪⎨⎧=++-=-+-=+-306015129101.2001.221113060129501.2001.221321321321321x x x x x x x x x x x x １） 第一个方程的两端乘12加到第二个方程的两端，第一个方程的两端乘－１加到第三个方程的两端，得⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--3060031110001.0001.00111321x x x２） 第二个方程的两端乘001.010-加到第三个方程的两端，得 ⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--60600311010001.0001.00111321x x x３） 从上述方程组的第三个方程依此求解，得()⎪⎩⎪⎨⎧==+-==+-=600300001.03100024011332321x x x x x x ３）高斯消去法的不足及其改进——高斯（全、列）主元素消去法在上例中，由于建模、计算等原因，系数2.001而产生0.0005的误差，实际求解的方程组为⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---306015129101.20005.22111321x x x ⎪⎩⎪⎨⎧===⇒70.4509.30142.2565321x x x注：数值稳定的算法高斯列主元素消去法就是在消元的每一步选取（列）主元素—一列中绝对值最大的元取做主元素，高斯列主元素消去法是数值稳定的方法。
Hilbert矩阵病态线性代数方程组的求解
实验一病态线性代数方程组的求解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分别计算。
a. n=5b. n=8c. n=10d. n=15取不同的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这称为矩阵的平衡问题，但是这样计算量比原问题本身要多。
数值分析实验报告--实验6--解线性方程组的迭代法
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。
数值求解Hilbert病态线性方程组
病态线性代数方程组的求解理论的分析表明，求解病态的线性代数方程组是困难的。
考虑方程组Hx = b 的求解，其中H 为Hilbert 矩阵，n n ij h H ⨯=)(，11-+=j i h ij ，n j i ,...,2,1,=1. 估计Hilbert 矩阵2-条件数与阶数的关系；2. 选择问题的不同维数，分别用Gauss 消去法，Jacobi 迭代，GS 迭代和SOR 迭代求解，比较结果；3. 讨论病态问题求解的算法。
解： 1、取Hilbert 矩阵阶数最高分别为n=20和n=100。
采用Hilbert 矩阵的2-条件数作为实验的比较对象，画出的曲线如下图所示：lg(())n cond H n从图中可以看出，在n ≤13之前，图像接近直线，在n ＞13之后，图像趋于平缓，在一定的围上下波动。
为了比较图像的线性部分，作出一条直线与已知曲线进行比较。
比较直线的关系式为：833.1519.1))(lg(-=n H cond n ，结果下图所示。
nl g (c o n d (H n ))lg(cond(Hn))~n 关系图从图2中可以看出，当n 较小时，n H cond n ~))(lg(之间近似满足线性关系。
当n 继续增大到100时，n H cond n ~))(lg(关系图下图所示：从图中可以看出，图像的走势符合在n=20时的猜想，在n 大于一定的值之后，图像趋于平缓，且在一定围震荡，同时又有一定上升趋势，但上升速度很慢。
2、选择不同的阶数n ，设方程组的精确解为xz=(1,1,…,1)T进行计算，用四种方法解x_Guass1、x_Jacobi1、x_GS1、x_SOR1对比表如下nl g (c o n d (H n ))lg(cond(Hn))~n关系图nl g (c o n d (H n ))lg(cond(Hn))~n 关系图Gauss消去法求解：选择问题的阶数为3~8时，用Gauss消去法求得的解与精确解一致，当阶数为9~14时，解开始出现偏差，而且n越大，偏差越大。
给定解计算求解病态方程组HX = b
给定解计算求解病态方程组HX = b计算实习题：考虑方程组HX = b 的求解，其中系数矩阵为Hilbert 矩阵：,,1(),,,1,2,.1i j n n i j H h h i j n i j ⨯===+-通过给定解计算出右端项的办法确定具体问题。
matlab 源程序：*高斯解法*Guass.mclear,clc;n=？;A=zeros(n,n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);endendfor i=1:nb(i)=0;for j=1:nb(i)=1/(i+j-1)+b(i);endenddet=1;x=zeros(n,1);for k=1:n-1for i=k+1:nif A(k,k)==0error('Errorr');endm=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(k,k);for k=n:-1:1for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);endAbX*雅可比解法*Jacobi.mclear,clc;n=?;A=zeros(n,n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);endendfor i=1:nb(i)=0;for j=1:nb(i)=1/(i+j-1)+b(i);endendAD=diag(diag(A));E=tril(A,-1);F=triu(A,1);DEFx=zeros(n,1);B=-(D^(-1))*(E+F);b=b';g=(D^(-1))*b;y=B*x+g;eps=0.000001;k=1;while (norm(y-x)>=eps)&(k<1000) x=y;y=B*x+g;k=k+1;yEnd*高斯-雅可比解法*GS.mclear,clc;n=6;A=zeros(n,n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);endendfor i=1:nb(i)=0;for j=1:nb(i)=1/(i+j-1)+b(i);endendAD=diag(diag(A));E=tril(A,-1);F=triu(A,1);DEFx=zeros(n,1);B=-((D+E)^(-1))*F;b=b';g=((D+E)^(-1))*b;y=B*x+g;yeps=0.000001;k=1;while (norm(y-x)>=eps)&(k<1000) x=y;y=B*x+g;k=k+1;%for k=1:10%y=B*x+g';%if abs(x-y)<eps% break;% end%x=y;endy'*松弛法*SOR.mclear,clc;n=?;A=zeros(n,n);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);endendfor i=1:nb(i)=0;for j=1:nb(i)=1/(i+j-1)+b(i);endendt=b/A;tx=t;x=zeros(n,1);D=diag(diag(A));E=tril(A,-1);F=triu(A,1);DEFw=1;B=-(((1/w)*D+E)^-1)*((1-(1/w))*D+F); b=b';g=(((1/w)*D+E)^-1)*b;y=B*x+g;eps=0.000001;k=1;while (norm(y-x)>=eps)&(k<1000) x=y;y=B*x+g;k=k+1;y'end。
Hilbert矩阵病态性分析
【六年级作文】美丽的清华育英，我爱你_550字清华育英，是我的母校。
她坐落在一个幽静的小镇上，周围是绿树掩映，鲜花盛开的田野。
清晨，鸟儿在校园里欢快地唱着歌，给整个校园带来了无尽的生机和活力。
走进校园，映入眼帘的是一片宽敞明亮的操场。
操场的边上是一排整齐的花坛和绿树成荫的小道。
如果你仔细倾听，还能听到操场上响亮的朗朗读书声。
这是老师和同学们在一起学习，讨论，一同度过快乐的校园时光。
校园的另一头是一座高大宽敞的教学大楼。
教学大楼里面有一间件丑小鸭变成白天鹅一样美丽多彩的图书馆。
步入图书馆，满目的书架上摆满了各种各样的书籍，凝聚着知识的智慧。
在这里，我可以尽情地沉浸在书海中，阅读名著，了解世界，获取知识，增长见识，更好地了解自己。
在清华育英读书，我们不仅能学到知识，还能领略到美妙的文化艺术之美。
校园里有一间优雅美丽的音乐室，每天都能听到琴声、笛声、歌声在这里回荡。
学校还设有美术室、舞蹈室，有专业老师指导，为同学们的艺术特长潜质展现平台。
校园里随处可见同学们的临摹绘画，跳舞练琴。
清华育英是一个有爱的大家庭。
在这里，老师们像家长一样关爱着我们，时时刻刻帮助我们，指导我们，让我们茁壮成长。
同学们则像兄弟姐妹一样相互扶持，同甘共苦，共同度过每一个难关。
在这个大家庭里，我们也有很多美好的记忆。
每当到了盛夏，学校会组织我们去郊外游山玩水，感受大自然的风景，感受生活的美好。
每到冬季，学校还会组织我们冰雪运动，享受冬日的快乐。
而每年一次的艺术节和运动会更是我们的喜爱，我们可以尽情地展现自我，施展才华。
美丽的清华育英，我爱你！你是一座美丽的学府，你如同一片绿洲，培育着一代又一代的学子。
在这里我们可以学到知识，丰富我们的头脑；在这里我们可以结交朋友，拓宽我们的视野；在这里我们还可以享受阳光下的快乐。
我爱清华育英，我为自己是一名清华育英的学子而骄傲！愿清华育英越办越好，越来越美，愿我们的友谊天长地久！。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
画出Hilbert矩阵2-条件数的对数和阶数的关系
n=200时
n=20时
从图中可以看出，
1）在n小于等于13之前，图像近似直线
log（cond（H））~1.519n-1.833
2）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势
3.结合计算结果，试讨论病态线性方程组的求解。
第1小题：
condition.m%第1小题程序
t1=20;%阶数n=20
x1=1:t1;
y1=1:t1;
fori=1:t1
H=hilb(i);
y1(i)=log(cond(H));
end
plot(x1,y1);
xlabel('阶数n');
ylabel('2-条件数的对数（log(cond(H)）');
fork=i:n
A(j,k)=A(j,k)-l(j,i)*A(i,k);
end
b(j)=b(j)-l(j,i)*b(i);
end
end
%回代过程
x(n)=b(n)/A(n,n);
fori=n-1:-1:1
c=A(i,:).*x;
x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);
end
Jacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数
function[x,n]=Jacobi(A,b,x0,e,m)
n=length(A);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=D\(L+U);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=(D-w*L)\((1-w)*D+w*U);
f=(D-w*L)\b*w;
x=B*x0+f;
n=1;
whilenorm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
ifn>m
disp('SOR超松弛迭代次数过多，迭代可能不收敛');
nC_G;
xCG(n)=normC_G;
xnCG(n)=nC_G;
end
Gauss.m%Gauss消去法
functionx=Gauss(A,b)
n=length(b);
l=zeros(n,n);x=zeros(1,n);
%消去过程
fori=1:n-1
forj=i+1:n
l(j,i)=A(j,i)/A(i,i);
f=D\b;
x=B*x0+f;
n=1;
whilenorm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
ifn>m
disp('Jacobi迭代次数过多，迭代可能不收敛');
break;
end
end
G_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数
[J,nJ]=Jacobi(A,b,x0,e,max);
[GS,nGS]=G_S(A,b,x0,e,max);
[S_R,nS_R]=SOR(A,b,x0,e,max,w);
[C_G,nC_G]=CG(A,b,x0,e,max);
normG=norm(G'-x);
xGauss(n)=normG;
normJ=norm(J-x);
第2小题：
solve.m%m第2小题主程序
N=4000;
xGauss=zeros(N,1);
xJacobi=zeros(N,1);
xnJ=zeros(N,1);
xGS=zeros(N,1);
xnGS=zeros(N,1);
xSOR=zeros(N,1);
xnSOR=zeros(N,1);
xCG=zeros(N,1);
nJ;
xJacobi(n)=normJ;
xnJ(n)=nJ;
normGS=norm(GS-x);
nGS;
xGS(n)=normGS;
xቤተ መጻሕፍቲ ባይዱGS(n)=nGS;
normS_R=norm(S_R-x);
nS_R;
xSOR(n)=normS_R;
xnSOR(n)=nS_R;
normC_G=norm(C_G-x);
function[x,n]=G_S(A,b,x0,e,m)
n=length(A);
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
B=(D-L)\U;
f=(D-L)\b;
x=B*x0+f;
n=1;
whilenorm(x-x0)>e
x0=x;
x=B*x0+f;
n=n+1;
ifn>m
disp('Gauss-Seidel迭代次数过多，迭代可能不收敛');
break;
end
end
SOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子
function[x,n]=SOR(A,b,x0,e,m,w)
n=length(A);
xnCG=zeros(N,1);
forn=1:N;
x=ones(n,1);
t=1.1;%初始值偏差
x0=t*x;%迭代初始值
e=1.0e-8;%给定的误差
A=hilb(n);
b=A*x;
max=100000000000;%可能最大的迭代次数
w=0.5;%SOR迭代的松弛因子
G=Gauss(A,b);
title('2-条件数的对数（log(cond(H)）与阶数n的关系图');
t2=200;%阶数n=200
x2=1:t2;
y2=1:t2;
fori=1:t2
H=hilb(i);
y2(i)=log(cond(H));
end
plot(x2,y2);
xlabel('阶数n');
ylabel('2-条件数的对数（log(cond(H)）');
（
理论分析表明，数值求解病态线性方程组很困难。考虑求解如下的线性方程组的求解Hx = b，期中H是Hilbert矩阵， ， ，i，j = 1,2，…，n
1.估计矩阵的2条件数和阶数的关系
2.对不同的n，取 ，分别用Gauss消去，Jacobi迭代，Gauss-seidel迭代，SOR迭代和共轭梯度法求解，比较结果。
break;
end
end
CG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数
function[x,n]=CG(A,b,x0,e,m)