分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理及Matlab程序

1、、实验目的及题目1.1实验目的：(1) 学会用高斯列主元消去法，LU分解法，Jacobi迭代法和Gauss-Seide迭代法解线性 方程组。(2) 学会用Matlab编写各种方法求解线性方程组的程序。1.2实验题目：1. 用列主元消去法解方程组：乂 +冷 +3& =42x + x? X3 + X4 = 12x - x? _ X33& 二 _3I Xi 2X2 3X3 X4 二 42. 用LU分解法解方程组Ax=b，其中f48-24 0 -12)(424A =024126 201224-2163. 分别用Jacobi迭代法和Gauss-Seide迭代法求解方程组:10x x? +

2、 2X3 = 118X2 _X3 +3X4 = _11|2为 一 x210x3 = 6l' x-i 3x2x3 11x4 = 25二、实验原理、程序框图、程序代码等2.1实验原理2.1.1高斯列主元消去法的原理Gauss消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式：+02X2 十川+bmXn = g1b22X2+H|+b2nXn 巳2+bnnXn F这个过程就是消元，然后再回代就好了。具体过程如下:对于k =1,2, III,n -1，若akO,依次计算mik 二 a / aka(k_1) =ai(k)_mkakjk) b(5=b(k)-mikbkk),

3、i,j=k + 1,川,n然后将其回代得到：xn二心“心)nXk =(bkk) - v a：1%)/ ak,k = n -1,n -2山1 j 土 1以上是高斯消去但是高斯消去法在消元的过程中有可能会出现akk）二0的情况，这时消元就无法进行了,即使主元数akk）=O,但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差 的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元 素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。2.1.2直接三角分解法（LU分解）的原理先将矩阵A直接分解为A = LU则求解方程组的问题就等价于求解两个三角形

4、方程组。 直接利用矩阵乘法，得到矩阵的三角分解计算公式为：Uii =au,i =1,2,|)(,nhi = aii /uii ,i = 2,1H,n*k AUj =akj _迟 IkmUmj, j =k,k +1,川，n,k = 2,3Hnm=1k 4lik =ik -、TimUmk)/Ukk,i 二 k 1,k 2H,n且k = n由上面的式子得到矩阵A的LU分解后，求解Ux=y的计算公式为工 F* yy =b -送1必=2,3,川nXn = y / UnnnXi =（yi - ' UijXj）/Uii,i 二 n-1,n-2,|,1 Ij十以上为LU分解法。2.1.3Jacobi迭

5、代法和Gauss-Seidel迭代法的原理（1）Jcaobi 迭代设线性方程组的系数矩阵A可逆且主对角元素ai1,a22,.,ann均不为零,令D 二 diagQi,a22 ,.,ann并将A分解成A 二 A- D D从而(1)可写成Dx 二 D - A x b令x = Bk f1其中B = I - DA,fi二Db以为迭代矩阵的迭代法（公式）xL)= Bx(k)+ f1称为雅可比Jacobi）迭代法,其分量形式为x(k1) = I bi -aHi =1,2，n,nZj 二j出k(k) aij Xj二 0,1,2，其中 x（0）=（x1（0用L.xfH为初始向量（2）Gauss-Seidel

6、迭代由雅可比迭代公式可知，在迭代的每一步计算过程中是用 xk的全部分量来计算xk1的（k 十）（k 卅）所有分量，显然在计算第i个分量Xi时，已经计算出的最新分量X1，,X.没有被利用。把矩阵A分解成A = D - L - U（6）其中D - diag ,a22,ann ，- L,-U分别为的主对角元除外的下三角和上三角部分，于 是,方程组（1）便可以写成D - L x 二 Ux b即x = B2 x f2其中f2 = D - L 'bB2 二 D 一 L，U ,以为迭代矩阵构成的迭代法（公式）Xk1 ®kf2(8)称为高斯一塞德尔迭代法，用分量表示的形式为(k 卅)1 豪(

7、M1)J( k) 1Xi=1 bi EajXj - Z ajXjaiiz. i = 1,2，n, k = 0,1,2,.2.2程序代码2.2.1高斯列主元的代码function Gauss(A,b)%A为系数矩阵，b为右端项矩阵m, n=size(A);n=len gth(b);for k=1: n-1pt,p=max(abs(A (k:n ,k);%找出列中绝对值最大的数p=p+k-1;if p>kt=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

8、);%开始消元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=0Ab=A,b;endendx=zeros( n,1);%开始回代x( n)=b( n)/A( n,n);for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endfor k=1:n fprintf('x%d=%fn',k,x(k);end2.2.2 LU 分解法的程序function LU(A,

9、b) %A 为系数矩阵， b 为右端项矩阵 m,n=size(A);%初始化矩阵 A ，b，L 和 Un=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)*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=yy(

10、1)=b(1);for k=2:n y(k)=b(k)-L(k,1:k-1)*y(1:k-1);endx=zeros(n,1); x(n)=y(n)/U(n,n);for k=n-1:-1:1 x(k)=(y(k)-U(k,k+1:n)*x(k+1:n)/U(k,k);end for k=1:nfprintf('x%d=%fn',k,x(k);end2.2.3Jacobi 迭代法的程序epe为精function Jacobi(A,b,eps) %A 为系数矩阵， b 为后端项矩阵， 度m,n=size(A);D=diag(diag(A);% 求矩阵 DL=tril(A)-D;%

11、 求矩阵 LU=triu(A)-D;% 求矩阵 Utemp=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:n fprintf('x%d=%fn',k,x(k);end2.2.4Gauss-Seidel 迭代程序function Gauss_Seidel(A,b,eps)%A为系数矩阵，b为后端项矩阵，epe为精度m,n=size(A);D=diag(diag(A);%求矩阵

12、 DL=D-tril(A);%求矩阵 LU=D-triu(A);%求矩阵 Utemp=1;x=zeros(m,1);k=0;while abs(max(x)-temp)>eps temp=max(abs(x);k=k+1;%记录循环次数x=i nv(D-L)*U*x+i nv(D-L)*b;%Gauss_Seidel 的迭代公式endfor k=1:n fprintf('x%d=%fn',k,x(k);end三、实验过程中需要记录的实验数据表格3.1 第一题(高斯列主元消去)的数据>> A=1 1 0 3;2 1 -1 1; 3 -1 -1 3;-1 2 3

13、-1;>> b=4;1;-3;4;>> Gauss(A,b)x1=-1.333333x2=2.333333x3=-0.333333x4=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 =48.0000 -24.00000-12.0000012.00

14、0012.00006.00000014.0000-1.000000012.9286x1=0.521179x2=1.005525x3=-0.375691x4=-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.467396x2=-2.358678x3=0.657604x4=2.8423973.4第三题用Gauss_Seide迭代的数据>> A=10 -1 2 0;0 8 -1 3;2 -1 10

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

16、到想要的结果。(2) 第二题中在 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 argume n得不至 U想要的结果。4.2 解决方案(1) 针对第一题中由于程序的第二行漏了一个分号导致输出了m和n的值，第8行中将 Ab 改为 A 问题就解决了。(2) 由于程序后面出现了矩阵y故输出的事矩阵y的值，但是我们要的事x的值，故只需 要将y改成x,或者直接把y去掉就解决了问题。(3) 在 function 文件中命名不能出现“ -”应该将其改为下划线“ _”，所以将 M 文件名 “Gauss-Seidel(A,b,eps”) 改成“ Gauss_Seidel(A,b,ep