数学建模与科学计算第五章.doc

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