[数学]数学建模中南大学数模课件

1、科学计算与数学建模中南大学数学科学与计算技术学院 小行星轨道方程计算问题2022/7/8第五章 小行星轨道方程计算问题 线性方程组求解的直接法小行星轨道方程问题 5.1线性方程组直接解法概述5.2直接解法5.3小行星轨道方程问题的模型求解5.42022/7/85.1.1 问题的引入一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，其单位为天文测量单位，在5个不同的时间对小行星作了5次观察，测得轨道上的5个点的坐标数据如下表：表 5.1.1 轨道上的5个点的坐标数据试确立小行星的轨道方程，并画出小行星的运动轨线图形。123455.7646.2866.7597

2、.1687.4080.6481.2021.8232.5263.360 5.1 小行星轨道方程问题 2022/7/85.1.2 模型的分析 由开普勒第一定律知，小行星轨道为一椭圆，设椭圆的一般方程为： ，需要确定系数 利用已知的数据，不妨设 欲确定系数 等价于求解一个线性方程组： 可写成矩阵的形式： 2022/7/8 其中， ，2022/7/8 5.1.3 模型的假设 假设：(1)小行星轨道方程满足开普勒第一定律；(2)以上所测得数据真实有效。 5.1.3 模型的建立 该问题的模型为： 可见，解答上述问题就是对线性方程进行求解。 2022/7/8 5.2 线性方程组直接解法概述 直接法：利用一系

3、列递推公式计算有限步，能直接得到方程组的精确解。当然，实际计算结果仍有误差，譬如舍入误差。舍入误差的积累有时甚至会严重影响解的精度 求解线性方程组最基本的一种直接法是消去法。这是一个众所周知的古老方法，但用在现代电子计算机上仍然十分有效。消去法的基本思想是，通过将一个方程乘或除以某个常数，以及将两个方程思想是，通过将一个方程乘或除以某个常数，以及将两个方程相加减这两种手续，逐步减少方程中的变元的数目，最终使每个方程仅含一个变元，从而得出所求的解。其中高斯消去法是广泛应用的方法，其求解过程分为消元过程和回代过程两个环节。消元过程将所给的方程组加工成上三角方程组。所归结的方程组再通过回代过程得出它

4、的解。高斯消去法由于添加了回代的过程，算法结构稍复杂，但这种算法的改进明显减少了计算量。 直接法比较适用于中小型方程组。对高阶方程组，即使系数矩阵是稀疏的，但在运算中很难保持稀疏性，因而有存储量大，程序复杂等不足。 2022/7/85.3.1 高斯消去法 消去法是一个古老的求解线性方程组的方法。由它改进得选主元法是目前计算机上常用的有效的求解低阶稠密矩阵线性方程组的方法。例 5.3.1 用 消去法解方程组5.3 直接解法2022/7/8解：第 1 步， 加到 ， 加到 ，得等价方程组： 第 2 步， 加到 得等价的方程组：2022/7/8 第 3 步，回代法求解 即可求得该方程组的解为： 用矩

5、阵法描述的约化过程即为： 这种求解过程称为具有回代的高斯消去法。 此例可见Gauss消去法的基本思想是：用矩阵A的初等行变换将系数矩阵化为具有简单形式的矩阵（如上三角阵，单位矩阵等），而三角形方程组是很容易回代求解的。2022/7/8 一般的，设有个未知数的线性方程组为： 2022/7/8 设 则 化为： 为方便， 则消去法为： 第1步： 计算 用 乘 的第一方程加到第 个方程中去 （即实行行的初等变换） ，消去第2个到第n个方程中的未知数 得与 等价方程组：2022/7/8 记为： 其中 式中元素 为进一步需要计算的元素，公式为：第 ， 步，继续上述过程消元。设第 步到第 步计算已完成，得到

6、与原方程组等价的方程组： 记为 ，下面进行第 步消元法：2022/7/8设 ，计算乘数 用 中第 个方程加到第 个方程 消去 中第 个方程 的未知数 得到与原方程组等价的方程组： 2022/7/8记为 其中 中元素计算公式为： 最后，重复上述过程，即 且设 共完成 步消元计算，得到与 等价的三角形方程组。2022/7/8 再用回代法求解 的解，计算公式为： 元素 称为约化的主元素。将 化为 的过程称为消元过程。 由消元过程和回代过程求解线性方程组的方法称为Gauss 消去法。 的求解过程 称为回代过程。2022/7/8定理（Gauss消去法）设 若约化的主元素 则可通过Gauss消去法（不进行

