Doolittle分解法的构造过程

Doolittle分解法的构造过程把矩阵A写成下列两个矩阵相乘的形式：A=LU，其中L为下三角矩阵，U为上三角矩阵。这样我们可以把线性方程组Ax=b写成Ax=（LU)x=L(Ux)=b令Ux=y,则原线性方程组Ax=b化为两个简单的三角方程组：Ux=y和Ly=b。于是可首先求解Ly=b得到向量y，然后求解Ux=y,从而求解线性方程组Ax=b的目的。设矩阵A的Doolittle分解为：为求出矩阵L和U，根据矩阵乘法规则，有用公式写法有：ijjjjiiiiijiauuullllULj000,,0,1,,,,213211LUuuuuuulllaaaaaaaaaAnnnnnnnnnnnn222112112121212222111211111nkjikkjikkjikijulula1,min1niulululanjuulululainkkkikkikijjnkkkjkkjkj,,2,,2,11111111111111111111iiikkijknkikkijkkijkjiijikkjiknkikkjikkjkijululululauulululajiji111111111,有时当于是可得Doolittle分解公式：u1j=a1jj=1,2,…,nli1=ai1/u11i=2,3,…,n2）Doolittle分解法算法1.输入变量个数n、系数矩阵A、常数项b2如果a11=0,则输出“LU分解失败”提示并终止，否则aj1aj1/a11（j=2,…,n）ijuulalijulauiiikkijkjiikkjikijijji/)(1111注：因为L和U中的三角零元素都不必存储，L的对角元素也因为都是1没有必要再记录在程序中。这样为节省计算机存储单元，用原来的n阶系数矩阵A把L和U贮存起来。算法中系数矩阵A的下三角存储L各元素而上三角存储U的元素（包括对角元）3）Dolittle分解法程序Clear[a,b,x];n=Input[“线性方程组阶数n=”];a=Input[系数矩阵A=]；b=Input[常数项b=]；t=1;eps1=0.00000001;If[Abs[a[[1,1]]]eps1,t=0;Break[],Do[a[[j,1]]=a[[j,1]]/a[[1,1]],{j,2,n}]];If[t==1,Do[a[[i,i]]=a[[i,i]]-Sum[a[[i,k]]*a[[k,i]],{k,1,i-1}];If[Abs[a[[i,i]]]eps1,t=0;Print[不能分解];Break[],nkjkkjkjkkikiikijkjijiikkjikijijiiikkiikiiiiaxabxnnkForaaaaaaaaaniijForLUaaaaaniFor1111111/)(1.61,...,1,.6/)(2.51.5,,2,1.50.41.3,,3,2.3分解，停止输出不能如果Do[a[[i,j]]=a[[i,j]]-Sum[a[[i,k]]*a[[k,j]],{k,1,i-1}];a[[j,i]]=a[[j,i]]-Sum[a[[j,k]]*a[[k,i]],{k,1,i-1}];a[[j,i]]=a[[j,i]]/a[[i,i]],{j,i+1,n}]],{i,2,n}]];If[t==1,Do[b[[i]]=b[[i]]-Sum[a[[i,k]]*b[[k]],{k,1,i-1}],{i,1,n}];Print[Ly=b的解y=,b];Do[b[[i]]=(b[[i]]-Sum[a[[i,k]]*b[[k]],{k,i+1,n}])/a[[i,i]],{i,n,1,-1}]]u=IdentityMatrix[n];Do[Do[u[[i,j]]=a[[i,j]],{i,1,j}],{j,1,n}]l=a-u+IdentityMatrix[n];Print[LU分解紧凑格式=,MatrixForm[a]]Print[矩阵L=,MatrixForm[l]]Print[矩阵U=,MatrixForm[u]]Print[解向量x=,b]说明本程序用于求线性方程组Ax=b的解。程序执行后，先通过键盘输入线性方程组阶数n、系数矩阵A、常数项b。程序即可给出Doolittle分解的紧凑格式、分解矩阵L和U、Ly=b的解y，和Ax=b的解x。如果出现主对角线元素aii=0给出Doolittle分解失效提示。