消元法实验报告13

西京学院数学软件实验任务书课程名称数学软件实验班级数0901学号0912020n4姓名王斌实验课题线性方程组直接三角分解法（ 分解，分解），平方根法（ 分解， 分解）实验目的熟悉线性方程组直接三角分解法（ 分解，分解），平方根法（ 分解， 分解）实验要求运用 等其中一种语言完成实验内容线性方程组直接三角分解法（ 分解，分解）线性方程组平方根法（ 分解， 分解）成绩教师一.线性方程组的直接三角分解法(Doolittle分解)设人为非奇异矩阵，且有分解式A=LU,其中L为单

位下三角矩阵，U为上三角矩阵，则称此分解为Doolittle

分解。此分解的目的是将Ax=b分解为两步，首先由Ly=b

解出y,再由Ux=y解出x。以下为实现此分解的程序：function[x,l,u]=Doolittle(A,b)clcclearallformatshortn=input('请输入矩阵的阶数：')；A=zeros(n,n);fori=1:nforj=1:nA(i,j)=input('请输入矩阵中的元素：')；endendAifdet(A)~=0fori=1:nb(i)=input('请输入b中的元素：')；endb=b'u=zeros(n,n);l=eye(n,n);u(1,:)=A(1,:);fori=2:nform=1l(i,m)=A(i,m)/u(1,1);fork=2:nforj=k:nu(k,j)=A(k,j)-sum(l(k,1:k-1)*u(1:k-1,j));fork=2:n-1l(k+1:n,k)=(A(k+1:n,k)-sum(l(k+1:n,1:k-1)*u(1:k-1,k)))/u(k,k);endendendendulendy=zeros(n,1);y(1)=b(1);fork=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);fork=n-1:-1:1x(k)=(y(k)-u(k,k+1:n)*x(k+1:n))/u(k,k);endendend此程序运行时，首先输入一个系数矩阵A,和列向量b

如果系数矩阵A是非奇异矩阵，则矩阵最后可求得A分

解后的矩阵L和矩阵U，以及列向量x。以下是运行结果:CominanclWindow请输入矩胜的阶数：3请输入矩除的元素：1请输入矩阵的元素：-1请输入矩除的元素：3请输入矩阵的元素：2请输入矩除的元素：-4请输入矩阵的元素：6请输入矩阵的元素：4请输入矩阵的元素：-9请输入矩阵的元素：2A=TOC\o"1-5"\h\z1 -1 32-4 64-9 2请输入b中的元素：1请输入b中的元素：41 -1 30-2 00 0 -10=0000 0 00000 1.0000 00000 2.5000 1.0000ans=-2.4000-1.00000.8000二. 分解法。设是（ ）2实对称矩阵，是非奇异的下三角矩阵，则称A二LLt为矩阵的 分解。其计算步骤为将分解后，首先由 求得，再由l求得解向量X以下为实现此分解的 程序function[x]=pingfg(A,b)clcclearallformatshortn=input('请输入矩阵的阶数：')；fori=1:nforj=1:nA(i,j)=input('请输入矩阵中的元素');endendAfori=1:nb(i)=input('请输入b中的元素：')；endb=b'[n,n]=size(A);L=zeros(n,n);L(1,1)=sqrt(A(1,1));fork=2:nL(k,1)=A(k,1)/L(1,1);endfork=2:n-1L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1)J2));fori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k)endendL(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1)J2));%y=zeros(n,1);fork=1:nj=1:k-1;y(k)=(b(k)-L(k,j)*y(j))/L(k,k);endx=zeros(n,1);U=L';fork=n:-1:1j=k+1:n;x(k)=(y(k)-U(k,j)*x(j))/U(k,k);end此程序运行时，首先输入一个系数矩阵A(A的阶数＞=2)且A为实对称正定矩阵，和列向量b矩阵最后可求得A分解后的非奇异的下三角矩阵L,以及列向量x。以下是运行结果：L=TOC\o"1-5"\h\z4 0 01 2 02-3 0ans=-111»I三.改进平方根分解法 分解因为用分解后求时需要开方，于是可以改因为用分解后求时需要开方，于是可以改进平法根分解法，使其分解时不需要开方。即将 分解为LDLt从而由 解得 再由DLtx 解得列向量x从而解的此线性方程组的解。以下是实现此过程的程序：functionimprovecholeskyclcclearallformatshortn=input('请输入矩阵的阶数：')；fori=1:nforj=1:nA(i,j)=input('请输入矩阵中的元素：');endendAfori=1:nb(i)=input('请输入b中的元素：')；endb=b'L=zeros(n,n);D=diag(n,0);S=L*D;fori=1:nL(i,i)=1;endfori=1:nforj=1:nif(eig(A)<=0)|(A(i,j)~=A(j,i))disp('wrong');break;endendendD(1,1)=A(1,1);fori=2:nforj=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endLDy=zeros(n,1);x=zeros(n,1);fori=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i);endfori=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));endx此程序运行时，首先输入一个系数矩阵A(A的阶数＞=2)且A为实对称正定矩阵，和列向量b矩阵最后可求得A分解后的非奇异的下三角矩阵L和对角矩阵D,以及列向量x。以下是运行结果：

CommandWindow请输入矩阵的阶数：3请输入矩阵的元素：2清输入矩阵的元素：-1清输入矩阵的元素：1清输入矩阵的元素：-1清输入矩阵的元素：2请输入矩阵的元素：1请输入矩阵的元素