幻灯片 1

Size: px

Start display at page:

Download "幻灯片 1"

罐柏
5 years ago
Views:

1 线性代数方程组浙江大学控制系

2 本章内容高斯消去法 LU 分解特殊矩阵和矩阵求逆误差分析条件数迭代方法 205/6/2 数值计算方法 2

3 误差分析和方程组条件数利用逆矩阵确定方程组是否病态的方法缩放系数矩阵 A, 使其每一行的最大元素为计算缩放后的逆矩阵, 如果 A - 中有元素的值大于几倍, 则方程组是病态的将逆矩阵与原矩阵相乘, 检查结果是否接近单位阵如果不接近单位阵, 则方程组是病态的求逆矩阵的逆矩阵, 与原系数矩阵对比如果不相等, 则方程组是病态的 205/6/2 数值计算方法 3

4 向量和矩阵的范数向量的范数三维欧氏空间中向量长度概念的推广定义 n 设对任意向量 R, 按一定的规则有一实数与之对应, 记为, 若满足, 0, 且 0当且仅当 0; ( 正定 ) 2, α α, α为任意实数 ( 齐次 ) 3, + y + y, 对任意, y R ( 三角不等式 ) 则称为向量的范数 n 205/6/2 数值计算方法 4

5 向量和矩阵的范数几种向量范数向量的 2 范数 ( 欧几里德范数 ) 向量的范数向量的范数 ( 极大值范数或一致向量范数 ) 向量的 p 范数 n / n i i + + {, } { } n ma n p i i n i i ma p / p n i n i 205/6/2 数值计算方法 5

6 向量和矩阵的范数矩阵的范数定义 : 对任意 n阶方阵 A, 按一定的规则有一实数与之对应, 记为 A 若 A 满足 : A 0 且, A 0当且仅当 A 0; ( 正定 ) 2 αa α A, α为任意实数 ; ( 齐次 ) 3 A + B A + B, 对任意 A, B两个 n阶方阵 ;( 三角不等式 ) 4 AB A B ;( 矩阵乘法不等式, 相容性条件 ) 则称 A 为矩阵 A的范数 205/6/2 数值计算方法 6

7 向量和矩阵的范数几种矩阵范数矩阵的 2 范数 ( 谱范数 ) 矩阵的范数 ( 列和范数 ) 矩阵的范数 ( 行和范数 ) ( T ) 2 AA A 2 λ ma A ma aij A n j n i ma i n n j a ij 205/6/2 数值计算方法 7

8 向量和矩阵的范数谱和谱半径 A R n n 的特征值为 λ,λ 2,,λ n, 称 A 的所有特征值的集合为 A 的谱 ρ( A) ma λ i 称为 A 的谱半径 i n A 为 A 的任意一种范数, 有 ρ(a) A ( 按行 ) 严格对角占优阵如果 A 满足条件 n a > a ( i, 2,, n) ii j j i ij 即 A 的每一行对角元素的绝对值都严格大于同行其他元素绝对值之和 205/6/2 数值计算方法 8

9 病态方程组和矩阵条件数例 : (,) T 若系数矩阵有微小扰动则若右端有微小扰动则 A 0 0 () 若同时对系数矩阵和右端有扰动, 则 b (50, 48.5) T (2) (2.97, 0.99) T (3) (48.5, ) T A A () b b (2) (3) /6/2 数值计算方法 9

10 病态方程组和矩阵条件数扰动方程 ( A+ A ) b+ b 如果方程组的系数或常数项有微小改变时, 解会发生很大的改变, 则称这种方程组为病态的扰动方程的解与原方程的解相对误差不大, 称为良态方程方程组的条件问题一般说来, 在用计算机解方程组时, 实际上解的都是扰动方程, 这是由于计算过程中不可避免地会产生舍入误差对于良态问题, 只要数值方法是稳定的, 就可以得到较好的结果而对于病态问题, 即使算法是稳定的, 其计算结果有时也会很坏 205/6/2 数值计算方法 0

11 扰动方程组的误差界 AB A B b 的扰动 A ( + ) b+ b A b A 的扰动 b A A A b A b ( A+ A)( + ) b A A ( + ) A A A + A A A b b 同时扰动 ( A+ A)( + ) b+ b A A A A b A + A b A A A A 205/6/2 数值计算方法

12 病态方程组和矩阵条件数条件数设 A R n n, 非奇异,Cond(A) A A - K( A) A A 称为谱条件数 2 2 条件数的性质对任何非奇异矩阵 A, 有 Cond(A) 对任何非奇异矩阵 A, 非零常数 c, 有 Cond(cA) Cond(A) 若 P 为正交矩阵, 则 K(P), 且 K(PA)K(AP)K(A) 若 AA T, λ K( A) λ ma min ( A) ( A) 205/6/2 数值计算方法 2

