平方根法和改进的平方根法解线性方程组（Matlab程序）

%定理2.2.3：对阵正定矩阵的楚列斯基（Cholesky）分解

%设A为n阶对阵正定矩阵，则存在一个可逆的下三角矩阵G，使得

%A=GG’，当限定G的对角元为正时，这种分解是唯一的

%--------A=GG’的分解算法-------

%参考教材：《数值分析》李乃成，梅立泉，科学出版社

clear;clc;

A=[9,18,9,-27;

18,45,0,-45;

9,0,126,9;

-27,-45,9,135];

b=[1 2 16 8]';

n=length(b);%方程个数n

G=zeros(n,n);

G(1,1)=sqrt(A(1,1));

G(2:n,1)=A(2:n,1)/G(1,1);

for j=2:n-1

G(j,j)=sqrt(A(j,j)-sum(G(j,1:j-1).^2));

for i=j+1:n

G(i,j)=(A(i,j)-sum(G(i,1:j-1).*G(j,1:j-1)))/G(j,j);

end

end

G(n,n)=sqrt(A(n,n)-sum((G(n,1:n-1)).^2));

G

%--------用A=GG’分解求解方程组Ax=b（平方根法）----

%Gy=b，G'x=y;

x=zeros(n,1);%未知向量

y=zeros(n,1);%中间向量

y(1)=b(1)/G(1,1);

for i=2:n

y(i)=(b(i)-sum(G(i,1:i-1)'.*y(1:i-1)))/G(i,i);

end

y

%---------由G'x=y求出x----------

%方法类似于U由x=y解出x

G=G';

x(n)=y(n)/G(n,n);

for i=n-1:-1:1

x(i)=(y(i)-sum(G(i,i+1:n)'.*x(i+1)))/G(i,i);

end

x