平方根法求解对称正定线性方程组的Matlab实现

数值分析实验二：用平方根法和改进的平方根法求线性方程组

1.平方根法

function [] = Square_Root(A,b)
tic
%使用CHOLESKY分解的平方根法
%#ok<*NOPRT> 
n=length(b);%矩阵维度
A(1,1)=sqrt(A(1,1));%计算l(1,1)
for i=2:n
    A(i,1)=A(i,1)/A(1,1);%L的第一列元素
end
for k=2:n
    Sigma1=0;
    for j=1:k-1
        Sigma1=Sigma1+A(k,j)*A(k,j);%lkk减去的求和式子
    end
    A(k,k)=sqrt(A(k,k)-Sigma1);%计算lkk
    if k+1<=n %判断lkk是否为最后一个对角元素
        for i=k+1:n
            Sigma2=0;
            for r=1:k-1
                Sigma2=Sigma2+A(i,r)*A(k,r);%lik减去的求和式子
            end
            A(i,k)=(A(i,k)-Sigma2)/A(k,k);%计算lik
        end
    end
end
L=tril(A);
U=triu(A');
y=L\b;x=U\y;
for i=1:n %将解横向展示
    value=x(i,1);
    fprintf('x%d=%.4f; ',i,value)
end
fprintf('\n')%将光标移到下一行
toc

2.改进的平方根法 

function [] = SR_fixed(A,b)
tic
%改进的平方根算法  Crout分解+LD'L分解
%#ok<*NOPRT>
n=length(b);
A(2,1)=A(2,1)/A(1,1); %计算l(2,1)
for i=2:n
    Sigma1=0;
    for r=1:i-1
        Sigma1=Sigma1+A(i,r)^2*A(r,r);
    end
    A(i,i)=A(i,i)-Sigma1;
    if i~=n %当求到d(n,n)时停止
        for j=1:i
            Sigma2=0;
            for r=1:j-1
                if r~=1 %避免遍历时出现r=1:0的情况
                    Sigma2=Sigma2+A(i,r)*A(r,r)*A(j,r);
                end
            end
            A(i+1,j)=(A(i+1,j)-Sigma2)/A(j,j);
        end
    end
end
D=zeros(n,n);L=tril(A);L_T=triu(A');
for k=1:n %提取对角元素and把L的对角元素改为1
    L(k,k)=1;
    L_T(k,k)=1;
    D(k,k)=A(k,k);
end
L,D,L_T %展示分解结果
y=L\b;x=D*L_T\y;
for i=1:n %将解横向展示
    value=x(i,1);
    fprintf('x%d=%.4f; ',i,value)
end
fprintf('\n') %将光标移到下一行
toc

注：tic/toc为程序执行计时语句 

将题设矩阵输入

>> A=[6,1,0;1,4,1;0,1,14]

A =

     6     1     0
     1     4     1
     0     1     14

>> b=[6,24,322]'

b =

     6
    24
   322

 调用两个函数计算可得结果

>> Square_Root(A,b)

L =

    2.4495         0         0
    0.4082    1.9579         0
         0    0.5108    3.7066


U =

    2.4495    0.4082         0
         0    1.9579    0.5108
         0         0    3.7066

x1=1.0000; x2=-0.0000; x3=23.0000; 
历时 0.022259 秒。
>> SR_fixed(A,b)

L =

    1.0000         0         0
    0.1667    1.0000         0
         0    0.2609    1.0000


D =

    6.0000         0         0
         0    3.8333         0
         0         0   13.7391


L_T =

    1.0000    0.1667         0
         0    1.0000    0.2609
         0         0    1.0000

x1=1.0000; x2=0.0000; x3=23.0000; 
历时 0.003965 秒。

另附运行截图