改进的平方根法及其程序实现

目录摘要 (2)0 引言 (3)1 预备知识 (3)1.1 TLL分解理论 (3)1.2 Cholesky分解法 (4)1.3 算法描述 (5)2改进的平方根法 (6)3TLDL分解算法描述 (7)4 应用举例 (8)5 程序实现 (10)5.1 程序码源 (10)5.2 实例计算 (12)6 结束语 (13)参考文献 (14)致谢 (15)改进的平方根法及其程序实现摘要: 针对对称正定方程组的解法,本文先对Cholesky分解法进行了分析研究,在此基础上给出了改进的平方根法(即TLDL分解法),此方法有效地避免了原平方根法开方运算所带来的误差和不便,并通过算法描述、实例计算,用C程序实现了TLDL分解,进一步提高了矩阵运算的速度和精度.关键词: 对称正定矩阵, 平方根法, TLDL分解, 算法Improved Methods of Square Root and Realizationof Its ProgramAbstract: Aiming at studying solutions of symmetric positive definition matrix in linear equations. Initially, the text has conducted a series of analyses and researches towards decomposition proposed by Cholesky. Then based on theses researches and analyses, it offers the improved methods of square–root (also called TLDL decom- position), which effectively avoid some errors and inconvenience brought by the process of extracting root. At the same time, it achieves the TLDL decomposition through the means of algorithm description, example calculation as well as applicat- ion of C program, further enhancing the speed and accuracy in matrix operation.Key words: Symmetric positive definition matrix, method of square root,TLDL decomposition, algorithm0 引言很多工程中的科学计算,例如应用有限元法解结构力学问题时,最后往往归结为求解系数矩阵为对称正定方程组解的问题.由于对称正定矩阵各阶顺序主子式以及全部的特征值均大于零, 这种特征也使得其三角分解具有更为简单的形式, 不同的分解也导出了一些不同的解法. 平方根法(即Cholesky 分解法), 就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法, 其计算量和存储量约为普通消去法的一半, 且无需选主元就能求得较为精确的数值解, 但由于在平方根法中含有多次开方运算, 因此给计算带来了许多不便, 而在原平方根法的基础上, 给出改进的平方根法(即T LDL 分解法), 成功避免了开方运算带来的的麻烦, 因此在各种工程计算中应用更为广泛.1 预备知识1.1 T LL 分解理论定理 1.1（矩阵的LU 分解定理） 设A 为n 阶矩阵, 如果A 的顺序主子式0(1,2,,1)i D i n ≠=- , 则A 可分解为一个单位下三角矩阵L 和一个上三角矩阵U 的乘积, 且这种分解是唯一的.定理1.2（对称正定矩阵的T LL 分解定理） 设A 为n 阶对称正定矩阵, 则必存在非奇异下的三角矩阵L , 使T A LL = , (1.1) 并且当L 的对角线元素均为正数时, 这种分解是唯一的.证明 因为A 对称正定, 则它的各阶顺序主子式均大于零. 由LU 分解定理可知, 矩阵A 存在唯一的Doolitle 分解, 即 1A LU = 其中, 1L 为单位下三角矩阵, U 为上三角矩阵. 记 D =diag 11(,,)nn u u ,1P D U -=, 则P 为单位上三角矩阵, 且1A L DP =. 由A 的对称性有T A A =, 可得 11T T T T A P D L L DP A ===由于T P 为单位下三角矩阵, T D D =, 从而由分解的唯一性可知, 1T P L =, 于是 11T A L DL =显然有110,(1,2,)k k k k k kk A L U L U u u k n ==⋅=>= ,这样0kk u >,(1,2,)k n = , 故可令 12D =diag则111122221111()()T TT A L D D L L D L D LL ===其中,121L L D =为非奇异的下三角矩阵. 而知限定L 的对角元素为正数时,这种分解是唯一的, 定理得证.1.2 Cholesky 分解法设A 为对称正定矩阵, 由上面T LL 分解定理可知, 则存在一个实的非奇异下三角矩阵L , 使T A LL =, 当限定L 的对角线元素为正数时, 这时A 的分解是唯一的.设11212212n n nn l l l L l l l ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦ (1.2) 其中, 0(1,2,,)ii l i n >= , 0()ij l i j =<.由T A LL =比较A 和T LL 的对应元素, 可求得L 的元素ij l 如下:由 21111a l =,1111i i a l l =, 得11l =,1111/,1,2,,i i l a l i n == .假设L 的第1k -列元素已经求得, 下面求L 的第k 列元素,,,ik l i k n = . 注意到 11k ik ij kj ik kk j a l l l l -==+∑,得计算公式, 对1,2,,j n = , 有121/2111(),(1,2,,)()/,(,1,,1)k kk kk kj j k ik ik ij kj kkj l a l k n l a l l l k n n -=-=⎧=-=⎪⎪⎨⎪=-=-⎪⎩∑∑ (1.3) 由此可逐行求得L 的全部元素, 从而由Ly b =及T L x y =得到求解方程组Ax b =的公式111()/,(1,2,,)()/,(,1,,1)k k k kj j kk j nk k jk j kk j k y b l y l k n x y l x l k n n -==+⎧=-=⎪⎪⎨⎪=-=-⎪⎩∑∑ (1.4) 上述方法称为Cholesky 分解法. 由于计算的对角线元素要作n 次开平方运算, 故Cholesky 分解法又称平方根法.1.3 算法描述步1 输入对称正定矩阵A 和右端向量b ; 步2 Cholesky 分解:11111111,/,2,,,i i d t a l a d i n ==== 对2,,k n = 计算: 11,k k kk kj kj j d a t l -==-∑11,/,1,,;k ik ik ij kj ik ik k j t a t l l t d i k n -==-==+∑步3 用向前消去法解下三角方程组Ly b =: 11y b =,对2,,k n = 计算 11k k k kj j j y b l y -==-∑;步4 解对角形方程组Dz y =:对1,,k n = 计算 /k k k z y d =;步5 用回代法解上三角形方程组T L x z =:n n x z =,对1,,1k n =- 计算 1nk k jkj j k x z lx =+=-∑.不难验证Cholesky 分解法的乘除计算总量约为32/6()n O n +, 为一般矩阵LU 分解计算量的一半. 虽然如此, 但其增加的n 次开平方运算是非常不利的,下面引出改进的平方根法—T LDL 分解法.2 改进的平方根法为了避免开方运算, 对A 作LU 分解, 即A LU =, 则111212122212111n n n n nn u u u l u u A l l u ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦/ 提出U 矩阵的对角元素 / 121112122212111111n n n n nn u u u l u u l l u ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦由A 对称正定, 可得0ii u >, 令111222nn n u d u d D u d ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦ 可证1212111n n Tu u u L ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦即T A LDL =.1121212212111111n n n n n d l l l d l A l l d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ (1.5) L 是对角元素为1的单位下三角矩阵.对矩阵A 作LU 分解, 共计算2n 个矩阵元素; 对称矩阵的T LDL 分解, 只需计算(1)2n n -个元素, 减少了近一半的工作量. 借助LU 分解计算公式, 容易得到T LDL 分解计算公式.设A 有LU 分解形式111211~212221211()1n n T n n n d d l d l l d d l A L DL LU l l d ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦其中~ij i ij i ji u d l d l ==.在分解中可套用LU 分解公式, 只要计算下三角矩阵L 和D 的对角元素k d .计算中只需要保存()ij L l =的元素, T L 的i 行j 列的元素用L 的ji l 表示. 由于对称正定矩阵的各阶主子式均大于零, 直接调用LU 分解公式可完成T LDL 分解计算.对于1,,k n = ,有11~~11k k k kk kk kj jk kk kj j jk j j d u a l u a l d l --====-=-∑∑121k k kk j kj j d a d l -==-∑ (1.6)11~~11()/()/k k ik ik ij jk kk ik j ij kj k j j l a l u u a d l l d --===-=-∑∑11()/,1,,k ik ik j ij kj k j l a d l l d i k n -==-=+∑ . (1.7)3 T LDL 分解算法描述步1 输入方程组阶数n 、系数矩阵A 和常数项b . 步2 FOR :1k = TO n ;121:k k kk j kj j d a d l -==-∑;FOR :1i k =+ TO n ; 11:()/k ik ik j ij kj k j l a d l l d -==-∑;步3 解方程组的步骤从略.由()T L DL x b =, 解方程组Ax b =可分为三步完成: (1)解方程组Lz b =, 其中T z DL x =. 11,1,2,,i i i ij j j z b l z i n -==-=∑ (1.8)(2)解方程组Dy z =, 其中T y L x =. /,1,2,,i i i y z d i n == (1.9)(3)解方程组T L x y =. 1,,1,,1ni i jij j i x y lx i n n =+=-=-∑ . (1.10)可以看出, 改进后的T LDL 分解乘除运算量约为32/6()n O n +, 计算过程也无需开方运算, 使得其运算更加简单易行, 因此改进后的T LDL 分解法相对于原平方根法具有更好的实效性和可行性.4 应用举例例 试用T LDL 分解法求解方程组123164844543842210x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦ 解 由T A LDL =, 可得Ax b =的方程组T LDL x b =, 令T DL x y =, 则Ly b =.计算步骤:(1)将A 直接分解为T A LDL =, 求出,L D ; (2)求解方程组Ly b =; (3)求解方程组1T L x D y -=. 现有1.16d =, 1214d l =, 2114l =, 1318d l =, 3112l =2122.154d d l =-=, 21312.3214d l d l l =--, 3232l =-, 39d =即由Ly b =可得1231411341013122y y y ⎡⎤⎢⎥-⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥-⎣⎦解得1234418y y y -⎡⎤⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦ 由1T L x D y -=有21311212323132311648145411842211l l d A l d l l l d ⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦⎣⎦111413122L ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎣⎦1649D ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦1231111424311221x x x ⎡⎤⎡⎤⎢⎥-⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦解得1232494x x x ⎡⎤⎢⎥⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦-⎢⎥⎣⎦5 程序实现5.1 程序码源应用TLDL 分解法解线性方程组, 为了便于计算, 编写如下程序, 通过计算机来实现线性方程组的求解, 简洁方便. 具体程序如下: Purpose ：T LDL 分解法解方程组 # include <stdio.h> # include <math.h># define MAX_N 20 //(x_i,y_i)的最大维数//int main() { int n; int i,j,k;int mi;double mx,tmp;static double a[MAX_N][MAX_N],b[MAX_N],x[MAX_N],y[MAX_N],z[MAX_N];static double L[MAX_N][MAX_N],d[MAX_N];//输入Ax=b的维数//printf("\n Input n value(dim of Ax=b):");scanf("%d",&n);if(n>MAX_N){printf("The input n is larger then MAX_N,please redefine the MAX_N.\n");return 1;}if(n<=0){printf("Please input a number between 1 and %d.\n", MAX_N);return 1;}//输入Ax=b的A矩阵//printf("Now input the matrix a(i,j),i,j=0,…, %d:\n", n-1);for(i=0;i<n;i++)for(j=0;j<n;j++)scanf("%lf",&a[i][j]);//输入b矩阵//printf("Now input the matrix b(i),i=0,…, %d:\n", n-1);for(i=0;i<n;i++)scanf("%lf",&b[i]);for(i=0;i<n;i++)L[i][i]=1; //u矩阵对角元素为1//for(k=0;k<n;k++){d[k]=a[k][k]; //计算d_i//for(j=0;j<=k-1;j++)d[k]-=(L[k][j] *d[j]*L[k][j]);for(i=k+1;i<n;i++) //计算L矩阵// {L[i][k]=a[i][k];for(j=0;j<=k-1;j++)L[i][k]-=(L[i][j] *d[j]*L[k][j]);L[i][k]/=d[k];}}for(i=0;i<n;i++) //求解L′z=b的z// {z[i]=b[i];for(j=0;j<=i-1;j++)z[i]-=L[i][j]*z[j];}for(i=0;i<n;i++) //求解Dy=z的y//y[i]=z[i]/d[i];for(i=n-1;i>=0;i--) //求解Ly=x// {x[i]=y[i];for(j=i+1;j<n;j++)x[i]-=L[j][i]*x[j];}printf("Solve … x_i=\n"); //输出//for(i=0;i<n;i++)printf("%f\n",x[i]);return 0;}5.2 实例计算用T LDL 分解法求解方程组:123164844543842210x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦ 程序输入输出Input n value(dim of Ax=b): 3Now input the matrix a(i,j),i,j=0, (2)16 4 8 4 5 -4 8 -4 22Now input the matrix b(i),i=0, (2)-4 3 10Solve … x_i=2.0000004.000000-2.2500006 结束语从矩阵分解角度看, T LDL 分解法与消元法本质上没有多大的区别, 其基本思想还是通过等价变换将线性方程组化为结构简单、易于求解的形式. 在实际问题中经常遇到对称正定矩阵方程组的求解问题, 对于这种具有特殊性质的系数矩阵, 采用改进的平方根法是一种行之有效的计算方法, 它是在普通消去法的基础上建立起来的, 算法思想虽有些复杂, 但其计算量约为普通消去法的一半, 且有效地避免了开方运算所带来的不便, 因此改进后的T LDL 分解法相对于原平方根法具有更好的实效性和可行性. 并通过C 程序求解使得这类问题变得更加简单易行, 在大量工程计算中有着广泛而重要的应用.参考文献[1] 刘元亮. Matlab 在对称正定矩阵的改进平方根分解法中的应用[J]. 湖南: 怀化学院学报. 2004,(02).[2] 郭丽杰, 周硕, 秦万广. 对称矩阵的改进Cholesky 分解在特征值问题中的应用[J]. 吉林: 东北电力学院学报. 2003,(02).[3] 杜廷松,沈艳君,覃太贵.数值分析及实验[M].北京:科学出版社.2007.[4] 肖筱南,赵来军,党林立.数值计算方法与上机实习指导[M].北京:北京大学出版社. 2004.[5] 马昌凤,林伟川.现代数值计算方法[M].北京: 科学出版社.2008.[6] 张韵华,奚梅成,陈长松.数值计算方法与算法[M].北京:科学出版社.2000.[7] 李庆扬,王能超,易大义.数值分析(4版)[M].武汉:华中科技大学出版社.2006.[8]严克明,欧志英,刘树群.数值计算基础[M].甘肃:甘肃人民出版社.2006.[9] Ferenc Szidsrovzky, Sidney Yokowitz. Principles and procedures of numerical analysis[M].上海:上海科学技术文献出版社.1982.[10] Rainer Kress. Numerical analysis[M].北京:世界图书出版社.2003.致谢本论文经过将近半年的努力即将告以尾声, 从论文的选题到最后完稿, 整个过程都经过了认真地考虑、仔细地查阅和细心地修改. 在此, 我首先要感谢我的导师郭晓斌老师, 不管是我的论文选题还是论文的撰写, 以及资料的查阅等各方面, 他都给了我莫大的帮助与启发, 尤其是在论文的几次修改过程中, 郭老师以他广博的专业学识、严谨的治学精神和他耐心的指导态度, 才使我的论文能顺利完成. 谨再次向郭老师致以崇高的敬意和诚挚的感谢. 也祝愿郭老师身体健康、工作顺利, 在学术研究上取得更加辉煌的成就, 为更多的学子导航.同时, 我也要对给予本论文参考文献的所有学术专家和老师致以真挚的谢意, 是他们出版的书籍与发表的学术论文给了我很大的启示与指导, 才将论文完成.其次, 在即将毕业之际, 我要感谢数信学院所有的老师对我的细心教育与培养, 让我在四年的学习生涯中不仅学到了扎实的专业知识, 而且他们的言传身教使我受益匪浅, 他们严谨的治学态度和耐心教导学生的博爱精神也是我永远学习的榜样, 并将积极影响着我今后的学习和工作.最后, 我还要再次感谢我的母校——西北师范大学给了我发展的平台, 感谢各位老师四年来对我的耐心教导和栽培.张登科 2011年5月。