13 病态方程组和矩阵条件数线性方程组 Ab 解的相对误差直接与 A 的条件数相关 A 的条件数 Cond(A) 相对大 (>>), 称 Ab 是病态方程组 / 坏条件, 或 A 是病态的 ; 当 A 的条件数 Cond(A) 相对小, 称 Ab 是良态方程组 / 好条件, 或 A 是良态的 0.99 A λ ma cond( A) 2 K( A) >> λmin 方程是病态方程 205/6/2 数值计算方法 3

14 矩阵条件数 Hilbert 矩阵是一个著名的病态矩阵, 它是一个对称正定矩阵, 当 n 3 时, Hilbert 矩阵是病态矩阵 n 越大, 条件数越大 Matlab 中 hilb() 构造 hilbert 矩阵, invhilb() 可求精确逆 A 2 3 n n n + 2 n n+ n+ 2 2n 205/6/2 数值计算方法 4

15 相对误差的事后估计 ( 近似解可靠性判别 ) 设 Ab,A 为非奇异矩阵, 为精确解, 为计算解, 残差向量 rb-a 则近似解的相对误差估计 r ( ) r ( ) cond( A) cond( A) b b Cond(A) 越小, 相对误差越小近似解的精度不仅依赖于残差向量 r, 也与矩阵 A 的条件数有关 205/6/2 数值计算方法 5

16 相对误差的事后估计 A A cond A A A ( ) 假定在求解 Ab 的过程中, 得到的解满足 r 0.00 绝对误差上界 A r 若给定 b 4, 相对误差上界 r ( ) 0.00 cond( A) b 4 205/6/2 数值计算方法 6

17 病态方程组的判别当 A 的行列式值相对小, 或 A 某些行 / 列近似线性相关, 方程组可能病态若用选主元消去法求解 Ab, 在 A 消去中出现小主元, 方程组可能病态当 A 元素数量级相差很大且无一定规则, 方程组可能病态估计条件数, 若条件数较大, 则方程组病态 205/6/2 数值计算方法 7

18 迭代求精技术设为精确解, 为得到的近似解, 则 + 残差向量 r b A A A + A A r 通过求解上述方程, 可以得到修正因子由于舍入误差的影响, 同样 2 + 不会是精确解, 可从 2 出发重复以上步骤对于 LU 分解, 只需计算残差向量再进行回代当 Ab 不过分病态时, 迭代求精是较成功的提高近似解精度的方法 205/6/2 数值计算方法 8

19 本章内容高斯消去法 LU 分解特殊矩阵和矩阵求逆误差分析条件数迭代方法 205/6/2 数值计算方法 9

迭代法的思想 a + a + a b 2 2 3 3 a + a + a b 2 22 2 23 3 2 a + a + a b 3 32 2 33 3 3 2 3 b a a a 2 2 3

20 迭代法的思想 a + a + a b a + a + a b a + a + a b b a a a b a a a b a a a 高斯 - 赛得尔方法 ( 异步迭代法 ) 雅可比方法 ( 同步迭代法 ) 205/6/2 数值计算方法 20

21 迭代法 A 任取初始向量 b G + f (0), 作 () (0) G + f ( k+ ) ( k) G + f 构造向量序列 { (k) } 求方程的近似解的方法, 称为一阶定常迭代法,G 为该迭代法的迭代矩阵如果对任意取初始近似 (0) ( ) *, 都有 lim k, 称迭代 k 法为收敛, 否则称迭代法为发散若迭代法收敛, 则称 (k) 为第 k 步迭代得到方程组的近似解 205/6/2 数值计算方法 2

22 迭代法问题构造迭代方法迭代的收敛性和收敛速度 205/6/2 数值计算方法 22

23 基本迭代方法 A 分裂为 AM-N 分裂阵 M 可选择的非奇异阵 Md 易于求解 M 选为 A 的某种近似 Ab MN+b M - N+M - b 迭代 : (0) 为初始向量 (k+) M - N+M - b 选取不同的 M 阵就得到不同迭代法 a 0 0 a2 a n a a 0 0 a a a a n A D L U 205/6/2 nn n 数值计算方法 n2 23

24 雅可比 (Jacobi) 迭代法 ( 同步迭代法 ) 设 A 为非奇异矩阵, 且 a ii 0, 选取 MD 和 ND-AL+U (k+) J (k) +f JD - (L+U) fd - b a a2... a n b a2 a22... a 2n 2 b a a... a b n n2 nn n n () 雅可比迭代法, 每迭代一次主要是计算一次矩阵乘向量, 即 J (k) (2) 计算过程中, 原始数据 A 始终不变 (3) 计算中需要两组工作单元来保存 (k) 及 (k+) a 0 a2... a n b a a 0... a b a a a... 0 b n 2 2 nn n n n2 n n a2 a n b 0... a ( ) a a k+ ( k) a ( k ) 2 a2n b + ( k) a22 a 22 + a 22 ( k+ ) ( k) n n a n bn a nn a nn n ( k+ ) ( k) i bi a ij j 205/6/2 数值计算方法 a ii j 24 j i

