第五章 小行星轨道方程计算问题——线性方程组求解的直接法.doc

第五章 小行星轨道方程计算问题 线性方程组求解的直接法5.1 小行星轨道方程问题5.1.1 问题的引入一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，其单位为天文测量单位，在5个不同的时间对小行星作了5次观察，测得轨道上的5个点的坐标数据如下表：表5.1.1 轨道上的5个点的坐标数据123455.7646.2866.7597.1687.4080.6481.2021.8232.5263.360试确立小行星的轨道方程，并画出小行星的运动轨线图形。5.1.2 模型的分析 由开普勒第一定律知，小行星轨道为一椭圆，设椭圆的一般方程为： ，需要确定系数。利用已知的数据，不妨设，欲确定系数等价于求解一个线性方程组： 可写成矩阵的形式： 其中，。5.1.3 模型的假设假设：(1)小行星轨道方程满足开普勒第一定律；(2)以上所测得数据真实有效。5.1.4 模型的建立该问题的模型为：，可见，解答上述问题就是对线性方程进行求解。5.2 线性方程组直接解法概述直接法：利用一系列递推公式计算有限步，能直接得到方程组的精确解。当然，实际计算结果仍有误差，譬如舍入误差，舍入误差的积累有时甚至会严重影响解的精度。求解线性方程组最基本的一种直接法是消去法。这是一个众所周知的古老方法，但用在现代电子计算机上仍然十分有效。消去法的基本思想是，通过将一个方程乘以或除以某个常数，以及将两个方程相加减这两种手续，逐步减少方程中的变元的数目，最终使每个方程仅含一个变元，从而得出所求的解。其中高斯消去法是广泛应用的方法，其求解过程分为消元过程和回代过程两个环节。消元过程将所给的方程组加工成上三角方程组。所归结的方程组再通过回代过程得出它的解。高斯消去法由于添加了回代的过程，算法结构稍复杂，但这种算法的改进明显减少了计算量。直接法比较适用于中小型方程组。对高阶方程组，即使系数矩阵是稀疏的，但在运算中很难保持稀疏性，因而有存储量大，程序复杂等不足。5.3 直接解法5.3.1 高斯消去法消去法是一个古老的求解线性方程组的方法。由它改进的选主元法是目前计算机上常用的有效的求解低阶稠密矩阵线性方程组的方法。例5.3.1 用消去法解方程组解 第1步，加到(5.3.2)，加到(5.3.3)，得等价方程组： 第2步，加到(5.3.5)得等价的方程组：第3步，回代法求解(5.3.6)即可求得该方程组的解为：用矩阵法描述的约化过程即为：.这种求解过程称为具有回代的消去法。 此例可见消去法的基本思想是：用矩阵的初等行变换将系数矩阵化为具有简单形式的矩阵（如上三角阵，单位矩阵等），而三角形方程组是很容易回代求解的。 一般的，设有个未知数的线性方程组为： （5.3.7）设则(5.3.7)化为为方便，。则消去法为：第步：计算，用乘(5.3.7)的第一方程加到第个方程中 (即实行行的初等变换) 。消去第2个到第个方程中的未知数得与(5.3.7)等价方程组： (5.3.8) 记为：其中(5.3.8)式中元素为进一步需要计算的元素，公式为：。 第步，继续上述过程消元。设第1步到第步计算已完成，得到与原方程组等价的方程组： (5.3.9)记为，下面进行第步消元法：设，计算乘数用乘 (5.3.9)中第个方程加到第个方程消去(5.3.9)中第个方程的未知数，得到与原方程组等价的方程组： (5.3.10)记为。其中中元素计算公式为： (5.3.11)最后，重复上述过程，即且设共完成步消元计算，得到与(5.3.7)等价的三角形方程组。 (5.3.12)再用回代法求解(5.3.12)的解，计算公式为： (5.3.13) 元素称为约化的主元素。将(5.3.7)化为(5.3.12)的过程称为消元过程。 由消元过程和回代过程求解线性方程组的方法称为消去法。(5.3.12)的求解过程(5.3.13)称为回代过程。定理5.3.1 （Gauss消去法）设若约化的主元素则可通过Gauss消去法（不进行两行的初等变换两行交换位置）将方程组化为等价的三角形方程组(5.3.12)。消元和求解的计算公式为：消元计算 回代计算5.3.2 矩阵的三角分解下面用矩阵理论进一步来分析消去法，设约化主元素由于对实行的初等变换相当于用初等矩阵左乘，于是，消去法第1步：,有：其中：(为初等三角矩阵)消去法第步消元过程：则有 (5.3.14)其中：利用递推公式(5.3.14)则有： (5.3.15)由(5.3.15)得： (5.3.16)其中为由乘数构成的下三角阵，为上三角矩阵，(5.3.16)表明，用矩阵理论来分析消去法，得到一个重要结果，即在条件下消去法实质上是将分解成两个三角矩阵的 显然，可由消去法及行列式性质可知，如果，则有。其中（的顺序主子式）反之，可用归纳法证明：如果的顺序主子式满足：则总结以上讨论，可得如下重要定理：定理5.3.2 （矩阵的三角分解）设，如果的顺序主子式有，则可分解为一个单位下三角矩阵与一个上三角矩阵的乘积，即且分解是唯一的。证明 现仅就来证明唯一性，存在性上面已证。假若(5.3.17)且对非奇异时考虑，为单位下三角阵，为上三角阵，由假设知存在（因为可逆故可逆），从而由(5.3.17)有。上式右端为上三角阵，左端为单位下三角阵，因此左右两端应为单位矩阵。故即分解是唯一的。称矩阵的三角分解（杜利特尔）分解。其中 在以上定理条件下，同样可有下面的三角分解：。其中为下三角矩阵，为单位上三角矩阵，称之为（克劳特）分解。如前例中系数矩阵的分解为：现设，若有分解，则(1)(2)，而求解这两个三角形方程组是很容易的。5.3.3 Gauss消去法的计算量定理5.3.3 设为阶非奇异矩阵，则用消去法解所需要的乘除法次数及加减法的次数分别为：但如果用（克莱姆）法则解，就需要计算个阶行列式，若行列式用子式展开，总共需要次乘法，如时消去法需要次乘除法，而克莱姆法则需要次乘法，由此可见，法则方程组的工作量太大，不便于使用。如果计算是在每秒作次乘除法计算机上进的，那么用消去法解阶方程组约需秒即可完成，而用法则大约需小时才能完成（大约相当于年）可见，法则完全不适于在计算机上求解高维方程组。5.3.4 Gauss主元素消去法用消去法解时，设非奇异，可能出现，这时必须进行行交换的Gauss消去法，但在实际计算中即使，但其绝对值很小时，用作除数，会导致中间结果矩阵的元素数量级严重增长和舍入误差的扩散，使最后结果不可靠。例5.3.2 设有方程组解 方程组的解方法一：用消去法求解（用具有舍入的3位浮点数进行运算）。 回代得解，与精确解比较，这个结果很差。方法二：用具有行交换的消去法（避免小主元）。回代得解：。 这个解对于具有舍入的3位浮点数进行计算，是一个很好的结果。方法一计算失败的原因，是用了一个绝对值很小的数作除数，乘数很大，引起约化中间结果数量误差很严重地增长，再舍入就使得计算结果不可靠了。这个例子告诉我们，在采用高斯消去法解方程组时，小主元可能导致计算失败，故在消去法中应避免采用绝对值很小的主元素。 对一般矩阵方程组，需要引进主元的技巧，即在高斯消去法的每一步应该选取系数矩阵或消元后的低阶矩阵中的绝对值最大的元素作为主元素，保持乘数以便减少计算过程中的舍入误差对计算解的影响。这个例子还告诉我们，对同一个数值问题，用不同的计算方法，得到的精度大不一样。一个计算方法，如果用此方法的计算过程中舍入误差得到控制，对计算结果影响较小，称此方法为数值稳定的；否则，如果用此计算方法的计算过程中舍入误差增长迅速，计算结果受舍入误差影响较大，称此方法为数值不稳定的。因此，我们解数值问题时，应选择和使用数值稳定的计算方法，否则，如果使用数值不稳定的计算方法去解数值计算问题，就可能导致计算失败。5.3.5 完全主元素消去法设有线性方程组，其中为非奇异矩阵。方程组的增广矩阵为第1步：首先在中选主元素，即选择使，再交换的第一行与第行元素，交换的第一列与第列元素（相当于交换未知数与），将调到位置（交换后的增广矩阵为，其元素仍记为，然后进行消元法计算。 第步：继续上述过程，设已完成第步到第步计算，约化为下述形式（为简单起见，仍记元素为元素为）：第步选主元区域 于是第步计算：对于按下述步骤从(1)计算到(3)：(1) 选主元素：选择使；(2) 如果，则交换第行与第行元素，如果，则交换的第列与第列元素；(3) 消元计算 ；(4) 回代求解。经过上面的过程，即从第步到第步完成选主元，交换两行，交换两列，消元计算。原方程组约化为：其中为未知数调换后的顺序。回代求解：5.3.6 列主元消去法 完全主元素消去法是解低阶稠密矩阵方程组的有效方法，但完全主元素方法在选主元时要花费一定的时间。现介绍一种在实际计算中常用的部分选主元（即列主元）消去法。列主元消去法，即是每次选主元时，仅依次按列选取绝对值最大的元素作为主元素，且仅交换两行，再进行消元计算。 设列主元消去法已经完成第步到第步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组： 第步选主元区域第步计算如下：对于按下述步骤从(1)计算到(4)：(1) 按列主元，即确定使；(2) 如果，则为非奇异，停止计算；(3) 如果，则交换第行第行元素；(4) 消元计算： ；(5) 回代计算：计算解在常数项内得到。例5.3.3 用列主元消去法求解方程组解 精确解为（舍入值）： 回代计算得到计算解：本例是具有舍入的位浮点数进行运算，所得的计算解还是比较准确的。 下面是完全主元素消去法框图（图5.3.1）：输入 选主元素 否是否交换行输出 否交换列且消元计算回代（当） ）调整求解输入计算解图5.3.1 完全主元素消去法框图用完全主元素消元法解，可用一整型数组开始记录未知数次序，即，最后记录调整后未知数的坐标。系数阵存在二维数组内，常数项存在内，解存在数组内。例5.3.4 若在计算过程中，只取位有效数字，试用列主元素法求解：解 第一步，选为主列元，将(5.3.18) (5.3.20)对调位置，第二步，选为列主元，不需换行， 由回代求解得：与原方程组准确解比较可知，本题用位有效数字计算的列主元素法，其精度较高。 下面用矩阵运算来描述列主元素法：记 是初等排列阵（由交换单位矩阵的第行与行所得）。则列主元素法为：其中的元素满足，由(5.3.21)得：简记为：，其中。考察时的 (5.3.22)其中则由排列阵性质（左乘矩阵是对矩阵进行行变化）:已知为单位下三角阵，其元素绝对值，记。由(5.3.22)得：。其中为排列阵，为单位下三角阵，为上三角阵。这表明，对应用列主元素法相当于对先进行一系列行变换后对再应用消去法。在实际计算中我们只能在计算过程中做关于行的变换。有结论：定理5.3.4 （列主元素三角分解定理） 若为非奇异性矩阵，则存在排列矩阵使。其中为单位下三角阵，为上三角阵。 存放在的下三角部分，存放在的上三角部分。由整数型数组记录可知的情况。5.3.7 GaussJordan消去法 消去法总是消去对角线下方的元素。现考虑一种修正，即消去对角线下方和上方的元素。这即为消去法。 设用消去法已完成步，于是化为等价方程组其中：在第步计算时，考虑对上述矩阵的第行上、下都进行消元计算(1) 按列选主元素，即定义使；(2) 换行（当）：交换第行与第行元素；(3) 计算乘数 (可存放在的单元中) （）；(4) 消元计算；(5) 计算主元素。上述过程全部执行完后有： 这表明用方法将约化为单位矩阵，计算解就在常数项位置得到，因此用不着回代求解。用方法解方程组的计算量大约需要次乘除法，要比消去法大些。但用方法求一个矩阵的逆矩阵还是比较合适的。定理5.3.5 （法求逆矩阵）设的逆矩阵，即求阶矩阵，使其中为单位矩阵，将按列写成，为列向量，为单位列向量。于是求解，等价于求解个方程组：，所以我们可以用法求解。例5.3.5 对求。解故：5.3.8 Gauss消去法的变形(1)直接三角分解 为求解将进行分解。即将原问题转化为求解两个三角形方程组：(1)，求；(2)，求。不选主元的三角分解法设，且有，其中为单位下三角阵，上三角阵。即 (5.3.23)下面说明：的元素可由步直接计算出来，其中第步定出的第行和的第列元素。由(5.3.23)有：(5.3.24)，得的第1行元素(5.3.25)，得的第1列元素到第列元素。由(5.3.23)利用矩阵乘法有： 故 (5.3.26)又由(5.3.23)有：故： (5.3.27)因此可得的第行和第列的全部元素。直接分解法约需乘除法，和消去法计算量基本相当。对计算和式，可采用双精度累加以提高精度。选主元的三角分解法在直接三角分解中，如果计算将要中断，或者当绝对值很小时，按分解公式计算可能引起舍入误差的积累。但当为非奇异时，可通过交换的行实现矩阵的分解。因此可采用与列主之消去法类似的方法将直接三角分解法修改为部分选主之的三角分解法。设已完成第步分解，这时有：第步分解需要到(5.3.26) (5.3.27)两式，为避免(5.3.27)中用小的数作除数，先引进量：由(5.3.26)及定义，易见则由(5.3.27)有：若，则我们可以用作为交换的第行与行（但我们将交换后的新元素仍记为及）于是有。控制了误差传播，再进行第步分解。对一般的非奇异矩阵求逆的方法：则：。5.3.9 平方根法 利用对称正定矩阵的三角分解得到的求解对称正定方程组的一种有效方法平方根法。设为对称阵，即，且的所有顺序主子式均非零，则知可以唯一分解为：。为利用的非奇异性，将再分解为：为对角矩阵，为单位上三角阵，则： (5.3.30)又，由分解的唯一性即得：，代入上面(5.3.30)中得：。定理5.3.6 （对称阵的三角分解）设为阶对称阵，且的所有顺序主子式均非零。则可唯一分解为：。其中为单位下三角阵，为对角阵。 当为对称正定矩阵时，则的所有顺序主子式，而设，则于是：从而：，其中为下三角阵。定理5.3.7 （对称正定矩阵的分解） 若为阶对称正定矩阵，则存在一个实的非奇异下三角矩阵，使。当限定的对角元素为正时，这种分解是唯一的。 下面来考虑计算元素的公式，由，由矩阵乘法及（时）有：。于是得解对称正定方程组的平方根法计算公式： 求解，即求解两个三角形方程组：(1), (2)， 由(5.3.31)知，则，从而。 这表明分解过程中元素的数量级不会增长太快且对角元素恒为整数。于是不选主元素的平方根是一个数值稳定的方法。当求出的第列元素时，的第行亦得出，所以平方根法大约需要次乘除法，约为一般直接分解法计算量的一半。由于为对称阵，因此在计算机中只需存储的下三角部分元素，共需个元素。可用一维数组存储，即：。 矩阵的元素一维数组的表示为：。的元素存放在的相应位置。 由公式(5.3.31)可见，用平方根法解对称正定方程组时，计算的元素需要进行开方计算，为避免开方计算，可对平方根法进行改进。例5.3.6 用平方根法求解：精确解为。(1)分解计算： 故，(2) 求解两个三角方程：解，得，代入，解得：。5.3.10 追赶法在一些实际问题中常有解三对角线性方程组Ax=f的问题，即： (5.3.35)其中满足条件：(1) (2) (3) 对于具有条件(2)的方程组(5.3.35)，我们介绍下面的追赶法求解。追赶法具有计算量少，方法简单，算法稳定的特点。定理5.3.8 设有三对角线性方程组，且满足条件(2)，则A为非奇异矩阵。证明 用归纳法证明。显然时，有：现假设定理对n-1阶的满足条件(2)的三对角矩阵成立，求证对满足条件(2)的阶三对角矩阵定理亦成立。由条件 ，则利用消元法的第步有：，显然，其中，且有故知满足条件(2)，利用归纳设知，故定理5.3.9 设满足条件(2)为三对角阵。则的所有顺序主子式都不为零。即：证明 由于是满足条件(2)的阶三对角阵。因此的任一个顺序主子式亦满足条件(2)的阶三对角矩阵。由上一个定理即知：根据这一结论以及三角分解定理知，