目录 2. 高斯消去法 2.. 顺序消去法 2..2 列主元消去法 2..3 全主元消去法 2..4 选主元消去法的应用三角形方程组和三角分解 2.2. 三角方程组的解法 Gauss 变换 Doolittle 分解选主元三角分解平方根

线性方程组的直接解法

目录 2. 高斯消去法 2.. 顺序消去法 2..2 列主元消去法 2..3 全主元消去法 2..4 选主元消去法的应用 2 2.2 三角形方程组和三角分解 2.2. 三角方程组的解法 2.2.2 Gauss 变换 2.2.3 Doolittle 分解 3 2.3 选主元三角分解 4 2.4 平方根法及改进的平方根法 2.4. 平方根法 2.4.2 改进的平方根法 5 2.5 追赶法

线性方程组的求解问题是一个古老的数学问题九章算术 : 详细记载了消元法 9 世纪初, 西方有了 Gauss 消去法求解大型线性方程组则是在 20 世纪计算机问世后才成为可能

线性方程组数值解法的分类直接法迭代法

直接法定义 : 在没有舍入误差的情况下经过有限次运算可求得精确解的方法举例 : 高斯消去法平方根法追赶法适用范围 : 低阶稠密矩阵方程组某些大型稀疏方程组 ( 如大型带状方程组 )

迭代法定义 : 采取逐次逼近的方法, 亦即从一个初始向量出发, 按照一定的计算格式, 构造一个无穷序列, 其极限才是方程组的精确解, 只经过有限次运算得不到精确解举例 : Jacobi 迭代 Gauss-Seidel 迭代超松弛迭代适用范围 : 大型稀疏方程组

2. 高斯消去法 2.. 顺序消去法 2..2 列主元消去法 2..3 全主元消去法 2..4 选主元消去法的应用 2 2.2 三角形方程组和三角分解 2.2. 三角方程组的解法 2.2.2 Gauss 变换 2.2.3 Doolittle 分解 3 2.3 选主元三角分解 4 2.4 平方根法及改进的平方根法 2.4. 平方根法 2.4.2 改进的平方根法 5 2.5 追赶法

2. 高斯消去法 2.. 顺序消去法 2..2 列主元消去法 2..3 全主元消去法 2..4 选主元消去法的应用 2 2.2 三角形方程组和三角分解 3 2.3 选主元三角分解 4 2.4 平方根法及改进的平方根法 5 2.5 追赶法

2.. 顺序消去法定义 ( 顺序消去法 ) 在逐步消元的过程中, 把系数矩阵约化成上三角矩阵, 从而将原方程组约化为容易求解的等价三角方程组, 再通过回代过程逐一求出各未知数

2.. 顺序消去法 a () x + a () 2 x 2 + a () 3 x 3 = b (), () a () 2 x + a () 22 x 2 + a () 23 x 3 = b () 2, (2) a () 3 x + a () 32 x 2 + a () 33 x 3 = b () 3. (3) a() 2 (2) + () a () ================== (3) + () a() 3 a () a (2) i j = a () i j + a () 设 a () 0 a () x + a () 2 x 2 + a () 3 x 3 = b (), j b (2) i = b () i + b () a() i a (2) 22 x 2 + a (2) 23 x 3 = b (2) 2, a (2) 32 x 2 + a (2) 33 x 3 = b (2) 3. a () a() i a (), i, j = 2, 3,, i = 2, 3.

2.. 顺序消去法 a () x + a () 2 x 2 + a () 3 x 3 = b (), (3) + (2) a (2) 22 ================== a(2) 32 a (2) 22 x 2 + a (2) 23 x 3 = b (2) 2, a (2) 32 x 2 + a (2) 33 x 3 = b (2) 3. a (3) 33 = a(2) 33 + a(2) b (3) 3 = b (2) 3 + b (2) 设 a (2) 22 0 a () x + a () 2 x 2 + a () 3 x 3 = b (), a (2) 22 x 2 + a (2) 23 x 3 = b (2) 2, 23 a(2) 32, 2 a (2) 22 a(2) 32 a (2) 22 a (3) 33 x 3 = b (3) 3..

2.. 顺序消去法一般情形考察 n 元线性方程组 A () X = b (), 其中 A () = a () a () 2 a () n a () 2 a () 22 a () 2n.. a () n a () n2 a () nn., X = x x 2. x n, b () = b () b () 2. b () n

2.. 顺序消去法若约化的主元 a (k) kk 0 (k =, 2,, n), 则经过 for k =, 2,, n- for i = k+,, n l ik = a (k) ik / a(k) kk, for j = k+,, n+ a (k+) i j = a (k) i j l ik a (k) k j end end end 顺序消元法可得 a () a () 2 a () n 0 a (2) 22 a (2) 2n.. 0 0 a (n) nn. x x 2. x n = b () b (2) 2. b (n) n

2.. 顺序消去法 x n = b (n) n / a (n) nn for i = n-, n-2,, n end x i = b(i) i a (i) i j x j j=i+ / ai ii 回代公式

2.. 顺序消去法若遇到 a (k) kk = 0, 则消去过程无法进行若 a (k) kk 不为零但很小, 尽管消去过程可以进行下去, 但用其做除数, 会引起计算结果的严重失真 0.0000x + 2x 2 = 2, x + x 2 = 3.

2..2 列主元消去法定义 ( 列主元消去法 ) 在消元过程中, 每次选主元时, 仅依次按列选取绝对值最大的元素作为主元素, 它只进行行交换, 而不产生未知数次序的调换列主元消去法能有效地避免顺序消元过程中的两个问题, 它是直接法中最常用的一种方式

