MATLAB的方程组解法.ppt

MATLAB的方程 组 解法 第六讲 王文健wwj527 MATLAB数据处理与应用 2011 2012学年选修课 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 2 主要内容 线性方程 组 的解法非线性方程 组 解法MATLAB统计分析 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 3 简介 线性方程 组 直接法 在没有舍入的情况下 通过有限步四则运算求得方程组的准确解迭代法 先给出一个解的初始值 然后按一定的法则逐步求出解的各个更准确的近似值的方法非线性方程 组 迭代法 不动点迭代法 Newton迭代法 BroydenBroyden迭代法 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 4 线性方程 组 的解法 线性方程的解法对于线性方程可以直接调用roots函数来求解例如 求解x3 5x2 8x 6 0P 15 86 roots P 求解2x9 43x7 x6 8x5 14x4 5x3 x2 10 x 12 0a 20431 814 51 1012 roots a TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 5 线性方程 组 的解法 线性方程组的解法直接法 利用符号运算 或 来求解直接法一般基于高斯消去法 主元素消去法 平方根法和追赶法等 在MATLAB中这些算法已被编成现成的库函数或运算符 可以直接调用进行求解 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 6 线性方程 组 的解法 线性方程组的解法建立系数矩阵和常数矩阵A 21 51 1 30 6 02 12 14 76 B 89 50 X A B TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 7 线性方程 组 的解法 线性方程组的解法利用矩阵的LU QR和Cholesky分解法求解 这三种方法对求解大型方程组非常有用优点是运算速度快 可以节省磁盘空间 节省内存LU分解法 Gauss消去分解 可以把任意矩阵分解为下三角矩阵的基本变换形式和上三角矩阵的乘积A LUL为下三角矩阵 U为上三角矩阵A X b变成L U X b X U L b TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 8 线性方程 组 的解法 LU分解法A 21 51 1 30 6 02 12 14 76 B 89 50 L U lu A X U L B TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 9 线性方程 组 的解法 Cholesky分解法如果A为对称正定矩阵 则Cholesky分解可将矩阵A分解为上三角矩阵和其转置的乘积 即A R R 其中R上三角矩阵A X b变成R R X bX R R b 举例 A 1648 45 4 8 422 B 28526 R chol A X R R B 正定矩阵设M是n阶实系数对称矩阵 如果对任何非零向量X x 1 x n 都有X MX 0 就称M正定 PositiveDefinite 正定矩阵在相合变换下可化为标准型 即单位矩阵 所有特征值大于零的对称矩阵 或厄米矩阵 也是正定矩阵 另一种定义 一种实对称矩阵 正定二次型f x1 x2 xn X AX的矩阵A A 称为正定矩阵 判定定理1 对称阵A为正定的充分必要条件是 A的特征值全为正 判定定理2 对称阵A为正定的充分必要条件是 A的各阶顺序主子式都为正 判定定理3 任意阵A为正定的充分必要条件是 A合同于单位阵 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 10 在A变换的作用下 向量 仅仅在尺度上变为原来的 倍 称 是A的一个特征向量 是对应的特征值 本征值 是 实验中 能测得出来的量 与之对应在量子力学理论中 很多量并不能得以测量 当然 其他理论领域也有这一现象 又称本征值 设A是向量空间的一个线性变换 如果空间中某一非零向量通过A变换后所奇异矩阵特征值得到的向量和X仅差一个常数因子 即AX kX 则称k为A的特征值 X称为A的属于特征值k的特征向量或特征矢量 eigenvector 如在求解薛定谔波动方程时 在波函数满足单值 有限 连续性和归一化条件下 势场中运动粒子的总能量 正 所必须取的特定值 这些值就是正的本征值 设M是n阶方阵 I是单位矩阵 如果存在一个数 使得M I是奇异矩阵 即不可逆矩阵 亦即行列式为零 那么 称为M的特征值 特征值的计算方法n阶方阵A的特征值 就是使齐次线性方程组 A I x 0有非零解的值 也就是满足方程组 A I 0的 都是矩阵A的特征值 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 11 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 12 线性方程 组 的解法 QR分解法对于任何长方矩阵A 都可以进行QR分解 其中Q为正交矩阵 R为上三角矩阵 即A QRA X b变成Q R X bX R Q b 举例 A 1648 45 4 8 422 B 28526 Q R qr A X R Q B 如果 AA E E为单位矩阵 A 表示 矩阵A的转置矩阵 或A A E 则n阶实矩阵A称为正交矩阵 若A为正交阵 则满足以下条件 1 A是正交矩阵2 AA E E为单位矩阵 3 A 是正交矩阵4 A的各行是单位向量且两两正交5 A的各列是单位向量且两两正交6 Ax Ay x y x y R正交矩阵通常用字母Q表示 举例 A r11r12r13 r21r22r23 r31r32r33 则有 r11 2 r12 2 r13 2 r21 2 r22 2 r23 2 r31 2 r32 2 r33 2 1r11 r12 r21 r22 r31 r32 0等性质正交方阵是欧氏空间中标准正交基到标准正交基的过渡矩阵 在矩阵论中 实数正交矩阵是方块矩阵Q 它的转置矩阵是它的逆矩阵 如果正交矩阵的行列式为 1 则我们称之为特殊正交矩阵 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 13 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 14 线性方程 组 的解法 迭代解法在MATLAB中迭代法非常适合求解大型系数矩阵的方程组Jacobi迭代法Gauss Serdel迭代法超松弛迭代法两步迭代法 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 15 线性方程 组 的解法 Jacobi迭代法线性方程组Ax b 如果系数矩阵A为非奇异矩阵 则A可写成A D L U的分解形式 其中D为对角矩阵 元素为A的对角元素 L和U分别为严格的下三角矩阵和上三角矩阵 故迭代形式为 xi 1 D 1 L U xi D 1b对应的Jacobi迭代公式的矩阵为 x k 1 Bx k f式中 B D 1 L U I D 1Af D 1b TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 16 线性方程 组 的解法 Jacobi迭代法举例 A 10 1 0 110 2 0 210 b 9 7 6 jacobi A b 0 0 1 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 17 线性方程 组 的解法 Gauss Seidel迭代法线性方程组Ax b Gauss Seidel迭代形式为 X k 1 GX k f式中 G为Gauss Seidel迭代矩阵 G D L 1b f D L 1b D为对角矩阵 L和U为严格的下三角和上三角矩阵编制Gauss Seidel迭代法的M文件 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 18 线性方程 组 的解法 Gauss Seidel迭代法举例 A 10 1 0 110 2 0 210 b 9 7 6 GS A b 0 0 0 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 19 线性方程 组 的解法 SOR迭代法线性方程组Ax b SOR迭代形式为 X k 1 B0X k f式中 B0为SOR迭代矩阵 B0 D wL 1 1 w D wU f w D wL 1bD为对角矩阵 L和U为严格的下三角和上三角矩阵编制SOR迭代法的M文件 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 20 线性方程 组 的解法 SOR迭代法举例 A 43 1 34 1 1 14 b 19 30 27 w 1 03 SOR A b 1 1 1 w TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 21 线性方程 组 的解法 双步迭代法线性方程组Ax b 如果A为对称正定矩阵 双步迭代法形式为 D L x k 1 2 Ux k b D U x k 1 Lx k 1 2 b式中 D为对角阵 L和U为严格的下三角和上三角矩阵编制SOR迭代法的M文件 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 22 线性方程 组 的解法 双步迭代法举例 A 10 122 111 13 2 1103 23 18 b 8 25 11 17 TS D A b 0 0 0 0 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 23 非线性方程的解法 对比法假定函数在 a b 区间上是单调函数 且f a f b 0 则在区间上方程f x 0至少有一个根思想 通过判断函数f x 的符号 逐步将有限区间缩小 使在足够小的区间内方程有且仅有一个根 近似根就是最终区间的重点位置举例 求解非线性方程xex 2 01 编制函数文件exam m2 命令窗口输入 DF exam 0 10 0 05 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 24 非线性方程的解法 迭代法不定点迭代法思想 非线性方程f x 0 如果在区间上连续 可将函数表示为x g x 迭代格式为 xk 1 g xk 如果limxk x 则迭代过程将收敛于方程的根举例 求解非线性方程ex 3x2 3 01 编制函数文件fun1 m2 命令窗口输入 x y n iterate 10 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 25 非线性方程的解法 迭代法切线迭代法 牛顿迭代法思想 迭代格式xk 1 xk f xk fg xk 给出导数举例 求解非线性方程ex 3x2 01 编制函数文件ff m2 编制导数文件dff m3 命令窗口输入 x n NT 2 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 26 非线性方程的解法 迭代法割线迭代法思想 当不存在解析的导函数时 割线法可近似来求解导函数 迭代格式xk 1 xk f xk xk xk 1 f xk f xk 1 举例 求解非线性方程ex 3x2 01 编制函数文件GX m2 编制函数文件ff m3 命令窗口输入 y GX 2 1 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 27 非线性方程的解法 求解非线性方程的MATALB函数求解非线性方程f x 0的函数fsolve调用格式 x fsolve fun x0 变量fun为方程表达式或函数文件名 x0为估计的方程根的初始值x fsolve fun x0 options 根据优化参数options进行最小化 默认值为0 如为1显示运算的中间步骤 x fval fsolve fun x0 返回值为目标函数在求得的根x对应的函数值 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 28 非线性方程的解法 求解非线性方程的MATALB函数求解非线性方程f x 0的函数fsolve调用格式 x fval exitflag fsolve 返回值为函数调用的退出标记exitflag 用于描述当前退出状态 x fval exitflag output fsolve 函数返回值为一个结构输出 用于描述优化的信息 x fval exitflag output jacobian fsolve 函数返回值为在点x处函数的Jacobian数值 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 29 非线性方程的解法 求解非线性方程的MATALB函数举例 求解f x xex 1 0的数值解命令窗口输入 xx 0 5 options 1 1 优化过程显示运算的中间步骤 xx fsolve x exp x 1 xx options x fval exitflag output Jacobian fsolve x exp x 1 xx TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 30 非线性方程的解法 求解非线性方程的MATALB函数fzero函数调用格式 x fzero fun x0 fun为求解的方程 x0为估计的方程的根或根可能存在的区间 如果x0为一个数值 此函数就在x0附近求出方程的根 如果x0为一个区间 且函数的根包含在该区间内则输出正常的方程根 否则输出错误信息x fzero fun x0 options 通过options指定的优化参数进行求解 x fval fzero fun x0 fval该函数返回值为方程根处的函数值 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 31 非线性方程的解法 求解非线性方程的MATALB函数fzero函数调用格式 x fval exitflag fzero exitflag返回函数fzero退出的函数值 x fval exitflag output fzero output返回优化的信息 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 32 非线性方程的解法 求解非线性方程的MATALB函数举例 求解f x 3x2 exp x 1 0的根命令窗口输入 x0 4 赋初始值fzero 3 x 2 exp x 1 x0 x1 410 区间未包含根 fzero 3 x 2 exp x 1 x1 x1 35 区间包含根fzero 3 x 2 exp x 1 x1 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 33 非线性方程组的解法 非线性方程组的求解求解更为复杂 繁琐 一般要采用迭代法进行求解不动点迭代法牛顿迭代法Broyden迭代法 TribologyResearchInstituteSOUTHWESTJIAOTONGUNIVERSITY 34 非线性方程组的解法 不动点迭代法思想 有n个未知数和n个方程的非线性方程组 F x 0 其迭代形式为