cholesky.asv

为配套王能超版本数值分析简明教程编写,方便学生教师演示!

function [L,d,x]=cholesky(A,b)
[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=diag(ones(1,n));
        for i=1:n
            for j=1:i-1
                s=0;
                for l=1:j-1
                   s=s+d(l)*L(i,l)*L(j,l);
                end
                L(i,j)=(A(i,j)-s)/d(j);
            end
            s=0;
            for l=1:i-1
                s=s+d(l)*L(i,l)^2;
            end
            d(i)=A(i,i)-s;
        end
    end
    d=diag(d);
else
    disp('此矩阵不为方阵');
end
y(1)=b(1);
for i=2:n
    y(i)=b(i)-L(i,1:i-1)*y(1:i-1)';
end
B=inv(d)*y';
C=L';
x(n)=B(n);
for i=n-1:-1:1
    x(i)=B(i)-C(i,i+1:n)*x(i+1:n)';
end