25 高斯 - 赛德尔 (Gauss-Seidel) 迭代法 ( 异步迭代法 ) 选取 MD-L 和 NM-AU (k+) G (k) +f G(D-L) - U f(d-l) - b b a a i n ( k+ ) ( k+ ) ( k) i i ij j ij j aii j j i+ ()G-S 迭代法每迭代一次主要是计算一次矩阵乘向量 (2) 计算 (k+) 的第 i 个分量 (k+) i 时, 利用已计算出的最新分量 (k+) j (j,,i ), 因此, 计算中只需要一组工作单元来保存 (k) 或 (k+) 205/6/2 数值计算方法 25

26 迭代终止准则近似相对百分误差 ε 对所有的 i, 近似相对百分误差小于预设的 ε s 时, 迭代终止其他的终止准则 a i ( k+ ) ( k) i i 00 ( k + ) i ( k+ ) ( k) 205/6/2 数值计算方法 26

27 迭代法例求解线性代数方程组 Jacobi 迭代 J ( k+ ) ( k) ( k) ( k+ ) ( k) ( k) ( k+ ) ( k) ( k) G-S 迭代 ( k+ ) ( k) ( k) ( k+ ) ( k+ ) ( k) ( k+ ) ( k+ ) ( k+ ) /6/2 数值计算方法 27

28 迭代法例 G-S 方法比 Jacobi 方法收敛快迭代次数 Jacobi 方法 G-S 方法 ε a 2 ε a 2 3 ε a 3 ε a 2 ε a 2 3 ε a /6/2 数值计算方法 28

29 迭代法的收敛性 G+f * G * +f 迭代 (k+) G (k) +f 引入误差向量 e (k) (k) - * 误差向量的递推公式 : e (k+) Ge (k) 于是 e (k) Ge (k-) G 2 e (k-2) G k e (0) (e (0) (0) - (*) ) ( k) k (0) k (0) k (0) e Ge G e q e 定理 : 设 G R n n, 则 G k 0( 零矩阵 )( 当 k 时 ) 的充要条件为 G 所有特征值满足 λ i <(i,2,,n) 或 G 的谱半径 ρ(g)< 205/6/2 数值计算方法 29

30 迭代法的收敛性一阶定常迭代法收敛性的基本定理设有方程组 G+f, 有迭代法 (k+) G (k) +f, 则对任选初始向量 (0), 迭代法收敛的充要条件是 ρ(g)< 迭代收敛的充分条件设有方程组 G+f 及一阶定常迭代法 (k+) G (k) +f, 如果有 G 的某种范数 G r q<, 则迭代法收敛 ; 误差估计 q < q q * ( k) ( k) ( k ) ( k) ( k ) 事后估计 * ( k ) k q q () (0) 事前估计当 q 时, 迭代法收敛缓慢 q 越小收敛速度越快 ; 可以事先估计保证误差 * (k) <ε 所需要的迭代次数 205/6/2 数值计算方法 30

31 迭代法的收敛性 matlab 中计算特征值 :eig() G 推论 : Jacobi 迭代法收敛的充要条件是 ρ(j)< (JD - (L+U)) G-S 迭代法收敛的充要条件是 ρ(g)< (G(D-L) - U) 例 : A J ρ(j).6> Jacobi 迭代不收敛 G ρ(g)< G-S 迭代收敛例 2: A J 0 ρ(j)< Jacobi 迭代收敛对于给定的线性方程组, 用雅可比法和 G-S ρ(g)> G-S 迭代不收敛法求解时可能都收敛或都不收敛, 也可能一 205/6/2 数值计算方法个收敛另一个不收敛

32 迭代法的收敛性如果 A 为 ( 按行或列 ) 严格对角占优阵, 则 Ab 的 Jacobi 方法和 G-S 方法都收敛实际上, 如果 A 为严格对角占优阵, 可以证明, J < 并且 G <, 因此,Jacobi 方法和 G-S 方法都收敛 G-S 迭代法比 Jacobi 迭代法收敛得快如果 A 为对称正定矩阵,G-S 迭代收敛例 : 判断解 Ab 的 Jacobi 方法和 G-S 方法的收敛性 0 2 A a > a 22 > a 33 > A 是严格对角占优矩阵, 两种方法都收敛 205/6/2 数值计算方法 32

