Ch3 线性方程组求解的数值方法

1. Ch3 线性方程组求解的数值方法 Numerical Solutions of Systems of Linear Equations

2. §3.0 Introduction 工程问题（电子网络，船体放样等） 自然科学（实验数据的最小二乘法曲线拟合） 线性方程组 解非线性方程组 解常/偏微分方程组（差分法，有限元法）

3. 3.0 Introduction 线性方程组的系数矩阵： （1）低价稠密矩阵（阶数<150）； （2）大型稀疏矩阵（阶数高且零元素较多）。 求解线性方程组的方法： （1）直接法: 经过有限步算术运算，求得精确解（假设计算过程没有舍入误差）。 如Gauss消去法，三角分解法。 （2）间接法(迭代法)： 通过迭代序列，逐步逼近方程组的解。如Jacobi迭代法，Gauss-Seidel迭代法。

4. 3.0 Introduction 【最简单的情形】三角形线性方程组： 简记为 ，A为上三角矩阵。 若所有aii0，可用“回代”过程得到方程组的解。

5. 3.0 Introduction function X=backsub(U,b) % 解上三角方程组--回代过程%Backsubstitution in Gauss elimination % Input--U is a n x n upper-trianglular matrix % --b is a n x 1 constant vctor % Output--X is the solution vector of UX=b % Usage: X=backsub(U,b) %Find the dimension of b and initialize X n=length(b); X=zeros(n,1); X(n)=b(n)/U(n,n); for p=n-1:-1:1 X(p)=(b(p)-U(p,p+1:n)*X(p+1:n))/U(p,p); end Backsub.m

6. 3.0 Introduction 【作业】（1）编写Matlab程序解下列上三角形方程组： （2）编写Matlab程序解下列下三角形方程组：

7. §3.1 Gauss消去法与LU分解法 直接法的基本思想： 将线性方程组化成与之等价的上三角形或下三角形，再用回代法求解。它的核心是矩阵分解。 核心：矩阵分解。

8. Gauss消去法（Gaussian elimination）： （对方程组的三种变换） （1）交换两个方程的次序； （2）用一个非零常数乘一个方程； （3）将一个方程的非零倍数加到另一个方程上去。 这等价与对增广矩阵进行三种行初等变换，将它化为行阶梯形。

9. §3.1 Gaussian Elimination and LU Decomposition 【例3.1】解下列线性方程组： 解：用行初等变换将方程组的增广矩阵化为行阶梯形： 解得x=(0, -1, 1)’.

10. §3.1 Gaussian Elimination and LU Decomposition 【注】上述过程可用矩阵表示为： LU=PA 其中

11. §3.1 Gaussian Elimination and LU Decomposition 【高斯消元法思路】 (1) 用消元法将A化为上三角阵 upper-triangular matrix ; (2) 回代求解 backward substitution 。

12. §3.1 Gaussian Elimination and LU Decomposition 记 Step 1：设，计算因子 将增广矩阵/* augmented matrix */ 第 i 行  mi1 第1行，得到 其中 Step k：设 ，计算因子 且计算 消元 共进行 ？ 步 n 1

13. §3.1 Gaussian Elimination and LU Decomposition What if ? What if ? Then we must find the smallest integer k  i with , and interchange the k-th row with the i-th row. 回代 What if we can’t find such k? 定理 若A的所有顺序主子式/* determinant of leading principal submatrices */均不为0，则高斯消元无需换行即可进行到底，得到唯一解。 No unique solution exists. No unique solution exists. 注：事实上，只要 A非奇异，即 A1存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯一解。

14. §3.1 Gaussian Elimination and LU Decomposition 高斯消元法的矩阵表达----矩阵分解 课本P40-42

15. §3.1 Gaussian Elimination and LU Decomposition 3.1.2 矩阵的LU分解（LU Decomposition） 若在Gauss消去法过程中，没有进行交换行的变换，即P=I, 则分解结果为 这样的分解式称为矩阵的LU分解。 其中 【Matlab】 >> [L,U,P]=lu(A) 【作业】P.67 题1, 2

