function [L,U,Y,X]=qw2014210705(A,b)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %切记!先将该文件重命名为doolittle.m, 然后存入到 bin文件中。。。。。。。。。 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%A为输入方阵
%然后在写入函数[L,U]=doolittle(A)即可获得B矩阵的Doolittle分解。 N=size(A);%获取A的行数和列数 if(N(1,1)~=N(1,2))%判断是否为方阵
disp('您输入的矩阵不是方阵,请重新输入'); return; end
n=N(1,1);%获取A的维数 Y=zeros(n,1); X=zeros(n,1); A1=A;
YY=0;%判断是否顺序主子式不等于0,如果存在任何一个顺序主子式为0,则该值变为1. for i=1:n
if(det(A1(1:i,1:i))==0) YY=1; end end if(YY==1)
disp('您输入的矩阵不能进行Doolittle分解') return; end L=eye(n); U=zeros(n);
for i=1:n%求解L的第一列,U的第一行 U(1,i)=A1(1,i); if(i>=2)
L(i,1)=A1(i,1)/U(1,1); end end
for k=2:n%以一行一列的方式,计算U和L阵中的其他元素 for j=k:n s=0; for r=1:k-1
s=s+L(k,r)*U(r,j); end
U(k,j)=A1(k,j)-s; end
for m=k+1:n s=0; for r=1:k-1
s=s+L(m,r)*U(r,k); end
L(m,k)=(A1(m,k)-s)/U(k,k); end end %求解过程 for i=1:n s=0; for j=1:n
s=s+Y(j,1)*L(i,j); end
Y(i,1)=(b(i,1)-s)/L(i,i); end nn=n; for i=1:n s=0; for j=1:n
s=s+X(j,1)*U(nn,j); end
X(nn,1)=(Y(nn,1)-s)/U(nn,nn); nn=nn-1; end
相关推荐: