专题一：MATLAB求解方程.ppt

MATLAB求解方程 制作 陈学明 基本概念 稀疏矩阵 sparsematrix 矩阵中非零元素的个数远远小于矩阵元素的总数 并且非零元素的分布没有规律 与之相区别的是 如果非零元素的分布存在规律 如上三角矩阵 下三角矩阵 对称矩阵 则称该矩阵为特殊矩阵 稠密矩阵 非0元素占所有元素比例较大的矩阵 若n阶矩阵A的行列式不为零 即 A 0 则称A为非奇异矩阵 否则称A为奇异矩阵 把矩阵A的行换成相应的列 得到的新矩阵称为A的转置矩阵 记作A AA E E为单位矩阵 或A A E 则n阶实矩阵A称为正交矩阵 基本概念 求矩阵A的秩rank A 求矩阵A的迹trace A 求矩阵A的行列式det A 求矩阵V的1范数norm V 1 求矩阵V的2范数norm V 或norm V 2 求矩阵V的 范数norm V inf 魔方矩阵 魔方矩阵是每行 每列及两条对角线上的元素和都相等的矩阵 对于n阶魔方阵 其元素由1 2 3 n2共n2个整数组成 magic n 生成一个n阶魔方阵 各行各列及两条对角线和为 1 2 3 n2 n 范得蒙矩阵 范得蒙 Vandermonde 矩阵最后一列全为1 倒数第二列为一个指定的向量 其他各列是其后列与倒数第二列的点乘积 vander V 生成以向量V为基础向量的范得蒙矩阵 例 A vander 1 2 3 5 A 11118421279311252551 希尔伯特矩阵 Hilbert矩阵的每个元素hij 1 i j 1 hilb n 生成n阶希尔伯特矩阵invhilb n 求n阶的希尔伯特矩阵的逆矩阵注意1 高阶Hilbert矩阵一般为病态矩阵 所以直接求逆可能出现错误结论 故应该采用invhilb n 注意2 由于Hilbert矩阵本身接近奇异 所以建议处理该矩阵时建议尽量采用符号运算工具箱 若采用数值解时应该验证结果的正确性 托普利兹矩阵 toeplitz toeplitz矩阵除第一行第一列外 其他每个元素都与左上角的元素相同 toeplitz x y 生成一个以x为第一列 y为第一行的托普利兹矩阵 这里x y均为向量 两者不必等长 toeplitz x 用向量x生成一个对称的托普利兹矩阵 例 T toeplitz 1 5 1 7 T 12345672123456321234543212345432123 帕斯卡矩阵 二次项 x y n展开后的系数随n的增大组成一个三角形表 称为杨辉三角形 由杨辉三角形表组成的矩阵称为帕斯卡 Pascal 矩阵 pascal n 生成一个n阶帕斯卡矩阵 矩阵分解 矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质 统称为简单矩阵 的组合 以便于理论分析或数值计算 通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积 求解线性方程组不可避免地要用到矩阵分解的概念 MATLAB中 线性方程组的求解主要用到三种基本的矩阵分解 即对称正定矩阵的cholesky分解 一般方程的gaussian消去法和矩阵的正交分解 这三种分解由函数chol lu和qr完成 矩阵的逆 求方阵A的逆可用inv A 说明1 对于Hilbert求逆时 不建议用inv 可用invhilb直接产生逆矩阵说明2 符号工具箱中也对符号矩阵定义了inv 函数 即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵说明3 对于奇异矩阵用数值方法的inv 函数 会给出错误的结果 但符号工具箱中的inv 会给出正确的结论 例A 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 det A B inv A A 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 A sym A inv A 矩阵的伪逆pinv A 若A不是一个方阵 或A是一个非满秩的方阵时 矩阵A没有逆矩阵 但可以找到一个与A的转置矩阵A 同型的矩阵B 使得 A B A A B A B B此时称矩阵B为矩阵A的伪逆 例求矩阵A的伪逆A 0 0 0 0 1 0 0 0 1 pinv A 矩阵分解 矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质 统称为简单矩阵 的组合 以便于理论分析或数值计算 通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积 求解线性方程组不可避免地要用到矩阵分解的概念 MATLAB中 线性方程组的求解主要用到三种基本的矩阵分解 即对称正定矩阵的cholesky分解 一般方程的gaussian消去法和矩阵的正交分解 这三种分解由函数chol lu和qr完成 矩阵分解之LU分解 矩阵的三角分解又称LU分解 它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积 即A LU Matlab使用函数lu实现LU分解 其格式为 L U lu A 其中U为上三角阵 L为下三角阵或其变换形式 满足LU A L U P lu A U为上三角阵 L为下三角阵 P为单位矩阵的行变换矩阵 满足LU PA 例 A 123 456 789 L U lu A L U P lu A 矩阵分解之Cholesky分解 对于正定矩阵A 则存在一个实的非奇异上三角阵R 满足R R A 称为Cholesky分解 Matlab使用函数chol实现Cholesky分解 其格式为 R chol A 若A非正定 则产生错误信息 R p chol A 不产生任何错误信息 若A为正定阵 则p 0 R与上相同 若A非正定 则p为正整数 R是有序的上三角阵 例 A pascal 4 R p chol A 矩阵分解之QR分解 将矩阵A分解成一个正交矩阵Q与一个上三角矩阵R的乘积A QR 称为QR分解 Matlab使用函数qr实现QR分解 其格式为 Q R qr A Q R E qr A 求得正交矩阵Q和上三角阵R E为单位矩阵的变换形式 R的对角线元素按大小降序排列 满足AE QR Q R qr A 0 例 A 123 456 789 101112 Q R qr A 代数方程的求解 一般的代数方程包括 线性方程非线性方程超越方程超越方程一般没有解析解 而只有数值解或近似解 只有特殊的超越方程才可以求出解析解来 matlab是获得数值解的一个最强大的工具 常用的命令有fsolve fzero等 但超越方程的解很难有精确的表达式 因此在matlab中常用eval 函数得到近似数值解 再用vpa 函数控制精度 求解线性方程 求解线性方程 对于Ax b A为m n矩阵 1 m n 恰定方程组 2 mn 超定方程 多用于曲线拟合 线性方程的组的求解 若系数矩阵的秩r n n为方程组中未知变量的个数 则有唯一解 若系数矩阵的秩r n 则可能有无穷解 求线性方程组的唯一解或特解 直接法 求线性方程组的通解线性方程组的无穷解 对应齐次方程组的通解 非齐次方程组的一个特解齐次线性方程组 常数项全部为零的线性方程组 直接法求线性方程组 求逆矩阵的思想 x inv A b或x A 1 b左除法 原理上是运用高斯消元法求解 但Matlab在实际执行过程中是通过LU分解法进行的 x A b符号矩阵法 此法最接近精确值 但运算速度慢 x sym A sym b 化为行最简形 C A b R rref C 则R的最后一列元素就是所求之解 如果无解 那么返回的是一个单位矩阵 例1 A 0 040 040 12 0 56 1 560 32 0 241 24 0 28 b 3 1 0 x11 inv A b 法1 求逆思想x12 A 1 b 法1 求逆思想x3 A b 法2 左除法x4 sym A sym b 法3 符号矩阵法C A b rref C 法4 化为行最简形 例 2 求下列方程组的解 A 5600015600015600015600015 B 10001 R A rank A 求秩X A B 求解 例3 求下列方程组的解 A 11 3 1 3 1 34 15 9 8 B 140 R A rank A 求秩X A B 由于系数矩阵不满秩 该解法可能存在误差X 00 0 53330 6000 一个特解近似值 A 11 3 1 3 1 34 15 9 8 B 140 C A B 构成增广矩阵R rref C 若用rref求解 则比较精确 LU分解求方程组的解 LU分解又称Gauss消去分解 可把任意方阵分解为下三角矩阵和上三角矩阵的乘积 即A LU L为下三角阵 U为上三角阵 则 A X b L U X b所以X U L b 这样可以大大提高运算速度 若采用 l u p lu A 则Ax b的解为x u l p b 利用LU分解求方程组的解例4 求下列方程组的一个特解 解 A 42 1 3 12 1130 B 2108 L U lu A X U L B 例5 利用LU分解法求解方程组 A 2426 49615 26918 6151840 先录入系数矩阵b 9232247 L U lu A 对A矩阵做LU分解y L b 求解方程组Ly bx U y 求解方程组Ux y得到方程组的最终解 故方程组的最终解为 x 0 5000 2 0000 3 0000 1 0000 或A 2426 49615 26918 6151840 先录入系数矩阵b 9232247 L U P lu A 对A矩阵做LU分解y L P b 求解方程组Ly Pbx U y 求解方程组Ux y得到方程组的最终解 Cholesky分解求方程组的解 若A为对称正定矩阵 则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积 即 A R R其中R为上三角阵 方程A X b变成R R X b所以X R R b QR分解求方程组的解 对于任何长方矩阵A 都可以进行QR分解 其中Q为正交矩阵 R为上三角矩阵的初等变换形式 即 A QR方程A X b变形成QRX b所以X R Q b 若采用 Q R E qr X 则Ax b的解为x E R Q b 解线性方程组的一般函数文件如下 function x y line solution A b m n size A y ifnorm b 0 b非零 非齐次方程组ifrank A rank A b 方程组相容 有解 ifrank A n 有唯一解x A b else 方程组有无穷多个解 基础解系disp 原方程组有有无穷个解 其齐次方程组的基础解系为y 特解为x y null A r x A b end else 方程组不相容 无解 给出最小二乘解disp 方程组的最小二乘法解是 x A b endelse 齐次方程组ifrank A n 列满秩x zero n 1 0解else 非0解 无穷多个解disp 方程组有无穷个解 基础解系为x x null A r endendreturn 如在MATLAB命令窗口 输入命令A 2 2 1 1 4 3 1 2 8 5 3 4 3 3 2 2 b 4 6 12 6 x y line solution A b 及 A 2 7 3 1 3 5 2 2 9 4 1 7 b 6 4 2 x y line solution A b 分别显示其求解结果 迭代法求解线性方程组 迭代法适用于解大型稀疏方程组 万阶以上的方程组 系数矩阵中零元素占很大比例 而非零元按某种模式分布 背景 电路分析 边值问题的数值解和数学物理方程常用迭代法 Jacobi迭代法Guass Seidel迭代法SOR迭代法 编写实现Jacobo迭代法的函数jacobi m如下 functions jacobi a b x0 err jacobi迭代法求解线性方程组 a为系数矩阵 b为 ax b右边的 矩阵b x0为迭代初值 err为迭代误差ifnargin 3 nargin判断输入变量个数的函数err 1 0e 6 elseifnargin 3errorreturnendD diag diag a 构造对角阵DD1 inv D 求对角阵D的逆矩阵L tril a 1 构造严格下三角阵U triu a 1 构造严格上三角阵B D1 L U 求出迭代矩阵 f D1 b s B x0 f n 1 n为迭代次数whilenorm s x0 errn n 1x0 s s B x0 f s endn A 9 1 1 110 1 1 115 b 7 8 13 x0 0 0 0 err 0 00005 s jacobi A b x0 err 结果 n10 77780 80000 866720 96300 96440 971930 99290 99350 995240 99870 99880 999150 99980 99980 999861 00001 00001 000071 00001 00001 0000 例 编写实现Seidel迭代法的函数seidel m如下 functions seidel a b x0 err seidel迭代法求解线性方程组 a为系数矩阵 b为 ax b右边的 矩阵b x0为迭代初值 err为迭代误差ifnargin 3err 1 0e 6 elseifnargin 3errorreturnendD diag diag a 构造对角阵DL tril a 1 构造严格下三角阵U triu a 1 构造严格上三角阵C inv D L B C U f C b n 1 n为迭代次数s B x0 fwhilenorm s x0 errn n 1x0 s s B x0 f s endn A 9 1 1 110 1 1 115 b 7 8 13 x0 0 0 0 err 0 00005 s seidel A b x0 err 结果 n10 77780 87780 977020 98390 99610 998730 99940 99980 999941 00001 00001 000051 00001 00001 0000 例 发现 seidel收敛速度比jacobi快