16. §3.1 Gaussian Elimination and LU Decomposition 例：单精度解方程组 /* 精确解为 和 */ 8个 8个 8个 3.1.3 选主元消去法/* Pivoting Strategies */ 用Gaussian Elimination计算： 小主元 /* Small pivot element */可能导致计算失败。

17. §3.1 Gaussian Elimination and LU Decomposition 每一步选绝对值最大的元素为主元素，保证 。 Step k:① 选取 全主元消去法/* Complete Pivoting */ ② If ik kthen 交换第 k 行与第 ik行; If jk kthen 交换第 k 列与第 jk列; ③ 消元 注：列交换改变了 xi的顺序，须记录交换次序，解完后再换回来。  列主元消去法/* Partial Pivoting, or maximal column pivoting */ 省去换列的步骤，每次仅选一列中最大的元。

18. §3.1 Gaussian Elimination and LU Decomposition 例： 注意：这两个方程组在数学上严格等价。 对每一行计算　　　　　 。为省时间，si只在初始时计算一次。以后每一步考虑子列 中　 最大的 aik为主元。 例：  注：列主元法没有全主元法稳定。   标度化列主元消去法 /* Scaled Partial Pivoting */ 注：稳定性介于列主元法和全主元法之间。

19. §3.1 Gaussian Elimination and LU Decomposition 实际应用中直接调用Gauss Elimination解3阶线性方程组的结果： 结合全主元消去后的结果：

20. §3.2 Cholesky分解 定理 设A是对称正定矩阵，则A有以下形式的LU分解： A=PTP （3.2.1） 其中P是一个上三角矩阵, 上式称为A的Cholesky分解(Choleskian decomposition). 【Matlab】 >> P=chol(A) 【作业】P.68 题4

21. §3.3 向量范数与矩阵范数 — 误差的度量 3.3.1向量范数(vector norms) 常用向量范数： 这些范数满足： 一般范数：

22. §3.3 向量范数与矩阵范数 如，设向量 x=(2, -4, -1)T, 则 x = 0.2000 0.4000 0.6000 0.8000 N1 = 2 N2 = 1.0954 N7 = 0.8153 Ninf = 0.8000 【Matlab】 x=(1:4)/5 N1=norm(x,1) N2=norm(x) N7=norm(x,7) Ninf=norm(x,inf)

23. §3.3 向量范数与矩阵范数 3.3.2矩阵范数(matrix norms) 定义矩阵A的范数定义为 命题：矩阵的常用范数： 其中，λ1,…, λn是 称为 的谱半径。 的特征值，

24. §3.3 向量范数与矩阵范数 例设矩阵 则 可以证明，矩阵范数满足：

25. (3.4.1) §3.4经典迭代法 －Jacobi迭代法与Gauss-Seidel迭代法 3.4.1．Jacobi迭代法 设有n阶方程组

26. 若系数矩阵非奇异，且(i = 1, 2,…, n)，将方程组 §3.4经典迭代法 （4.1）改写成

27. （3.4.2） （3.4.3） §3.4经典迭代法 然后写成迭代格式 （3.4.2）式也可以简单地写为

28. A = B (3.4.6) §3.4经典迭代法 写成矩阵形式： U L D Jacobi 迭代阵

29. B §3.4经典迭代法 3.4.2 Gauss-Seidel迭代法 只存一组向量即可。 … … … … 写成矩阵形式： (3.4.7) Gauss-Seidel 迭代阵

30. (3.5.1) 推论：若迭代矩阵Ｂ的某种范数 ，则（3.5.1） §3.5 迭代法的分析 ——迭代法的收敛性、收敛速度与数值稳定性 3.5.1 收敛性 定理3.1对任意初始向量X(0)及常向量F，迭代格式(3.5.1) 收敛的充分必要条件是迭代矩阵B的谱半径(B) < 1。 确定的迭代法对任意初值X(0)均收敛于方程组 X = BX + F 的唯一解x*。

