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

M A T L A B-平方根法和改进平方根法求解线性方程组例题与程序
（2）设对称正定阵系数阵线方程组
123456784240240
0022121
32064114183562002161433232181
22
4
103943344111422025310114215006
3
3
4
21945x x x x x x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢
⎥⎢⎥
⎢⎥⎢⎥⎢⎥
----⎢⎥⎢
⎥⎢⎥
----⎢⎥⎢⎥⎢⎥
=
⎢⎥⎢⎥⎢⎥
----⎢⎥⎢
⎥⎢⎥
----⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢---⎢⎥⎢
⎥⎢--⎢⎥⎢⎢⎥⎣
⎦⎣⎦
⎣⎦⎥⎥
⎥⎥ (1,1,0,2,1,1,0,2)T x *=--
二、数学原理 1、平方根法
解n 阶线性方程组Ax=b 的choleskly 方法也叫做平方根法，这里对系数矩
阵A 是有要求的，需要A 是对称正定矩阵，根据数值分析的相关理论，如果A 对称正定，那么系数矩阵就可以被分解为的T A=L L •形式，其中L 是下三角矩阵，将其代入Ax=b 中，可得：T LL x=b 进行如下分解：
T L x
L b
y y ⎧=⎨
=⎩ 那么就可先计算y,再计算x ，由于L 是下三角矩阵，是T L 上三角矩阵，这样的计算比直接使用A 计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，
那么对于对称正定矩阵A 进行Cholesky 分解，我再描述一下过程吧： 如果你对原理很清楚那么这一段可以直接跳过的。

设T A=L L •，即
1112111112112122221222221212....................................n n n n n n nn n n nn nn a a a l l l l a
a a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
其中,,1,2,...,ij ji a a i j n ==
第1步，由矩阵乘法，2
1111
1111,i i a l a l l ==g 故求得
1
11111
,2,3,...i i a l l i n a ==
= 一般的，设矩阵L 的前k-1列元素已经求出 第k 步，由矩阵乘法得
1
1
221
1
k k kk km
kk
ik im km ik kk
m m a l l a l l l l --===+=+∑∑,
于是
1
1(2,3,...,n)1(),1,2,...kk k ik ik im km m kk l k l a l l i k k n l -=⎧=⎪⎪=⎨⎪=-=++⎪⎩
∑ 2、改进平方根法
在平方根的基础上，为了避免开方运算,所以用T
LDL A =
计算；其中，
1
1122.......
.
.
n d D D D d ⎤⎤⎡⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎢⎢⎥⎣
⎦⎣
⎣
；
得
1121212212111111n n n n n d l l l d l A l l d ⎡⎤⎡⎤
⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦
⎣⎦L L M MO O O M L
按行计算的L 元素及D 对元素公式 对于n i ,,2,1Λ=
1
1
(1,21)
j ij ij ik jk k t a t l j i -==-=-∑…，.
/(1,2,)ij ij j l t d j ==…，i-1.
1
1i i ii ik ik
k d a t l -==-∑
计算出LD T =的第i 行元素(1,2,i-1)ij t j =…，后，存放在A 的第i 行相置，然后再计算L 的第i 行元素，存放在A 的第i 行.D 的对角元素存放在A 的相应位置.
对称正定矩阵A 按T LDL 分解和按T LL 分解计算量差不多，但T LDL 分解不需要开放计算。

求解b Ly =,
y x DL T
=的计算公式分别如下公式。

1111
,i i i ik h
k y b y b l y ===⎧⎪⎨
=-⎪⎩∑ ()2,....,i n =
1/,
/n n n n i i i ki k
k i x y d x y d l x ===⎧⎪⎨
=-⎪⎩∑ ()1,....,2,1i n =-
三、程序设计 1、平方根法
function [x]=pfpf(A,b)
%楚列斯基分解 求解正定矩阵的线性代数方程 A=LL ’ 先求LY=b 再用L ’X=Y 即可以求出解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
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); %D为n*n的主对角矩阵
S=L*D;
for i=1:n %L的主对角元素均为1
L(i,i)=1;
end
for i=1:n
for j=1:n %验证A是否为对称正定矩阵
if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时，输出wrong
disp('wrong');break;end
end
end
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编程，使我受益匪浅。