1、克劳特(Crout)(LU)分解法求解线性方程组的matlab实现

1、克劳特（Crout）（LU)分解法求解线性方程组
function [x，L,U］=Crout（A，b)
%Crout分解法求解线性方程组
％系数矩阵:A
N=size（A);
n=N（1)；
L=zeros(n，n）; ％下三角矩阵
U=eye(n,n）; %上三角矩阵
L（1:n，1)=A(1：n，1）; %L的第一列
U（1,1:n)=A（1，1:n）/L(1，1); ％U的第一行？
for k=2:n
for i=k:n
L（i，k)=A(i,k）-L（i,1：（k—1)）＊U(1:(k—1),k）；
％L的第k列
end
for j=（k+1)：n
U(k,j）=（A(k，j)-L（k,1：（k—1）)＊U（1：(k—1),j))/（L（k,k));
%U的第k行
end
end
%y=inv(L）＊b；
％x=inv（U)*y;
y=SolveDownTriangle（L，b）;
x=SolveUpTriangle(U，y）；%求解线性方程组的解x
％x=U\(L\b);
function x=SolveUpTriangle（A,b）
%求解上三角矩阵Ax=b的解
N=size(A)；
n=N（1）；
for i=n：—1：1
if(i<n）
s=A(i,(i+1）:n）＊x(（i+1）：n，1）；
else
s=0；
end
x（i,1)=(b（i）—s）/A（i，i）;
end
function x=SolveDownTriangle（A,b）
%求解下三角矩阵Ax=b的解
N=size(A）；
n=N(1);
for i=1：n
if(i〉n）
s=A(i,1:（i—1))＊x（1:(i—1），1);
else
s=0;
end
x(i,1）=（b（i)-s)/A(i,i）;
end
％求解线性方程组的解
clc
clear
A=[12 -3 3;—16 3 —1;1 1 1];
b=[15;—13；6]；
％x=A\b
［x，L,U］=Crout(A,b)
解:
x =
1
2
3
L =
12.0000 0 0
—16.0000 —1.0000 0
1.0000 1。

2500 4。

5000
U =
1.0000 —0。

2500 0。

2500
0 1.0000 -3.0000
0 0 1。

0000
列主元LU分解
function ［L，U,x］=lux(A,b)
％LU 分解法解线性方程组(列主元LU分解）
［n，n]=size(A)；
p=eye（n);%p记录了选择主元时候所进行的行变换
for k=1:n-1
［r,m］=max(abs(A(k:n,k))); ％选列主元
m=m+k-1;
if（A（m，k）～=0)
if（m～=k)
A([k m］，：)=A（[m k],:）；
p([k m］）=p(［m k])；
end
for i=k+1:n
A(i,k）=A（i，k)/A(k,k);
j=k+1:n;
A（i,j）=A(i,j）—A（i，k)＊A(k，j); end
end
end
L=tril(A,—1)+eye(n，n）;
U=triu（A);
%解下三角矩阵 Ly=b
newb=p＊b；
y=zeros（n，1);
for k=1:n
j=1：k-1;
y(k）=(newb(k)—L（k,j)*y（j））/L(k,k)；
end
%解上三角方程组 Ux=y
x=zeros（n,1）；
for k=n：—1：1
j=k+1：n;
x(k)=(y(k)-U(k，j）＊x(j））/U(k，k);
end。