%就是下面的这个程序,请解释一下每一步的意思
function hl=zhjLU(A)
[n n] =size(A); RA=rank(A);
if RA~=n
disp('请注意:因为A的n阶行列式hl等于零,所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);
return
end
if RA==n
for p=1:n
h(p)=det(A(1:p, 1:p));
end
hl=h(1:n);
for i=1:n
if h(1,i)==0
disp('请注意:因为A的r阶主子式等于零,所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:'), hl;RA
return
end
end
if h(1,i)~=0
disp('请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:')
for j=1:n
U(1,j)=A(1,j);
end
for k=2:n
for i=2:n
for j=2:n
L(1,1)=1;L(i,i)=1;
if i>j
L(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);
else
U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);
end
end
end
end
hl;RA,U,L
end
end
你的代码我帮你解释了,果然是个复杂的活,如下
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循环,主要是判断各阶主子式是否为0
if 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对角线上的元素为1
if 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,下三角函数L
end %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 %并且退出程序
else
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存储.
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对角线上的元素为1
if 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,下三角函数L
end %RA~=n结束标志