7、两行的初等变换两行交换位置）将方程组化为等价的三角形方程组 。消元和求解的计算公式为： 1、消元计算 2、回代计算2022/7/85.3.2 矩阵的三角分解 下面用矩阵理论进一步来分析 Gauss消去法，设约化主元素 由于对 实行的初等变换相当于用初等矩阵左乘于是， 消去法第1步： , 则有： 其中： （ 为初等三角矩阵） 2022/7/8Gauss消去法第k步消元过程： 则有其中： 2022/7/8利用递推公式则有：由 得 : 2022/7/8其中L为由乘数构成的下三角阵,U为上三角矩阵， 表明，用矩阵理论来分析Gauss消去法，得到一个重要结果，即在 条件下Gauss消去法实质上是A将分解

8、成两个三角矩阵的2022/7/8显然，可由Gauss消去法及行列式性质可知，如果 则有 其中 为顺序主子式反之，可用归纳法证明：如果A的顺序主子式 满足：则 2022/7/8定理 5.3.2 （矩阵的三角分解）设 ，如果 的顺序主子式有 ，则 可分解为一个单位下三角矩阵与一个上三角矩阵的乘积，即 且分解是唯一的。证明 现仅就 来证明唯一性，存在性上面已证。假若 且对 非奇异时考虑， 为单位下三角阵， 为上三角阵，由假设知 存在（因为 可逆 故 可逆），从而由 有 ，上式右端为上三角阵，左边为单位下三角阵，因此左右两端应为单位矩阵。故 即分解是唯一的。 称矩阵的三角分解 （杜利特尔）分解。202

9、2/7/8 其中在以上定理条件下，同样可有下面的三角分解： ，其中L为下三角矩阵，U为单位上三角矩阵，称之为Crout（克劳特）分解。如前例中系数矩阵A的分解为： 现设 ，若如分解 ，则 而求解这两个三角形方程组是很容易的。2022/7/8定理 5.3.3 设A为n阶非奇异矩阵，则用Gauss消去法解 所需要的乘除法次数及加减法的次数分别为：5.3.3 Gauss消去法的计算量但如果用 （克莱姆）法则解 ，就需要计算 个n阶行列式，若行列式用子式展开，总共需要 次乘法，如 时 消去法需要430乘除法，而克莱姆法则却需要39916800次乘法，由此可见， 法则方程组的工作量太大，不便于使用。如果

10、计算是在每秒作 次乘除法计算机上进的，那么用 消去法解20阶方程组约需0.03秒即可完成，而用 法则大约需 小时才能完成（大约相当于 年）可见， 法则完全不适于在计算机上求解高维方程组。2022/7/8用Gauss消去法解 时，设A非奇异，但可能出现 ，这时必须进行行交换的Gauss消去法，但在实际计算中即使 ，但其绝对值很小时，用 作除数，会导致中间结果矩阵 的元素数量级严重增长和舍入误差的扩散，使最后结果不可靠。例5.3.2 设有方程组5.3.4 Gauss主元素消去法2022/7/8解 方程组得解方法一：用Gauss消去法求解。（用具有舍入的3位浮点数进行运算） 回代得解 ， ，与精确解

11、比较，这结果很差。2022/7/8方法二：用具有行交换的Gauss消去法（避免小主元） 回代得解：这个解对于具有舍入的3位浮点数进行计算，是一个很好的结果。方法一计算失败的原因，是用了一个绝对值很小的数作除数，乘数很大，引起约化中间结果数量误差很严重增长，再舍入就使得计算结果不可靠了。2022/7/8这个例子告诉我们，在采用高斯消去法解方程组时，小主元可能导致计算失败，故在消去法中应避免采用绝对值很小的主元素。 对一般矩阵方程组，需要引进主元的技巧，即在高斯消去法的每一步应该选取系数矩阵或消元后的低阶矩阵中的绝对值最大的元素作为主元素，保持乘数 以便减少计算过程中的舍入误差对计算解的影响。20

12、22/7/8 这个例子还告诉我们，对同一个数值问题，用不同的计算方法，得到的精度大不一样，一个计算方法，如果用此方法的计算过程中舍入误差得到控制，对计算结果影响较小，称此方法为数值稳定的；否则，如果用此计算方法的计算过程中舍入误差增长迅速，计算结果受舍入误差影响较大，称此方法为数值不稳定。因此，我们解数值问题时，应选择和使用数值稳定的计算方法，否则，如果使用数值不稳定的计算方法去解数值计算问题，就可能导致计算失败。2022/7/8设有线性方程组，其中为非奇异矩阵。方程组得增广矩阵为 第 步 ：首先在 中选主元素，即选择 使 ，再交换 的第一行与第 行元素，交换 的第一列与第 列元素（相当于交换

