第五章小行星轨道方程计算问题

第五章小行星轨道方程计算问题——线性方程组求解的直接法5.1小行星轨道方程问题问题的引入一天文学家要确定一颗小行星绕太阳运转的轨道，他在轨道平面内成立以太阳为原点的直角坐标系，其单位为天文丈量单位，在5个不同的时间对小行星作了5次察看，测得轨道上的5个点的坐标数据以下表：表轨道上的5个点的坐标数据x123455.7646.2866.7597.1687.408y0.6481.2021.8232.5263.360试确定小行星的轨道方程，并画出小行星的运动轨线图形。模型的剖析由开普勒第必定律知，小行星轨道为一椭圆，设椭圆的一般方程为：a1x22a2xya3y22a4x2a5y10，需要确定系数ai,i1,2,3,4,5;利用已知的数据，不如设xiyii1,2,3,4,5;欲确定系数ai等价于求解一个线性方程组：a1x122a2x1y1a3y122a4x12a5y110a1x222a2x2y2a3y222a4x22a5y210a1x322a2x3y3a3y322a4x32a5y310a1x422a2x4y4a3y422a4x42a5y410a1x522a2x5y5a3y522a4x52a5y510可写成矩阵的形式：AXb此中，x22xyy22x2ya11111111x22xyy22x22y2a212222Ax322x3y3y322x32y3，Xa3，b1x422x4y4y422x42y4a41x522x5y5y522x52y5a51模型的假定假定：小行星轨道方程知足开普勒第必定律；以上所测得数据真切有效。模型的成立该问题的模型为：33.22377.47010.419911.52801.2960a1139.513815.11151.444812.57202.4040a2145.684124.64333.323313.51803.6460a3151.380236.21276.380714.33605.0520a4154.878549.781811.289614.81606.7200a51可见，解答上述问题就是对线性方程进行求解。5.2线性方程组直接解法概括直接法：利用一系列递推公式计算有限步，能直接获得方程组的精准解。自然，实质计算结果仍有偏差，比如舍入偏差。舍入偏差的累积有时甚至会严重影响解的精度。求解线性方程组最基本的一种直接法是消去法。这是一个尽人皆知的古老方法，但用在现代电子计算机上仍旧十分有效。消去法的基本思想是，经过将一个方程乘或除以某个常数，以及将两个方程相加减这两种手续，逐渐减少方程中的变元的数目，最后使每个方程仅含一个变元，进而得出所求的解。此中高斯消去法是宽泛应用的方法，其求解过程分为消元过程和回代过程两个环节。消元过程将所给的方程组加工成上三角方程组。所归纳的方程组再经过回代过程得出它的解。高斯消去法因为增添了回代的过程，算法构造稍复杂，但这类算法的改良明显减少了计算量。直接法比较合用于中小型方程组。对高阶方程组，即便系数矩阵是稀少的，但在运算中很难保持稀少性，因此有储存量大，程序复杂等不足。5.3直接解法高斯消去法Gauss消去法是一个古老的求解线性方程组的方法。由它改良得选主元法是当前计算机上常用的有效的求解低阶浓密矩战线性方程组的方法。例用Gauss消去法解方程组2x12x22x31(5.3.1)3x12x24x31(5.3.2)2x13x29x35(5.3.3)3）2解第1步，（）（），1)加到(5.3.3)，得等价方程加到（5.3.1(22组：2x12x22x31x2x31(5.3.4)2x28x32(5.3.5)第2步，2加到（）得等价的方程组：2x12x22x31x2x3110x30(5.3.6)第3步，回代法求解（）即可求得该方程组的解为：x0,x1,x1.3212用矩阵法描绘的约化过程即为：2221(1)A,b3241/231395/2r1(2)r2r2r(1)rr1233

222122210111(2)0111.0282r22r3r300100这类求解过程称为拥有回代的Gauss消去法。此例可见Gauss消去法的基本思想是：用矩阵的初等行变换将系数矩阵A化为拥有简单形式的矩阵（如上三角阵，单位矩阵等），而三角形方程组是很简单回代求解的。一般的，设有n个未知数的线性方程组为：a11x1a12x2a1nxnb1a21x1a22x2a2nxnb2(5.3.7)an1x1an2x2annxnbn设A（aij)nn，X(x1,x2,xn)T，b(b1,b2,bn)T，则化为：AXb，为方便，AA(1)(aij(1))nn，bb(1)(b1(1),b2(1),bn(1))T，detA0。则消去法为：（1），a(1)mi1i1(i2,3,4n),，用mi1乘的第一方程加到第1步：a110计算a(1)11第i个方程中去（i2,3,n)。（即推行行的初等变换）Rmi1R1Ri(i2,3,n)消去第2个到第n个方程中的未知数x1,得与等价方程组：a11(1)a12(1)ain(1)x1b1(1)a22(2)a2(2)nx2b2(2)(5.3.8)a(2)a(2)xnb(2)n2nnn记为：A(2)Xb(2)此中（）式中元素aij(2)为进一步需要计算的元素，公式为：aij(2)aij(1)mi1a1(1)j，(i,j2,3n)，bi(2)bi(1)mi1b1(1),i2,3,n。第k，k1,2,n1步，持续上述过程消元。设第1步到第k1步计算已完成，获得与原方程组等价的方程组：a11(1)a12(1)a1(1)nx1b1(1)a22(2)a2n(2)b2(2)x2akk(k)akn(k)x3bk(k)(5.3.9)ank(k)ann(k)xnbn(k)记为A(K)Xb((k)，下边进行第k步消元法：(k)0，计算乘数mikaik(k)(k)(kk1,n)，用mik（）中第k个设akkakk方程加到第（ik1，，n)，（）中第i（ik1，，n)i个方程消去个方程的未知数xk，获得与原方程组等价的方程组：a11(1)a12(1)a1(1)nx1b1(1)(2)(2)a22a2nx2b2(2)akk(k)ak(k,k)1akn(k)bk(k)(k1)(k1)ak1,k1ak1,nbk(k1)1(k1)an,k1an(k,n1)xnbn(k1)记为A(k1)Xb(k1)此中A(k1),b(k1）中元素计算公式为：a(k1)a(k)ma(k)(i,jk1,n)ijijikkjbi(k1)bi(k)mikbk(k)(ik1.n)（）A(k1)与Ak前k行元素相同，b(k1)与b(k)前k个元素相同最后，重复上述过程，即k1,2,3n1;且设akk(k)0(k1,2,n1)，共完成n1步消元计算，获得与（）等价的三角形方程组。a(1)a(1)a(1)x1b(1)11121n1a22(2)a2(2)nx2b2(2)(5.3.12)ann(n)xnbn(n)再用回代法求解（）的解，计算公式为：xnb(n)ann(n)n(bi(i)aij(i)xj)ji1aii(i),(in1,n2,1)xi元素akk(k)称为约化的主元素。将（）化为（）的过程称为消元过程。由消元过程和回代过程求解线性方程组的方法称为Gauss消去法。（）的求解过程（）称为回代过程。定理（Gauss消去法）设AXb,ARnn。若约化的主元素akk(k)0,k1,2n则可经过Gauss消去法（不进行两行的初等变换—两行互换地点）将方程组化为等价的三角形方程组（）。消元和求解的计算公式为：10消元计算aik(k)mikakk(k)(ik1,n)aij(k1)aij(k)mikakj(k)(i,jk1,n),k1,2n1bi(k1)bi(k)mikbk(k)(ik1,n)20回代计算b(n)xnann(n)n(bi(i)aij(i)xj)ji1xiaii(i),(in1,n2,1)矩阵的三角分解下边用矩阵理论进一步来剖析Gauss消去法，设约化主元素a(k)0(k1,2n1)，A，于kk因为对A推行的初等变换相当于用初等矩阵左乘是，Gauss消去法第步：A(1)X(1)A(2)Xb(2),则有：1bL1A(1)A(2)L1b(1)b(2)此中：100L1m21100Gauss消去法第k步消元过程：（L1为初等三角矩阵）mn101A(k)Xb(k)A(k1)Xb(k1),则有LkA(k)A(k1),Lkb(k)b(k1),k1,2,n1(5.3.14此中：1第k列1第k列11Lk1Lk11mk1,k1mk1,k1mn,k1mn,k1利用递推公式则有：Ln1Ln2L2L1A(1)A(n)U,Ln1Ln2L2L1b(1)b(n)(5.3.15由(5.3.15)得：A=(L11L21Ln11)ULU此中1a11(1)a12(1)a1(1)nm211a22(2)a2(2)nLm31m321,Umn1mn21ann(n)L为由乘数组成的下三角阵，U为上三角矩阵，(5.3.16)表示，用矩阵理论来分析Gauss消去法，获得一个重要结果，即在akk(k)0,(k1,2n)条件下Gauss消去法实质上是将A分解成两个三角矩阵的ALU.明显，可由Gauss消去法及队列式性质可知，假如a(i)0(i1,2k)，ii则有det（A1)a11(1)0,det(Ai)a11(1)a22(2)aiii( )0(i2,3k)a11a1i此中A1（）（A的次序主子式）a11,Aiai1aii反之，可用概括法证明：假如A的次序主子式Ai知足：detAi()i0(1,k，2则,a(i)0。ii总结以上议论，可得以下重要定理：定理（矩阵的三角分解）设ARnn，假如A的次序主子式有detA(i)i0(1,，n-则A可分解为一个单位下三角矩阵与一个上三角矩阵的乘积，即ALU，且分解斯独一的。证明现仅就det(Ai)0(1,1)in来证明独一性，存在性上边已证。倘若A11且对A且对A非奇怪时考虑，L1,L为单位下三角阵，U1,ULULU（）为上三角阵，由假定知U11存在（因为detA0,L1可逆ALU11，故U1可逆），从而由（）有L1L1UU11，上式右端为上三角阵，左侧为单位下三角阵，因此左右两头应为单位矩阵。故L1L,U1U,即分解是独一的。称矩阵的三角分解ALU为Doolittle（杜利特尔）分解。此中1u11u12u1nLl211,Uu22u2nln1ln21unn在以上定理条件下，相同可有下边的三角分解：ALU，此中L为下三角矩阵，U为单位上三角矩阵，称之为Crout（克劳特）分解。如前例中系数矩阵A的分解为：2221222A3243111LU1392101212现设AXb，若如分解ALU，则AXbLUXb(1)LYbUXYX而求解这两个三角形方程组是很简单的。Guass消去法的计算量定理设A为n阶非奇怪矩阵，则用Guass消去法解AXb所需要的乘除法次数及加减法的次数分别为：（）MDn3n2nn31333(2)ASn(n1)(2n5)/6n33但假如用Gramer（克莱姆）法例解AXb，就需要计算n1个n阶队列式，若队列式用子式睁开，总合需要(n1)!次乘法，如n10时Gauss消去法需要430乘除法，而克莱姆法例却需要39916800次乘法，因而可知，Gramer法例方程组的工作量太大，不便于使用。假如计算是在每秒作105次乘除法计算机长进的，那么用Gauss消去法解20阶方程组约需0.03秒即可达成，而用Gramer法例大概需1.31011小时才能达成（大概相当于107年）可见，Gramer法例完整不适于在计算机上求解高维方程组。Gauss主元素消去法用Gauss消去法解AXb时，设A非奇怪，但可能出现akkk0，akkk0。但其绝对值很小时，用akkk作除数，会致使中间结果矩阵Ak的元素数目级严重增加和舍入偏差的扩散，使最后结果不行靠。0.0001x1x21例设有方程组x1x22解方程组得解X*(0.99989999,1.00010001)方法一：用Gauss消去法求解。（用拥有舍入的3位浮点数进行运算）0.00010011100000.00010011A,b1m210100001000012回代得解x21.00，x10.00，与精准解比较，这结果很差。方法二：用拥有行互换的Gauss消去法（防止小主元）A,br1r2112m20.0001001120.0001001101.001.00回代得解：x21.00,x11.00这个解对于拥有舍入的3位浮点数进行计算，是一个很好的结果。方法一计算失败的原由，是用了一个绝对值很小的数作除数，乘数很大，惹起约化中间结果数目偏差很严重增加，再舍入就使得计算结果不行靠了。这个例子告诉我们，在采纳高斯消去法解方程组时，小主元可能致使计算失败，故在消去法中应防止采纳绝对值很小的主元素。对一般矩阵方程组，需要引进主元的技巧，即在高斯消去法的每一步应当选用系数矩阵或消元后的低阶矩阵中的绝对值最大的元素作为主元素，保持乘数mik1以便减少计算过程中的舍入偏差对计算解的影响。这个例子还告诉我们，对同一个数值问题，用不同的计算方法，获得的精度大不相同，一个计算方法，假如用此方法的计算过程中舍入偏差获得控制，对计算结果影响较小，称此方法为数值稳固的；不然，假如用此计算方法的计算过程中舍入偏差增加快速，计算结果受舍入偏差影响较大，称此方法为数值不稳固。所以，我们解数值问题时，应选择和使用数值稳固的计算方法，不然，假如使用数值不稳固的计算方法去解数值计算问题，便可能致使计算失败。完整主元素消去法设有线性方程组AXb，此中A为非奇怪矩阵。方程组得增广矩阵为a11a12.a1nb1a21a22a2nb2A，bai1j1an1an2annbn第1步（k1）：第一在A中选主元素，即选择，，max1inaij0，再交i1j1使ai1j11jn换（A，b）的第一行与第i1行元素，互换A的第一列与第j1列元素（相当于互换未知数x1与xj1），将aj1（11，）地点（互换后的增广矩阵为A，b，其元素仍记i1调到为aij，bi，而后进行消元法计算。第k步：持续上述过程，设已达成第1步到第k1步计算，A，b约化为下述形式（为简单起见，仍记（）（k）元素为bi）：Ak元素为aij，ba11a12an1b1a22an2b2A,bakkaknbkankannbn第k步选主元地区于是第k步计算：对于k=1，2，，n-1按下述步骤从1计算到31选主元素：选择ik,jk使aik,jkaij0maxkjn(2)假如ikk，则互换A,b第k行与第ik行元素，假如jkk，则互换A的第k列与第jk列元素。（3）消元计算mikaik,(ik1,.........,n)akkaijaijmikak(ji,jk1,.....n...bibimikbk(ik1,.....n...回代求解。经过上边的过程，即从第1步到第n-1步达成选主元，互换两行，互换两列，消元计算，原方程组约化为：a11a12a1ny1b1a22a2ny2b2annynbn此中y1，y2.......，yn为未知数x1，x2.....xn调动后的次序。回代求解：bnynannn(biaijyj)ji1yi,(in1,........,2,1)aii列主元消去法完整主元素消去法是解低阶浓密矩阵方程组的有效方法，但完整主元素方法在选主元时要花销必定的时间。现介绍一种在实质计算中常用的部分选主元，（即列主元）消去法。列主元消去法，即是每次选主元时，仅挨次按列选用绝对值最大的元素作为主元素，且仅互换两行，再进行消元计算。设列主元消去法已经达成第1步到第k-1步的按列选主元，互换两行，消元计算获得与原方程组等价的方程组：a11(1)a12(1)a1n(1)b1(1)(2)(2)b2(2)a22a2nA(k)akk(k)akn(k)bkbk(k)bn(k)ank(k)ann(k)第k步选主元地区第k步计算以下：对于k=1，2，，n-1按下述步骤从1计算到41按列主元，即确定ik使aik,kmaxaik0kin2假如aik,k0，则A为非奇怪，停止计算。3假如ikk，则互换A,b第ik行第k行元素消元计算：aik

mik

aik

a,k(k

i

k1,....n..

)aij

aij

mik

*akj

(i,j

k1,...n..bi

bi

mik

*bk

(i

k1,...n...bnbn/ann5回代计算：n(in1,......2,1)bibiaijbj/aiij=i+1计算解x1,x2......xn在常数项内bn获得0.729x10.81x20.9x30.6867例用列主元消去法求解方程组x1x2x30.83381.331x11.21x21.1x31.000解精准解为（舍入值）：x*(0.2245,0.2841,0.3279)T0.72900.81000.90000.6867A,b1.0001.0001.0000.8338(r1r3,m210.7513,m310.5477)1.3311.2101.1001.0001.3311.2101.1001.00000.090900.17360.0825000.14730.29750.1390r3r2,m320.61711.3311.2101.1001.00000.14730.29750.1390000.010000.003280回代计算获得计算解：x0.2246,0.2812,0.3280T本例是拥有舍入的4位浮点数进行运算，所得的计算解仍是比较正确的。下边是完整主元素消去法框图（图）输入n，A，b，cIZII1,2,,n选主元素k=1，2，，n-1aik,jkmaxaijkinkjnc0aik,jk是输回代（当）ann0）c0出detA0bnannbn否bibiaijxjaiiin1,.......2,1ikkstop调整求解否k（）ibjkxIZG互换行i，，12.......n否互换列且输入计算解IZ(K)IZ(kj)xii1,2.....n消元计算stopaikmikaikakkaijaijmikakjbibimikbkik1,.......n(jk1,......n)n1图完整主元素消去法框图用完整主元素消元法解AXb，可用一整型数组IZn开始记录未知数x1，x2，.........xn序次，即IZn1，2......n，最后记录调整后未知数的足标。系数阵A存在二维数组An，n内，常数项b存在bn内，解存在数组Xn内。例若在计算过程中，只取3位有效数字，试用列主元素法求解：0.50x11.10x23.10x36.002.00x14.50x20.36x30.025.00x10.96x26.50x30.96解第一步，选5.00为主列元，将对换地点，5.00x10.96x26.50x30.962.00x14.50x20.36x30.020.50x11.10x23.10x36.002.000.40,m310.50m210.105.005.00m214.12x22.24x30.364*。*m311.00x22.45x35.905.320第二步，选4.12为列主元，不需换行，m321.000.2434.12*m322.99x35.99**由***回代求解得：x32.00，x21.00，x12.60与原方程组正确解x12.60，x21.00，x32.00比较，可知，本题用3位有效数字计算的列元素法是相当精准的。大批实践表示：列主元法为解线性方程组的精准方法。下边用矩阵运算来描绘列主元素法记Ik,ik是初等排排阵（由互换单位矩阵I的第k行与ik行所得）。则列主元素法为：L1I1,iA(1)A(2),L1I1,ib(1)b(2)11LkIk,ikA(k)A(k1),LkIk,ikb(k)b(k1)此中Lk的元素知足mik1k1，2，，n-1，由得：Ln1In1,in1L2I2,iL1I1,iAA(n)U21简记为：PAU，Pbb(n)，此中PLn1In1,in1L2I2,i2L1I1,i。下边观察n4时的1P。（）L3I3,iL2I2,iL1I1,iAUA4231L3I3，i3L2I3，i3I3，i3I2，i2L1I2，i2I3，i3I3，i3I2，i2I1，i1AL3L2L1PA(5.3.22)此中L1I3，i3I2，i2L1I2，i2I3，i3L2I3，i3L2I3，i3L3L3PI3，i3I2，i2I1，i1则由排排阵性质（左乘矩阵是对矩阵进行行变化。）已知Lk(k1,2,3)为单位下三角阵，其元素绝对值1，记L1L3L2L1。由（）得：PALU。此中P为排排阵，L为单位下三角阵，U为上三角阵。这表示，对AXb应用列主元素法相当于对A|b先进行一系列行变换后对PAXPb再应用消去法。再实质计算中我们只好在计算过程中，做对于Gauss行的变换。有结论：定理（列主元素三角分解定理）若A为非奇怪性矩阵，则存在摆列矩阵P使PA=LU。此中L为单位下三角阵，U为上三角阵。L寄存在A的下三角部分，U寄存在A的上三角部分。由整数型数组IP(n)记录可知P的状况。GaussJordan消去法Gaus消s去法老是消去对角线下方的元素。现考虑一种修正，即消去对角线下方和上方的元素。这即为GaussJordanGJ消去法。设用GJ消去法已达成k-1步，于是Axb化为等价方程组Akxbk此中：1a1,ka1,n|b11Akbk1ak1,kak1,n|ak,kak,n|bkan,kan,n|bn在第k步计算时，考虑对上述矩阵的第k行上、下都进行消元计算k1,2,,n1、按列选主元素，即定义ik使aik,kmaxaikkin2、换行（当ik≠k）。互换[A,b]第k行与第ik行元素。3、计算乘数mikaik/akki1n,ik,mkk1/akk(mik可寄存在aik的单元中)（mik1）aijmikakj（——且）4、消元计算aiji1n,jk1nikbimikbk（—且）bii1n,ikakjmkk（k,k—）5、计算主元素akjj1nbkmkkbk上述过程所有履行完后有：1?b1A,bA(kb(k1?1)1)b21?bn这表示用GJ方法将A约化为单位矩阵，计算解就在常数项地点获得，因此用不着回代求解。用GJ方法接方程组的计算量大概需要n3/2次乘除法，要比Gauss消去法大些。但用GJ方法求一个矩阵的逆矩阵仍是比较适合的。定理（GJ法求逆矩阵）设A的逆矩阵A1，即求n阶矩阵X，使AXIn此中In为单位矩阵，将X按列写成Xx1,x2,,xnIne1,e2,,en，xi为列向量，ei为单位列向量。于是求解AXIn，等价于求解n个方程组：Axieii1,2,n,，所以我们能够用GJ法求解AXIn。182求A1。例对A649104345解182100AI364910010r16r21821004345001012610182100r4rr(1)023401131649100104345001r28r110184980012610r2(2)r3001821r3(1)r3100952818r3(2)r20101032001821r3(18)r1故：9528181821Gauss消去法的变形直接三角分解为求解AXb将A进行L、U分解。即A将原问题转变为求解两个LU，三角形方程组（1）Lyb求y(2)Uxy求x。不选主元的三角分解法设detA0且有ALU此中L为单位下三角阵，U上三角阵。即1u11u12u1nl211u22u2nAl31l321（）ln1ln2ln31unn下边说明：L、U的元素可由n步直接计算出来，此中第r步定出U的第r行和L的第r列元素。由（）有：a1iu1i(i1,2,得U的第1行元素ai1(i2,3,得L的第1列元素到第r1列元ai1li1u11li1u11素。由利用矩阵乘法有：nr1arilrkukilrkukiuri(rk时lrk0)k1k1r1故uriarilrkuki,ir,r1,,n(5.3.26k1nr1又由有：airlikukrlikukrlirurrk1k1故：r1airlikukr(ir,r1,,n)k1(5.3.27)lirurr所以可得U的第r行和第r列的所有元素。L、U确定后，求解Lyb和Uxy的计算公式为：y1b1i1(5.3.28)yibilikyk(i2,3,,n)k1xnynunnnxiyiuikxk(in1,n2,,1)(5.3.29)ki1uii直接分解法约需n33乘除法，和Gauss消去法计算量基真相当。对计算和式，可采纳双精度累加以提升精度。选主元的三角分解法在直接三角分解中，假如urr0计算将要中止，或许当urr绝对值很小时，按分解公式计算可能惹起舍入偏差的累积。但当A为非奇怪时，可经过互换A的行实现矩阵PA的LU分解。所以可采纳与列主之消去法近似的方法将直接三角分解法改正为部分选主之的三角分解法。设已达成第r1步分解，这时有：u11u12u1nl21u22u2nl31l32u33Aur1,r1ur1,nlr,r1arrarnln1ln2ln,r1anrann第r步分解需要到（）（）两式，为防止（）中urr用小的数r1urr作除数，先引进量：Siairlikukr(ir,r1,,n)，由（）及Si定k1义，易见urrSr则由（）有：lirSi1,,n)。Sr(ir若SirmaxSi，则我们能够用Sir作为urr互换A的第r行与ir行（但我们rin将互换后的新元素仍记为lij及aij）于是有lir1(ir1,,n)。控制了偏差流传，再进行第r步分解。对一般的非奇怪矩阵A求逆的方法：PALU则：A1U1L1P平方根法利用对称正定矩阵的三角分解而获得的求解对称正定方程组的一种有效方法——平方根法设A为对称阵，即ATA，且A的所有次序主子式均非零，则知A可以独一分解为：ALU为利用U的非奇怪性，将U再分解为：u111u12u11u1nu11u22DU0U1un1,nun1,n1unn1D为对角矩阵，U0为单位上三角阵，则：ALULDU0（）又TTTT（）AA（）U0DL，由分解的独一性即得：U0L代入上边中得：ALDLT。定理（对称阵的三角分解）设A为n阶对称阵，且A的所有次序主子式均非零。则A可独一分解为：ALDLT。此中L为单位下三角阵，D为对角阵。当A为对称正定矩阵时，则A的所有次序主子式Di0，而ALDLT。设Ddiag(d1,d2,,dn),，则d1D10,Di2,3,,n)于是：diDi10(id1d2Ddnd1d1d2d2D1/2D1/2dndn进而：ALDLTLD1/2D1/2T1/2)(LD1/2TTL（LD）L1L1此中L1LD1/2为下三角阵。定理（对称正定矩阵的Cholesky分解）若A为n阶对称正定矩阵，则存在一个实的非奇怪下三角矩阵L，使LLT。当限制L的对角元素为正时，这类分解是独一的。下边来考虑计算L元素的公式，由l11l11l21ln1l21l22l22ln2A，由矩阵乘法及ljk0（jk时）ln1ln2nnlnnlii0(i1,2,,n)nj1ljjlij。于是得解对称正定方程组AXb的平方根法有：aijlikljklikljkk1k1计算公式：j1l11a11,ljj(ajjljk2)1/20,(j1,2,,n)(5.3.31)k1j1j1,2,,nl21a21a11,lij(aijk1likljk)ljj(ij1,,n)(5.3.32)求解AXb即求解两个三角形方程组：(1)Lyb求y(2)LTxy求xi1yi(bilikyk)lii(i1,2,,n)(5.3.33)1nxi(yilkixk)lii(in,n1,,1)(5.3.34)ki1j，则lj2由知ajjljk2(j1,2,,n)kajmjaax，j从j而k11jnmaxljk2maxajjM。j.k1jn这表示分解过程中元素ljk的数目级不会增加太快且对角元素ljj恒为整数。于是不选主元素的平方根是一个数值稳固的方法。当求出L的第j列元素时，LT的第j行亦得出，所以平方根法大概需要n36次乘除法，约为一般直接LU分解法计算量的一半。因为A为对称阵，所以在计算机中只要储存A的下三角部分元素，共需n(n1)2个元素。可用一维数组存储，即：A(n(n1)2)a11,a21,a22,,an1,an2,,ann。矩阵A的元素aij一维数组的表示为：A（i(i1)j)。L的元素寄存在A的相2应地点。由公式可见，用平方根法解对称正定方程组时，计算L的元素lii需要进行开方计算，为防止开方计算，可对平方根法进行改良。675x19例用平方根法求解：7138x210精准解为x*(1,1,2)T。586x39ALLT,l11a1162.4495,l21a21l11762.85771分解计算：l21249l22a22132.19856l31a31l11562.0412,l32(a32l31l21)l220.9856l33a33l312l3220.9285故，Ll11002.449500l21l2202.85772.19850l31l32l332.04120.98560.92852求解两个三角方程：解Lyb得y(3.6741,0.2273,1.8570)T代入LTxy解得：x（1.0,1.0,2.0)T。追赶法在一些实质问题中常有解三对角线性方程组Ax=f的问题，即：b1c1x1f1a2b2c2x2f2aibicixi(5.3.35)fian1bn1cn1xnfnanbn此中A知足条件：(1).|b1||c1|0(2).|bi||ai||ci|aici0,i2,...,n1(5.3.36)(3).|bn||an|0对于拥有条件2的方程组，我们介绍下边的追赶法求解。追赶法具有计算量少，方法简单，算法稳固的特色。定理设有三对角线性方程组Ax=f，且A知足条件2，则A为非奇怪矩阵。证明用概括法证明。明显n=2时，有：b1c1b1b2c1a20detAb2a2不然b1a2，由条件(）b1a21矛盾c1b2b2c1现假定定理对n-1阶的知足条件2的三对角矩阵成立，求证对知足条件2的n阶三对角矩阵定理亦成立。由条件b10，则利用消元法的第1步有：b1c10000b2c1a2c200b1Aa3b3c30A1an1bn1cn1ancn2c2a3b3c3c1a2，且明显，detAb1detB，此中B,2b2an1bn1cn1b1anbn有c1c1|0，故知B知足条件（），利||ba||b||||a||b||a||c22b1b1用概括设知detB0，故detA0。定理设Axf，A（2）为三对角阵。则A的所有次序主子式都知足条件不为零。即：detAk0k1,2,........n。证明因为A是知足条件（2）的n阶三对角阵。所以A的任一个次序主子式亦知足条件（2）的n阶三对角矩阵。由上一个定理即知：detAk0,(k1,...n)。依据这一结论以及三角分解定理知，这类矩阵A可进行三角分解：A=LU。在这里特其他有：b1c1111a2b2c2r2212r3313(5.3.37)an1bn1cn1n1anbnrnn1此中待定系数i,i,i（）可由利用矩阵乘法例则立刻得出：1.b11,c111c1112.airi,biirii1aii1i,(i2,...n)(5.3.38)3.ciii,(i2,...n)由1b10,|b1||c1|0。1c1。故0|1|1。b1下边概括证明：|i||ci|0，进而：0|i|1。i1,...n1上边已经考证对i1成立。现设对i1成立。求证对i成立。由假定0|i1|1，再由及假定条件2有: