LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序.docx

《LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序.docx》由会员分享，可在线阅读，更多相关《LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序.docx（10页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、一、实验目的及题目1.1 实验目的：（1） 学会用高斯列主元消去法，LU 分解法，Jacobi 迭代法和 Gauss-Seidel 迭代法解线性方程组。（2） 学会用 Matlab 编写各种方法求解线性方程组的程序。1.2 实验题目：1. 用列主元消去法解方程组：x + x + 3x = 4 1242x + x - x + x = 112342 x - x - x + 3x = -31234-x1+ 2x2+ 3x - x = 4342. 用 LU 分解法解方程组 Ax = b, 其中 48-240-12 4 -24241212A = 4 ， b = 06202 -2 -66216 -2 3.

2、 分别用 Jacobi 迭代法和 Gauss-Seidel 迭代法求解方程组：10x - x + 2x= -118123 x - x + 3x = -112 234x - x12+10x = 63-x1+ 3x2- x +11x34= 25二、实验原理、程序框图、程序代码等2.1 实验原理2.1.1 高斯列主元消去法的原理+ b x = g+ b x = g1n n12n n2Gauss 消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式：b x + b x + 11 112 2b x +22 2b x = gnn nn这个过程就是消元，然后再回代就好了。具体过程如下

3、：对于k = 1,2, n -1 ，若a(k ) 0, 依次计算kk南昌航空大学数学与信息科学学院实验报告m = a(k ) / a(k )ikikkk第9页a(k +1) = a(k ) - ma( k ), nijijik kjb(k +1) = b(k ) - mb( k ) , i, j = k +1,然后将其回代得到：iix = b(n) / a( n)ik k,1 nnnnnx = (b(k ) - a(k ) x ) / a(k ) , k = n -1,n - 2,以上是高斯消去。kkkjjkkj =k +1但是高斯消去法在消元的过程中有可能会出现a( k ) = 0 的情况，

4、这时消元就无法进行了，kk即使主元数a( k ) 0, 但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差kk的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。2.1.2 直接三角分解法（LU 分解）的原理先将矩阵 A 直接分解为 A = LU 则求解方程组的问题就等价于求解两个三角形方程组。直接利用矩阵乘法，得到矩阵的三角分解计算公式为：, n2, nu = a , i = 1,2,1i1il = a / u, i =i1i111, nk -1u kj= a -kjm=1l ukm

5、 mj, j = k, k +1,n, k = 2,3,l = (a - k -1 l u) / u , i = k +1,k + 2, ikikim mkkkm=1, n且k n由上面的式子得到矩阵 A 的 LU 分解后，求解 Ux=y 的计算公式为 y = b 11n y = b - i-1 ly , i = 2,3,iix = yij jj =1/ u nnnn,1iix = ( y - n u x ) / u, i = n -1,n - 2,以上为 LU 分解法。ij jiij =i+12.1.3 Jacobi 迭代法和 Gauss-Seidel 迭代法的原理（1） Jcaobi 迭代

6、设线性方程组Ax = b(1)11的系数矩阵 A 可逆且主对角元素a,a ,.,a22nn均不为零,令并将 A 分解成D = diag ( a ,a1122,.,a)nn从而(1)可写成A = (A - D)+ DDx = (D - A)x + b(2)令x = B x + f11其中 B = I - D-1 A, f = D-1b .(3)111以 B 为迭代矩阵的迭代法(公式)x(k +1) = B x(k ) + f11称为雅可比(Jacobi)迭代法,其分量形式为(4)1 nx ( k +1 ) =iaiib - aii jj =1j ix ( k )ji = 1,2,.n,k = 0

7、,1,2,.(5)()12n其中 x(0 ) = x(0 ) ,x(0 ) ,.x(0 ) T 为初始向量.（2） Gauss-Seidel 迭代由雅可比迭代公式可知,在迭代的每一步计算过程中是用 x (k )的全部分量来计算 x (k +1)的所有分量,显然在计算第 i 个分量 x (k +1)时,已经计算出的最新分量 x (k +1) ,.,x(k +1)没有被利i用。把矩阵 A 分解成1i-1A = D - L - U(6)其中 D = diag ( a11 ,a22 ,.,ann ), - L,-U 分别为 A 的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成(D -

8、L)x = Ux + b即x = B2其中x + f2B = (D - L)-1U ,22以 B 为迭代矩阵构成的迭代法(公式)f = (D - L)-1 b2(7)x(k +1) = B x(k ) + f22(8)称为高斯塞德尔迭代法,用分量表示的形式为1 i-1nx( k +1 ) =iab - aiijx( k +1 ) - ajijx( k )jiij =1j =i+1i = 1,2,n,k = 0,1,2,.2.2 程序代码2.2.1 高斯列主元的代码function Gauss(A,b)%A 为系数矩阵，b 为右端项矩阵m,n=size(A); n=length(b); for

