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

解线性n 阶方程组直接法—Cholesky 方法解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步，由矩阵乘法，211111111,i i a l a l l == 故求得111111,2,3,...i i a l l i n a === 一般的，设矩阵L 的前k-1列元素已经求出第k 步，由矩阵乘法得112211k k kk km kkik im km ik kk m m a l l a l l l l --===+=+∑∑, 于是11(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 -=⎧=⎪⎪=⎨⎪=-=++⎪⎩∑ 注意到21kkk km m a l ==∑，于是有21max km kk ii i nl a a ≤≤≤≤ 这充分说明分解过程中元素2km l 的平方不会超过系数矩阵A 的最大对角元，因而分解过程中舍入误差的放大收到了控制，用平方根法解对称正定方程组时可以不考虑选主元问题。

线性方程组的四种数值解法(电子科技大学物理电子学院，四川 成都 610054)摘要：本文介绍了四种求解线性方程组的数值解法: 雅克比迭代法、高斯赛德尔迭代法、高斯消去法和改进的平方根法的基本原理和算法流程，通过求解具体方程，对四种求解方法进行了对比。

对于雅克比迭代法和高斯赛德尔迭代法，研究了两种算法对求解同一方程组的迭代效率差异，结果表明高斯赛德尔迭代法达到同样精度所需迭代次数较少。

对于高斯消去法，通过选择列主元的方法提高算法的准确度，计算结果表明高斯消去法计算精确，且运算复杂度也不是很高。

对于改进的平方根法，其运算复杂度低，但对于给定的方程组有着严苛的要求。

关键词：雅克比迭代法；高斯赛德尔迭代法；高斯消去法；改进的平方根法；线性方程组引言线性方程组的求解在日常生活和科研中有着极其重要的应用，但在实际运算中，当矩阵的维数较高时，用初等方法求解的计算复杂度随维数的增长非常快，因此，用数值方法求解线性方程组的重要性便显现出来。

经典的求解线性方程组的方法一般分为两类：直接法和迭代法。

前者例如高斯消去法,改进的平方根法等，后者的例子包括雅克比迭代法，高斯赛德尔迭代法等。

这些方法的计算复杂度在可以接受的范围内，因此被广泛采用。

一般来说，直接法对于阶数比较低的方程组比较有效；而后者对于比较大的方程组更有效。

在实际计算中，几十万甚至几百万个未知数的方程组并不少见。

在这些情况下，迭代法有无可比拟的优势。

另外，使用迭代法可以根据不同的精度要求选择终止时间，因此比较灵活。

在问题特别大的时候，计算机内存可能无法容纳被操作的矩阵，这给直接法带来很大的挑战。

而对于迭代法，则可以将矩阵的某一部分读入内存进行操作，然后再操作另外部分。

本文使用上述四种算法求解对应的方程组，验证各种算法的精确度和计算速度。

1 算法介绍1.1 雅克比迭代法 1.1.1 算法理论设线性方程组(1)b Ax的系数矩阵A 可逆且主对角元素 均不为零,令并将A 分解成 (2)从而(1)可写成令其中. (3)以B 1为迭代矩阵的迭代法(公式)(4)称为雅克比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为(5)其中为初始向量.1.1.2 算法描述 1给定迭代初始向量X 0以及误差要求delta 2根据雅克比迭代公式计算出下一组向量 3判断X 是否满足误差要求，即||X k+1 – X k || < delta4若误差满足要求，则停止迭代返回结果；若否，则返回第二步进行下一轮迭代1.2 高斯赛德尔迭代法nna ,...,a ,a 2211()nna ,...,a ,a diag D 2211=()D D A A +-=()b x A D Dx +-=11f x B x +=b D f ,A D I B 1111--=-=()()111f x B x k k +=+⎩⎨⎧[],...,,k ,n ,...,i x a ba xnij j )k (j j i iii)k (i21021111==∑-=≠=+()()()()()Tn x ,...x ,x x 002010=1.2.1 算法理论由雅克比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i 个分量时,已经计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯—塞德尔（Gauss-Seidel ）迭代法.把矩阵A 分解成(6)其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成即其中(7)以为迭代矩阵构成的迭代法(公式)(8)称为高斯—塞德尔迭代法(公式),用变量表示的形式为(9)1.2.2 算法描述 1给定迭代初始向量X 0以及误差要求delta2根据高斯赛德尔迭代公式计算出下一组向量()k x ()1+k x ()1+k ix ()()1111+-+k i k x ,...,x 1+k()1+k x()1+k jx U L D A --=()nna ,...,a ,a diag D 2211=U ,L --A ()b Ux x L D +=-22f x B x +=()()b L D f ,U L D B 1212---=-=2B ()()221f x B x k k +=+⎩⎨⎧[],...,,k ,n ,,i x a x a b a xi j n i j )k (j ij )k (j ij i ii)k (i21021111111==∑∑--=-=+=++3判断X 是否满足误差要求，即||X k+1 – X k || < delta4若误差满足要求，则停止迭代返回结果；若否，则返回第二步进行下一轮迭代1.3 高斯消去法 1.3.1 算法理论下面三种变换称为初等行变换：1.对调两行；2.以数k ≠0乘某一行中的所有元素；3.把某一行所有元素的k 倍加到另一行对应的元素上去。

实验名称：改进平方根法学院：___数学学院______________________班级姓名：学号：实验日期 2015 年 05 月 26 日自评成绩：97一、实验目的(1)熟练掌握改进平方根法和共轭梯度法的迭代过程(2)尝试使用自己熟悉的计算机语言解决数学中的问题(3)通过上机实验来巩固课本中所学的知识二、实验内容与结果题目1：改进平方根法源程序1#include<iostream>using namespace std;int main(){double a[100][100],l[100][100],u[100][100],b[10],y[10],x[10];int i,j,k,n;cout<<"请输入矩阵的行数: ";cin>>n;cout<<"请输入右端项： "<<endl;for(j=0;j<n;j++){cin>>b[j];}cout<<"请输入矩阵： "<<endl;for(i=0;i<n;i++){for(j=0;j<n;j++){cin>>a[i][j];}}for(j=0;j<n;j++){u[0][j]=a[0][j];}for(j=1;j<n;j++){l[j][0]=u[0][j]/u[0][0];}for(i=1;i<n-1;i++){double s=0;for(k=0;k<i-1;k++){s=s+l[i][k]*u[k][i];}u[i][i]=a[i][i]-s;for(j=i+1;j<n;j++){double s=0;for(k=0;k<i-1;k++){s=s+l[i][k]*u[k][j];}u[i][j]=a[i][j]-s;l[j][i]=u[i][j]/u[i][i];}}double s=0;for(k=0;k<n-1;k++){s=s+l[n-1][k]*u[k][n-1];}u[n-1][n-1]=a[n-1][n-1]-s;y[0]=b[0];for(i=1;i<n;i++){double s=0;for(k=0;k<i-1;k++){s=s+l[i][k]*y[k];}y[i]=b[i]-s;}x[n-1]=y[n-1]/u[n-1][n-1];for(i=n-2;i>=0;i--){double s=0;for(k=i+1;k<n;k++){s=s+u[i][k]*x[k];}x[i]=(y[i]-s)/u[i][i];}cout<<"输出结果："<<endl;for(i=0;i<n;i++){cout<<"x("<<i<<")"<<"="<<x[i]<<endl;}return 0;}运行结果1。

浅析线性方程组的平方根解法在求解线性方程组时,直接解法有顺序高斯消元法、列主元高斯消元法、全主元高斯消元法、高斯约当消元法、消元形式的追赶法、LU 分解法、矩阵形式的追赶法，当我们遇到对称正定线性方程组时，我们就要用到平方根法(对称LLT 分解法)来求解，为了熟悉和熟练运用平方根法求解线性方程组，下面对运用平方根法求解线性方程组进行解析。

一、运用平方根法求解线性方程组涉及到的定理及定义我们在运用平方根法求解线性方程组时，要判定线性方程组Ax=b 的系数矩阵A 是否是对称正定矩阵，那么我们就要了解正定矩阵的性质和如下定理及定义：1、由线性代数知，正定矩阵具有如下性质：1) 正定矩阵A 是非奇异的2) 正定矩阵A 的任一主子矩阵也必为正定矩阵 3) 正定矩阵A 的主对角元素均为正数 4) 正定矩阵 A 的特征值均大于零 5) 正定矩阵A 的行列式必为正数定义一 线性方程组Ax=b 的系数矩阵A 是对称正定矩阵，那么Ax=b 是对称正定线性方程组。

