选主元doolittle分解法求解n元线性方程组 MATLAB实现

算法参考数值分析第四版 颜庆津 P23 （选主元的doolittle分解法算法）
运行结果截图：输入矩阵A，b。输出QA=LU的分解矩阵和Ax=b的解。


A=[1 8 2 3;-6 -3 8
1;2 4 4 2;10 5 -5 6]
b=[12;40;-50;80]





doolittle.m文件函数内容：
function [L,U,x] = doolittle(A,b)%用选主元的doolittle分解法求解线性方程组
z=size(A);
n=z(1);%b(1),b(2)分别是A的行和列.这里只处理n*n的非奇异矩阵
%错误检查
if z(1)~=z(2)%非方阵错误
error('MATLAB:Crout:Input Matrix should be a Square matrix. See Crout.');
end
if n~=rank(A)%非满秩矩阵错误
error('MATLAB:Crout:Input Matrix should be FULL RANK. See Crout.');
end

%初始化
M=1:n;%它的第k个元素M(k)记录第k个主元素所在的行号
S=1:n;%存储第一步中的中间量
L=zeros(n,n);
U=zeros(n,n);
x=1:n;
%U中的主对角线元素均为1
for i=1:n
L(i,i)=1;
end
%第一步 做分解QA=LU，Q是置换矩阵，即实现选主元素
for k=1:n
%第一小步，计算中间量S(i)
for i=k:n
sum=0;
for t=1:k-1
sum=sum+L(i,t)*U(t,k);
end
S(i)=A(i,k)-sum;
end
%第二小步。选行号ik
smax=S(k);%临时存储最大的S数组元素
ik=k;%临时存储最大元素的下标
for i=k+1:n
if S(i)>smax;
ik=i;
end
end
M(k)=ik;
%第三小步
if ik~=k;
for t=1:k-1
temp=L(k,t);L(k,t)=L(ik,t);L(ik,t)=temp;
end
for t=k:n
temp=A(k,t);A(k,t)=A(ik,t);A(ik,t)=temp;
end
temp=S(k);S(k)=S(ik);S(ik)=temp;
end
%第四小步
U(k,k)=S(k);
for j=k+1:n
sum=0;
for t=1:k-1
sum=sum+L(k,t)*U(t,j);
end
U(k,j)=A(k,j)-sum;
end
for i=k+1:n
L(i,k)=S(i)/U(k,k);
end
end
%第二步 求Qb
for k=1:n-1
t=M(k);
temp=b(k);b(k)=b(t);b(t)=temp;
end
%第三步 求解Ly=Qb和Ux=y
y=1:n;
y(1)=b(1);
for i=2:n
sum=0;
for t=1:i-1
sum=sum+L(i,t)*y(t);
end
y(i)=b(i)-sum;
end
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
sum=0;
for t=i+1:n
sum=sum+U(i,t)*x(t);
end
x(i)=(y(i)-sum)/U(i,i);
end

end %end of function