2..2 列主元消去法例用列主元消去法求解 2x + x 2 + 2x 3 = 5, 5x x 2 + x 3 = 8, x 3x 2 4x 3 = 4.

2..2 列主元消去法解 : 增广矩阵为 2 2 5 5 8 3 4 4 r r 2 ======= 5 8 2 2 5 3 4 4 r 2 2 5 r ======== r 3 5 r 2.8 4.2 5.6 5 8.4.6.8 r 2 r 3 ======= 5 8 2.8 4.2 5.6.4.6.8 r 3 + 2 r 2 ======== 5 8 2.8 4.2 5.6 0.5

2..2 列主元消去法列主元消去法 for k =, 2,, n- find i k k,, n s.t. a (k) i k,k = max k i n a (k) ik ; interchage the k, i k th rows in [A (k), b (k) ] ; for i = k+,, n l ik = a (k) ik / a(k) kk ; for j = k+,, n+ a (k+) i j = a (k) i j l ik a (k) k j ; end end end

2..3 全主元消去法定义 ( 全主元消去法 ) 全主元消去法选主元的范围更大, 对于 ( ) A () b () 来说, 在整个系数矩阵中选主元, 即将绝对值最大的元素经过行列变换使其置于 a () 的位置, 然后进行消元过程得到 ( ) A (2) b (2) 接下来在该矩阵中划掉第一行第一列后剩余的 n 阶子系数矩阵中选主元, 并通过行列交换置其于 a (2) 22 的位置, 然后进行消元 ;

2..3 全主元消去法例用全主元消去法求解 2x + x 2 + 2x 3 = 5, 5x x 2 + x 3 = 8, x 3x 2 4x 3 = 4.

2..3 全主元消去法解 : 增广矩阵为 x x 2 x 3 2 2 5 5 8 3 4 4 r r 2 ======= x x 2 x 3 5 8 r 2 2 5 r 2 2 5 ======== 3 4 4 r 3 5 r 2.8 4.2 5.6 x x 2 x 3 5 8.4.6.8 x x 3 x 2 r 2 r 5 8 3 ====== 4.2 2.8 5.6 c 2 c 3.6.4.8 r 3 + 8 2 r 2 ========= x x 3 x 2 5 8 2.8 4.2 5.6 /3 /3

2..4 选主元消去法的应用求逆矩阵应用一 : 求逆矩阵 ( ) A E Gauss Jordan Elimination =================== ( ) E A

2..4 选主元消去法的应用求逆矩阵例求矩阵 3 2 A = 23 2 2 的逆矩阵.

2..4 选主元消去法的应用求逆矩阵 3 2 0 0 23 0 0 2 2 0 0 r 2 + 23 r, r 3 + 23 r =================== r ( 23) r r 2 ======= r + 478 226 r 2, r 3 + 522 226 r 2 ======================== r (2.26) r + 365 09 r 3, r 2 + 673 09 r 3 ======================== r (.09) 23 0 0 3 2 0 0 2 2 0 0 0.478 0.044 0 0.044 0 0 2.26.522 0.478 0 0.522 2.044 0 0.044 0 0.365 0.2 0.057 0 0 0.673 0.442 0.2 0 0 0.09 0.673 0.365 0 0 0.452 0.88 0.358 0 0 0.886 0.452 0.660 0 0 0.660 0.358 0.98

2..4 选主元消去法的应用求行列式应用二 : 求行列式设有矩阵 a a 2 a n a 2 a 22 a 2n A =... a n a n2 a nn 用主元消去法将其化为上三角矩阵, 并设对角元素为 b, b 22,, b nn, 故 A 的行列式为 det(a) = ( ) m b b 22 b nn, 其中 m 为所施行的行列交换的次数

2. 高斯消去法 2 2.2 三角形方程组和三角分解 2.2. 三角方程组的解法 2.2.2 Gauss 变换 2.2.3 Doolittle 分解 3 2.3 选主元三角分解 4 2.4 平方根法及改进的平方根法 5 2.5 追赶法

2.2. 三角方程组的解法下三角形方程组考察 Ly = b () 其中 b = (b,, b n ) T R n 已知,y = (y,, y n ) T R n 未知, 而且 l ii 0, i =, 2,, n. L = l l 2 l 22 l 3 l 32 l 33,...... l n l n2 l n3 l nn

2.2. 三角方程组的解法由方程组 () 的第一个方程 l y = b 得 2 由方程组 () 的第二个方程 y = b l. l 2 y + l 22 y 2 = b 2 得 y 2 = b 2 l 2 y l 22.

2.2. 三角方程组的解法 3 一般地, 若已求得 y,, y i, 则由方程组 () 的第 i 个方程 l i y + l i2 y 2 + + l i,i y i + l ii y i = b i 得 y i = b i i j= l i jy j l ii.

2.2. 三角方程组的解法解下三角形方程组 : 前代法 for j = :n- b(j) = b(j) / L(j, j) b(j+: n) = b(j+: n) - b(j)l(j+:n, j) end b(n) = b(n) / L(n, n) 该算法所需加减乘除的次数 : n (2i ) = n 2, 即该算法的运算量为 n 2. i=

2.2. 三角方程组的解法前代法 (matlab 代码 :fs.m) function b = fs(l, b, n) 2 for j = :n- 3 b(j) = b(j) / L(j,j); 4 b(j+: n) = b(j+: n) - b(j) * L(j+:n, j); 5 end 6 b(n) = b(n) / L(n,n); 7 end

2.2. 三角方程组的解法上三角形方程组考察 Ux = y (2) 其中 y = (y,, y n ) T R n 已知,x = (x,, x n ) T R n 未知, 而 U = 且 u ii 0, i =, 2,, n. u u 2 u 3 u n u 22 u 23 u n u 33 u 3n.... u nn,

