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
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??? endend
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=yx=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