6月7日更新:
LUDecomp[m_]:=Block[{l,u,n},
n=Dimensions[m][[1]];l=IdentityMatrix[n];u=ConstantArray[0,{n,n}];u[[1]]=m[[1]];
Do[l[[r+1;;n,r]]=(m[[r+1;;n,r]]-l[[r+1;;n,;;r-1]].u[[;;r-1,r]])/u[[r,r]];
u[[r+1,r+1;;n]]=m[[r+1,r+1;;n]]-l[[r+1,;;r]].u[[;;r,r+1;;n]];,{r,1,n-1}];{l,u}]
提高了计算效率几百倍,原本for循环的效率太低了
附老外的算法(真的短....):
LU[mat_]:=Module[{n,L,U},
n=Length@mat;L=IdentityMatrix[n];U=mat;
Do[L[[k;;n,k]]=U[[k;;n,k]]/U[[k,k]];
U[[(k+1);;n,k;;n]]-=L[[(k+1);;n,{k}]].U[[{k},k;;n]],{k,1,n-1}];{L,U}]
—————————原内容—————————
虽然有自带的函数,但是那个函数会进行行变换,所以我又自己写了一个
LUDecomp[m_]:=Block[{lu,n},
n=Dimensions[m][[1]];
lu=ConstantArray[0,{n,n}];
For[r=1,r<=n,r++,
For[i=r,i<=n,i++,lu[[r,i]]=m[[r,i]]-Sum[lu[[r,k]]*lu[[k,i]],{k,1,r-1}]];
For[i=r+1,i<=n,i++,lu[[i,r]]=(m[[i,r]]-Sum[lu[[i,k]]*lu[[k,r]],{k,1,r-1}])/lu[[r,r]]];
];
lu
]
返回结果中L和U是在一起的,L是结果的下三角部分,U是上三角部分
LowerTriangularize[%,-1]+IdentityMatrix@Dimensions@%
UpperTriangularize[%%]