数值模拟导论-第三讲 求解线性系统的基本方法 雅克比·怀特 感谢Deepak Ramaswamy, Michal Rewienski,Karen Veroy and Jacob White 解的存在和唯一性 高斯消元法 LU分解 对角元和等比级数 逐步逼近法 适应条件 摘要 应用范围 SMA-HPC ?2003 MIT ·无电源或刚性支承 ·对角和严格对角占优矩阵 ·n×n方阵 线性方程 SMA-HPC ?2003 MIT 要求加权变量x，使得矩阵M各列的加权和等于右边的b。 线性方程 SMA-HPC ?2003 MIT 疑问解答 ·给定Mx=b －这个方程是否有解？ －解是否唯一？ ·看是否有解? 存在一组变量x1,…..xn，使得： 由此看出：只有当b在由M各列组成 的向量空间内，解才存在。 11 2 2 ... NN x MxM xM b+++ = nullnull null 11 2 2 ... NN xM xM x M b+++ = nullnull null 线性方程 SMA-HPC ?2003 MIT 疑问解答续 ·解是否唯一？ 假设存在非零的变量，使得 又如果，则M(x+y)=b。 由此可见：当且仅当矩阵M的各列向 量线性无关，方程解唯一。 11 2 2 ... NN yM yM y M b+++ = nullnull null 11 2 2 ... NN yM yM y M b+ ++ = nullnull null 1 ,, n yy… 线性方程 SMA-HPC ?2003 MIT 疑问解答续 方阵 ·给定，其中M为N*N方阵 －如果对任意给定的b，方程均有解，那么对每一 个b所对应的解都是唯一的。 因为对任意的b都有解，那么矩阵M列向量所组成的 向量空间必定是N维向量空间。又因为矩阵M有N 列，所以矩阵M的列向量必定线性无关。 各列向量相互线性不相关的方阵称之为非奇 异阵。 高斯消元基础 SMA-HPC ?2003 MIT 重要工具 用高斯消元法解线性方程 Mxb= ·一种“直接”的方法 有限步求准确解（不计 舍入误差）。 ·求解增广矩阵的精确解 ·计算所消耗的机时。 高斯消元法基础 SMA-HPC ?2003 MIT 举例说明 3*3方阵 高斯消元法基础 SMA-HPC ?2003 MIT 举例说明 解题思路 用矩阵的第一行消去第二和第三行的x 1 1 x 高斯消元法基础 SMA-HPC ?2003 MIT 举例说明 矩阵形式 乘子 对角元 高斯消元法基础 SMA-HPC ?2003 MIT 举例说明 消去第三行中的X2 乘子 对角元 高斯消元基础 SMA-HPC ?2003 MIT 举例说明 符号简化表示 高斯消元法基础 SMA-HPC ?2003 MIT 举例说明 消去第三行的x2 乘子 对角元 高斯消元法基础 SMA-HPC ?2003 MIT 举例说明 形成的三角阵 高斯消元法基础 SMA-HPC ?2003 MIT 举例说明 右边向量的变化 高斯消元法基础 SMA-HPC ?2003 MIT 举例说明 组合各部分 LU分解基础 SMA-HPC ?2003 MIT 解方程 Mxb= 第一步 第二步消元过程 解方程 Lyb= 第三步回代过程 解方程 Ux y= 回顾线性代数我们可以知道，一个矩阵A可以用高斯消去法分解成一个下三角阵和一个 上三角阵乘积的形式。高斯消去法的基本思想是用第一行消去除第一行外所有行中的 x1，然后用第二行消去除第一，二行外所有行中的x2，依此类推，从而将系统转化成为 上三角矩阵。在这里我们需要了解更多的关于高斯消去的知识。 LU分解基础 SMA-HPC ?2003 MIT 解方程 Mxb= 第一步 第二步消元过程 解方程 Lyb= 第三步回代过程 解方程 Ux y= 回顾线性代数我们可以知道，一个矩阵A可以用高斯消去法分解成一个下三角阵和一个 上三角阵乘积的形式。高斯消去法的基本思想是用第一行消去除第一行外所有行中的 x1，然后用第二行消去除第一，二行外所有行中的x2，依此类推，从而将系统转化成为 上三角矩阵。在这里我们需要了解更多的关于高斯消去的知识。 LU分解基础 SMA-HPC ?2003 MIT 解三角矩阵 矩阵特点 第一个方程只有y1一个未知量第二 个方程只有y1和y2两个未知量 LU分解基础 SMA-HPC ?2003 MIT 解三角矩阵 算法 解三角矩阵方程是直接的但又是费时的。y1可以用一个除法运算直接求得，y2可以用一 乘法，一个减法，一个除法求得。一旦求得yk－1，yk可以由k－1步乘法，k－2步加法，一 减法和一个除法求得。算法的运算数量是： N除法＋0加/减＋1加/减＋、、、、N-1加/减＋0乘法＋1乘法＋………N-1乘法 注释：求y1 求y2 求yn求y1 求y2 求yn ＝（N-1）（N-2）加/减＋（N-1）（N-2）乘＋N除 ＝N 2 步 LU分解基础 SMA-HPC ?2003 MIT 分解 图片 上图便是LU分解的图形表示。第一步，用第一个方程消去第二到第四方程中的x1。这一 过程我们用除第一行外的各行分别减去第一行乘以某个缩放因子，从而使系数a21， a31，a41变为零。再用缩放因子（又称之为乘子）代替这些零位。对于第二行，乘子是 a21/a11，因为第二行减去第一行乘以a21/a11，a21位刚好为零。由于在消去过程中 a22，a23，a24的值也会随之改变，因此我们将他们变成蓝色。同样在消去a31和a41的 过程中，a31和a41也被他们的乘子所代替。在这一过程中第三行其余的位置的值也会随 之改变，因此也将他们变为蓝色。 用同样的方法处理第二行。计算消去第三行和第四行中x2的乘子，并且用这些乘子 代替出现的零。注意在消去过程中改变的量，将他们改为绿色。最后一步，便是用第三 行消去第四行中的x3，更新第四行的各个位置的值，并且将a44变为粉红色。 我们可以看到乘子在代替矩阵中的零的位置之后，在消去过程中他们的值并没有改 变。 LU分解基础 SMA-HPC ?2003 MIT 分解 算法 for i=1到n-1 { 每一行 for j=i+1到n { 每一要消去的目标行 ji ji ii M M M = M ii 为对角元 for k=i+1到n { 对角元后面的元素 jk jk ji jk M MMM← ? 乘子 } } } SMA-HPC ?2003 MIT 分解 对角元为零 LU分解基础 第i步 第i行 第j行 如果Mii＝0怎么办？不能得到 ji ii M M 做一下简单变形（部分对角元） 0 0 ii ii If M find M j i = ≠ > 第j行和第i行交换。 LU分解基础 SMA-HPC ?2003 MIT 分解 零对角元 两重要定理 1）只有当M为非奇异矩阵时列主元法（交换行）才 能有效。 2）对严格对角占优矩阵LU分解将不会产生零对角 元。 定理：在对严格对角占优的矩阵进行高斯消元时 不会产生零对角元。 证明：1）求出第一步消元后的矩阵。 2）考察（n－1）×（n－1）的次矩阵。 仍然是完全对角占优矩阵。 第一步 第一步消元后的第二行 由此得出 LU分解基础 SMA-HPC ?2003 MIT 数值问题 小对角元 特例 我们能够解释这种现象吗？ 为了求出精确解，我们进行第一步消元： 出 第二步回代： 求出： 在浮点运算中 LU分解基础 SMA-HPC ?2003 MIT 数值问题 小对角元 浮点运算的一个性质 双精度数 主要问题： 结论： 避免大小数之间相加减 符号位 有效数位8 LU分解基础 SMA-HPC ?2003 MIT 数值问题 小对角元 回过头来看这个例子 LU分解基础 SMA-HPC ?2003 MIT 数值问题 小对角元 列主元法减小误差 如果 max ii ji ji M M > < 那么第i行和最大的 ji M 所在行交换 可以看出乘子变小了这一项被圆整 列主元法是怎么起作用的应看下面： 求得 我们知道不用列主元法时是 ，圆整为 右边的数值3被圆整时忽略掉了，而现在他还保 留着。继续回代过程： LU分解基础 SMA-HPC ?2003 MIT 数值问题 小对角元 如果在LU分解过程中，矩阵是对角占优或采 用了列主元法来减小舍入误差那么存在以下的特 点： 1）乘子总是最小的。 2）LU因子各位上的最大值将小于等于原始矩阵各 位最大值的 (1) 2 n? 倍 列主元法是怎么起作用的应看下面： 求得 我们知道不用列主元法时是 ，圆整为 右边的数值3被圆整时忽略掉了，而现在他还保 留着。继续回代过程： 逐步逼近法 SMA-HPC ?2003 MIT 实例 数据表 用一个N次的多项式来近似表示数据： () 2 01 2 ... n n ft t t tα αα α=++ ++ 多项式插值 逐步逼近法 SMA-HPC ?2003 MIT 实例 矩阵形式 逐步逼近法 SMA-HPC ?2003 MIT 实例 1 x 在我们用一个高次的多项式函数来逼近曲线t时。当数据需要一个100次的多项式 来逼近时，我们需要大量的系数来逼近高阶系统，而不是用一两个系数。这里主 要是考虑到问题的精度，不久我们将会学习到。 逐步逼近法 SMA-HPC ?2003 MIT 扰动分析 误差范数 1 x 向量的范数 二范数 2 2 1 n i i xx = = ∑ 一范数 1 1 n i i xx = = ∑ 无穷范数 max i i xx ∞ = 正方形 单位园 逐步求解法 SMA-HPC ?2003 MIT 扰动分析 误差范数 1 x 矩阵的范数 向量导致的误差范数 1 max max xx Ax A Ax x = == 由矩阵A引起误差的范数的最大值等于x被A放大的最大倍数 简单范数的计算 1 1 max n ij j i AA = == ∑ 最大列和 1 max n ij i j AA ∞ = == ∑ 最大行和 2 A = 不易计算 逼近求解法 SMA-HPC ?2003 MIT 扰动分析 1 x 扰动计算 ( )( ) null LU x M MxxMxMxMxMxbδδ δδδδ++=+++= nullnullnullnullnullnullnullnullnullnull 实际值 分解舍入误差的扰动误差 因为Mx-b=0 所以 ( ) ( ) 1 M xMxx xMMxxδ δδδ δδ ? =? + ? =? + 取范数 1 1 xMM Mxx M δ δδ ? ≤+ 其中与误差相关的量为 1 xM MM xx M δδ δ ? ≤ + nullnullnullnullnull 条件数 由线性代数学我们知道解x的变化范围与原始的变化范围之间有一个与矩阵 M有关的倍数，这个倍数是 他被习惯的称之为条件数，由于条件数是基于范数，因此不存在它的准确 值。而是奇异矩阵M的一个比值系数。 奇异矩阵超出本课程的考虑范围，如果想了解请咨询Trefethen或Bau. 逐步逼近法 SMA-HPC ?2003 MIT 扰动分析 1 x 更清晰的几何逼近法 12 M MM ?? = ?? nullnull 解Mx=b我们会得到 11 2 2 xM xM b+ = null null 正交列逼近 当用向量来逼近时很难确定M1和与之对应的M2的值 逐步逼近法 SMA-HPC ?2003 MIT 几何分析 多项是插值 1 x 多项式值得级数接近 线性 问题行的缩放比例是否减小等比级数？ 行的缩放比例是否减小了条件数？ 定理：如果使用浮点运算，当进行行缩放时在某种意义上不会减少等比级数。 没有舍入误差 条件数 解的存在和唯一性 高斯消元法 LU分解 对角元和等比级数 逐步逼近法 适应条件 总结