2.2. 三角方程组的解法由方程组 (2) 的第 n 个方程 u nn x n = y n 得 2 由方程组 (2) 的第 n 个方程 x n = y n u nn. u n,n x n + u n,n x n = y n 得 x n = y n u n,n x n u nn.

2.2. 三角方程组的解法 3 一般地, 若已求得 x n,, x i+, 则由方程组 (2) 的第 i 个方程 u ii x i + u i,i+ x i+ + + u i,n x n = y i 得 x i = y i n j=i+ u i j x j u ii.

2.2. 三角方程组的解法解上三角形方程组 : 回代法 for j = n-:-:2 y(j) = y(j) / U(j, j) y(:j-) = y(:j-) - y(j)u(:j-, j) end y() = y() / U(, ) 该算法的运算量也为 n 2.

2.2. 三角方程组的解法回代法 (matlab 代码 :bs.m) function y = bs(u, y, n) 2 for j = n:-:2 3 y(j) = y(j) / U(j,j); 4 y(:j-) = y(:j-) - y(j) * U(:j-, j); 5 end 6 y() = y() / U(,); 7 end

2.2. 三角方程组的解法一般线性方程组考察 Ax = b (3) 其中 A R n n, x, b R n 若 A = LU, 即 A 能分解成一个下三角阵 L 和一个上三角阵 U 的乘积, 则 (3) 的解 x 可由以下两步得到 : 用前代法解 Ly = b 得 y; 2 用回代法解 Ux = y 得 x

2.2.2 Gauss 变换定义 ( 矩阵三角分解 ) 将矩阵 A 分解为一个下三角阵 L 和一个上三角阵 U 的乘积, 最自然的做法是通过一系列初等变换, 逐步将 A 约化为上三角阵, 并且保证这些初等变换的乘积是一个下三角阵

2.2.2 Gauss 变换方式一 :Gauss 变换定义 (Gauss 变换 ( 矩阵 )) L k =... l k+,k.... l n,k I l k e T k l k = (0,, 0, l k+,k,, l nk ) T Gauss 向量

2.2.2 Gauss 变换方式一 :Gauss 变换对于 x = (x,, x n ) T R n, L k x = (x,, x k, x k+ l k+,k x k,, x n l nk x k ) T. 取便有 l ik = x i x k, i = k +,, n, x k 0 L k x = (x,, x k, 0,, 0) T.

2.2.2 Gauss 变换方式一 :Gauss 变换性质 ( L k ) L k 的逆为 L k = I + l k e T k 证明. e T k l k = 0, (I + l k e T k )(I l ke T k ) = I l k e T k l k e T k = I.

2.2.2 Gauss 变换方式一 :Gauss 变换性质 (2 L k ) L k A = (I l k e T k )A = A l k(e T k A), e T k A 为 A 的第 k 行类似地,Ae k 为 A 的第 k 列例 2 3 2 3 2 3 4 4 5 7 = 2 3 2 2 5

2.2.2 Gauss 变换方式一 :Gauss 变换性质 (3 L k ) 若 j < k, 则 L j L k = (I l j e T j )(I l ke T k ) = I l je T j l k e T k. 证明. 因为当 j < k 时, 有 e T j l k = 0 例 2 3 4 = 2 3 4

