对称矩阵的平方根法

数值计算方法

09医软（1）班

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

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

秦强：利用c++实现算法，并写出程序 梅世友：调试并运行程序

马蕾：根据解方程组的方法写出流程图 乔琼：负责word的排版（解题流程，程序截图等）

一．实验名称

平方根法解对称方程组 二．实验目的

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

三．理论分析推导和程序实现 （一）、理论分析推导

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

1u/uu11u12u1nu111u22u2nu22= unnunn1211u23/u22u1n/u11u2n/u22 u(n1,n)/u(n1,n1)1或写成

~u U=D

其中D为对角矩阵，U1为单位上三角矩阵。于是 A=LU=LD（2.12） 又

A=A=uDL(AL分别表示A.L的转置矩阵) 由分解的唯一性，得u=L。由（2.12）知 A=LDL 且这种分解是唯一的。

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

~u

~~

AkLkUk（k=1,2,…，n)，

其中

a11a12a1kak1ak2akk1l2111lk1lk2lk,k11Ak=，

Lk

=

，

u11u12u1ku22u2kUk=

ukk由于A的各阶主子式皆大于零，所以 deAk=detLkdetUk=

u>0.

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

因为A为对称正定，所以进一步分解为

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

uii>0（i=1,2，…，n），于是可将D

u11u22D= unnu11 =u22unnu11u22 unn =DD

1212于是

A=LDL=LDDL =（LD121211212）LD=LL。 其中L1=LD2为下三角矩阵。当限定L1的对角元素为正时，这种分解也是唯一的。 由

l11l21l22 A=L1L1=ln1ln2lnnl11l21ln1l22ln2， lnn用比较法可导出计算L1元素的计算公式，对i=1,2，…，n， lii=

aiilik2k1i1，

i11 lji=(ajiljklki),ji1,i2,,n.

liik1这一方法称为平方根法，又称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(\"\\\\*Cholesky分解法解方程*\\n\");

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

printf(\"请输入方程组的系数:\\n\"); //数组a[][]的输入 for(i=0;ifor(j=0;jscanf(\"%f\ } }

printf(\"\\n请输入方程组的y:\\n\"); //数组b[]的输入 for(i=0;iprintf(\"\\n方阵A[][]为：\\n\"); //2.数组A[][]的输出 for(i=0;ifor(j=0;jprintf(\"%f \ }

printf(\"\\n\"); }

printf(\"\\n方程组y为:\\n\"); //数组b[]的输出 for(i=0;i//3.判断是否可能用平方根法 for(i=0;ifor(j=0;jif(a[i][j]==a[j][i]) n=n+1; } }

if(n!=size*size) {

printf(\"次线性方程组的系数不是对称矩阵，不可以用平方根法，请按键退出\\n\");

return 0; } else {

printf(\"次线性方程组的系数是对称矩阵，可以用平方根法，算法如下：\\n\\n\");

//4.平方根法的计算过程 计算Lii

l[0][0]=sqrt(a[0][0]); for(i=1;il[i][0]=a[i][0]/l[0][0];

printf(\"l[%d][0]=%f\\n\ }

for(i=1;il[i][i]=a[i][i];

printf(\"l[%d][%d]=%f\\n\

for(k=0;kl[i][i]=l[i][i]-l[i][k]*l[i][k];

}

l[i][i]=sqrt(l[i][i]);

printf(\"l[%d][%d]=%f\\n\\n\

输出

Lji

for(j=i+1;jl[j][i]=a[j][i];

printf(\"l[%d][%d]=%f\\n\ for(k=0;kl[j][i]=l[j][i]-l[j][k]*l[i][k]; }

l[j][i]=l[j][i]/l[i][i];

printf(\"l[%d][%d]=%f\\n\\n\ } }

//5.回代过程 求出

yi

y[0]=b[0]/l[0][0]; for(i=1;iy[i]=b[i];

for(k=0;ky[i]=y[i]-l[i][k]*y[k]; }

y[i]=y[i]/l[i][i];

printf(\"y[%d]=%f\\n\ }

printf(\"\\n\"); for(i=0;ifor(j=i+1;jl[i][j]=0; } }

//6.L矩阵的输出

printf(\"L矩阵为：\\n\\n\"); for(i=0;ifor(j=0;jprintf(\"%f \ }

printf(\"\\n\"); }

printf(\"\\n\"); 求出xi

x[size-1]=y[size-1]/l[size-1][size-1]; for(i=size-2;i>=0;i--) {

x[i]=y[i];

for(k=i+1;kx[i]=x[i]-l[k][i]*x[k]; }

x[i]=x[i]/l[i][i];

printf(\"x[%d]=%f\\n\ }

printf(\"\\n\");

printf(\"计算的x值为：\\n\"); 输出xi

//x[]的输出过程 for(i=0;i{

printf(\"x[%d]=%f \ }

printf(\"\\n\"); } }

四、算法的源程序截图

接上图

（程序参见课后习题3 P211）

五．实验流程图

开始 输入数据n及方程组的系数和b【】的值 N 矩阵是否对称 Y 输出矩阵A并根据A求出矩阵L1的值 i1lii=√aii-k1lik, lji=（aii-ljklki）/lii 2i1k1输出矩阵L1，由L1y=b，L1x=y解出xi，yi T输出 结束