13、未知数 与 ），将 调到 位置（交换后的增广矩阵为 ，其元素仍记为 ，然后进行消元法计算。 5.3.5 完全主元素消去法2022/7/8第 步：继续上述过程，设已完成第 步到第 步计算， 约化为下述形式（为简单起见，仍记 元素为 元素为 ）：域 于是第 步计算：对于 按下述步骤从 计算到第 步主元区域2022/7/8（1）选主元素：选择 使 （2）如果 ，则交换 第 行与第 行元素，如果 ，则交换 的第 列与第 列元素。（3）消元计算2022/7/8（4）回代求解。经过上面的过程，即从第 步到第 步完成选主元，交换两行，交换两列，消元计算，原方程组约化为：其中 为未知数 调换后的顺序。回代求解

14、：2022/7/8 5.3.6 列主元消去法 完全主元素消去法是解低阶稠密矩阵方程组的有效方法，但完全主元素方法在选主元时要花费一定的时间。现介绍一种在实际计算中常用的部分选主元，（即列主元）消去法。列主元消去法，即是每次选主元时，仅依次按列选取绝对值最大的元素作为主元素，且仅交换两行，再进行消元计算。 设列主元消去法已经完成第 步到第 步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组： 2022/7/8 第 步选主元区域2022/7/8第 步计算如下：对于 , 按下述步骤从 计算到 按列主元，即确定 使 如果 ，则 为非奇异，停止计算。 如果 ，则交换 第 行第 行元素 消元计算

15、：2022/7/8（5）回代计算： 计算解 在常数项内 得到例 5.3.3 用列主元消去法求解方程组解 精确解为（舍入值）：2022/7/8回代计算得到计算解： 本例是具有舍入的位浮点数进行运算，所得的计算解还是比较准确的。 下面是完全主元素消去法框图（图 ） 2022/7/8输入 选主元素 否是否交换行输出 否交换列且消元计算回代（当 ）调整求解输入计算解图5.3.1 完全主元素消去法框图2022/7/8用完全主元素消元法解 ，可用一整型数组 开始记录未知数 次序，即 ，最后记录调整后未知数的足标。系数阵 存在二维数组 内，常数项 存在 内，解存在数组 内。例 5.3.4 若在计算过程中，只

16、取3位有效数字，试用列主元素法求解：2022/7/8解 第一步，选 为主列元，将 对调位置， 第二步，选 为列主元，不需换行， 2022/7/8由 回代求解得： 与原方程组准确解 比较，可知，本题用3位有效数字计算的列元素法是相当精确的。大量实践表明：列主元法为解线性方程组的精确方法。下面用矩阵运算来描述列主元素法记 是初等排列阵（由交换单位矩阵 的第 行与 行所得）。则列主元素法为：其中 的元素满足 ，由 得： 2022/7/8 简记为： ，其中 。 下面考察 时的 其中， 则由排列阵性质（左乘矩阵是对矩阵进行行变化。）2022/7/8 已知 为单位下三角阵，其元素绝对值 ，记 。由 得：

17、。其中 为排列阵， 为单位下三角阵， 为上三角阵。这表明，对 应用列主元素法相当于对 先进行一系列行变换后对 再应用 消去法。再实际计算中我们只能在计算过程中，做关于行的变换。有结论：定理5.3.4 （列主元素三角分解定理） 若 为非奇异性矩阵，则存在排列矩阵 使 。其中 为单位下三角阵， 为上三角阵。 存放在 的下三角部分， 存放在 的上三角部分。由整数型数组记录可知 的情况。2022/7/8 5.3.7 Gauss-Jordan 消去法 消去法总是消去对角线下方的元素。现考虑一种修正，即消去对角线下方和上方的元素。这即为 消去法。 设用 消去法已完成 步，于是 化为等价方程组 其中： 20

18、22/7/8在第 步计算时，考虑对上述矩阵的第 行上、下都进行消元计算1、按列选主元素，即定义 使2、换行（当ikk）。交换A,b第k行与第ik行元素。3、计算乘数 ( 可存放在 的单元中) （ ） 4、消元计算计算主元素 2022/7/85、计算主元素 上述过程全部执行完后有： 这表明 用方法将 约化为单位矩阵，计算解就在常数项位置得到，因此用不着回代求解。用 方法接方程组的计算量大约需要 次乘除法，要比 消去法大些。但用 方法求一个矩阵的逆矩阵还是比较合适的。 2022/7/8定理 （G-J法求逆矩阵） 设 的逆矩阵 ，即求 阶矩阵 ，使 其中 为单位矩阵，将 按列写成 ， 为列向量， 为

