求直接三角分解法的matlab程序，每一步都需要解释，谢谢！

Question

%就是下面的这个程序，请解释一下每一步的意思function hl=zhjLU(A)[n n] =size(A); RA=rank(A); if RA~=ndisp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);returnendif RA==n for p=1:nh(p)=det(A(1:p, 1:p));endhl=h(1:n);for i=1:nif h(1,i)==0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'), hl;RAreturnendend if h(1,i)~=0 disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendhl;RA,U,Lendend

Answer 1

你的代码我帮你解释了，果然是个复杂的活，如下function hl=zhjLU(A) %函数名为zhiLU，输入矩阵A,来进行LU分解，返回值hl为A的各阶主子式的行列式[n n] =size(A);%返回矩阵A的维数
RA=rank(A); %返回矩阵A的秩if RA~=n %判断矩阵A的秩RA是否不等于A的维数n，当不等于n时，即小于n时执行if中的语句disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);%输出disp中的字符串，RA的值,hl的值(此时它的值为矩阵A的行列式)return %并且退出程序end % if end 语句块，end标志if语句结束if RA==n %判别A的秩是否等于n，当等于n时，即满秩时执行下面的语句for p=1:n%从1到n循环，主要是想获得1至n阶主子式的行列式h(p)=det(A(1:p, 1:p));%A(1:p, 1:p)为A的前p行前p列，即A的p阶主子式，从而h(p)为A的p阶主子式的值end %for语句结束标志hl=h(1:n);%把A的各阶主子式的行列式值赋给hl存储.for i=1:n %从1到n循环，主要是判断各阶主子式是否为0if h(1,i)==0 %判段第i阶主子式的行列式是否为0，若为0输出下面的语句disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'), hl;RA%输出disp中的字符串，各阶主子式的行列式的值，以及矩阵A的秩return %当有一个为0时退出程序end %if语句结束标志end %for语句结束标志if h(1,i)~=0 %如果执行到这一句，说明上边的for循环都执行了且没有return出去，则此时i的值为n，判断第n阶行列式是否为0，即还是判断A的行列式是否不为0，若不为0则输出下面语句disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')%输出disp中的字符串for j=1:n%对1到n循环U(1,j)=A(1,j);%把A的第一行中的各个值赋给U中第一行的各个变量end%for循环结束标志for k=2:n%从第2列到第n列进行循环，即依次给LU的第k列赋值for i=2:n%从2到n循环 ,与j的循环是相互联系的，即当i>j时对下三角矩阵L赋值，当j>=i时对上三角U赋值for j=2:n%从2到n循环 ,与i的循环是相互联系的，即当i>j时对下三角矩阵L赋值，当j>=i时对上三角U赋值L(1,1)=1;L(i,i)=1; %始终保持L对角线上的元素为1if i>j %当i>j时，此时就是要对L进行赋值了，因为当i<j时,L(i,j)=0,不用管L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);%根据LU分解的关系知,L(i,1)处的值=A(i,1)/U(1,1)，从而对L(i,1)进行赋值，即对L的第一列赋值L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k);%还是根据LU分解的关系L(i,k)处的值等于(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k)，对L(i,k)赋值，即对L的第k列赋值else %否则,即i<j时，此时就是要对U进行赋值了，因为当i>j时,U(i,j)=0,不用管U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);%根据LU分解的关系知U(k,j)的值等于A(k,j)-L(k,1:k-1)*U(1:k-1,j)进而对U的第k列赋值end %if结束标志end %for j=2:n结束标志end %for i=2:n结束标志end %for k=2:n结束标志hl;RA,U,L %输出矩阵的秩RA，上三角函数U,下三角函数Lend %if h(1,i)~=0结束标志end %RA==n结束标志 %%%%%%%%%%%%%%%%%%%%%%%%以上代码中有很多多余的，当判断A的秩为n之后其他的主子式的行列式都不为0了，判断是多余的，故我进行了改进,如下%%%%%%%%%%%%%%%%%%%%%%% function hl=LUfenJie(A) %函数名为zhiLU，输入矩阵A,来进行LU分解，返回值hl为A的各阶主子式的行列式[n n] =size(A);%返回矩阵A的维数
RA=rank(A); %返回矩阵A的秩if RA~=n %判断矩阵A的秩RA是否不等于A的维数n，当不等于n时，即小于n时执行if中的语句disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);%输出disp中的字符串，RA的值,hl的值(此时它的值为矩阵A的行列式)return %并且退出程序elsefor p=1:n%从1到n循环，主要是想获得1至n阶主子式的行列式h(p)=det(A(1:p, 1:p));%A(1:p, 1:p)为A的前p行前p列，即A的p阶主子式，从而h(p)为A的p阶主子式的值end %for语句结束标志hl=h(1:n);%把A的各阶主子式的行列式值赋给hl存储.
disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')%输出disp中的字符串for j=1:n%对1到n循环U(1,j)=A(1,j);%把A的第一行中的各个值赋给U中第一行的各个变量end%for循环结束标志for k=2:n%从第2列到第n列进行循环，即依次给LU的第k列赋值for i=2:n%从2到n循环 ,与j的循环是相互联系的，即当i>j时对下三角矩阵L赋值，当j>=i时对上三角U赋值for j=2:n%从2到n循环 ,与i的循环是相互联系的，即当i>j时对下三角矩阵L赋值，当j>=i时对上三角U赋值L(1,1)=1;L(i,i)=1; %始终保持L对角线上的元素为1if i>j %当i>j时，此时就是要对L进行赋值了，因为当i<j时,L(i,j)=0,不用管L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);%根据LU分解的关系知,L(i,1)处的值=A(i,1)/U(1,1)，从而对L(i,1)进行赋值，即对L的第一列赋值L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k);%还是根据LU分解的关系L(i,k)处的值等于(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k)，对L(i,k)赋值，即对L的第k列赋值else %否则,即i<j时，此时就是要对U进行赋值了，因为当i>j时,U(i,j)=0,不用管U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);%根据LU分解的关系知U(k,j)的值等于A(k,j)-L(k,1:k-1)*U(1:k-1,j)进而对U的第k列赋值end %if结束标志end %for j=2:n结束标志end %for i=2:n结束标志end %for k=2:n结束标志hl;RA,U,L %输出矩阵的秩RA，上三角函数U,下三角函数Lend %RA~=n结束标志