9、k=1:n-1pt,p=max(abs(A(k:n,k);%找出列中绝对值最大的数p=p+k-1; if pkt=A(k,:);A(k,:)=A(p,:);A(p,:)=t;%交换行使之变到主元位置上t=b(k);b(k)=b(p);b(p)=t;endm=A(k+1:n,k)/A(k,k);%开始消元A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1); if flag=0end endAb=A,b;x=zeros(n,1);%开始回代x(n)=b(n)/A(n

10、,n);for k=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);end for k=1:nfprintf(x%d=%fn,k,x(k);end2.2.2 LU 分解法的程序function LU(A,b)%A 为系数矩阵，b 为右端项矩阵m,n=size(A);%初始化矩阵 A，b，L 和 U n=length(b);L=eye(n,n);U=zeros(n,n);U(1,1:n)=A(1,1:n);%开始进行 LU 分解L(2:n,1)=A(2:n,1)/U(1,1);for k=2:nU(k,k:n)=A(k,k:n)-L(k,1:k-1)*

11、U(1:k-1,k:n);L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k)/U(k,k);endL%输出 L 矩阵U%输出 U 矩阵y=zeros(n,1);%开始解方程组 Ux=y y(1)=b(1);for k=2:ny(k)=b(k)-L(k,1:k-1)*y(1:k-1);end x=zeros(n,1);x(n)=y(n)/U(n,n); for k=n-1:-1:1x(k)=(y(k)-U(k,k+1:n)*x(k+1:n)/U(k,k);endfor k=1:nfprintf(x%d=%fn,k,x(k);end2.2.3 Jaco

12、bi 迭代法的程序function Jacobi(A,b,eps)%A 为系数矩阵，b 为后端项矩阵，epe 为精度m,n=size(A);D=diag(diag(A);%求矩阵 D L=D-tril(A);%求矩阵 LU=D-triu(A);%求矩阵 U temp=1;x=zeros(m,1); k=0;while abs(max(x)-temp)eps temp=max(abs(x);k=k+1;%记录循环次数x=inv(D)*(L+U)*x+inv(D)*b;%雅克比迭代公式endfor k=1:nfprintf(x%d=%fn,k,x(k);end2.2.4 Gauss-Seidel

13、迭代程序function Gauss_Seidel(A,b,eps)%A 为系数矩阵，b 为后端项矩阵，epe 为精度m,n=size(A);D=diag(diag(A);%求矩阵 D L=D-tril(A);%求矩阵 LU=D-triu(A);%求矩阵 U temp=1;x=zeros(m,1); k=0;while abs(max(x)-temp)eps temp=max(abs(x);k=k+1;%记录循环次数x=inv(D-L)*U*x+inv(D-L)*b;%Gauss_Seidel 的迭代公式endfor k=1:nfprintf(x%d=%fn,k,x(k);end三、实验过程中

14、需要记录的实验数据表格3.1 第一题（高斯列主元消去）的数据 A=1 1 0 3;2 1 -1 1; 3 -1 -1 3;-1 2 3 -1; b=4;1;-3;4; Gauss(A,b) x1=-1.333333 x2=2.333333 x3=-0.333333 x4=1.0000003.2 第二题（LU 分解法）数据 A=48 -24 0 -12;-24 24 12 12;0 6 20 2;-6 6 2 16; b=4; 4;-2;-2; LU(A,b)L =1.0000000-0.50001.00000000.50001.00000-0.12500.2500-0.07141.0000U

15、=48.0000-24.00000-12.0000012.000012.00006.00000014.0000-1.000000012.9286x1=0.521179 x2=1.005525 x3=-0.375691 x4=-0.2596693.3 第三题 Jacobi 迭代法的数据 A=10 -1 2 0;0 8 -1 3;2 -1 10 0;-1 3 -1 11;b=-11;-11;6;25;Jacobi(A,b,0.00005) x1=-1.467396 x2=-2.358678 x3=0.657604 x4=2.8423973.4 第三题用 Gauss_Seidel 迭代的数据 A=1

16、0 -1 2 0;0 8 -1 3;2 -1 10 0;-1 3 -1 11; b=-11;-11;6;25; Gauss_Seidel(A,b,0.00005) x1=-1.467357x2=-2.358740 x3=0.657597 x4=2.842405四、实验中存在的问题及解决方案4.1 存在的问题（1）第一题中在 matlab 中输入 Gauss(A,b)（数据省略）得到 m =4n =4Undefined function or variable Ab.Error in = Gauss at 8ap,p=max(abs(Ab(k:n,k); 没有得到想要的结果。（ 2 ）第二题中在

17、matlab 中输入 y=LU(A,b) 得到 y =4.00006.0000-5.0000-3.3571 不是方程组的解。（3）第三题中在用高斯赛德尔方法时在 matlab 中输入 Gauss-Seidel(A,b,eps)结果程序报错 Error using = GaussToo many output arguments.得不到想要的结果。4.2 解决方案（1） 针对第一题中由于程序的第二行漏了一个分号导致输出了 m 和 n 的值，第 8 行中将Ab 改为 A 问题就解决了。（2） 由于程序后面出现了矩阵y 故输出的事矩阵 y 的值，但是我们要的事x 的值，故只需要将 y 改成 x，或者

18、直接把 y 去掉就解决了问题。（3） 在 function 文件中命名不能出现“-”应该将其改为下划线“_”，所以将 M 文件名 “Gauss-Seidel(A,b,eps)”改成“Gauss_Seidel(A,b,eps)”就解决问题了。五、心得体会本次试验涉及到了用高斯列主元消去法，LU 分解法，Jacobi 迭代法以及 Gauss-seidel 迭代法等四种方法。需要对这些方法的原理都要掌握才能写出程序，由于理论知识的欠缺，我花了很大一部分时间在看懂实验的原理上，看懂了实验原理之后就开始根据原理编写程序， 程序中还是出现了很多的低级错误导致调试很久才能运行。通过这次试验使我深刻的体会到理论知识的重要性，没有理论知识的支撑是写不出程序来的。写程序时还会犯很多低级的错误，以后一定要加强理论知识的学习，减少编程时低级错误的产生。