一、实验目的及题目1.1 实验目的：（1）学会用高斯列主元消去法，LU 分解法，Jacobi 迭代法和Gauss-Seidel 迭代法解线性方程组。

（2）学会用Matlab 编写各种方法求解线性方程组的程序。

1.2 实验题目：1. 用列主元消去法解方程组：1241234123412343421233234x x x x x x x x x x x x x x x ++=⎧⎪+-+=⎪⎨--+=-⎪⎪-++-=⎩2. 用LU 分解法解方程组,Ax b =其中4824012242412120620266216A --⎛⎫ ⎪- ⎪= ⎪ ⎪-⎝⎭，4422b ⎛⎫ ⎪ ⎪= ⎪- ⎪-⎝⎭3. 分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解方程组：1232341231234102118311210631125x x x x x x x x x x x x x -+=-⎧⎪-+=-⎪⎨-+=⎪⎪-+-+=⎩二、实验原理、程序框图、程序代码等2.1实验原理2.1.1高斯列主元消去法的原理Gauss 消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式：1111221122222n n n n nn n nb x b x b x g b x b x g b x g +++=⎧⎪++=⎪⎨⎪⎪=⎩这个过程就是消元，然后再回代就好了。

具体过程如下： 对于1,2,,1k n =-，若()0,k kka ≠依次计算()()(1)()()(1)()()/,,1,,k k ik ik kk k k k ij ij ik kjk k k i i ik k m a a a a m a b b m b i j k n++==-=-=+然后将其回代得到：()()()()()1/()/,1,2,,1n n n n nn n k k k k k kj j kk j k x b a x b a x a k n n =+⎧=⎪⎨=-=--⎪⎩∑以上是高斯消去。

但是高斯消去法在消元的过程中有可能会出现()0k kka =的情况，这时消元就无法进行了，即使主元数()0,k kka ≠但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散。

因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。

然后换行使之变到主元位置上，再进行销元计算。

即高斯列主元消去法。

2.1.2直接三角分解法（LU 分解）的原理先将矩阵A 直接分解为A LU =则求解方程组的问题就等价于求解两个三角形方程组。

直接利用矩阵乘法，得到矩阵的三角分解计算公式为：1111111111,1,2,,/,2,,,,,1,,,2,3,()/,1,2,,i i i i k kj kj km mj m k ik ik im mk kkm u a i n l a u i nu a l u j k k n k nl a l u u i k k n k n-=-===⎧⎨==⎩⎧=-=+⎪⎪=⎨⎪=-=++≠⎪⎩∑∑且由上面的式子得到矩阵A 的LU 分解后，求解Ux=y 的计算公式为11111,2,3,/()/,1,2,,1i i i ij j j n n nn n i i ij j ii j i y b y b l y i nx y u x y u x u i n n -==+=⎧⎪⎨=-=⎪⎩=⎧⎪⎨=-=--⎪⎩∑∑以上为LU 分解法。

2.1.3Jacobi 迭代法和Gauss-Seidel 迭代法的原理 （1）Jcaobi 迭代设线性方程组b Ax = (1)的系数矩阵A 可逆且主对角元素nn a ,...,a ,a 2211均不为零,令 ()nn a ,...,a ,a diag D 2211=并将A 分解成()D D A A +-= (2) 从而(1)可写成()b x A D Dx +-= 令11f x B x +=其中b D f ,A D I B 1111--=-=. (3) 以1B 为迭代矩阵的迭代法(公式)()()111f x B x k k +=+ (4)称为雅可比(Jacobi)迭代法,其分量形式为⎩⎨⎧[],...,,k ,n ,...,i x a ba xnij j )k (j j i iii)k (i21021111==∑-=≠=+ (5)其中()()()()()Tn x ,...x ,x x 002010=为初始向量. （2）Gauss-Seidel 迭代由雅可比迭代公式可知,在迭代的每一步计算过程中是用()k x 的全部分量来计算()1+k x的所有分量,显然在计算第i 个分量()1+k i x 时,已经计算出的最新分量()()1111+-+k i k x ,...,x 没有被利用。

把矩阵A 分解成U L D A --= (6)其中()nna ,...,a ,a diag D 2211=,U ,L --分别为A 的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成()b Ux x L D +=- 即22f x B x +=其中()()b L D f ,U L D B 1212---=-= (7)以2B 为迭代矩阵构成的迭代法(公式)()()221f x B x k k +=+ (8)称为高斯—塞德尔迭代法,用分量表示的形式为⎩⎨⎧[],...,,k ,n ,,i x a x a b a xi j n i j )k (j ij )k (j ij i ii)k (i21021111111==∑∑--=-=+=++2.2程序代码2.2.1高斯列主元的代码function Gauss(A,b) %A 为系数矩阵，b 为右端项矩阵 [m,n]=size(A); n=length(b); for k=1:n-1[pt,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); %开始消元 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:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);endfor k=1:nfprintf('x[%d]=%f\n',k,x(k));end2.2.2 LU分解法的程序function LU(A,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(1)=b(1);for k=2:ny(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:1x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);endfor k=1:nfprintf('x[%d]=%f\n',k,x(k));end2.2.3 Jacobi迭代法的程序function Jacobi(A,b,eps) %A为系数矩阵，b为后端项矩阵，epe为精度[m,n]=size(A);D=diag(diag(A)); %求矩阵DL=tril(A)-D; %求矩阵LU=triu(A)-D; %求矩阵Utemp=1;x=zeros(m,1);k=0;while abs(max(x)-temp)>epstemp=max(abs(x));k=k+1; %记录循环次数x=-inv(D)*(L+U)*x+inv(D)*b; %雅克比迭代公式endfor k=1:nfprintf('x[%d]=%f\n',k,x(k));end2.2.4Gauss-Seidel迭代程序function Gauss_Seidel(A,b,eps) %A为系数矩阵，b为后端项矩阵，epe为精度[m,n]=size(A);D=diag(diag(A)); %求矩阵DL=D-tril(A); %求矩阵LU=D-triu(A); %求矩阵Utemp=1;x=zeros(m,1);k=0;while abs(max(x)-temp)>epstemp=max(abs(x));k=k+1; %记录循环次数x=inv(D-L)*U*x+inv(D-L)*b; %Gauss_Seidel的迭代公式endfor k=1:nfprintf('x[%d]=%f\n',k,x(k));end三、实验过程中需要记录的实验数据表格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)x[1]=-1.333333x[2]=2.333333x[3]=-0.333333x[4]=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.0000 0 0 0-0.5000 1.0000 0 00 0.5000 1.0000 0-0.1250 0.2500 -0.0714 1.0000U =48.0000 -24.0000 0 -12.00000 12.0000 12.0000 6.00000 0 14.0000 -1.00000 0 0 12.9286x[1]=0.521179x[2]=1.005525x[3]=-0.375691x[4]=-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)x[1]=-1.467396x[2]=-2.358678x[3]=0.657604x[4]=2.8423973.4第三题用Gauss_Seidel迭代的数据>> A=[10 -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)x[1]=-1.467357x[2]=-2.358740x[3]=0.657597x[4]=2.842405四、实验中存在的问题及解决方案4.1存在的问题（1）第一题中在matlab中输入>> Gauss(A,b)（数据省略）得到m =4 n =4 ??? Undefined function or variable "Ab".Error in ==> Gauss at 8[ap,p]=max(abs(Ab(k:n,k)));没有得到想要的结果。