LU和QR分解法解线性方程组

LU 和 QR 法解线性方程组 一 问题描述 求解方程组 2011612 6384 10278 5124 4 3 2 1 x x x x 3 7 7 2 要求 1 编写用三角 LU 分解法求解线性方程组 2 编写用正交三角 QR 分解法求解线性方程组 2 问题分析 求解线性方程组 Ax b 其实质就是把它的系数矩阵 A 通过各种变换成一个下三角或 上三角矩阵 从而简化方程组的求解 因此 在求解线性方程组的过程中 把系数矩阵 A 变换成上三角或下三角矩阵显得尤为重要 然而矩阵 A 的变换通常有两种分解方法 LU 分解法和 QR 分解法 1 LU 分解法 将 A 分解为一个下三角矩阵 L 和一个上三角矩阵 U 即 A LU 其中 L U 1 0 01 001 21 21 nn ll l nn n n u uu uuu 00 00 0 222 11211 2 QR 分解法 将 A 分解为一个正交矩阵 Q 和一个上三角矩阵 R 即 A QR 三 实验原理 1 LU 分解法 解 Ax b 的问题就等价于要求解两个三角形方程组 Ly b 求 y Ux y 求 x 设 A 为非奇异矩阵 且有分解式 A LU L 为单位下三角阵 U 为上三角阵 L U 的元素可以有 n 步直接计算定出 用直接三角分解法解 Ax b 要求 A 的所有顺序 主子式都不为零 的计算公式 i 2 3 n 2 1 niau lili 11 ual ilil 计算 U 的第 r 行 L 的第 r 列元素 i 2 3 n i r r 1 n 1 1 r k kirkriri ulau i r 1 n 且 r n rr r k krikirir uulal 1 1 求解 Ly b Ux y 的计算公式 3 2 1 1 11 niylby by i k kikii 1 2 1 1 nniuxuyx uyx ii n ik kikii nnnn 2 QR 分解法 4 实验步骤 1 LU 分解法 1 将矩阵 A 保存进计算机中 再定义 2 个空矩阵 L U 以便保存求出的三角矩阵的值 利用公式 将矩阵 A 分解为 LU L 为单位下三角阵 U 为上三角阵 2 可知计算方法有三层循环 先通过公式 计算出 U 矩阵的第一行元素 和 L 矩阵的第一列元素 li u il l 再根据公式 和 和上次的出的值 求出矩阵其余的元素 每次都要三次循环 求 下一个元素需要上一个结果 3 先由公式 Ly b 3 2 1 1 11 niylby by i k kikii 求出 y 因为 L 为下三角矩阵 所以由第一行开始求 y 4 再由公式 Ux y 1 2 1 1 nniuxuyx uyx ii n ik kikii nnnn 求出 x 因为 U 为上三角矩阵 所以由最后一行开始求 x 2 QR 分解法 5 程序流程图 1 LU 分解法 开始 输入系数矩阵 A 常数项b及n det 1 K 1 n 1 1 调选列主元子程序 i k 1 n 1 aik aik akk 即求乘数 mik j k 1 n 1 aik aij aik akj j bi bi aik bk i det akkdet k bn bn am 即求出xn i n 1 1 1 s 0 j i 1 n 1 s s aij bj j bi bi s aij 输出b及det 结束 i det amndet 回 代 过 程 消 元 过 程 由主程序转来 d akk p k i k 1 n 1 aik d d aik p i i d 0 P k j k n 1 t apj apj akj akj t j t bp bp bk bk t det det 返回主程序 返回主程序 输出失败标志 A 0 det 0 结束 选 列 主 元 2 QR 分解法 开始 输入数据 A b 矩阵A的转置 Householder变换 转置 矩阵相乘 输出上三角阵 R H3 H2 H1 A 输入正交阵Q H3 H2 H1 T 解上三角阵 输出结果 结束 6 实验结果 1 LU 分解法 2 QR 分解法 7 实验总结 为了求解线性方程组 我们通常需要一定的解法 其中一种解法就是通过矩阵的三角 分解来实现的 属于求解线性方程组的直接法 在不考虑舍入误差下 直接法可以用有限 的运算得到精确解 因此主要适用于求解中小型稠密的线性方程组 1 三角分解法 三角分解法是将 A 矩阵分解成一个上三角形矩阵 U 和一个下三角形矩阵 L 这样的 分解法又称为 LU 分解法 它的用途主要在简化一个大矩阵的行列式值的计算过程 求反 矩阵和求解联立方程组 不过要注意这种分解法所得到的上下三角形矩阵并非唯一 还可 找到数个不同 的一对上下三角形矩阵 此两三角形矩阵相乘也会得到原矩阵 2 QR 分解法 QR 分解法是将矩阵分解成一个正规正交矩阵 Q 与上三角形矩阵 R 所以称为 QR 分解 法 在编写这两个程序过程中 起初遇到不少麻烦 虽然课上老师反复重复着 算法不 难的 It s so easy 但是当自己实际操作时 感觉并不是那么容易 毕竟是要把实际的数 学问题转化为计算机能够识别的编程算法 所以在编写程序之前我们仔细认真的把所求解 的问题逐一进行详细的分析 最终转化为程序段 每当遇到问题时 大家或许有些郁闷 但最终还是静下心来反复仔细的琢磨 一一排除了错误 最终完成了本次实验 回头一想原来编个程序其实也没有想象的那么复杂 只要思路清晰 逐步分析 就可以 慢慢搞定了 通过这次实验 让我们认知到团队的作用力度是不容忽视的 以后不管干任何时都要 注重团队合作 遇到不懂得不明白的大家一起讨论 越讨论越清晰 愈接近最优答案 这 样不管干什么都能事半功倍 庆幸自己有这么个团队 也明白大家一起劳动的果实最珍贵 附 源代码 LR 分解法 include void input A 输入矩阵 A void input b 输入矩阵 b void output x float x 4 输出方程组的根 void main int i k r float A 4 4 4 2 1 5 8 7 2 10 4 8 3 6 12 6 11 20 b 4 2 7 7 3 x 4 u 4 4 l 4 4 y 4 给定的方程组 input A input b 显示矩阵 A b printf 矩阵 A 4 4 n for i 0 i 4 i for k 0 k 4 k printf 10f A i k printf n for i 0 i 4 i u 0 i A 0 i for i 0 i 4 i l i 0 A i 0 u 0 0 l i i 1 for r 1 r 4 r 计算矩阵 L 和 U for i r i 4 i float t 0 for k 0 k r k t t l r k u k i u r i A r i t for i r i 4 i float t1 0 for k 0 k r k t1 l i k u k r l i r A i r t1 u r r y 0 b 0 for i 1 i 4 i float t 0 for k 0 k 0 i float t 0 for k i 1 k 4 k t u i k x k x i y i t u i i printf n output x x 输入矩阵 A void input A float A 4 4 int i j printf input matrix A 4 4 n for i 0 i 4 i for j 0 j 4 j scanf f 输入矩阵 b void input b float b 4 int i printf input matrix b 4 n for i 0 i 4 i scanf f 输出方程的根 x void output x float x 4 int i printf 方程组的根为 n for i 0 i 4 i printf x d 13f i 1 x i printf n QR 分解法 include include define N 4 void matrix time double A N double B N double C N int n 两个矩阵相乘 结果存在矩阵 C N void matrix vec double A N double B N double C N int n 矩阵和向量相 乘 结果存在向量 C N double vec value double A int n 求向量的模 void vec time double a double H N int n 两个向量相乘得一个矩阵 void householder double a double H N int n int m 求解 Householder 矩阵函数 void matrix turn double A N int n 求矩阵装置 void solution double A N double b int n 求解上三角方程的解 void QR double A N double b int n void main double A 4 4 4 2 1 5 8 7 2 10 4 8 3 6 12 6 11 20 double b 4 2 7 7 3 int i j int n 4 printf 待求解的方程组系数矩阵 A 为 n for i 0 i n i 显示矩阵 A for j 0 j n j printf 10f A i j printf n QR A b n void matrix time double A N double B N double C N int n 两个矩阵相乘 结 果存在矩阵 C N int i j k for i 0 i n i for j 0 j n j C i j 0 for k 0 k n k C i j C i j A i k B k j void matrix vec double A N double B N double C N int n 矩阵和向量相乘 结果存 在向量 C N int i j for i 0 i n i C i 0 for j 0 j n j C i C i A i j B j double vec value double A int n 求向量的模 double Sum 0 int i for i 0 i n i Sum Sum A i A i Sum sqrt Sum return Sum void vec time double a double H N int n 两个向量相乘得一个矩阵 int i j for i 0 i n i for j 0 j n j H i j a i a j void householder double a double H N int n int m 计算 Householder 矩阵 int i j double first 存放首元素 double vec v 存放向量的模 double a1 N 0 for i 0 i n i for j 0 j n j if i j H i j 1 else H i j 0 first a m 取矩阵 A 转置的第 m 行 取矩阵 A 的第 m 列数 vec v vec value 计算 m 列元素所构成向量的模 a m a m vec v w 的分子部分 vec v vec value w 的分母部分 vec v vec v vec v for i m i n i for j m j n j H i j a i a j 2 vec v void matrix turn double A N int n 求矩阵的转置 double a N N 0 int i j for i 0 i n i for j 0 j n j a i j A i j for i 0 i n i for j 0 j 0 i for j n 1 j i j sum A i j x j x i b i sum A i i sum 0 printf 矩阵的 QR 分解求解结果为 n for i 0 i n i printf X d 10f n i 1 x i void QR double A N double b int n int i j k t double H1 N N 0 H2 N N 0 H3 N N 0 H1 H2 存放相邻两次的 Householder 矩阵 H3 为中间变量矩阵 double temp N N 0 double Q N N 0 R N N 0 Q1 N N 0 double b1 N 0 for i 0 i n i 将矩阵 A 临时存放在 temp 中 for j 0 j n j temp i j A i j for i 0 i n i H2 i i 1 单位阵 matrix turn temp n 矩阵 A 的转置 方便后续取 A 的某一列元素 for i 0 i n 1 i Householder 的一次变换 householder temp i H1 n i 得到 Householder 矩阵 H1 matrix time H1 A Q n 矩阵 H1 和 A 相乘放在 Q 中 for k 0 k