Jacobi迭代法分量形式matlab,LU分解法列主元高斯法Jacobi迭代法Gauss-Seidel法的原理及Matlab程序...

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

南昌航空大学数学与信息科学学院实验报告 一、实验目的及题目 1.1 实验目的： (1)学会用高斯列主元消去法，LU分解法，Jacobi迭代法和Gauss-Seidel迭代法解线性方程组。 (2)学会用Matlab编写各种方法求解线性方程组的程序。 1.2 实验题目： 1. 用列主元消去法解方程组： 2. 用LU分解法解方程组其中 ， 3. 分别用Jacobi迭代法和Gauss-Seidel迭代法求解方程组： 二、实验原理、程序框图、程序代码等 2.1实验原理 2.1.1高斯列主元消去法的原理 Gauss消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式： 这个过程就是消元，然后再回代就好了。具体过程如下： 对于，若依次计算 然后将其回代得到： 以上是高斯消去。 但是高斯消去法在消元的过程中有可能会出现的情况，这时消元就无法进行了，即使主元数但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。 2.1.2直接三角分解法(LU分解)的原理 先将矩阵A直接分解为则求解方程组的问题就等价于求解两个三角形方程组。 直接利用矩阵乘法，得到矩阵的三角分解计算公式为： 由上面的式子得到矩阵A的LU分解后，求解Ux=y的计算公式为 以上为LU分解法。 2.1.3Jacobi迭代法和Gauss-Seidel迭代法的原理 (1)Jcaobi迭代 设线性方程组 (1) 的系数矩阵A可逆且主对角元素均不为零,令 并将A分解成 (2) 从而(1)可写成 令 其中. (3) 以为迭代矩阵的迭代法(公式) (4) 称为雅可比(Jacobi)迭代法,其分量形式为 (5) 其中为初始向量. (2)Gauss-Seidel迭代 由雅可比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i个分量时,已经计算出的最新分量没有被利用。 把矩阵A分解成 (6) 其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 即 其中 (7) 以为迭代矩阵构成的迭代法(公式) (8) 称为高斯—塞德尔迭代法,用分量表示的形式为 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>k t=A(k,:);A(k,:)=A(p,:);A(p,:)=t; %交换行使之变到主元位置上 t=b(k);b(k)=b(p);b(p)=t; end m=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~=0 Ab=[A,b]; end end x=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); end for k=1:n fprintf( x[%d]=%f\n ,k,x(k)); end 2.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:n U(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); end L %输出L矩阵 U %输出U矩阵 y=zeros(n,1); %开始解方程组Ux=y y(1)=b(1); for k=2:n y(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:1 x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k); end for k=1:n fprintf( x[%d]=%f\n ,k,x(k)); end 2.2.3 Jacobi迭代法的程序 function Jacobi(A,b,eps) %A为系数矩阵，b为后端项矩阵，epe为精度 [m,n]=size(A); D=diag(diag(A)); %求矩阵D L=tril(A)-D; %求矩阵L U=triu(A)-D; %求矩阵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; %雅克比迭代公式 end for k=1:n fprintf( x[%d]=%f\n ,k,x(k)); end 2.2.4 Gauss-Seidel迭代程序 functi