定义二 如果方阵A 满足A=AT ，那么A 是对称阵。

2.1.4 平方根法和改进的平方根法如果A 是n 阶对称矩阵，由定理2还可得如下分解定理：定理2 若A 为n 阶对称矩阵，且A 的各阶顺序主子式都不为零，则A 可惟一分解为：A ＝LDLT ，其中L 为单位下三角阵，D 为对角阵。

证明 因为A 的各阶顺序主子式都不为零，所以A 可惟一分解为：A ＝LU 因为 ，所以可将 U 分解为：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=nn u u u U 2211⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛11122211112 u u u u u u n nn n 1DU =其中 D 为对角矩阵，U1为单位上三角阵．于是：A ＝LDU1＝L(DU1)因为A 为对称矩阵，所以，A ＝AT ＝U1TDTLT ＝U1T(DLT)，由 A 的 LU 分解的惟一性即得：L ＝U1T ，即U1＝LT ，故A ＝LDLT 。

M A T L A B-平方根法和改进平方根法求解线性方程组例题与程序（2）设对称正定阵系数阵线方程组12345678424024000221213206411418356200216143323218122410394334411142202531011421500633421945x 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 xL by 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 aa a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦其中,,1,2,...,ij ji a a i j n ==第1步，由矩阵乘法，211111111,i i a l a l l ==g 故求得111111,2,3,...i i a l l i n a === 一般的，设矩阵L 的前k-1列元素已经求出 第k 步，由矩阵乘法得112211k k kk kmkkik im km ik kkm m a l l a l l l l --===+=+∑∑, 于是11(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、改进平方根法在平方根的基础上，为了避免开方运算,所以用TLDL A =计算；其中，11122.........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Λ=11(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.11i i ii ik ikk 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 分解不需要开放计算。

数值分析思考题11、讨论绝对误差（限）、相对误差（限）与有效数字之间的关系。

2、相对误差在什么情况下可以用下式代替？3、查阅何谓问题的“病态性”，并区分与“数值稳定性”的不同点。

4、取，计算，下列方法中哪种最好？为什么？（1）(33-，（2）(27-，（3）()313+，（4）()611，（5）99-数值实验数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。

求解方法大致可分为直接法和迭代法两大类。

直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。

当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。

如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。

Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。

对正定对称矩阵，采用平方根方法无需选主元。

方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。

数值计算方法上机题目11、实验1. 病态问题实验目的：算法有“优”与“劣”之分，问题也有“好”和“坏”之别。

所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。

希望读者通过本实验对此有一个初步的体会。

数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。

病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。

问题提出：考虑一个高次的代数多项式re x xex x*****-==141.≈)61∏=-=---=201)()20)...(2)(1()(k k x x x x x p （E1-1）显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。

现考虑该多项式方程的一个扰动0)(19=+xx p ε (E1-2)其中ε是一个非常小的数。

（2）设对称正定阵系数阵线方程组12345678424024000221213206411418356200216143323218122410394334411142202531011421500633421945x 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 xL by 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 aa a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦其中,,1,2,...,ij ji a a i j n ==第1步，由矩阵乘法，211111111,i i a l a l l ==故求得111111,2,3,...i i a l l i n a === 一般的，设矩阵L 的前k-1列元素已经求出 第k 步，由矩阵乘法得112211k k kk kmkkik im km ik kkm m a l l a l l l l --===+=+∑∑, 于是11(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、改进平方根法在平方根的基础上，为了避免开方运算,所以用TLDL A =计算；其中，11122.........n d D D D d ⎤⎤⎡⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎢⎢⎥⎣⎦⎣⎣；得1121212212111111n n n n n d l l l d l A l l d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦按行计算的L 元素及D 对元素公式 对于n i ,,2,1 =11(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.11i i ii ik ikk d a t l -==-∑计算出LD T =的第i 行元素(1,2,i-1)ij t j =…，后，存放在A 的第i 行相置，然后再计算L 的第i 行元素，存放在A 的第行.D 的对角元素存放在A 的相应位置.对称正定矩阵A 按T LDL 分解和按T LL 分解计算量差不多，但TLDL 分解不需要开放计算。

目录摘要 (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月。

浅析线性方程组的平方根解法在求解线性方程组时,直接解法有顺序高斯消元法、列主元高斯消元法、全主元高斯消元法、高斯约当消元法、消元形式的追赶法、LU 分解法、矩阵形式的追赶法，当我们遇到对称正定线性方程组时，我们就要用到平方根法(对称LLT 分解法)来求解，为了熟悉和熟练运用平方根法求解线性方程组，下面对运用平方根法求解线性方程组进行解析。

一、运用平方根法求解线性方程组涉及到的定理及定义我们在运用平方根法求解线性方程组时，要判定线性方程组Ax=b 的系数矩阵A 是否是对称正定矩阵，那么我们就要了解正定矩阵的性质和如下定理及定义：1、由线性代数知，正定矩阵具有如下性质：1) 正定矩阵A 是非奇异的2) 正定矩阵A 的任一主子矩阵也必为正定矩阵 3) 正定矩阵A 的主对角元素均为正数 4) 正定矩阵 A 的特征值均大于零 5) 正定矩阵A 的行列式必为正数定义一 线性方程组Ax=b 的系数矩阵A 是对称正定矩阵，那么Ax=b 是对称正定线性方程组。

定义二 如果方阵A 满足A=AT ，那么A 是对称阵。

2.1.4 平方根法和改进的平方根法如果A 是n 阶对称矩阵，由定理2还可得如下分解定理：定理2 若A 为n 阶对称矩阵，且A 的各阶顺序主子式都不为零，则A 可惟一分解为：A ＝LDLT ，其中L 为单位下三角阵，D 为对角阵。

证明 因为A 的各阶顺序主子式都不为零，所以A 可惟一分解为：A ＝LU 因为 ，所以可将 U 分解为：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=nn u u u U O 2211⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛11122211112M O ΛΛu u u u u u n nn n 1DU =其中 D 为对角矩阵，U1为单位上三角阵．于是：A ＝LDU1＝L(DU1)因为A 为对称矩阵，所以，A ＝AT ＝U1TDTLT ＝U1T(DLT)，由 A 的 LU 分解的惟一性即得：L ＝U1T ，即U1＝LT ，故A ＝LDLT 。

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

设，即其中第1步，由矩阵乘法，故求得一般的，设矩阵L的前k-l列元素已经求出第k步，由矩阵乘法得于是2、改进平方根法在平方根的基础上，为了避免开方运算，所以用计算；其中，；得按行计算的元素及对元素公式对于、、计算出的第行元素后，存放在的第行相置，然后再计算的第行元素，存放在的第行、的对角元素存放在的相应位置、对称正定矩阵按分解和按分解计算量差不多，但分解不需要开放计算。

求解, 的计算公式分别如下公式。

3、程序设计1、平方根法function[x]二pfpf(A,b)%楚列斯基分解求解正定矩阵的线性代数方程A二LL'先求LY二b再用L' X二Y即可以求出解X[n,n]=size(A) ;L(1, l)=sqrt(A(l, 1)) ；for k=2：nL(k, l)=A(k, 1)/L(1,1) ； endfor k=2： n-1 L(k, k) =sqrt (A(k, k)_sum(L(k, 1：k-1) > 2)) ； for i=k+l:n L(i,k) = (A(i,k)-stun(L(i, l：kT)、*L(k, l：kT)))/L(k,k)； endendL (n, n)=sqrt (A(n, n) -sum(L(n, 1: n-1)、2)) ； %解下三角方程组Ly二b相应的递推公式如下，求出y矩阵y二zeros (n, 1) ；%先生成方程组的因变量的位置，给定y的初始值for k=l： n j=l： 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+l：n； x(k) = (y (k)-U(k,j)*x(j))/U(k,k)；end » A二[4,2,-4,0,2,4,0,02,2,-1,- 2,1,3,2,01,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,32,1,-8,- 1,22,4,-10,-34,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、148160、152810、91202、01852、改进平方根法function[x]=improvecholesky(A,b,n)%用改进平方根法求解Ax=bL=zeros(n,n)； %L为n*n矩阵D=diag(n, 0) ； %D 为n*n 的主对角矩阵S=L*D； for i=l: n %L 的主对角元素均为 1 L(i, i) = l； endfor i=l: n for j=l：n %验证 A 是否为对称正定矩阵if (eig(A)<=0) (A(i, j)~=A(j,i))%A的特征值小于0或A非对称时，输出wrongdisp(wrong) ； break； endendendD(1,1)=A(1, 1) ；%将A 分解使得A=LDLTfor i二2：n for j=l：i~l S(i,j)=A(i,j)-sum(S(i,l：j~ l)*L(j,l： j-1))； L(i,l：i-l)=S(i,l：i-l)/D(l：i-l,l：i-l); end D(i, i)=A(i, i)-sum(S(i, 1： i-l)*L(i, 1： i-1)) ； endy=zeros (n, 1) ； % x, y 为阶矩阵x=zeros (n, 1); for i=l: n y(i) = (b(i)-sum(L(i, 1： i-l)*D(l： i-1,1： i-l)*y (1： i-1)))/D(i, i) ；%通过LDy=b 解得y 的值endfor i=n：-1:1 x(i)=y(i)-sum(L(i+l：n, i)*x(i+l：n)) ；%通过LTx=y 解得x 的值end» A二[4,2,-4,0,2,4,0,02,2,-1,-2,1,3,2,01,14,1,-8,- 3,5,6 0,-2,1,6,-1,-4,-3,32,1,-8,-1,22,4,-10,-34,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、148160、152810、91202、01854、结果分析和讨论平方根法和改进平方根法求解线性方程组的解为x二(121、1481, -140、1127,29、7515,-60、1528,10、9120,-26、7963,5、4259, -2、0185) To与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差，所以在求解的过程中还是应该选择主元以此消除误差，提高精度。