对称矩阵的平方根法

数值计算方法

09医软（1）班

本组实验同学：刘康康秦强梅世友马蕾乔琼

任务分配：刘康康：分析平方根法解对称矩阵的思想，并写出推导过程

秦强：利用c++实现算法，并写出程序

梅世友：调试并运行程序

马蕾：根据解方程组的方法写出流程图

乔琼：负责word的排版（解题流程，程序截图等）

一．实验名称

平方根法解对称方程组

二．实验目的

理解平方根法解对称方程组的基本思想，原理以及公式的推导过程，和用c++实现算法

三．理论分析推导和程序实现

（一）、理论分析推导

设A为对称矩阵，且A的所有顺序主子式均不为0，由定理4（P178）可知，A可以唯一地分解为A=LU，其中L为单位下三角矩阵，U为上三角矩阵，为了利用A的对称性，再将U进行分解，即

⎥⎥⎥⎥⎦⎤⎢⎢⎢

⎢⎣⎡nn n n u u u u u u 22211211=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡nn

u u u 22

11⎥⎥

⎥

⎥⎥⎥⎦

⎤⎢

⎢

⎢⎢⎢⎢⎣

⎡

---1///1//1)1,1(),1(22

222

2311

111

12

n n n n n n u u u u u u u u u

u

或写成

U=D u ~

其中D 为对角矩阵，U1为单位上三角矩阵。于是 A=LU=LD u ~

（2.12） 又

A=A '=u '~

D L '(A 'L '分别表示A.L 的转置矩阵) 由分解的唯一性，得u '~

=L 。由（2.12）知 A=LD L ' 且这种分解是唯一的。

如果A 是对称正定的，利用代数知识可以证明ii u （i=1，2，…，n ）皆大于零。事实上，将A=LU 按分块形式写出，有

k k k U L A =（k=1,2,…，n)，

其中

k

A =⎥⎥⎥⎦

⎤

⎢⎢⎢⎣⎡kk k k k a a a a a a 2111211，

k

L =

⎥⎥⎥⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎢⎢⎣⎡-11111,2121k k k k l l l l ，

k U =⎥

⎥⎥

⎥⎦

⎤⎢⎢⎢⎢⎣⎡kk

k

k u u u u u u 22211211

由于A 的各阶主子式皆大于零，所以 de k A =det k L det k U =

∏=k

i ii

u 1

>0.

上式对k=1,2，…，n 成立，故有

ii u >0（i=1,2，…，n ）。

因为A 为对称正定，所以ii u >0（i=1,2，…，n ），于是可将D

进一步分解为

D=⎥

⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎣⎡nn

u u u 2211 =⎥⎥⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎢⎣

⎡nn u u u

22

11⎥⎥⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎢⎣

⎡nn u u u

22

11

=D 2

1

D 2

1

于是

A=LD L '=L D 2

1D 2

1L ' =（L D

2

1）'⎪⎪⎭

⎫ ⎝⎛21LD =L L '。 其中L1=L D 21

为下三角矩阵。当限定L1的对角元素为正时，这种分解也是唯一的。 由

A=L11L '=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n l l l n 2n 1n 222111l l l ⎥

⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎣⎡n

l l l n 2n 221n 2111l l l ， 用比较法可导出计算L1元素的计算公式，对i=1,2，…，n ， ii l =

∑-=-1

1

2

i k ik ii l a ，

ji l =.,,2,1),(1

1

1

n i i j l l a l i k ki jk ji ii ++=-∑-=

这一方法称为平方根法，又称Cholesky 分解。 （二）程序的实现

1.输入方程组的系数和b 的值并判断矩阵是否为对称矩阵 #include "stdio.h" #include "math.h" #define N 20 int main() { int i,j,k; int size; int n=0;

float a[N][N],l[N][N]; float b[N],x[N],y[N]; //float S,H,M,L;

printf("\t\t\t\t*Cholesky 分解法解方程*\n");

printf("请输入方阵A的n："); scanf("%d",&size);

printf("\n");

printf("请输入方程组的系数:\n"); //数组a[][]的输入

for(i=0;i

{

for(j=0;j

{

scanf("%f",&a[i][j]);

}

}

printf("\n请输入方程组的y:\n"); //数组b[]的输入

for(i=0;i

scanf("%f",&b[i]);

printf("\n方阵A[][]为：\n");

//2.数组A[][]的输出

for(i=0;i

{

for(j=0;j

{