19、单位列向量。于是求解 ，等价于求解 个方程组： ，所以我们可以用 法求解 。例 5.3.5 对 求 。2022/7/8解 2022/7/8 故： 2022/7/8 5.3.8 直接三角分解为求解 将 进行 分解。即 将原问题转化为求解两个三角形方程组 求 求 。 不选主元的三角分解法设 且有 其中 为单位下三角阵， 上三角阵。即 2022/7/8下面说明： 的元素可由 步直接计算出来，其中第 步定出 的第 行和 的第 列元素。由 有： 得 的第1行元素 得 的第1列元素到第 列元素。由 利用矩阵乘法有： 故 2022/7/8 又由 有： 故： 因此可得 的第 行和第 列的全部元素。2022/7

20、/8 直接分解法约需 乘除法，和 消去法计算量基本相当。对计算 和式，可采用双精度累加以提高精度。（B）选主元的三角分解法 在直接三角分解中，如果 计算将要中断，或者当 绝对值很小时，按分解公式计算可能引起舍入误差的积累。但当 为非奇异时，可通过交换 的行实现矩阵 的 分解。因此可采用与列主之消去法类似的方法将直接三角分解法修改为部分选主之的三角分解法。 设已完成第 步分解，这时有： 2022/7/8 2022/7/8第 步分解需要到 两式，为避免 中 用小的数 作除数，先引进量： 由 及 定义，易见 则由 有： 若 ，则我们可以用 作为 交换 的第 行与 行（但我们将交换后的新元素仍记为 及

21、 ）于是有 。控制了误差传播，再进行第 步分解。 对一般的非奇异矩阵 求逆的方法： 则：2022/7/8 5.3.9 平方根法 利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法平方根法设 为对称阵，即 ，且 的所有顺序主子式均非零，则知 可以唯一分解为：为利用 的非奇异性，将 再分解为：2022/7/8 为对角矩阵， 为单位上三角阵，则： 又 ，由分解的唯一性即得： 代入上面 中得： 。 定理 5.3.6 （对称阵的三角分解）设 为 阶对称阵，且 的所有顺序主子式均非零。则 可唯一分解为： 。其中 为单位下三角阵， 为对角阵。 当 为对称正定矩阵时，则 的所有顺序主子式 ，而

22、 设 ，则 于是：2022/7/8从而： 其中 为下三角阵。2022/7/8定理 5.3.7 （对称正定矩阵的 分解） 若 为 阶对称正定矩阵，则存在一个实的非奇异下三角矩阵 ，使 。当限定 的对角元素为正时，这种分解是唯一的。 下面来考虑计算 元素的公式，由 由矩阵乘法及 （ 时） 2022/7/8 有： 。于是得解对称正定方程组 的平方根法计算公式： 求解 即求解两个三角形方程组:(1) (2)2022/7/8由 知 ，则 ，从而 。 这表明分解过程中元素 的数量级不会增长太快且对角元素 恒为整数。于是不选主元素的平方根是一个数值稳定的方法。 当求出 的第 列元素时， 的第 行亦得出，所以

23、平方根法大约需要 次乘除法，约为一般直接 分解法计算量的一半。由于 为对称阵，因此在计算机中只需存储 的下三角部分元素，共需 个元素。可用一维数组存储，即： 。2022/7/8 矩阵 的元素 一维数组的表示为： 。 的元素存放在 的相应位置。 由公式 可见，用平方根法解对称正定方程组时，计算 的元素 需要进行开方计算，为避免开方计算，可对平方根法进行改进。例 5.3.6 用平方根法求解： 精确解为 。2022/7/8(1)分解计算：故，(2)求解两个三角方程：解 得 代入 解得： 。2022/7/8 5.3.10 追赶法在一些实际问题中常有解三对角线性方程组Ax=f的问题，即：2022/7/8

24、 其中 满足条件： 对于具有条件 的方程组 ，我们介绍下面的追赶法求解。追赶法具有计算量少，方法简单，算法稳定的特点。2022/7/8定理5.3.8 设有三对角线性方程组Ax=f，且A 满足条件 ，则A为非奇异矩阵。证明 用归纳法证明。显然n=2时，有：现假设定理对n-1阶的满足条件 的三对角矩阵成立，求证对满足条件 的 阶三对角矩阵定理亦成立。由条件 ，则利用消元法的第 步有：2022/7/8 显然 ，其中 2022/7/8且有 故知 满足条件 ，利用归纳设知 ，故定理 5.3.9 设 满足条件 为三对角阵。则 的所有顺序主子式都不为零。即：证明 由于A是满足条件（2）的n阶三对角阵。因此A的任一个顺序主子式亦满足条件（2）的n阶三对角矩阵。由上一个定理即知： 根据这一结论以及三角分解定理知，这种矩阵A可进行三角分解：A=LU。在这里特别的有：