病态线性方程组的求解

1、安徽工业大学数理科学与工程学院病态线性方程组的求解专业数学与应用数学班级数*班学号*姓名*指导教师*二。一五年五月一、设计目的：为了更加透彻的熟悉数值分析课程，学习各种数学软件的使用，锻炼自己 对知识的实际运用能力。二、引言：用直接法求解AX=F线性方程组，对于系数矩阵A对角占优是很有效的。方 程阶数不高时，人们经常使用；而当方程组阶数大时，由于积累误差，导致结果 失真。为了克服误差积累问题，通常用迭代法。它具有可达到所要求的精度和 对计算内存要求不大的优点，对求解大型线性方程组，迭代法计算时间远比直 接法少，所以在实际计算中，迭代法也被人们广泛使用。然而迭代法要研究迭 代格式的收敛性，如Ja

2、cobi迭代对系数矩阵为病态矩阵不收敛，为此我们提供 一种修改的Jacobi迭代，并给出一些数值例子来说明有较好的效果。三、解线性方程组迭代法的描述设线性方程组AX=P这里A:a j ,X:x i,F:f i,1 i,j n,为了更广泛 地应用，对A只限制实的非奇异矩阵，那么，若给定初值x，我们熟知的有： Jacobi 迭代：xk 1)=(fajXj(k)/a ii j 二 i, 1 - i - n四、求解病态线性方程组的另一种迭代解法设线性方程BX=F,这里系数矩阵B是病态的，指的是矩阵条件数是较大的。 条件数越大，就越难求得准确解，为此，我们将方程的两端同加DX项，那 么 相应的Jacob

3、i迭代有：X(k = (D A)TF (D - H)X(k)(1)这里，A为 B 的对角阵，即 A: b ii,H:b j j #i 14i,j n,记 M=(D + A)-1(D - H) , AX(k) = X(k4l) - X(k)，那么有：X(k) = M ：X (k,)=M k.：X(0)由此可见(1)式收敛，迭代矩阵M的谱半径应满足p(M )1 ,谱半径若用M的特征值判断，则较为繁锁；若 B不可约，根据线性代数对角占优简单迭 代必收敛的性质，为简单起见，我们取D为对角线阵，即D: d i,为保证收敛就得取di=Sign2 |bij |, bii ,符号Signa,b的含义是与b同号

4、，数值取a,这是充分条件，实际计算中有时可放宽处理，比如可取di=Sign Max|bij |, Di 或者 di =Sign Max|bij |, bii j # i ,因而相应 ijj的Jacobi迭代修改为：xi(k+1) =(f i-Z bijXj(k) +dixi(k)/(b ii +di)(2)j二i下面列举普通Jacobi迭代不收敛，解不出正确结果，而用修改的Jacobi迭代可求出满意结果的例子：例1 :2 -11-2A= 111 F= -3_11-2_ 0精确解为X=L1 -1 -1，普通Jacobi迭代不收敛，取di=Sign Max aij |, aii j # i ,1

5、J 3用(2)式迭代 40 步，X=L 1.000035 -1.000053 -1.000003例2 :方程系数为Hilbert阵H: hj =1/（i j 一1） F:fi 八 hj 1 i,j n精确解X=1 11T,普通Jacobi迭代发散，取dj=Sign | hj |,hii 用修改的（2）式迭代，n=1200时迭代10870步结果摘录如下：系数为Hilbert阵1200阶的结果x(i) iterate 10870 steps valuasion14910.9999691.0002780.9995930.9998561.000061.0001411.0001491.0001241.0

6、000861.0000491.0000021.0000021.0000021.0000021.0000021.0000021.0000021.0000021.0000021.0000018910.9999980.9999980.9999980.9999980.99999811910.9999980.9999980.9999980.9999980.9999981.0000041.0000041.0000041.0000041.0000041.0000041.0000041.0000041.0000041.000004iteration = 10870 max error = 4.07000000

7、0000462e-04五、结论与问题以上数值试验表明该迭代算法对求解病态线性方程组是有效的，具优点 是可达到预定的精度。求解病态方程组大条件数的系数矩阵，要获得较正确的解是很困难的。 本文的算法迭代虽能保证收敛，但与精确解的误差到底怎样，还应将解代入 原方程，衡量、检验求得解的可信程度；另外，如何克服求解病态线性方程 组收敛缓慢的问题，还有待于进一步工作。附录例1的matlab计算程序:clearm=40;n=3;A=2 -1 1;1 1 1;1 1 -2;F=-2 -3 0;x=zeros(n,m);for i=1:nx(i,1)=0;endfor i=1:nif i=1d(i)=max(a

8、bs(A(i,2:n)*A(i,i)/abs(A(i,i);elseif i=nd(i)=max(abs(A(i,1:n-1)*A(i,i)/abs(A(i,i);elsed(i)=max(max(abs(A(i,1:i-1),max(abs(A(i,i-1:n)*A(i ,i)/abs(A(i,i);endendfor k=2:mfor i=1:nx(i,k)=(F(i)-A(i,:)*x(:,k-1)+A(i,i)*x(i,k-1)+d(i)*x(i,k-1)/(A(i,i)+d(i);endendx(：,m)例2的matlab计算程序：clearn=1200;m=10870;for i=1:nfor j=1:nH(i,j)=1/(i+j-1);endendfor i=1:nf(i)=sum(H(i,:);endfor i=1:nd(i)=sum(abs(H(i,:)*H(i,i)/abs(