平方根法

function [x]=pingfg(A,b)
%乔累斯基分解
[n,n]=size(A);
L=zeros(n,n);%实际上不用为 L 申请空间,使用 A 即可
L(1,1)=sqrt(A(1,1));
for k=2:n
    L(k,1)=A(k,1)/L(1,1);
end
for k=2:n-1
    L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).^2));
    for i=k+1:n
         L(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k);
    end
end
L(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).^2));
%解下三角方程组 Ly=b
y=zeros(n,1);
for k=1:n
   j=1:k-1;
   y(k)=(b(k)-L(k,j)*y(j))/L(k,k);
end
%解上三角方程组  L'x=y
x=zeros(n,1);
U=L';
for k=n:-1:1
    j=k+1:n;
    x(k)=(y(k)-U(k,j)*x(j))/U(k,k);
end

posted on 2008-12-13 21:24  秋波渡  阅读(6081)  评论(0编辑  收藏  举报