数值线性代数Cholesky分解法解线性方程组MATLAB实现

算法思想

如何利用电子计算机来快速、有效地求解线性方程组是数值线性代数研究的核心问题，而且也是目前人在继续研究的重大课题之一。
Cholesky分解法又叫做平方根法，是求解对称正定线性方程组最常用的方法之一。A是一个对称正定的矩阵，则存在一个对角元均为正数的下三角阵L,使得A=LL’，称为Cholesky分解，而后我们可以由下面三步求解：
(1)计算A的Cholesky分解：A=LL’；
(2)求解Ly=b得y；
(3)求解L’x=y得x。

矩阵创建

function [A,b]=creatMaxtrix(n)
    T=zeros(n-1,n-1);
    for i=1:n-1
        T(i,i)=2;
    end
    for i=1:n-2
        T(i,i+1)=-1;
        T(i+1,i)=-1;
    end
    I=eye(n-1,n-1);
    k=(n-1)^2;
    A=zeros(k,k);
    b=ones(k,1);
    j=1;
    for i=1:n-1
        A(j:j+n-2,j:j+n-2)=T+2*I;
        j=j+n-1;
    end
    j=1;
    for i=1:n-2
        A(j:j+n-2,j+n-1:j+2*n-3)=-I;
        A(j+n-1:j+2*n-3,j:j+n-2)=-I;
        j=j+n-1;
    end
end

该矩阵是一个(n-1)² 阶的矩阵，其具有这样几个特点，
(1)A是块三对角阵，共有五条对角线上有非零元素；
(2)A是不可约对角占优的；
(3)A是对称正定的，而且是稀疏的。
矩阵的形式如下

可通过spy来观察矩阵的分布，一个n=11时的矩阵分布情况如下

Cholesky分解

function []=Cholesky(n)
[A,b]=creatMaxtrix(n);
[n,~]=size(b);
x=A\b;
%Cholesky分解
for k=1:n
    A(k,k)=sqrt(A(k,k));
    A(k+1:n,k)=A(k+1:n,k)/A(k,k);
    for j=k+1:n
        A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);
    end
end
%计算时将y,x存在b的存储单元里
%解下三角形方程组：前代法
L=tril(A,0);
for j=1:n-1
    b(j)=b(j)/L(j,j);
    b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);
end
b(n)=b(n)/L(n,n);
%解上三角形方程组：回代法
U=L';
for j=n:-1:2
    b(j)=b(j)/U(j,j);
    b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);
end
b(1)=b(1)/U(1,1);
%用二范数来衡量误差
norm(x-b)

总结

运行结果如下

>> Cholesky(11)
ans =
   1.9304e-14

与Gauss消去法对比，函数中只有对于矩阵的分解部分变动较大。运行结果表明算法的准确性，但是像前一篇文章中所说，由于实际问题的复杂，往往需要考虑更多。