matlab迭代法求某数平方根,MATLAB平方根法和改进平方根法求解线性方程组例题与程序要点.doc...

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

(2)设对称正定阵系数阵线方程组

数学原理

平方根法

解n阶线性方程组Ax=b的choleskly方法也叫做平方根法，这里对系数矩阵A是有要求的，需要A是对称正定矩阵，根据数值分析的相关理论，如果A对称正定，那么系数矩阵就可以被分解为的形式，其中L是下三角矩阵，将其代入Ax=b中，可得：

进行如下分解：

那么就可先计算y,再计算x，由于L是下三角矩阵，是上三角矩阵，这样的计算比直接使用A计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，

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

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

设，即

其中

第1步，由矩阵乘法，故求得

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

第k步，由矩阵乘法得

于是

改进平方根法

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

得

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

对于

.

.

计算出的第行元素后，存放在的第行相置，然后再计算的第行元素，存放在的第行.的对角元素存放在的相应位置.

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

程序设计

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