《MATLAB解方程》PPT课件.ppt

第7章MATLAB解方程与函数极值7 1线性方程组求解7 2非线性方程数值求解7 3常微分方程初值问题的数值解法7 4函数极值 7 1线性方程组求解7 1 1直接解法1 利用左除运算符的直接解法对于线性方程组Ax b 可以利用左除运算符 求解 x A b 例7 1用直接解法求解下列线性方程组 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 x A b 2 利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积 常见的矩阵分解有LU分解 QR分解 Cholesky分解 以及Schur分解 Hessenberg分解 奇异分解等 1 LU分解矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式 线性代数中已经证明 只要方阵A是非奇异的 LU分解总是可以进行的 MATLAB提供的lu函数用于对矩阵进行LU分解 其调用格式为 L U lu X 产生一个上三角阵U和一个变换形式的下三角阵L 行交换 使之满足X LU 注意 这里的矩阵X必须是方阵 L U P lu X 产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P 使之满足PX LU 当然矩阵X同样必须是方阵 实现LU分解后 线性方程组Ax b的解x U L b 或x U L Pb 这样可以大大提高运算速度 例7 2用LU分解求解例7 1中的线性方程组 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 L U lu A x U L b 或采用LU分解的第2种格式 命令如下 L U P lu A x U L P b 2 QR分解对矩阵X进行QR分解 就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式 QR分解只能对方阵进行 MATLAB的函数qr可用于对矩阵进行QR分解 其调用格式为 Q R qr X 产生一个一个正交矩阵Q和一个上三角矩阵R 使之满足X QR Q R E qr X 产生一个一个正交矩阵Q 一个上三角矩阵R以及一个置换矩阵E 使之满足XE QR 实现QR分解后 线性方程组Ax b的解x R Q b 或x E R Q b 例7 3用QR分解求解例7 1中的线性方程组 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 Q R qr A x R Q b 或采用QR分解的第2种格式 命令如下 Q R E qr A x E R Q b 3 Cholesky分解如果矩阵X是对称正定的 则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积 设上三角矩阵为R 则下三角矩阵为其转置 即X R R MATLAB函数chol X 用于对矩阵X进行Cholesky分解 其调用格式为 R chol X 产生一个上三角阵R 使R R X 若X为非对称正定 则输出一个出错信息 R p chol X 这个命令格式将不输出出错信息 当X为对称正定的 则p 0 R与上述格式得到的结果相同 否则p为一个正整数 如果X为满秩矩阵 则R为一个阶数为q p 1的上三角阵 且满足R R X 1 q 1 q 实现Cholesky分解后 线性方程组Ax b变成R Rx b 所以x R R b 例7 4用Cholesky分解求解例7 1中的线性方程组 命令如下 A 2 1 5 1 1 5 0 7 0 2 1 1 1 6 1 4 b 13 9 6 0 R chol A Errorusing cholMatrixmustbepositivedefinite命令执行时 出现错误信息 说明A为非正定矩阵 7 1 2迭代解法迭代解法非常适合求解大型系数矩阵的方程组 在数值分析中 迭代解法主要包括Jacobi迭代法 Gauss Serdel迭代法 超松弛迭代法和两步迭代法 1 Jacobi迭代法对于线性方程组Ax b 如果A为非奇异方阵 即aii 0 i 1 2 n 则可将A分解为A D L U 其中D为对角阵 其元素为A的对角元素 L与U为A的下三角阵和上三角阵 于是Ax b化为 x D 1 L U x D 1b与之对应的迭代公式为 x k 1 D 1 L U x k D 1b这就是Jacobi迭代公式 如果序列 x k 1 收敛于x 则x必是方程Ax b的解 Jacobi迭代法的MATLAB函数文件Jacobi m如下 function y n jacobi A b x0 eps ifnargin 3eps 1 0e 6 elseifnargin epsx0 y y B x0 f n n 1 end 例7 5用Jacobi迭代法求解下列线性方程组 设迭代初值为0 迭代精度为10 6 在命令中调用函数文件Jacobi m 命令如下 A 10 1 0 1 10 2 0 2 10 b 9 7 6 x n jacobi A b 0 0 0 1 0e 6 2 Gauss Serdel迭代法在Jacobi迭代过程中 计算时 已经得到 不必再用 即原来的迭代公式Dx k 1 L U x k b可以改进为Dx k 1 Lx k 1 Ux k b 于是得到 x k 1 D L 1Ux k D L 1b该式即为Gauss Serdel迭代公式 和Jacobi迭代相比 Gauss Serdel迭代用新分量代替旧分量 精度会高些 Gauss Serdel迭代法的MATLAB函数文件gauseidel m如下 function y n gauseidel A b x0 eps ifnargin 3eps 1 0e 6 elseifnargin epsx0 y y G x0 f n n 1 end 例7 6用Gauss Serdel迭代法求解下列线性方程组 设迭代初值为0 迭代精度为10 6 在命令中调用函数文件gauseidel m 命令如下 A 10 1 0 1 10 2 0 2 10 b 9 7 6 x n gauseidel A b 0 0 0 1 0e 6 例7 7分别用Jacobi迭代和Gauss Serdel迭代法求解下列线性方程组 看是否收敛 命令如下 a 1 2 2 1 1 1 2 2 1 b 9 7 6 x n jacobi a b 0 0 0 x n gauseidel a b 0 0 0 7 2非线性方程数值求解7 2 1单变量非线性方程求解在MATLAB中提供了一个fzero函数 可以用来求单变量非线性方程的根 该函数的调用格式为 z fzero fname x0 tol trace 其中fname是待求根的函数文件名 x0为搜索的起点 一个函数可能有多个根 但fzero函数只给出离x0最近的那个根 tol控制结果的相对精度 缺省时取tol eps trace 指定迭代信息是否在运算中显示 为1时显示 为0时不显示 缺省时取trace 0 例7 8求f x x 10 x 2 0在x0 0 5附近的根 步骤如下 1 建立函数文件funx m functionfx funx x fx x 10 x 2 2 调用fzero函数求根 z fzero funx 0 5 z 0 3758 7 2 2非线性方程组的求解对于非线性方程组F X 0 用fsolve函数求其数值解 fsolve函数的调用格式为 X fsolve fun X0 option 其中X为返回的解 fun是用于定义需求解的非线性方程组的函数文件名 X0是求根过程的初值 option为最优化工具箱的选项设定 最优化工具箱提供了20多个选项 用户可以使用optimset命令将它们显示出来 如果想改变其中某个选项 则可以调用optimset 函数来完成 例如 Display选项决定函数调用时中间结果的显示方式 其中 off 为不显示 iter 表示每步都显示 final 只显示最终结果 optimset Display off 将设定Display选项为 off 例7 9求下列非线性方程组在 0 5 0 5 附近的数值解 1 建立函数文件myfun m functionq myfun p x p 1 y p 2 q 1 x 0 6 sin x 0 3 cos y q 2 y 0 6 cos x 0 3 sin y 2 在给定的初值x0 0 5 y0 0 5下 调用fsolve函数求方程的根 x fsolve myfun 0 5 0 5 optimset Display off x 0 63540 3734 将求得的解代回原方程 可以检验结果是否正确 命令如下 q myfun x q 1 0e 009 0 23750 2957可见得到了较高精度的结果 7 3常微分方程初值问题的数值解法7 3 1龙格 库塔法简介7 3 2龙格 库塔法的实现基于龙格 库塔法 MATLAB提供了求常微分方程数值解的函数 一般调用格式为 t y ode23 fname tspan y0 t y ode45 fname tspan y0 其中fname是定义f t y 的函数文件名 该函数文件必须返回一个列向量 tspan形式为 t0 tf 表示求解区间 y0是初始状态列向量 t和y分别给出时间向量和相应的状态向量 例7 10设有初值问题 试求其数值解 并与精确解相比较 精确解为y t 1 建立函数文件funt m functionyp funt t y yp y 2 t 2 4 t 1 2 求解微分方程 t0 0 tf 10 y0 2 t y ode23 funt t0 tf y0 求数值解y1 sqrt t 1 1 求精确解t y y1 y为数值解 y1为精确值 显然两者近似 例7 11求解著名的VanderPol方程 例7 12有Lorenz模型的状态方程 试绘制系统相平面图 7 4函数极值MATLAB提供了基于单纯形算法求解函数极值的函数fmin和fmins 它们分别用于单变量函数和多变量函数的最小值 其调用格式为 x fmin fname x1 x2 x fmins fname x0 这两个函数的调用格式相似 其中fmin函数用于求单变量函数的最小值点 fname是被最小化的目标函数名 x1和x2限定自变量的取值范围 fmins函数用于求多变量函数的最小值点 x0是求解的初始值向量 MATLAB没有专门提供求函数最大值的函数 但只要注意到 f x 在区间 a b 上的最小值就是f x 在 a b 的最大值 所以fmin f x1 x2 返回函数f x 在区间 x1 x2 上的最大值 例7 13求f x x3 2x 5在 0 5 内的最小值点 1 建立函数文件mymin m functionfx mymin x fx x 3 2 x 5 2 调用fmin函数求最小值点 x fmin mymin 0 5 x 0 8165