2.2.2 Gauss 变换 } {{ } A L } {{ } L A L 2 L 3 } {{ } L 2 L A } {{ } L 3 L 2 L A

2.2.2 Gauss 变换 L A = A = 2 4 3 2 0 4 3 3 8 7 9 5 6 7 9 8 2 0 4 3 3 8 7 9 5 6 7 9 8 = 2 0 3 5 5 4 6 8

2.2.2 Gauss 变换 L 2 L A = L A = 3 4 2 0 3 5 5 4 6 8 2 0 3 5 5 4 6 8 = 2 0 2 2 2 4

2.2.2 Gauss 变换 L 2 L A = 2 0 2 2 2 4 L 3 L 2 L A = 2 0 2 2 2 4 = 2 0 2 2 2 = U.

2.2.2 Gauss 变换 L = 2 4 3, L2 = 2 0 4 3 3 8 7 9 5 6 7 9 8 } {{ } A L L 2 L 3 = = 3 4 2 4 3 3 4 2 4 3 3 4 } {{ } L, L3 = 2 0 2 2 2 } {{ } U

2.2.2 Gauss 变换对于一般矩阵 A, 记 A () = A, 则 A (2) = L A () A (3) = L 2 A (2) = L 2 L A (). A (k) = L k L A () = A (k) A (k) 2 A (k) 22 其中 A (k) 是 k 阶上三角阵,A (k) 22 为 A (k) 22 = a (k) kk a (k) kk.. a (k) nk a (k) nn

2.2.2 Gauss 变换若 a (k) kk 0, 则可确定一个 Gauss 变换 L k, 使得 L k A (k) 第 k 列的最后 n k 个元素为 0 取 L k = I l k e T k, 则其中 A (k+) 是 k 阶上三角阵 l k = (0,, 0, l k+,k,, l nk ) T, l ik = a (k) ik /a(k) kk, i = k +,, n A (k+) = L k A (k) = A (k+) A (k+) 2 A (k+) 22 如此进行 n 步, 最终所得矩阵 A (n) 即为所要求的上三角形式

2.2.2 Gauss 变换上述步骤可描述为 L n L A = A (n). 令 L = (L n L ) 及 U = A (n), 则 A = LU, 其中 L 是一个单位上三角阵事实上, 由于 e T j l i = 0 ( j < i), 则 L = L L 2 L n = (I + l e T )(I + l 2e T 2 ) (I + l n e T n ) = I + l e T + l 2e T 2 + l n e T n.

2.2.2 Gauss 变换 A (k+) = L k A (k) = (I l k e T k )A(k) = A (k) l k e T k A(k). A (k+) 与 A (k) 的前 k 行元素相同因为 e T k A(k) 是 A (k) 的第 k 行,l k 的前 k 个分量为 0 A (k+) 中第 k + 到 n 行元素需要更新 a (k+) i j a (k+) ik = 0, i = k +,, n, = a (k) ik l ika (k) k j, i, j = k +,, n

2.2.2 Gauss 变换存储方式 : A (k) 中第 k + n 行元素在计算出 A (k+) 以后不再有用, 故 A (k) 中相应位置上的元素可用新值更新由于 A (k+) 中的第 k 列对角线以下的元素均为 0, 无需存储, 故这些位置可用于存储 l k 中的非 0 元

2.2.2 Gauss 变换算法 ( 计算三角分解 :Gauss 消去法 ) for k = :n- A(k+:n, k) = A(k+:n, k) / A(k, k) A(k+:n, k+: n) = A(k+:n, k+: n) - A(k+:n, k) A(k, k+: n) end 运算次数 : n [(n k) + 2(n k) 2 ] = k= n(n ) 2 + = 2 3 n3 + O(n 2 ) n(n )(2n ) 3

2.2.2 Gauss 变换 LU 分解的 matlab 代码 function A = GsLU(A) [m n] = size(a); for k = :n- A(k+:n, k) = A(k+:n, k) / A(k, k); A(k+:n, k+: n) = A(k+:n, k+: n) - A(k+:n, k) * A(k, k+: n); end end

2.2.2 Gauss 变换 a (k) kk (k =,, n ) 称作主元三角分解的条件当且仅当主元 a (k) kk (k =,, n ) 均不为 0 时, 算法才能进行到底

2.2.2 Gauss 变换定理主元 a (i) ii (i =,, k) 均不为 0 A 的 i 阶顺序主子阵 A i (i =,, k) 非奇异证明当 k = 时,a () 0 A 非奇异 2 假设定理直到 k 成立, 下证 : 若 A,, A k 非奇异, 则 A k 非奇异 a (k) kk 0

2.2.2 Gauss 变换证明由归纳假设,a (i) ii 0(i =,, k ), 故 A (k) = L k L A = A (k) A (k) 2 A (k) 22 其中 A (k) 是 k 阶上三角阵, 而 A (k) 的 k 阶顺序主子阵形如 A (k) a (k) kk

2.2.2 Gauss 变换证明. 记 L,, L k 的 k 阶顺序主子阵为 (L ) k,, (L k ) k, 则 A (k) (L k ) k (L ) k A k = a (k) kk = det A k = a (k) kk det A(k) = A k 非奇异当且仅当 a (k) kk 0

2.2.2 Gauss 变换定理 ( 矩阵三角分解的条件 ) A 的各阶顺序主子阵均非奇异 = 存在唯一的单位下三角阵 L 和上三角阵 U, 使得 A = LU.

2.2.3 Doolittle 分解方式二 :Doolittle 分解 a a 2 a n a 2 a 22 a 2n...... a n a n2 a nn 图 : Doolittle 分解 l 2 =..... l n l n2 u u 2 u n u 22 u 2n.... u nn

2.2.3 Doolittle 分解方式二 :Doolittle 分解图 : Doolittle 分解运算次序 : 先行后列, 先 U 后 L u u 2 u 3 u,n u n l 2 u 22 u 23 u 2,n u 2n l 3 l 32 u 33 u 3,n u 2n.......... l n, l n,2 l n,3 u n,n u n,n l n l n2 l n3 l n,n u nn

2.2.3 Doolittle 分解方式三 :Crout 分解 a a 2 a n a 2 a 22 a 2n...... a n a n2 a nn 图 : Crout 分解 l l 2 l 22 =..... l n l n2 l nn u 2 u n u 2n....

2.2.3 Doolittle 分解方式三 :Crout 分解图 : Crout 分解运算次序 : 先列后行, 先 L 后 U l u 2 u 3 u,n u n l 2 l 22 u 23 u 2,n u 2n l 3 l 32 l 33 u 3,n u 2n.......... l n, l n,2 l n,3 l n,n u n,n l n l n2 l n3 l n,n l nn

2.2.3 Doolittle 分解方式二 :Doolittle 分解 n a i j = l ik u k j i, j =, 2,, n k= l ii =, l i j = 0( j > i), u i j = 0(i > j) = 先算 U 的第 k 行 ( j k) n k k a k j = l kr u r j = l kr u r j = l kr u r j + l kk u k j r= r= r= k = u k j = a k j l kr u r j r= 再算 L 的第 k 列 (i > k) n k k a ik = l ir u rk = l ir u rk = l ir u rk + l ik u kk r= r= r= = l ik = a ik k r= l kru r j u kk

2.2.3 Doolittle 分解方式二 :Doolittle 分解先行后列 for k = : n for j = k,, n % 计算第 k 行 u k j = a k j k r= l lru r j end for i = k+,, n % 计算第 k 列 l ik = (a ik k r= l iru rk )/u kk end end

例利用 Doolittle 分解求解线性方程组 2 3 4 x 2 3 4 2 3 2 0 0 3 4 4 9 3 x 2 x 3 x 4 = 5 0 7

计算 U 的第一行,L 的第一列, 得 u =, u 2 = 2, u 3 = 3, u 4 = 4, l 2 = a 2 u = 3, l 3 = a 3 u = 2, l 3 = a 3 u = 4.

2 计算 U 的第二行,L 的第二列, 得 u 22 = a 22 l 2 u 2 = 2, u 23 = a 23 l 2 u 3 = 3, u 24 = a 24 l 2 u 4 =, l 32 = a 32 l 3 u 2 u 22 = 3, l 42 = a 42 l 4 u 2 u 22 = 3.

3 计算 U 的第三行,L 的第三列, 得 u 33 = a 33 l 3 u 3 l 32 u 23 = 3, u 34 = a 34 l 3 u 4 l 32 u 24 = 2, l 43 = a 43 l 4 u 3 l 42 u 23 u 33 = 2. 4 计算 U 的第四行, 得 u 44 = a 44 l 4 u 4 l 42 u 24 l 43 u 34 = 4,

2.2.3 Doolittle 分解方式二 :Doolittle 分解 2 3 4 2 3 4 3 4 2 3 2 0 0 3 4 4 9 3 = 3 2 3 4 3 2 2 3 3 2 4

2.2.3 Doolittle 分解方式二 :Doolittle 分解 2 3 4 2 3 4 3 4 2 3 2 0 0 3 4 4 9 3 = 3 2 3 4 3 2 2 3 3 2 4 y 2 3 2 3 4 3 2 y 2 y 3 y 4 = 5 0 7 求得 Y = ( 2,, 7, 6) T.

2.2.3 Doolittle 分解方式二 :Doolittle 分解 2 3 4 2 3 4 3 4 2 3 2 0 0 3 4 4 9 3 = 3 2 3 4 3 2 2 3 3 2 4 2 3 4 x 2 2 3 3 2 4 x 2 x 3 x 4 = 7 6 求得 X = (, 2, 3, 4) T.

2.2.3 Doolittle 分解对称矩阵的三角分解定理若 A 为 n 阶对称矩阵, 且 A 的各阶顺序主子式都不为 0, 则 A 可惟一分解为 A = LDL T, 其中 L 为单位下三角阵,D 为对角阵证明 : 由 A = l 2..... l n l n2 u u 2 u n u 22 u 2n.... u nn

2.2.3 Doolittle 分解对称矩阵的三角分解因为 u ii 0(i =,, n), 故 U 可分解为 u u 22 U =... u 2 u u n u u 2n u 22.... = DU u nn 其中 D 为对角矩阵,U 为单位上三角阵于是 A = LDU = L(DU ), 因为 A 为对称阵, 故 A = A T = U T DT L = U T (DLT ). 由 A 的 LU 分解的惟一性即得 L = U T.

2.3 选主元三角分解例求的 LU 分解 A = [ 0 ] A 非奇异, 但 Gauss 消去第一步就会失效

2.3 选主元三角分解例求的 LU 分解 A = [ 0 20 ] [ L = 设 ɛ machine 0 6, 则 [ L = 0 0 20 0 0 20 LŨ = ] [ ] 0 20, U = 0 0 20 ] [ ] 0 20, Ũ = 0 0 20 [ 0 20 0 ]

2.3 选主元三角分解例设求 A = [ 0 20 Ax = b ] [, b = 0 ], 近似解 : [ 0 x = ] 精确解 : x [ ]

2.3 选主元三角分解交换方程次序 [ ] [ ] [ x2 0 0 20 = x 近似解 ˆL = [ 0 ], Û = ˆx [ [ ] 0 + 0 20 ] ]

2.3 选主元三角分解图 : 置换矩阵 I pq =... 0.. 0... p q = p q ( ) e e p e q e p+ e q e p e q+ e n

2.3 选主元三角分解具体步骤假定消去过程已经进行到了 k 步, 即已经确定了 k 个 Gauss 变换 L,, L k R n n 和 2(k ) 个初等置换矩阵 P,, P k, Q,, Q k R n n 使得 A (k) = L k P k L P A Q Q k A (k) A (k) 2 = 0 A (k) 22

2.3 选主元三角分解具体步骤 2 第 k 步 : 选 a (k) pq = max { a (k) : k i, j n} 若 a (k) pq = 0, 则 R(A) = k, 停止! i j 若 a (k) pq 0, 交换 A (k) 的第 k p 行及第 k q 列, 记交换后的 A (k) 22 为 ã (k) kk ã (k) kn Ã (k) 22 =. ã (k) nk ã (k) nn.

2.3 选主元三角分解具体步骤 2 第 k 步 : 选 a (k) pq = max { a (k) : k i, j n} 计算 Gauss 变换 i j L k = I l k e T k l k = (0,, 0, l k+,k,, l n,k ) T l i,k = ã (k) ik /ã(k) kk, i = k +,, n. 则 A (k+) = L k P k A (k ) Q k A (k+) A (k+) 2 = 0 A (k+) 22 其中 A (k+) 为 k 阶上三角阵,P k = I kp, Q k = I kq R n n

2.3 选主元三角分解全主元 Gauss 消去法的矩阵表述存在置换矩阵 P k, Q k 和初等下三角阵 L k, k =,, n, 使得为上三角阵令 L n P n L P A Q Q r = U Q = Q Q n P = P n P L = P(L n P n L P ) 则有 PAQ = LU 以下说明 L 为一个单位下三角阵

2.3 选主元三角分解 } {{ } P =I 4 A = A () = 3 0 3 9 2 3 2 3 2 2 3 3 3 0 3 9 2 3 2 3 2 2 3 3 } {{ } A (), } {{ } Q =I 3 = 3 3 2 3 2 9 3 2 2 3 0 3 } {{ } Ã () 0.2308 0.0769 0.2308 } {{ } L 3 3 2 3 2 9 3 2 2 3 0 3 } {{ } Ã () = 3 3 2.3077 8.5385 0.7692 2.7692.8462.923 9.3077 2.5385 0.7692 } {{ } A (2)

2.3 选主元三角分解 } {{ } P 2 =I 23 A (2) = 3 3 2.3077 8.5385 0.7692 2.7692.8462.923 9.3077 2.5385 0.7692 3 3 2.3 8.54 0.77 2.77.85.92 9.3 2.54 0.77 } {{ } A (2) } {{ } Q 2 =I 24, = 3 2 3.92.85 2.77 0.77 8.54.3 0.77 2.54 9.3 } {{ } Ã (2) 0.0645 0.0645 } {{ } L 2 3 2 3.92.85 2.77 0.77 8.54.3 0.77 2.54 9.3 } {{ } Ã (2) = 3 2.92 3.85 5.77.92.85 2.77 8.42.3 2.42 9.3 } {{ } A (3)

2.3 选主元三角分解 A (3) = 3 2 3.92.85 2.77 8.42.3 2.42 9.3, } {{ } P 3 =I 34 3 2 3.92.85 2.77 8.42.3 2.42 9.3 } {{ } A (3) 0.2 } {{ } L 3 3 2.92 5.77 3.85.92 2.77.85 9.3 2.42.3 8.42 } {{ } Q 3 =I 34 } {{ } Ã (3) = = 3 2.92 5.77 3.85.92 2.77.85 9.3 2.42.3 8.42 } {{ } Ã (3) 3 2 3.92 2.77.85 9.3 2.42 8.2 } {{ } A (4)

2.3 选主元三角分解将上述步骤合并, 即得 L 3 P 3 L 2 P 2 L P A Q Q 2 Q 3 = A (4) = A Q = (L 3 P 3 L 2 P 2 L P ) A (4) = P L P 2 L 2 P 3 L 3 A (4) = P 3 P 2 P } {{ } P A Q }{{} Q = P 3 P 2 L P 2 L2 P 3 L3 } {{ } L A (4) }{{} U

2.3 选主元三角分解观察 L = P 3 (P 2 L P 2 L 2 )P 3L 3, 设 P 2 = I 23, P 3 = I 34 以及 L = L P 3 = P 2 = 2 3 4 3 2 4 3 4 @ 4 2 @ 3, L2 = P 2 = P 3 = 由此可以看出 L 是一个单位下三角阵 @ 3 @ 4 3 2 4 3 4 @ 4 2 @ 3 L 3 = # 4 L2 === 3 2 @ 3 4 @ 4 L3 === 3 2 @ 3 4 @ 4 # 4

2.3 选主元三角分解定理 L 是一个单位下三角阵 = P n P (L n P n L P ) = P n P 2 L P 2L2 P nln

2.3 选主元三角分解定理存在置换矩阵 P, Q 以及单位下三角阵 L 和上三角阵 U, 使得 PAQ = LU, 且 l i j,u 的非零对角元的个数正好等于 A 的秩

2.3 选主元三角分解算法.2.( 全主元三角分解 ) for k = :n- 确定 p, q(k<=p, q<=n) 使得 A(p,q) = max{ A(i,j) : i=k:n, j=k:n} A(k,:n) <-> A(p,:n) A(:n,k) <-> A(:n,q) u(k)=p v(k)=q if A(k,k)!=0 A(k+:n,k) = A(k+:n,k)/ A(k,k) A(k+:n,k+: n) = A(k+:n,k+: n) - A(k+:n,k)* A(k,k+: n) else stop end end

2.3 选主元三角分解例对如下矩阵进行全主元三角分解 2 A = 4 3 3 8 7 9 5 6 7 9 8,

2.3 选主元三角分解 2 4 3 3 A =, 8 7 9 5 6 7 9 8 2 8 7 9 5 4 3 3 4 3 3 = 8 7 9 5 2 6 7 9 8 6 7 9 8 } {{ }} {{ }} {{ } P A (0) Ã (0) 8 7 9 5 8 7 9 5 4 3 3 3 3 2 2 2 2 = 2 0 3 5 5 4 4 4 4 3 7 9 7 6 7 9 8 4 4 4 4 } {{ }} {{ }} {{ } L Ã A ()

2.3 选主元三角分解 2 4 3 3 A =, 8 7 9 5 6 7 9 8 2 8 7 9 5 4 3 3 4 3 3 = 8 7 9 5 2 6 7 9 8 6 7 9 8 } {{ }} {{ }} {{ } P A (0) Ã (0) 8 7 9 5 8 7 9 5 4 3 3 3 3 2 2 2 2 = 2 0 3 5 5 4 4 4 4 3 7 9 7 6 7 9 8 4 4 4 4 } {{ }} {{ }} {{ }

2.3 选主元三角分解 8 7 9 5 3 3 A () 2 2 2 =, 3 5 5 4 4 4 7 9 7 4 4 4 8 7 9 5 8 7 9 5 3 3 7 9 7 2 2 2 4 4 4 = 3 5 5 3 5 5 4 4 4 4 4 4 7 9 7 3 3 4 4 4 2 2 2 } {{ }} {{ }} {{ } P 2 A () Ã () 3 7 2 7 } {{ } L 2 8 7 9 5 7 4 9 4 7 4 = 3 5 5 2 4 4 4 4 7 7 3 3 6 2 2 2 2 7 7 } {{ }} {{ } 8 7 9 5 7 4 9 4 7 4

2.3 选主元三角分解 8 7 9 5 3 3 A () 2 2 2 =, 3 5 5 4 4 4 7 9 7 4 4 4 8 7 9 5 8 7 9 5 3 3 7 9 7 2 2 2 4 4 4 = 3 5 5 3 5 5 4 4 4 4 4 4 7 9 7 3 3 4 4 4 2 2 2 } {{ }} {{ }} {{ } P 2 A () Ã () 3 7 2 7 } {{ } L 2 8 7 9 5 7 4 9 4 7 4 = 3 5 5 4 4 4 3 3 2 2 2 } {{ } Ã () 8 7 9 5 7 4 9 4 2 7 7 4 4 7 6 2 7 7 } {{ } A (2)

2.3 选主元三角分解 0 } {{ } P 3 3 7 2 7 } {{ } L 8 7 9 5 7 4 9 4 A (2) =, 2 4 7 7 6 2 7 7 8 7 9 5 8 7 9 5 7 4 9 4 2 7 7 4 4 7 6 2 7 7 } {{ } A (2) 8 7 9 5 7 4 9 4 7 4 6 7 2 7 2 7 4 7 } {{ } Ã (2) 7 4 = = 7 4 9 4 7 4 6 7 2 7 2 7 4 7 } {{ } Ã (2) 8 7 9 5 7 4 9 4 7 4 6 7 2 7 2 } {{ } A () 3

2.3 选主元三角分解算法.2.( 列主元三角分解 ) for k = :n- 确定 p(k<=p<=n) 使得 A(p,k) = max{ A(i,j) : i=k:n, j=k:n} A(k,:n) <-> A(p,:n) u(k)=p if A(k,k)!=0 A(k+:n,k) = A(k+:n,k)/ A(k,k) A(k+:n,k+: n) = A(k+:n,k+: n) - A(k+:n,k)* A(k,k+: n) else stop end end

2. 高斯消去法 2 2.2 三角形方程组和三角分解 3 2.3 选主元三角分解 4 2.4 平方根法及改进的平方根法 2.4. 平方根法 2.4.2 改进的平方根法 5 2.5 追赶法

2.4. 平方根法对称正定矩阵适用对象 : 对称正定矩阵方程组定义 ( 对称正定矩阵 ) 设 A 是 n 阶实对称矩阵, 若 0 x R n, 恒有 x T Ax > 0, 则称 A 为对称正定矩阵

2.4. 平方根法对称正定矩阵的性质性质若 A 对称正定, 则 A 非奇异 2 任一主子矩阵 A k, k =, 2,, n 必正定 3 a ii > 0, i =, 2,, n 4 λ i (A) > 0, i =, 2,, n 5 det(a) > 0

2.4. 平方根法定理对称矩阵 A 正定 A 的各阶顺序主子式 det (A i ) > 0, i =, 2,, n.

2.4. 平方根法定理 (Cholesky 分解 ) 对称矩阵 A 正定 = 存在惟一的主对角元皆正的下三角阵 L, 使得 A = LL T.

2.4. 平方根法证明. A 对称正定 A 各阶顺序主子式 > 0 = A = LD L T, L : 单位下三角阵 L T 非奇异 = 存在惟一的向量 x R n 使得 L T x = e i 于是所以 0 < x T Ax = x T LD L T x = ( L T x) T D L T x = e T i De i = d i. A = LD L T = LD /2 D /2 L T = (D /2 L) T D /2 L LL T

2.4. 平方根法 a a 2 a n a 2 a 22 a 2n =... a n a n2 a nn l l 2 l 22..... l n l n2 l nn l l 2 l n l 22 l n2.... l nn

2.4. 平方根法方式二 :Doolittle 分解图 : 平方根法运算次序 l l 2 l 22 l 3 l 32 l 33........ l n, l n,2 l n,n l n l n2 l n,n l nn

2.4. 平方根法方式二 :Doolittle 分解由矩阵乘法 a i j = n l ik u k j = k= 注意到 l i j = 0( j > i), 知计算第 j 行时当 i = j 时, k= a j j = n l ik l jk (i, j =, 2,, n), k= j j l 2 jk = l j j = a j j l jk l jk k= 当 i > j 时, j j a i j = l ik l jk = l i j l j j = a i j l ik l jk = l i j = a i j j k= l ikl jk l j j k= k= /2

2.4. 平方根法方式二 :Doolittle 分解 for j =, 2,, n j /2 l j j = a j j l 2 jk k= for i = j+, j+2,, n j l i j = a i j l ik l jk /l j j end end k= 平方根法

2.4. 平方根法例用平方根法求解 2 2 8 4 4 6 x x 2 x 3 = 0 2 3, L = l l 2 l 22 l 3 l 32 l 33

2.4. 平方根法例用平方根法求解 2 2 8 4 4 6 x x 2 x 3 = 0 2 3, L = l l 2 l 22 l 3 l 32 l 33 解验证 A 的对称正定性 : a = > 0, 2 2 8 = 8 4 > 0, 2 2 8 4 4 6 = 6 > 0

2.4. 平方根法例用平方根法求解 2 2 8 4 4 6 解 x x 2 x 3 = 0 2 3, L = l l 2 l 22 l 3 l 32 l 33 l = a =, l 2 = a 2 = 2 l l 3 = a 3 = l l 22 = a 22 l 2 2 = 2 l 32 = a 32 l 3 l 2 l 22 = l 33 = a 33 l 2 3 l2 32 = 2

2.4. 平方根法例用平方根法求解解 2 2 8 4 4 6 x x 2 x 3 = 0 2 3, L = 2 2 2 求解 LY = b, 得 3 求解 L T X = Y, 得 Y = (0,, 2) T. X = (,, ) T.

2. 高斯消去法 2 2.2 三角形方程组和三角分解 3 2.3 选主元三角分解 4 2.4 平方根法及改进的平方根法 2.4. 平方根法 2.4.2 改进的平方根法 5 2.5 追赶法

2.4.2 改进的平方根法平方根法的局限计算 l i j 时需要用到开方运算只能求解对称正定线性方程组而在很多工程问题中, 经常得到的是一个系数矩阵对称但不一定正定的线性方程组为了避免开方运算和求解这类方程组, 可采用改进平方根法

2.4.2 改进的平方根法 A = LDL T = l 2.... l n l n2 d d 2... d n l 2 l n l n2....

2.4.2 改进的平方根法存储方式及运算次序 a a 2 a 22 a 3 a 32 a 33........ a n, a n,2 a n,n a n a n2 a n,n a nn d l 2 d 2 l 3 l 32 d 3........ l n, l n,2 d n l n l n2 l n,n d n

2.4.2 改进的平方根法由矩阵乘法 a i j = 注意到 l i j = 0( j > i), 知计算第 j 列时当 i = j 时, n l ik d k l jk k= j j j a j j = d k l 2 jk = d k l 2 jk + d j = d j = a j j d k l 2 jk k= k= k= 当 i > j 时, a i j = j j l ik d k l jk = l i j d j l j j = a i j l ik d k l jk = l i j = a i j j k= l ikd k l jk d j k= k=

2.4.2 改进的平方根法改进平方根法 for j =, 2,, n j d j = a j j l 2 jk d k end end k= for i = j+, j+2,, n j l i j = a i j l ik d k l jk /d j k=

2.4.2 改进的平方根法等价于 find Y s.t. LY = b; find X s.t. DL T X = Y. AX = b 即 for i =, 2,, n i y i = b i l ik y k k= end for i = n, n-,, x i = y n i l ki x k d i k=i+ end

2.4.2 改进的平方根法例用改进平方根法求解 5 4 0 4 6 4 4 6 4 0 4 5 x x 2 x 3 x 4 = 2 2

2.4.2 改进的平方根法例用改进平方根法求解 d = a = 5 5 4 0 4 6 4 4 6 4 0 4 5 x x 2 x 3 x 4 = 2 2 l 2 = a 2 d = 4/5 l 3 = a 3 d = /5 l 4 = a 4 d = 0 d 2 = a 22 d l 2 2 = 2.8 l 32 = a 32 l 3 d l 2 d 2 =.4286 l 42 = a 4 l 4 d l 2 d 2 =

2.4.2 改进的平方根法例用改进平方根法求解 5 4 0 4 6 4 4 6 4 0 4 5 x x 2 x 3 x 4 = 2 2 d 3 = a 33 d l 2 3 d 2l 2 32 = 2.485 l 43 = a 43 l 4 d l 3 l 42 d 2 l 32 d 3 = 4 d 4 = a 44 d l 2 4 d 2l 2 42 d 3l 2 43 = 5

2.4.2 改进的平方根法 5 求 LY = b 得 Y = (2, 0.6, 0.7428, 0.83333) T 6 求 DL T X = Y 得 X = (.00002,.00003,.00004,.00002) T

2.5 追赶法适用范围 : 三对角线性方程组 b c a 2 b 2 c 2 a 3 b 3 c 3......... a n b n c n a n b n x x 2 x 3. x n x n = d d 2 d 3. d n d n 系数矩阵 A 是三对角矩阵, 它常常是按行严格对角占优的, 即 b > c > 0, b i a i + c i, a i 0, c i 0, i = 2, 3,, n, b n > a n > 0.

2.5 追赶法 b c a 2 b 2 c 2 a 3 b 3 c 3 = l 2 l 3 u c u 2 c 2 u 3 c 3.................. a n b n l n u n () u = b

2.5 追赶法 b c a 2 b 2 c 2 a 3 b 3 c 3 = l 2 l 3 u c u 2 c 2 u 3 c 3.................. a n b n l n u n (2) l 2 u = a 2 = l 2 = a 2 u l 2 c + u 2 = b 2 = u 2 = b 2 l 2 c

2.5 追赶法 b c a 2 b 2 c 2 a 3 b 3 c 3 = l 2 l 3 u c u 2 c 2 u 3 c 3.................. a n b n l n u n (3) l 3 u 2 = a 3 = l 3 = a 3 u 2 l 3 c 2 + u 3 = b 3 = u 3 = b 3 l 3 c 2

2.5 追赶法 b c u c a 2 b 2 c 2 l 2 u 2 c 2 a 3 b 3 c 3 = l 3 u 3 c 3.................. a n b n l n u n u = b, 递推关系 l i = a i, u i i = 2, 3,, n u i = b i l i c i

2.5 追赶法例求解微分方程 u xx = f, u(a) = u(b) = 0. x (a, b), 解网格剖分 : a = x 0 < x < < x n = b, h = b a n 2 差分离散 : u xx (x i ) u i+ 2u i + u i h 2, i =, 2,, n

2.5 追赶法解 3 生成矩阵 : 2 2......... 2 2 u u 2. u n 2 u n = f f 2. f n 2 f n

2.5 追赶法解 4 用追赶法求解线性方程组 2 2 2 2 2 /2 = 2/3 3/4 4/5 2 3/2 4/3 5/4 6/5

2.5 追赶法存储方式系数矩阵与右端项的存储用四个 n 维向量 a, b, c, d 分别来存储三条对角线上的元素及右端项的值 2 l 与 u 的存储 l 的各元素存储在 a 对应的元素位置,u 的各元素存储在 b 对应的元素位置上 3 未知量 x 的存储 x 的各元素存储在 d 对应的元素位置

目录 2. 高斯消去法 2.. 顺序消去法 2..2 列主元消去法 2..3 全主元消去法 2..4 选主元消去法的应用 三角形方程组和三角分解 2.2. 三角方程组的解法 Gauss 变换 Doolittle 分解 选主元三角分解 平方根

目录 2. 高斯消去法 2.. 顺序消去法 2..2 列主元消去法 2..3 全主元消去法 2..4 选主元消去法的应用三角形方程组和三角分解 2.2. 三角方程组的解法 Gauss 变换 Doolittle 分解选主元三角分解平方根