西京学院数学软件实验任务书西京学院数学软件实验任务书课程名称数学软件实验班级数 0901学号0912020114姓名王斌实验课题线性方程组直接三角分解法（Doolittle 分解，Grout 分解） ，平方根法（Cholesky 分解，LDLT 分解）实验目的熟悉线性方程组直接三角分解法（Doolittle 分解，Grout 分解） ，平方根法（Cholesky 分解，LDLT 分解）实验要求运用 Matlab/C/C+/Java/Maple/Mathematica 等其中一种语言完成实验内容线性方程组直接三角分解法（Doolittle 分解，Grout 分解）线性方程组平方根法（Cholesky 分解，LDLT 分解）成绩教师一线性方程组的直接三角分解法（Doolittle 分解）设为非奇异矩阵，且有分解式=LU，其中 L 为单位下三角矩阵，U 为上三角矩阵，则称此分解为Doolittle 分解。此分解的目的是将 Ax=b 分解为两步，首先由 Ly=b 解出 y，再由 Ux=y 解出 x。以下为实现此分解的程序：functionx,l,u=Doolittle(A,b)clcclear allformat shortn=input(请输入矩阵的阶数: ); A=zeros(n,n);for i=1:nfor j=1:nA(i,j)=input(请输入矩阵中的元素： );endendAif det(A)=0for i=1:nb(i)=input(请输入b中的元素： );endb=bu=zeros(n,n);l=eye(n,n);u(1,:)=A(1,:);for i=2:nfor m=1l(i,m)=A(i,m)/u(1,1);for k=2:nfor j=k:nu(k,j)=A(k,j)-sum(l(k,1:k-1)*u(1:k-1,j);for k=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);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);endendend此程序运行时，首先输入一个系数矩阵 A，和列向量b 如果系数矩阵 A 是非奇异矩阵，则矩阵最后可求得 A分解后的矩阵 L 和矩阵 U，以及列向量 x。以下是运行结果：二Cholesky 分解法。设 A 是 n（n=2）阶实对称矩阵，L 是非奇异的下三角矩阵，则称为矩阵 A 的 Cholesky 分解。其计TALL算步骤为将 A 分解后，首先由 Ly=b,求得 y，再由x=y,求TL得解向量 x。以下为实现此分解的 Matlab 程序:function x=pingfg(A,b)clcclear allformat shortn=input(请输入矩阵的阶数:); for i=1:nfor j=1:nA(i,j)=input(请输入矩阵中的元素);endendAfor i=1:nb(i)=input(请输入b中的元素：); endb=bn,n=size(A);L=zeros(n,n);L(1,1)=sqrt(A(1,1);for k=2:nL(k,1)=A(k,1)/L(1,1);endfor k=2:n-1L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).2);for i=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).2);%y=zeros(n,1);for k=1:nj=1:k-1;y(k)=(b(k)-L(k,j)*y(j)/L(k,k);endx=zeros(n,1);U=L;for k=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。以下是运行结果：三改进平方根分解法(LDLT 分解)因为用 Cholesky 分解后求 x 时需要开方，于是可以改进平法根分解法，使其分解时不需要开方。即将 A 分解为从而由 Ly=b 解得 y,再由=y 解得列向量TLDLTDL xx。从而解的此线性方程组的解。以下是实现此过程的Matlab 程序：function improvecholesky clcclear allformat shortn=input(请输入矩阵的阶数: ); for i=1:nfor j=1:nA(i,j)=input(请输入矩阵中的元素: );endendAfor i=1:nb(i)=input(请输入b中的元素: ); endb=bL=zeros(n,n); D=diag(n,0); S=L*D;for i=1:n L(i,i)=1;endfor i=1:n for j=1:n if (eig(A)=2）且 A 为实对称正定矩阵，和列向量 b 矩阵最后可求得 A 分解后的非奇异的下三角矩阵 L 和对角矩阵 D，以及列向量 x。以下是运行结果：