数值计算方法LU分解法实验

LU分解求解线性方程组【实验目的】1、 实现LU算法，并利用该算法求解线性方程组【实验内容】1.编写LU分解法程序LU_Decom.m,调用该程序，求解方程组2x+2x+3x+4x=12x+4x+9x+16x=1 ,,,a<! 1 2 3 4 的根。4x+8x+24x+63x=36x1+116x2+51x3+100x：=-29【算法流程说明】一、LU分解，得到L矩阵与U矩阵的值计算U11，U12，..q.，…U的值；计算L,L，…L，…L的值；1121i1 n1计算U的第2行和L的第2列，U的第3行和L的第三列，U的第r行和L的第r列，……U的第n行和L的第n列3.1计算U「『U,++「•••，U门，…U『的值3.1.13.2计算Lr3.1.13.2计算Lr+1,r3.2.1L=irUr,j=^r,j-,-1'r,Hjk=1,L，…,L，…L的值r+2,r i,r n,rair--/",'k=1JUrr二、 回代过程求y值求炉勺值:七=。(1)求y,y...,y，…y的值2 3kn2.1 y(k)=b(k)-如y(j)*L(k,j)j=1三、 回代过程求x值求x的值求x1，x2,…,x^，…尤］的值2.1 x(k)=y(k)-—x(j)*U(k,j)j=k+1【C语言算法】//LU分解，得到L矩阵与U矩阵的值For(r=2;r<=n;r++)(〃计算U的第r行各元素的值For(j=r;j<=n;j++)(For(k=1;k<=r-1;k++)(a(r,j)=a(r,j)-L(r,k)*U(k,j);}U(r,j)=a(r,j);}〃计算L的第r列各元素的值For(i=r+1;r<=n;r++)(For(k=1;k<=r-1;k++)(a(i,r)=a(i,r)-L(i,k)*U(k,r);}L(i,r)=a(i,r)/U(r,r);}}//回代求y值y(1)=b(1);For(k=2;k<=n;k++)(For(j=1;j<=k-1;j++)(b(k)=b(k)-y(j)*L(k,j);}y(k)=b(k);}//回代求x值x(n)=y(n)/U(n,n);For(k=n-1;k>=1;k--)(For(j=n;j>=k+1;j++)(b(k)=b(k)-y(j)*L(k,j);}x(k)=b(k)/U(k,k);【实验要求】请同学们参考上次实验的gauss_original.m文件，复习Matlab程序的格式和语法，对照给大家的C语言算法写出LU分解方法解多元一次线性方程组的程序，并调用该程序计算出方程组的解。本次实验要求写实现代码。实验报告请写明【实验目的】、【实验内容】、【源代码及注释】和【实验心得】function[L,U,x,y]=LU_Decom(A,b)%求矩阵A的LU分解，其中%A为要分解的矩阵；L为单位下三角矩阵；U为上三角矩阵；flag为指标变量，=1时表示计算成功，=0表示计算失败[n,m]=size(A);nb=length(b);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;endifn~=nberror('ThecolumnsofAmustequalthelengthofb!');endL=eye(n);U=zeros(n);y=zeros(n,1);x=zeros(n,1);%LU分解forj=1:nU(1,j)=A(1,j);L(j,1)=A(j,1)/U(1,1);endforr=2:nforj=r:nfork=1:r-1A(r,j)=A(r,j)-L(r,k)*U(k,j);endU(r,j)=A(r,j);endfori=r+1:nfork=1:r-1A(i,r)=A(i,r)-L(i,k)*U(k,r);endL(i,r)=A(i,r)/U(r,r);endend%由下三角矩阵L得到变量y的解y(1)=b(1);fori=2:nforj=1:i-1b(i)=b(i)-L(i,j)*y(j);endy(i)=b(i)；end%由上三角矩阵U得到变量x的解x