迭代法的收敛性例 : 用 G-S 方法求解 3 + 7 + 3 76 2 3 + 5 + 3 28 2 3 2 + 3 5 2 3 不收敛! Iteration ε a 2 ε a 2 3 ε a 3 2 3 4 5 6 2.000-96.5-995.0-2049 2.03640 5-2.05790 5 0.7 09.83 09.90 09.89 09.90.0990 0.80000 4.

33 迭代法的收敛性例 : 用 G-S 方法求解不收敛! Iteration ε a 2 ε a 2 3 ε a [ A] [ A] 不是所有的矩阵都可以通过重排变为严格对角占优阵 /6/2 数值计算方法

34 逐次超松驰迭代法 (SOR,Successive Over-Relaation) 高斯 - 赛得尔方法的基础上进行修改 ( ) ( ) ( ) ( ) ω ( ) ( ) ( ) ( b a a ) i n k+ k k+ k k+ k i ω i + ω i i + i ij j ij j aii j j i ω,sor 方法即为 G-S 方法 0<ω<, 结果为当前迭代结果和上一次迭代结果的加权平均, 称为低松弛方法用于使得非收敛方程组收敛或者克服振荡加速收敛 <ω<2, 超松弛方法隐含假设 : 新值沿正确方向向真实解移动, 但是移动的速度慢用于加速已知是收敛的方程组的收敛速度根据经验确定 ω 值 b a a i n ( k+ ) ( k+ ) ( k) i i ij j ij j aii j j i+ 0<ω<2, 称为松驰因子 G-S 方法得到的解 205/6/2 数值计算方法 34

35 SOR 例松弛因子满足 (k) * <0 5 的迭代次数 /6/2 数值计算方法最佳松弛因子精确解为 * [ ] T, 初始向量 (0) [ ] T

36 SOR 收敛性 SOR 的矩阵形式 ω ( b a a ) i n ( k+ ) ( k) ( k+ ) ( k) i i + i ij j ij j aii j j i M(D- ω L)/ω ( ) ( ) [ ] ( k ) ( k) D ωl + ω D+ ωu + ωb ( k + ) ( k) Gω + f Gω ( D ωl) ( ω) D+ ωu ω ω f ( D L) b [ ] SOR 法收敛的充要条件是 ρ(gω)< 如果 A 为对称正定阵, 且 0<ω<2, 则解 Ab 的 SOR 法收敛 205/6/2 数值计算方法 36

37 MATLAB 中的函数矩阵分析线性方程组函数描述函数描述 cond 计算矩阵条件数右除 / 和左除 \ help slash norm 计算矩阵或向量的范数 lu lu 分解 inv 矩阵求逆 chol cholesky 分解 pinv 求矩阵伪逆 det 计算行列式的值 rank 求秩 eig,eigs 矩阵特征值 205/6/2 数值计算方法 37

38 第三章总结各种方法方法稳定性精度应用范围编程难度备注图解法差受限比数值方法耗时 Cramer 法则受舍入误差影响受限方程数多于 3 个时计算复杂列主元高斯消去 LU 分解受舍入误差影响一般, 适用于方程组系数矩阵为低阶稠密矩阵带状矩阵中等, 计算量较小, 存储量较大实现消去法 ; 进行矩阵求逆计算迭代法 (Jacobi 和 G-S) 可能不收敛优秀收敛时, 适用于大型稀疏线性方程组较简单, 计算量有时较大, 存储量较小 205/6/2 数值计算方法 38

39 第三章总结重要内容直接法高斯消去 LU 分解迭代 Jacobi Gauss-Seidel SOR 具体算法问题及改进范数条件数病态方程组 205/6/2 数值计算方法 39

目录三种迭代格式 Jacobi 迭代法 Gauss-Seidel 迭代超松弛迭代法 Jacobi 与 G-S 迭代的收敛性分析收敛的充分必要条件收敛的充分条件及误差估计收敛速度平均

目录三种迭代格式 Jacobi 迭代法 Gauss-Seidel 迭代超松弛迭代法 Jacobi 与 G-S 迭代的收敛性分析收敛的充分必要条件收敛的充分条件及误差估计收敛速度平均线性方程组的古典迭代解法目录 1 3.1 三种迭代格式 3.1.1 Jacobi 迭代法 3.1.2 Gauss-Seidel 迭代 3.1.3 超松弛迭代法 2 3.2 Jacobi 与 G-S 迭代的收敛性分析 3.2.1 收敛的充分必要条件 3.2.2 收敛的充分条件及误差估计 3 3.3 收敛速度 3.3.1 平均收敛速度和渐近收敛速度 3.3.2 模型问题 3.3.3 Jacobi 和