31. §3.5迭代法的分析 定义1：如果矩阵的每一行中，不在主对角线上的所有元素绝对值之和小于主对角线上元素的绝对值，即 则称矩阵A (按行) 严格对角占优。 定理3.2：若A是严格对角占优，则A是非奇异的。 定理3：若线性方程组AX = b的系数矩阵A按行严格对角占优，则Jacobi迭代法和Gauss-Seidel迭代法对任意给定初值均收敛。

32. §3.5迭代法的分析 3.5.2 迭代法的误差估计与收敛速度 定理3.4设X*是方程组AX = b的同解方程X = BX + F的准确解，若迭代公式中迭代矩阵B的某种范数， 则有 (1) (2)

33. §3.6 超松驰法(SOR)及分块迭代法 迭代 （3.6.1） 加速 （3.6.2） 其中ω称为松弛因子． 定理3.５SOR方法收敛的必要条件是

34. §3.7 条件数、病态方程组与算法的稳定性 3.7.1病态方程组与数值稳定的算法 定义 对线性方程组Ax=b, 如果解x关于问题（矩阵A和b）的“微小”变化（即舍入误差）不敏感，则称此方程组是“好的”，或称为“良态的”，否则，称之为“坏的”，或称为“病态的”。 定义 对求解线性方程组Ax=b的一个(迭代)算法 ，如果关于初始值的误差，不随计算逐步放大，则称算法是数值稳定的，否则，称之为是数值不稳定的。 数值不稳定的算法是不可用的！

35. §3.7条件数、病态方程组与算法的稳定性 例 （病态方程组）考虑线性方程组 已知其精确解是 u =(1,1,1,1)T. (1)对右边常数向量b作“微小”扰动，得方程 它的解是 v = (9.2, -12.6, 4.5, -1.1)T. (2)对系数矩阵A作“微小”扰动，得方程组 注意：上述两种情况的解都是精确解！ 它的解是 w = (-81, 137, -34, 22)T.

36. §3.7条件数、病态方程组与算法的稳定性 一般地，对线性方程组 Au=b （1） (1) 考虑对此线性方程组的第一类扰动方程组 设其解为v=u+u,即设 两方程相减得 于是 由此得到 (2) （2）式右端给出了解的相对误差界。

37. §3.7条件数、病态方程组与算法的稳定性 (2) 考虑对线性方程组（1）的第二类扰动方程组 设其解为v=u+u,即设 两方程相减得 于是 由此得到 (3) 进一步得到 (4) （4）式右端给出了解的相对误差界。

38. §3.7条件数、病态方程组与算法的稳定性 3.7.2条件数与方程组病态性 定义 设A是一个非奇异矩阵，||.||是与某种向量范数相容的矩阵范数，正数 cond(A):=||A-1|| ||A|| 称为矩阵A的条件数（conditional number）。 Remark: 由上面分析可知，系数矩阵的条件数很大(>>1)（此时矩阵也称为“病态矩阵”）的线性方程组Ax=b是病态的；反之，系数矩阵的条件数很小的线性方程组Ax=b是良态的。 Example:著名的Hilbert 矩阵是病态的。 【Matlab】 >> A=hilb(n); >> c=cond(A)

39. §3.7条件数、病态方程组与算法的稳定性 定理3.6(条件数的性质)设A是一个非奇异矩阵. (1) cond(A)≥1, cond(A)=cond(A-1), con(aA)=cond(A), a ≠0. (2) 如果 Q是正交矩阵，则 cond2(Q)=1, cond2(A)= cond2(QA)= cond2(AQ), 其中cond2(A)=||A-1||2 ||A||2.

40. §3.8 稀疏矩阵的计算 例 （比较稀疏矩阵算法与普通矩阵算法：机器时间与占用空间）考虑线性方程组

41. §3.8 稀疏矩阵的计算 % sparse and full matrix n=100; row=2:n; col=1:n-1; val=ones(1,n-1); offd=sparse(row, col,val,n,n); a=sparse(1:n,1:n.4*ones(1,n),n,n); b=full(a); rh=[1:n]; tic; x=a/rh; t1=toc tic; y=b/rh; t2=toc 【实验要求】逐步增大n的值，观察发生的现象。