2线性方程组的直接解法(12).ppt

1 现代数值计算 第二章线性方程组的数值解法 第二章线性方程组的数值解法 2 0引言 2 1Gauss消去法 2 2矩阵的三角分解 2 3QR分解和奇异值分解 在自然科学和工程技术中很多问题的解决常常归结为解线性代数方程组 例如 电学中的网络问题用最小二乘法求实验数据的曲线拟合问题解非线性方程组问题用差分法或者有限元方法解常微分方程偏微分方程边值问题等都导致求解线性代数方程组 2 0引言 这些方程组的系数矩阵大致分为两种一种是低阶稠密矩阵 例如 阶数大约为 150 另一种是大型稀疏矩阵 即矩阵阶数高且零元素较多 2 0引言 设有线性方程组Ax b 其中为非奇异阵 关于线性方程组的数值解法一般有两类 直接法与迭代法 2 0引言 1 直接法就是经过有限步算术运算 可求得方程组精确解的方法 若计算过程中没有舍入误差 但实际计算中由于舍入误差的存在和影响 这种方法也只能求得线性方程组的近似解 本章将阐述这类算法中最基本的高斯消去法及其某些变形 这类方法是解低阶稠密矩阵方程组的有效方法 近十几年来直接法在求解具有较大型稀疏矩阵方程组方面取得了较大进展 2 0引言 2 迭代法就是用某种极限过程去逐步逼近线性方程组精确解的方法 迭代法具有需要计算机的存贮单元较少 程序设计简单 原始系数矩阵在计算过程中始终不变等优点 但存在收敛性及收敛速度问题 迭代法是解大型稀疏矩阵方程组 尤其是由微分方程离散后得到的大型方程组 的重要方法 第6章介绍迭代法解线性方程组 2 0引言 高斯 Gauss 消去法是解线性方程组最常用的方法之一它的基本思想是通过逐步消元 行的初等变换 把方程组化为系数矩阵为三角形矩阵的同解方程组 然后用回代法解此三角形方程组 简单形式 得原方程组的解 例如 直接法的基础 2 1Gauss消去法 下面讨论一般的解n阶方程组的高斯消去法 1 消去过程将原方程组记为A 1 x b 1 其中A 1 aij 1 n n aij n n b 1 b 1 第一次消元 2 1Gauss消去法 1 消去过程 1 第一次消元 其中 2 1Gauss消去法 注 若a11 1 0 则在第1列中至少有一个元素不为0 可交换行后再消元 2 第k次消元 2 1Gauss消去法 1 消去过程 2 第k次消元 注 为减少计算量 令 则 2 1Gauss消去法 1 消去过程 3 当k n 1时得完成第n 1次消元后得到与原方程组等价的三角形方程组A n x b n 注 当det A 0时 显然有aii i 0 i 1 n 称aii i 为主元素 2 1Gauss消去法 2 回代过程求解三角形方程组A n x b n 得到求解公式注 求解过程称为回代过程 2 1Gauss消去法 3 算法设计令i k 1 k 2 n 2 1Gauss消去法 3 算法设计functionx gauss a m size a n m 1 fork 1 n 1fori k 1 na i a i a k a i k a k k endendforj n 1 1a j n 1 a j n 1 a j j 1 n a j 1 n n 1 a j j endx a 1 n n 1 2 1Gauss消去法 3 算法设计i k 1 k 2 n 2 1Gauss消去法 3 算法设计在Excel中 2 1Gauss消去法 3 算法设计在Excel中 2 1Gauss消去法 4 Gauss消去法的计算量以乘除法的次数为主 1 消元过程 第k步时 n k n k n k 1 n k n k 2 共有 2 1Gauss消去法 4 Gauss消去法的计算量 1 消元过程 2 回代过程 求xi中 乘n i次 除1次 共n i 1次 i 1 n 1 共有 2 1Gauss消去法 4 Gauss消去法的计算量 1 消元过程 2 回代过程 3 总次数为注 当n 20时约为2670次 比克莱姆法则9 7 1020次大大减少 2 1Gauss消去法 5 说明 1 若消元过程中出现akk k 0 则无法继续 2 若akk k 0 但较小 则小主元做除数将产生大误差 3 每次消元都选择绝对值最大者作主元 称为高斯主元消去法 2 1Gauss消去法 5 说明 4 通常第k步选择 第k列主对角元以下元素绝对值最大者作主元 该行与第k行对调 称为列主元消去法 2 1Gauss消去法 5 说明 例1 用舍入三位有效数字求解线性方程组 准确解是x1 10 0 x2 1 00 解 1 不选主元的Gauss消去法计算结果 x1 10 0 x2 1 01 此解无效 2 按列选主元的Gauss消去法计算结果 x1 10 0 x2 1 00 2 1Gauss消去法 5 说明 例2 求解线性方程组和 准确解是x1 1 0 x2 1 0 2 1Gauss消去法 5 说明 5 行标度化当线性方程组的系数矩阵的元素大小相差很大时 则可能引进因丢失有效数字而产生的误差 并且舍入误差的影响严重 即使用Gauss主元消去法得到的解也不可靠 为避免这一问题 可将系数矩阵的行元素按比例增减以使元素的变化减小 2 1Gauss消去法 5 说明 5 行标度化如用每行元素的最大模除该行各元素 使它们的模都不大于1 这叫行标度化 其目的是要找到真正的主元 消元过程仍是对原方程组进行的 只不过在Gauss列主元消去法的算法中 按列选主元ck时 应修改为其中 2 1Gauss消去法 6 列主元消去法算法 2 1Gauss消去法 6 列主元消去法算法functionx gaussa a m size a n m 1 x zeros n 1 fork 1 n 1 c i max abs a k n k q i k 1 ifq kd a q a q a k a k dendfori k 1 na i a i a k a i k a k k endendforj n 1 1x j a j n 1 a j j 1 n x j 1 n a j j end 2 1Gauss消去法 6 列主元消去法算法在Excel中 2 1Gauss消去法 2 2 1LU分解和LDU分解1 原理若A LU 则Ax b LUx b 若其中L U为三角形矩阵 则方程组易解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解1 原理定义 1 称为单位下三角阵 2 设L为单位下三角阵 U为上三角阵 若A LU 则称A可三角 LU 分解 并称A LU为A的三角分解 杜利特尔分解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解1 原理定理 1 A aij n n有唯一LU分解 A的顺序主子式 k 0 k 1 2 n 1 2 若A aij n n可逆 则存在置换阵P 使PA的n个顺序主子式全不等于0注 由Ax b PAx Pb知 经适当行交换后 A总可三角分解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解2 LU分解设A的各阶顺序主子式不为0 且A LU 即 1 主对角线 含 上边 当列标m 行标k时 lkm 0j k k 1 nk 1 2 n 2 2矩阵的三角分解 2 2 1LU分解和LDU分解2 LU分解设A的各阶顺序主子式不为0 且A LU 即 2 主对角线下边 当行标m 列标k时 umk 0i k 1 k 2 nk 1 2 n 1欲求lik和ukj 2 2矩阵的三角分解 2 2 1LU分解和LDU分解2 LU分解 1 主对角线 含 上边 当列标m 行标k时 lkm 0j k k 1 nk 1 2 n 2 主对角线下边 当行标m 列标k时 umk 0i k 1 k 2 nk 1 2 n 1设k 1 a1j l11u1j u1j a1jj 1 2 nai1 li1u11 li1 ai1 u11i 2 3 n 欲求lik和ukj 2 2矩阵的三角分解 2 2 1LU分解和LDU分解2 LU分解一般地 最后 2 2矩阵的三角分解 2 2 1LU分解和LDU分解3 求解方程组下三角Ly b 上三角Ux y 2 2矩阵的三角分解 2 2 1LU分解和LDU分解3 求解方程组 例3 用杜丽特尔法解方程组解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解3 求解方程组 例3 用杜丽特尔法解方程组解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解3 求解方程组 例3 用杜丽特尔法解方程组解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解3 求解方程组 例3 用杜丽特尔法解方程组解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解3 求解方程组 例3 用杜丽特尔法解方程组解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解3 求解方程组 例3 用杜丽特尔法解方程组解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解4 紧凑格式 2 2矩阵的三角分解 2 2 1LU分解和LDU分解5 代码functionx lua a b m size a n m 1 x zeros n 1 y zeros n 1 fork 1 na k k n a k k n a k 1 k 1 a 1 k 1 k n a k 1 n k a k 1 n k a k 1 n 1 k 1 a 1 k 1 k a k k endU triu a 0 L eye n tril a 1 fori 1 ny i b i L i 1 i 1 y 1 i 1 endfori n 1 1x i y i U i i 1 n x i 1 n U i i end 2 2矩阵的三角分解 2 2 1LU分解和LDU分解5 代码紧凑格式 functionx luj a m size a n m 1 x zeros n 1 fork 1 na k k n 1 a k k n 1 a k 1 k 1 a 1 k 1 k n 1 a k 1 n k a k 1 n k a k 1 n 1 k 1 a 1 k 1 k a k k endU triu a 0 forj n 1 1x j a j n 1 U j j 1 n x j 1 n U j j end 2 2矩阵的三角分解 2 2矩阵的三角分解 2 2 1LU分解和LDU分解6 Excel实现 2 2矩阵的三角分解 2 2 1LU分解和LDU分解6 Excel实现 2 2 1LU分解和LDU分解7 LDU分解设A LU 令D diag u11 u22 unn 则U1 D 1U仍是一个单位上三角阵故矩阵的三角分解除了杜丽特尔分解 还有如下 1 克劳特分解 A LU L为下三角阵 U为单位上三角阵 2 LDU分解 A LDU L为单位下三角阵 D为对角阵 U为单位上三角阵上述分解均具有唯一性 2 2矩阵的三角分解 2 2 2乔列斯基分解1 原理定理 若A正定 则A可唯一分解为A LDLT其中L为单位下三角 D为元素全是正数的对角阵证明 设A LDU L为单位下三角阵 D为对角阵 U为单位上三角阵 由对称性LDU UTDLT 由唯一性知L UT 所以A LDLT 2 2矩阵的三角分解 2 2 2乔列斯基分解1 原理定理 若A正定 则A可唯一分解为A LDLT其中L为单位下三角 D为元素全是正数的对角阵证明 设A LDLT D diag d1 d2 dn 对于ei 0 0 1 0 0 T 存在xi 0 使得LTxi ei另外 xiTAxi xiT LDLT xi LTxi TD LTxi eiTDei di由于A正定 有di xiTAxi 0 2 2矩阵的三角分解 2 2 2乔列斯基分解1 原理定理 若A正定 则A可唯一分解为A LDLT其中L为单位下三角 D为元素全是正数的对角阵乔列斯基分解 若A正定 则A可唯一分解为A GGT其中G为下三角阵记则有A LDLT LD1 2D1 2LT LD1 2 LD1 2 T GGT 2 2矩阵的三角分解 2 2 2乔列斯基分解2 算法设A GGT则有 从而 2 2矩阵的三角分解 2 2 2乔列斯基分解2 算法设A GGT则有 注 由上式看出 所以即算法是稳定的 2 2矩阵的三角分解 2 2 2乔列斯基分解2 算法Ax b GGTx b Gy b GTx y求解公式 上述方法称为 平方根法 注 不能选主元作行交换 破坏对称性 2 2矩阵的三角分解 2 2 2乔列斯基分解 例4 用平方根法解方程组解 2 2矩阵的三角分解 2 2 2乔列斯基分解3 代码functionx choleskey a b m size a n m 1 x zeros n 1 y zeros n 1 G zeros m forj 1 nG j j sqrt a j j G j 1 j 1 G j 1 j 1 G j 1 n j a j 1 n j G j 1 n 1 j 1 G j 1 j 1 G j j endfori 1 ny i b i G i 1 i 1 y 1 i 1 G i i endfori n 1 1x i y i G i 1 n i x i 1 n G i i end 2 2矩阵的三角分解 2 2 2乔列斯基分解4 Excel中的计算 2 2矩阵的三角分解 2 2 3追赶法在实际问题中 经常遇到以下形式的方程组这种方程组的系数矩阵称为三对角矩阵 以下针对该方程组的特点提供一种简便有效的算法 追赶法 2 2矩阵的三角分解 2 2 3追赶法作克劳特分解 A LU 易知L U具如下形式 比较第i行的元素得 2 2矩阵的三角分解 2 2 3追赶法作克劳特分解 A LU 易知L U具如下形式 计算过程 l1 u1 l2 u2 ln 2 2矩阵的三角分解 2 2 3追赶法有了A的LU分解后 线性方程组Ax d等价于两个简单的方程组 Ly d Ux y Ly d即 2 2矩阵的三角分解 2 2 3追赶法有了A的LU分解后 线性方程组Ax d等价于两个简单的方程组 Ly d Ux y Ux y即 2 2矩阵的三角分解 2 2 3追赶法分解公式是计算y的公式是计算过程 l1 u1 y1 l2 u2 y2 ln yn 追 2 2矩阵的三角分解 2 2 3追赶法分解公式是计算y的公式是计算过程 2 2矩阵的三角分解 2 2 3追赶法分解公式是计算x的公式是 赶 2 2矩阵的三角分解 2 2 3追赶法 2 2矩阵的三角分解 2 2 3追赶法 例5 用追赶法解解 追 2 2矩阵的三角分解 2 2 3追赶法 例5 用追赶法解解 赶 2 2矩阵的三角分解 2 2 3追赶法 2 2矩阵的三角分解 2 2 3追赶法functionx chase a b c d m size b n m 2 L linspace 0 0 n U L x L y L L 1 b 1 y 1 d 1 L 1 fork 1 n 1U k c k L k L k 1 b k 1 a k U k y k 1 d k 1 a k y k L k 1 endx n y n fori n 1 1 1x i y i U i x i 1 end 2 2矩阵的三角分解 2 2 3追赶法 2 2矩阵的三角分解 2 2 3追赶法functiond chase a b c d m size b n m 2 d 1 d 1 b 1 fork 1 n 1c k c k b k b k 1 b k 1 a k c k d k 1 d k 1 a k d k b k 1 endfori n 1 1 1d i d i c i d i 1 end 2 2矩阵的三角分解 2 2 3追赶法 2 2矩阵的三角分解 定理 若A为对角占优 diagonallydominant 的三对角阵 且满足 b1 c1 0 bi ai ci aici 0 i 2 3 n 1 bn an 0 则追赶法可解以A为系数矩阵的方程组注 1 如果A是严格对角占优阵 则不要求三对角线上的所有元素非零 2 2矩阵的三角分解 注 1 如果A是严格对角占优阵 则不要求三对角线上的所有元素非零 2 根据不等式可知 分解过程中 矩阵元素不会过分增大 算法保证稳定 3 运算量为O 6n 2 2矩阵的三角分解 2 3QR分解和奇异值分解 矩阵分解是将矩阵分解为数个具有特殊性质的矩阵因子的乘积 除了三角分解以外 还有QR分解 QRFactorization 和奇异值分解法 SingularValueDecompostion 2 3QR分解和奇异值分解 2 3 1正交矩阵 定义2 3 1 若矩阵Q Rn n 且满足QQT QTQ I则称矩阵Q为正交矩阵正交矩阵Q有如下性质 1 Q 1 QT 2 det Q 1 3 Qx的长度与x的长度相等 下面介绍几类特殊的正交矩阵 2 3QR分解和奇异值分解 2 3 1正交矩阵1 置换矩阵将单位矩阵任意两行 列 交换得到的矩阵 称为置换矩阵 譬如 将单位矩阵的第i行和第j列交换 得到置换矩阵Pij 任意个置换矩阵的乘积仍然是置换矩阵 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 对于某个角度 记s sin c cos 那么 是一个正交矩阵 记w x y T为二维平面中的一个向量 用极坐标表示为w rsin rcos T 那么即Gw表示将向量w顺时针旋转 角所得到的向量 如图2 1所示 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 推广到n n的情形 形如的矩阵称为Givens矩阵或Givens变换 或称 平面 旋转矩阵 旋转变换 其中 为旋转的角度 显然 G i j 也是正交矩阵 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 若x Rn n y G i j x 则y的分量为如果要使yj 0只要选择 满足即可 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 例6 用Givens变换将海森伯格 Hdssenberg 型矩阵化为上三角矩阵 解 首先消去A 2 1 构造Givens变换G 1 2 其中 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 例6 用Givens变换将海森伯格 Hdssenberg 型矩阵化为上三角矩阵 解 从而G 1 2 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 例6 用Givens变换将海森伯格 Hdssenberg 型矩阵化为上三角矩阵 解 A1 G 1 2 A 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 例6 用Givens变换将海森伯格 Hdssenberg 型矩阵化为上三角矩阵 解 其次消去A1 3 2 构造Givens变换G 2 3 其中 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 例6 用Givens变换将海森伯格 Hdssenberg 型矩阵化为上三角矩阵 解 从而得A2 G 2 3 A1 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 例6 用Givens变换将海森伯格 Hdssenberg 型矩阵化为上三角矩阵 解 最后消去A2 4 3 构造Givens变换G 3 4 其中 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 例6 用Givens变换将海森伯格 Hdssenberg 型矩阵化为上三角矩阵 解 从而 上三角矩阵R为A3 G 3 4 A2 2 3QR分解和奇异值分解 2 3 1正交矩阵2 旋转矩阵 Givens变换 例6 用Givens变换将海森伯格 Hdssenberg 型矩阵化为上三角矩阵 2 3QR分解和奇异值分解 2 3 1正交矩阵3 反射矩阵 Householder变换 设w Rn 且 w 2 1 则P I 2wwT称为Householder变换或者Householder矩阵Householder矩阵有如下性质 1 PT P 即P是对称矩阵 2 PPT P2 I 2wwT 2wwT 4w wTw wT I 即P是正交阵 3 如图2 2 设w是R3上的一个单位向量 并设S为过原点且与w垂直的平面 则一切v Rn可以分解成v v1 v2 其中v1 S v2 S 2 3QR分解和奇异值分解 2 3 1正交矩阵3 反射矩阵 Householder变换 Householder矩阵有如下性质 1 PT P 即P是对称矩阵 2 PPT P2 I 2wwT 2wwT 4w wTw wT I 即P是正交阵 3 如图2 2 设w是R3上的一个单位向量 并设S为过原点且与w垂直的平面 则一切v Rn可以分解成v v1 v2 其中v1 S v2 S 不难验证Pv1 v2 Pv2 v2 所以Pv v1 v2这样 v经变换后的象Pv是v关于S对称的向量 2 3QR分解和奇异值分解 2 3 1正交矩阵3 反射矩阵 Householder变换 Householder矩阵有如下性质 1 PT P 即P是对称矩阵 2 PPT P2 I 2wwT 2wwT 4w wTw wT I 即P是正交阵 3 如图2 2 设w是R3上的一个单位向量 并设S为过原点且与w垂直的平面 则一切v Rn可以分解成v v1 v2 其中v1 S v2 S 所以 Householder变换又称镜面反射变换Householder矩阵也称初等反射矩阵 2 3QR分解和奇异值分解 2 3 1正交矩阵3 反射矩阵 Householder变换 一个重要的应用是对x 0 求Householder矩阵P 使得Px ke1 由正交矩阵的性质可知 Px 2 ke1 2 x 2 即k x 2 由上面所讨论的P的构造 有 令 u x ke1 w u u 2设x x1 xn T 为了使x ke1计算时不损失有效数位取k sign x1 x 2 则u x1 sign x1 x 2 x2 xn T从而P I uuT 其中 u 22 1 x 2 x 2 x1 1 e1 1 0 0 T 2 3QR分解和奇异值分解 2 3 1正交矩阵3 反射矩阵 Householder变换 例7 已知x 3 5 1 1 T 求Householder矩阵P 使得Px 6e1 其中 x 2 6 解 取k 6 u x ke1 9 5 1 1 T u 2 108 1 54则 2 3QR分解和奇异值分解 2 3 2QR分解本节给出正交三角分解 又称QR分解 的存在性定理和唯一性定理定理2 3 4设A Rn n 则存在正交阵P 使得PA R 其中R为上三角阵 2 3QR分解和奇异值分解 2 3 2QR分解定理2 3 4设A Rn n 则存在正交阵P 使得PA R 其中R为上三角阵 证明 首先考虑A的第一列a1 a11 a21 an1 T 可找到Householder矩阵P1 使得P1a1的元素除了第1个以外都为0 同理 找到P2使得P2P1A的第二列对角元以下元素为0 而第一列对角元以下元素与P1A一样是0 依次这样下去 可以得到Pn 1Pn 2 P1A R 其中R为上三角形矩阵 P Pn 1Pn 2 P1为正交阵 给出构造性证明 2 3QR分解和奇异值分解 2 3 2QR分解定理2 3 4设A Rn n 则存在正交阵P 使得PA R 其中R为上三角阵 定理2 3 5设A Rn n 且A非奇异 则存在正交阵Q与上三角阵R 使得A有如下分解A QR且当R的对角元均为正时 分解是唯一的 该定理保证了A可分解为A QR 若A非奇异 则R也非奇异 如果不规定R的对角元为正 则分解不是唯一的 2 3QR分解和奇异值分解 2 3 2QR分解 例8 用Householder变换作矩阵A的QR分解解 找Householder矩阵P1 R3 3 使则有 2 3QR分解和奇异值分解 2 3 2QR分解 例8 用Householder变换作矩阵A的QR分解解 再找 R2 2 使 1 44949 3 44949 T 0 T 得且 2 3QR分解和奇异值分解 2 3 2QR分解 例8 用Householder变换作矩阵A的QR分解解 这是一个下三角矩阵 但对角元皆为负数 只要令D I R P2P1A就是对角元为正的上三角矩阵 使得A QR 其中 2 3QR分解和奇异值分解 2 3 2QR分解 例8 用Householder变换作矩阵A的QR分解 2 3QR分解和奇异值分解 2 3 2QR分解QR分解是计算特征值的有力工具 也是用于其它矩阵计算问题 包括解方程组Ax b 这只要令y QTb 再解上三角形组Rx y 这个计算过程是稳定的 也不必选主元 但是计算量比高斯消去法将近大一倍 2 3QR分解和奇异值分解 2 3 2QR分解Matlab以qr函数来执行QR分解法 其语法为 Q R qr A 其中Q正交矩阵 而R为上三角矩阵 此外 原矩阵A不必为方阵 如果矩阵A Rm n 则Q Rm m且R Rm n 2 3QR分解和奇异值分解 2 3 3奇异值分解奇异值分解时线性代数中一种重要的矩阵分解 在信号和图像处理 统计和数据压缩等领域有重要应用 定义2 3 7设A Cm n AHA的n个特征值的非负平方根叫做A的奇异值 记为 i A AH表示矩阵A的共轭转置 如果把AHA的特征值记为 i A 则 i A 2 3QR分解和奇异值分解 2 3 3奇异值分解关于矩阵的奇异值分解 singularvaluedecomposition SVD 我们有如下定理 定理2 3 9 奇异值分解 设A Cm n 则存在酉阵U Cm n和V Cm n