MATLAB-平方根法和改进平方根法求解线性方程组例题与程序

（2）设对称正定阵系数阵线方程组

42224102214302004021812243344121411681344115310163340x10x620256x32033x423103x5914x622142x71521945x80

x(1,1,0,2,1,1,0,2)T

二、数学原理 1、平方根法

解n阶线性方程组Ax=b的choleskly方法也叫做平方根法，这里对系数矩阵A是有要求的，需要A是对称正定矩阵，根据数值分析的相关理论，如果A对称正定，那么系数

TTLLA=L•L矩阵就可以被分解为的形式，其中L是下三角矩阵，将其代入Ax=b中，可得：x=b

进行如下分解：

yLTxLyb

那么就可先计算y,再计算x，由于L是下三角矩阵，是L上三角矩阵，这样的计算比

T

直接使用A计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，

那么对于对称正定矩阵A进行Cholesky分解，我再描述一下过程吧：

如果你对原理很清楚那么这一段可以直接跳过的。

设A=L•L，即

Ta11a12a21a22an1an2...a1nl11l11l21...ln1l...a2nll...l212222n2.....................annln1ln2...lnnlnn

其中

aijaji,i,j1,2,...,n

第1步，由矩阵乘法，

2a11l11,ai1li1l11故求得

l11a11,li1ai1,i2,3,...na11

一般的，设矩阵L的前k-1列元素已经求出

第k步，由矩阵乘法得

akkll, 　aiklimlkmliklkk2km2kkm1m1k1k1

于是

k12lkkakklkmm1(k2,3,...,n)k1l1(all)　　,　ik1,k2,...n ikimkmiklm1kk

2、改进平方根法

在平方根的基础上，为了避免开方运算,所以用ALDL计算；其中，

d1D...d1dn...d1dn...11D2D2dn；

T得

1l121Aln1ln2

d1d211l121dnl1nl2n1

按行计算的L元素及D对元素公式

对于i1,2,,n

j1tijaijtikljk(j1,2…，i1)k1.

lijtij/dj(j1,2,…，i-1).

diaiitiklikk1i1

计算出TLD的第i行元素tij(j1,2,…，i-1)后，存放在A的第i行相置，然后再计算L的第i行元素，存放在A的第i行.D的对角元素存放在A的相应位置.

TTTLLLDLLDL对称正定矩阵A按分解和按分解计算量差不多，但分解不需要开放计算。求解Lyb, DLTxy的计算公式分别如下公式。

y1b1,i1 yibilikyhk1 i2,....,n

xnyn/dn,nxiyi/dilkixkki1 in1,....,2,1

三、程序设计 1、平方根法

function [x]=pfpf(A,b)

%楚列斯基分解 求解正定矩阵的线性代数方程 A=LL’ 先求LY=b 再用L’可以求出解X

[n,n]=size(A);

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

for k=2:n

L(k,1)=A(k,1)/L(1,1);

end

X=Y 即

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矩阵

y=zeros(n,1);%先生成方程组的因变量的位置，给定y的初始值

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矩阵

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

>> A=[4,2,-4,0,2,4,0,0

2,2,-1,-2,1,3,2,0

-4,-1,14,1,-8,-3,5,6

0,-2,1,6,-1,-4,-3,3

2,1,-8,-1,22,4,-10,-3

4,3,-3,-4,4,11,1,-4

0,2,5,-3,-10,1,14,2

0,0,6,3,-3,-4,2,19];

>> b=[0;-6;20;23;9;-22;-15;45];

>> x=pfpf(A,b)

x =

121.1481

-140.1127

29.7515

-60.1528

10.9120

-26.7963

5.4259

-2.0185

2、改进平方根法

function [x]=improvecholesky(A,b,n) %用改进平方根法求解Ax=b

L=zeros(n,n); %L为n*n矩阵

D=diag(n,0); S=L*D;

for i=1:n L(i,i)=1;

end

for i=1:n for j=1:n if (eig(A)<=0)|(A(i,j)~=A(j,i)) disp('wrong');break;end

end

end

%D为n*n的主对角矩阵

%L的主对角元素均为1

%验证A是否为对称正定矩阵

%A的特征值小于0或A非对称时，输出wrong

D(1,1)=A(1,1); %将A分解使得A=LDLT

for i=2:n

for j=1:i-1

S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');

L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);

end

D(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');

end

y=zeros(n,1); % x，y为n*1阶矩阵

x=zeros(n,1);

for i=1:n

y(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i);

%通过 LDy=b解得y的值

end

for i=n:-1:1

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

%通过LTx=y解得x的值

end

>> A=[4,2,-4,0,2,4,0,0

2,2,-1,-2,1,3,2,0

-4,-1,14,1,-8,-3,5,6

0,-2,1,6,-1,-4,-3,3

2,1,-8,-1,22,4,-10,-3

4,3,-3,-4,4,11,1,-4

0,2,5,-3,-10,1,14,2

0,0,6,3,-3,-4,2,19];

>> b=[0;-6;20;23;9;-22;-15;45];

>> n=8;

>> x=improvecholesky(A,b,n)

x =

121.1481

-140.1127

29.7515

-60.1528

10.9120

-26.7963

5.4259

-2.0185

四、结果分析和讨论

平方根法和改进平方根法求解线性方程组的解为x=（121.1481，-140.1127，

29.7515，-60.1528，10.9120，-26.7963，5.4259，-2.0185）T。与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差，所以在求解的过程中还是应该选择主元以此消除误差，提高精度。

五、完成题目的体会与收获

对称正定矩阵的平方根法及改进平方根法是目前解决这类问题的最有效的方法之一，合理利用的话，能够产生很好的求解效果。改进平方根法较平方根法，因为不用进行开方运算，所以具有一定的求解优势。通过求解此题，学会了平方根法和改进平方根法matlab编程，使我受益匪浅。