pfg.m

来自「为配套王能超版本数值分析简明教程编写,方便学生教师演示!」· M 代码 · 共 55 行

function [L,x]=pfg(A,b)%进行平方根分解法的条件为A为正定矩阵
[n,m]=size(A);
k=1;
if n==m
    for i=1:n;
        for j=1:n;
            if A(i,j)~=A(j,i)
                k=0;
                disp('此矩阵不是对称矩阵');
                break;
            end
        end
        if k==0
            break;
        end
    end
    if k~=0
        for i=1:n;%判断是否为正定矩阵
            B=A(1:i,1:i);
            if det(B)<=0
                disp('此矩阵的顺序主子式不全大于零,不为正定据阵');
                k=0;
            end 
        end
    end
    if k~=0
        L(1,1)=A(1,1)^(1/2);
        for i=1:n
            for j=1:i-1
                s=0;
                for l=1:j-1
                   s=s+L(i,l)*L(j,l);
                end
                L(i,j)=(A(i,j)-s)/L(j,j); 
                s=0;
                for l=1:i-1
                    s=s+L(i,l)^2;
                end
                L(i,i)=(A(i,i)-s)^(1/2);
            end
        end
        y(1)=b(1)/L(1,1);
        for i=2:n
            y(i)=(b(i)-L(i,1:i-1)*y(1:i-1)')/L(i,i);
        end
        M=L';
        x(n)=y(n)/M(n,n);
        for i=n-1:-1:1
            x(i)=(y(i)-M(i,i+1:n)*x(i+1:n)')/M(i,i);
        end
    end
else
    disp('此矩阵不为方阵');
end