线性方程组求解的数值方法.ppt

1、Ch3 线性方程组求解的数值方法,Numerical Solutions of Systems of Linear Equations,3.0 Introduction,线性方程组,工程问题（电子网络，船体放样等）,自然科学（实验数据的最小二乘法曲线拟合）,解非线性方程组,解常/偏微分方程组（差分法，有限元法）,3.0 Introduction,线性方程组的系数矩阵：,（1）低价稠密矩阵（阶数150）；,（2）大型稀疏矩阵（阶数高且零元素较多）。,求解线性方程组的方法：,（1）直接法: 经过有限步算术运算，求得精确解（假设计算过程没有舍入误差）。 如Gauss消去法，三角分解法。,（2）间接法

2、(迭代法)： 通过迭代序列，逐步逼近方程组的解。如Jacobi迭代法，Gauss-Seidel迭代法。,3.0 Introduction,【最简单的情形】三角形线性方程组：,简记为,，A为上三角矩阵。,若所有aii0，可用“回代”过程得到方程组的解。,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 so

3、lution 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,3.0 Introduction,3.0 Introduction,【作业】（1）编写Matlab程序解下列上三角形方程组：,（2）编写Matlab程序解下列下三角形方程组：,3.1 Gauss消去

4、法与LU分解法,直接法的基本思想：,将线性方程组化成与之等价的上三角形或下三角形，再用回代法求解。它的核心是矩阵分解。,核心：矩阵分解。,Gauss消去法（Gaussian elimination）：,（对方程组的三种变换）,（1）交换两个方程的次序；,（2）用一个非零常数乘一个方程；,（3）将一个方程的非零倍数加到另一个方程上去。,这等价与对增广矩阵进行三种行初等变换，将它化为行阶梯形。,【例3.1】解下列线性方程组：,解：用行初等变换将方程组的增广矩阵化为行阶梯形：,解得 x=(0, -1, 1).,3.1 Gaussian Elimination and LU Decomposition

5、,【注】上述过程可用矩阵表示为：,LU=PA,其中,3.1 Gaussian Elimination and LU Decomposition,【高斯消元法思路】,(1) 用消元法将A化为上三角阵 upper-triangular matrix ;,(2) 回代求解 backward substitution 。,3.1 Gaussian Elimination and LU Decomposition,3.1 Gaussian Elimination and LU Decomposition,消元,记,其中,Step k：设 ，计算因子 且计算,共进行 ？ 步,n 1,回代,No uniqu

6、e solution exists.,What if we cant find such k ?,No unique solution exists.,定理,若A的所有顺序主子式 /* determinant of leading principal submatrices */ 均不为0，则高斯消元无需换行即可进行到底，得到唯一解。,注：事实上，只要 A 非奇异，即 A1 存在，则可通过逐次消元及行交换，将方程组化为三角形方程组，求出唯一解。,3.1 Gaussian Elimination and LU Decomposition,高斯消元法的矩阵表达-矩阵分解,课本P40-42,3.1

7、Gaussian Elimination and LU Decomposition,3.1.2 矩阵的LU分解（LU Decomposition）,若在Gauss消去法过程中，没有进行交换行的变换，即P=I, 则分解结果为,3.1 Gaussian Elimination and LU Decomposition,这样的分解式称为矩阵的LU分解。,【Matlab】, L,U,P=lu(A),其中,【作业】P.67 题1, 2,3.1.3 选主元消去法 /* Pivoting Strategies */,例：单精度解方程组,用Gaussian Elimination计算：,8个,小主元 /* S

8、mall pivot element */ 可能导致计算失败。,3.1 Gaussian Elimination and LU Decomposition, 全主元消去法 /* Complete Pivoting */,每一步选绝对值最大的元素为主元素，保证 。,Step k: 选取, If ik k then 交换第 k 行与第 ik 行; If jk k then 交换第 k 列与第 jk 列;, 消元,注：列交换改变了 xi 的顺序，须记录交换次序，解完后再换回来。, 列主元消去法 /* Partial Pivoting, or maximal column pivoting */,省去

9、换列的步骤，每次仅选一列中最大的元。,3.1 Gaussian Elimination and LU Decomposition,例：,注：列主元法没有全主元法稳定。,例：, 标度化列主元消去法 /* Scaled Partial Pivoting */,对每一行计算 。为省时间，si 只在初始时计算一次。以后每一步考虑子列 中 最大的 aik 为主元。,注：稳定性介于列主元法和全主元法之间。,3.1 Gaussian Elimination and LU Decomposition,实际应用中直接调用Gauss Elimination 解3阶线性方程组的结果：,结合全主元消去后的结果：,3.

10、1 Gaussian Elimination and LU Decomposition,3.2 Cholesky分解,【Matlab】,设A是对称正定矩阵，则A有以下形式的LU分解：,A=PTP （3.2.1）,其中P是一个上三角矩阵, 上式称为A的Cholesky分解(Choleskian decomposition)., P=chol(A),定理,【作业】P.68 题4,3.3 向量范数与矩阵范数, 误差的度量,3.3.1 向量范数 (vector norms),常用向量范数：,这些范数满足：,一般范数：,如，设向量 x=(2, -4, -1)T, 则,3.3 向量范数与矩阵范数,【Mat

