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[%%]