11、lab】,x=(1:4)/5 N1=norm(x,1) N2=norm(x) N7=norm(x,7) Ninf=norm(x,inf),x = 0.2000 0.4000 0.6000 0.8000 N1 = 2 N2 = 1.0954 N7 = 0.8153 Ninf = 0.8000,定义 矩阵A的范数定义为,3.3 向量范数与矩阵范数,3.3.2 矩阵范数 (matrix norms),命题：矩阵的常用范数：,其中，1, n是,的特征值，,称为 的谱半径。,例 设矩阵,3.3 向量范数与矩阵范数,则,可以证明，矩阵范数满足：,3.4 经典迭代法,Jacobi迭代法与Gauss-Seid

12、el迭代法,3.4.1Jacobi迭代法,设有n阶方程组,(3.4.1),若系数矩阵非奇异，且 (i = 1, 2, n)，将方程组,（4.1）改写成,3.4 经典迭代法,然后写成迭代格式,（3.4.2）,（3.4.2）式也可以简单地写为,（3.4.3）,3.4 经典迭代法,写成矩阵形式：,L,U,D,Jacobi 迭代阵,3.4 经典迭代法, ,只存一组向量即可。,写成矩阵形式：,Gauss-Seidel 迭代阵,3.4.2 Gauss-Seidel迭代法,(3.4.7),3.4 经典迭代法,定理3.1 对任意初始向量X(0)及常向量F，迭代格式(3.5.1),(3.5.1),推论：若迭代矩

13、阵的某种范数 ，则（3.5.1）,收敛的充分必要条件是迭代矩阵B的谱半径(B) 1。,确定的迭代法对任意初值X(0)均收敛于方程组,X = BX + F,3.5 迭代法的分析,迭代法的收敛性、收敛速度与数值稳定性,3.5.1 收敛性,的唯一解x*。,定义1：如果矩阵的每一行中，不在主对角线上的所有元素绝对值之和小于主对角线上元素的绝对值，即,则称矩阵A (按行) 严格对角占优。,定理3.2：若A是严格对角占优，则A是非奇异的。,3.5 迭代法的分析,定理3：若线性方程组AX = b的系数矩阵A按行严格对角占优，则Jacobi迭代法和Gauss-Seidel迭代法对任意给定初值均收敛。,3.5.

14、2 迭代法的误差估计与收敛速度,定理3.4 设X*是方程组AX = b的同解方程X = BX + F的准确解，若迭代公式中迭代矩阵B的某种范数，,(1),(2),则有,3.5 迭代法的分析,迭代,加速,（3.6.2）,3.6 超松驰法(SOR)及分块迭代法,定理3. SOR方法收敛的必要条件是,其中称为松弛因子,（3.6.1）,3.7 条件数、病态方程组与算法的稳定性,3.7.1 病态方程组与数值稳定的算法,定义 对线性方程组Ax=b, 如果解x关于问题（矩阵A和b）的“微小”变化（即舍入误差）不敏感，则称此方程组是“好的”，或称为“良态的”，否则，称之为“坏的”，或称为“病态的”。,定义 对

15、求解线性方程组Ax=b的一个(迭代)算法 ，如果关于初始值的误差，不随计算逐步放大，则称算法是数值稳定的，否则，称之为是数值不稳定的。,数值不稳定的算法是不可用的！,例 （病态方程组）考虑线性方程组,3.7条件数、病态方程组与算法的稳定性,已知其精确解是 u =(1,1,1,1)T.,(1)对右边常数向量b作“微小”扰动，得方程,(2)对系数矩阵A作“微小”扰动，得方程组,它的解是 v = (9.2, -12.6, 4.5, -1.1)T.,它的解是 w = (-81, 137, -34, 22)T.,注意：上述两种情况的解都是精确解！,一般地，对线性方程组 Au=b （1）,两方程相减得,于

16、是,(2),由此得到,设其解为v=u+u,即设,(1) 考虑对此线性方程组的第一类扰动方程组,（2）式右端给出了解的相对误差界。,3.7条件数、病态方程组与算法的稳定性,两方程相减得,于是,(3),由此得到,设其解为v=u+u,即设,(2) 考虑对线性方程组（1）的第二类扰动方程组,（4）式右端给出了解的相对误差界。,进一步得到,(4),3.7条件数、病态方程组与算法的稳定性,3.7.2 条件数与方程组病态性,定义 设A 是一个非奇异矩阵，|.|是与某种向量范数相容的矩阵范数，正数 cond(A):=|A-1| |A| 称为矩阵A的 条件数（conditional number）。,Remar

17、k: 由上面分析可知，系数矩阵的条件数很大(1)（此时矩阵也称为“病态矩阵”）的线性方程组Ax=b是病态的；反之，系数矩阵的条件数很小的线性方程组Ax=b是良态的。, A=hilb(n); c=cond(A),Example: 著名的Hilbert 矩阵是病态的。,【Matlab】,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.,3.7条件数、病态方程组与算法的稳定性,3.8 稀疏矩