代数方程和常微分方程求解.ppt

第8章代数方程和常微分方程求解 代数方程是未知数和常数进行有限次代数运算所组成的方程 它包括有理方程和无理方程 代数方程的解称为的根或零点 其求解一般是通过代数几何来进行 微分方程是含有一个或是多个导数的方程 只有一个自变量及其导数的微分方程称为常微分方程 包含两个以上自变量及其偏导数的微分方程称为偏微分方程 工程上许多物理规律 设计过程的模拟和评价 凡是涉及质量和能量运动设计分析的问题 都最终归结到微分方程 8 1代数方程求解 8 1 1代数方程图解法符号绘图函数fplot 和ezplot 也可以用于图解法求代数方程的根 它适用于求解维数较少的一维方程或二维方程组 对于一维方程图解 其解就是函数曲线与x轴交点所对应的变量数值 如果有多个交点 则表示该方程有多个解 如果没有交点 则表示该方程没有解 例如 在例5 3使用符号绘图函数绘制代数方程的图形 图5 3左图 中可见 函数在区间 5 5 内与x轴有3个交点 因此该代数方程该区间内有3个实根 对于二维方程组图解 其解就是两条函数曲线的交点所对应的坐标数值 如果只有1个交点 或切点 则表示该方程组有1个解 如果有2个交点 则表示该方程组有2个解 如果没有交点 则表示该方程没有解 例8 1用图解法求解二维联立方程 a 2 b 2 定义横轴区间ezplot x 2 y 2 1 69 a b axis equal 控制坐标轴比例相等holdon gridon ezplot 2 4 x 3 y 1 5 a b line a b 0 0 line 0 0 b a xlabel bfx ylabel bfy title bf二维代数方程组的图解法 gtext bff 1 x 2 y 2 1 3 2 gtext bff 2 2 4x 3 y 1 5 8 1 2代数方程的解析解求非线性方程或方程组解析解的函数调用格式 X solve fun x 其中 fun是符号方程的函数表达式 x是自变量 X是解析解 应当指出 函数solve fun x 也可以用于求线性方组的解析解 例8 2求非线性解方程组解析解 二维非线性方程组的解析解symsabxy f1 x 2 x y a f2 y 2 x y b disp 二维非线性方程组的解析解 X Y solve f1 f2 x y M文件运行结果 二维非线性方程组的解析解 x a a b 1 2 a a b 1 2 Y 1 a b 1 2 b 1 a b 1 2 b 8 1 3线性方程组的数值解最简便方法是使用矩阵左除或是矩阵求逆的方法 求解线性方程组AX b X A bX inv A b其中 A是方程组的系数矩阵 b是常数向量 X是解析解 例8 3求线性方程组的数值解 线性方程组的数值解AA 1 1 1 1 1 2 1 4 2 3 1 5 3 1 2 11 bb 5 2 2 0 线性方程组常数向量disp 采用矩阵左除求出线性方程组的解 xx AA bbdisp 采用矩阵求逆求出线性方程组的解 zx inv AA bbdisp 计算残量 r AA zx bbdisp 计算残量的模 R norm r M文件运行结果 采用矩阵左除或矩阵求逆求出线性方程组的解 xx zx 1 00002 00003 0000 1 0000计算残量 r 1 0e 014 0 08880 2220 0 44410 1776计算残量的模 R 5 3475e 015 8 1 4非线性方程的数值解1 一维非线性方程对于一维非线性方程求解 可以看作是单变量的极小化问题 通过不断缩小搜索区间来逼近一维问题的真解 因此 可以使用一维非线性方程优化解函数来求解 其调用格式是 x fx flag fzero fun x0 其中 输入参数中 fun是非线性方程的函数表达式 x0是根的初值 输出参数中 x是非线性方程的数值解 fx是数值解的函数值 返回参数flag 0时 表示求解成功 否则求解出现问题 函数fzero所使用的算法为二分法 secant法和逆二次插值法的组合 例8 4求解一维非线性方程 求解单变量x非线性方程x0 0 1 解的初值 xz fz flag fzero atan x exp x x0 disp 求解成功性判断参数 flagdisp 非线性方程的解 xzdisp 非线性方程解的函数值 fzM文件运行结果 求解成功性判断参数 flag 1非线性方程的解 xz 0 6066非线性方程解的函数值 fz 1 1102e 016 2 多维非线性方程组对于多维非线性方程组使用多维非线性方程组优化解函数求解的格式 x fval flag fsolve fun x0 其中 输入参数中fun是非线性方程组的向量函数表达式 x0是根的初值 输出参数中x是非线性方程 组 的数值解 fval是数值解的函数值 返回参数flag 0时 表示求解成功 函数fsolve的作用是从根的初值x0开始 以逐渐减少误差的算法 搜索出满足多维非线性方程组fun的实根x和对应的函数值fval 例8 5求解二维非线性方程组x0 2 3 解的初值 定义非线性方程组表达式f和向量xfun inline 2 x 1 x 2 2 exp x 1 x 1 3 x 1 x 2 exp x 2 x xn fval flag fsolve fun x0 disp 求解成功性判断参数 flagifflag 0disp 方程组的解成功 elseifflag 0disp 方程组的解不成功 enddisp 非线性方程组的解 xndisp 非线性方程组解的函数值 fval M文件运行结果 求解成功性判断参数 flag 1方程组的解成功非线性方程组的解 xn 0 99781 2755非线性方程组解的函数值 fval 1 0e 006 0 1945 0 3372两个非线性方程解的函数值非常接近于0 说明其解的误差很小 且判断参数flag 0 解成功 8 2常微分方程求解 求解微分方程必须事先对自变量的某些值规定出函数或是导数的值 若在自变量为零的点上 给出初始条件 称为初值问题 最普遍的自变量是 时间 例如 弹性系统的自由振动 若以时间为零来限定位移和速度 这是一个初值问题 若在自变量为非零的点上 给出边界条件 称为边值问题 最普遍的自变量是 位移 例如 描述梁弯曲变形的微分方程 边界条件总是规定在梁的两端 8 2 1常微分方程的解析解在MATLAB中 用大写的字母D表示导数 Dy表示一阶导数 D2y表示二阶导 表示微分方程 在MATLAB中 由函数dsolve进行常微分方程 组 的求解问题 其调用格式是 y dsolve e1 e2 en c1 c2 cn v1 v2 vn 它用于求解常微分方程组e1 e2 en在初值条件c1 c2 cn下的特解 如果不给出初值条件 则求出解常微分方程组的通解 v1 v2 vn是求解变量 例9 1 求常微分方程的通解在此问题中 微分方程的MATLAB表达式为 Dy y 1 y 2 0没有给出初值条件 只能求出解常微分方程组的通解 调用求解微分方程组的语句为 y dsolve Dy y 1 y 2 0 x 运算结果 y 1 1 exp 2 x C1 1 2 1 1 exp 2 x C1 1 2 即该微分方程的解析解为 求常微分方程组的通解 没有给出初值条件 只能求出解常微分方程组的通解 调用求解微分方程组的语句为 x y dsolve Dx y exp z Dy exp 2 z z 运算结果 x cos z C2 sin z C1 1 5 exp 2 z 1 2 exp z y sin z C2 cos z C1 2 5 exp 2 z 1 2 exp z 其中 C1与C2是积分常数 即该微分方程的解析解为 求常微分方程组当和条件下的特解 在此问题中 两个微分方程的MATLAB表达式为 e1 Dx 2 x Dy 10 cos t e2 Dx Dy 2 y 4 exp 2 t 初值条件表达式为 C1 x 0 2C2 y 0 0 系统函数dsolve默认自变量为t 所以调用求解微分方程组的语句为 x y dsolve Dx 2 x Dy 10 cos t Dx Dy 2 y 4 exp 2 t x 0 2 y 0 0 运算结果 x 2 exp t sin t 4 cos t 3 sin t 2 exp 2 t y 2 exp t cos t sin t 2 cos t 即该微分方程的解析解为 8 2 2常微分方程的数值解MATLAB中用于求解一阶微分方程组初值问题数值解的常用函数是ode45 它是基于四阶五级Runge Kutta变步长算法 采用误差作为监测指标的闭环控制 既保证有较高的运算速度和中等计算精度 也保证预期的数值稳定性 是大部分场合求解非刚性常微分方程数值解的首选方法 求解非刚性常微分方程数值解的常用函数还有ode23 它是基于二阶三级Runge Kutta的单步算法 计算精度较低 函数ode45的调用格式是 t x ode45 fun t0 tf x0 其中 fun是描述微分方程的符号函数或是用inline 函数描述的微分方程 t0表示起始时刻 tf表示终止时刻 x0是初值问题的初始状态变量 例9 2试求洛伦茨 lorenz 系统动态模型状态方程的数值解 初值 终止时刻 其中 系数 描述洛伦茨 lorenz 系统动态模型状态方程3维非线性常微分方程组的函数文件functionxdot lorenzeq t x flag beta rho sigma 输入参数中 beta rho sigma是方程组有关项的系数 变量flag用于占位xdot beta x 1 x 2 x 3 rho x 2 rho x 3 x 1 x 2 sigma x 2 x 3 求解洛伦茨 lorenz 系统3维非线性常微分方程组tf 100 x0 0 0 1e 10 tf是终止时刻 beta 8 3 rho 10 sigma 28 给各系数赋值 t x ode45 lorenzeq 0 tf x0 beta rho sigma subplot 1 2 1 plot t x 绘制2维线图gridontitle bf状态变量的时间响应图 subplot 1 2 2 plot3 x 1 x 2 x 3 绘制3维线图axis 1042 2020 2025 title bf状态变量的相空间三维图 8 2 3高阶微分方程的降阶求解n阶微分方程一般形式 已知输出变量各阶导数初值选择一组状态变量可以把n阶微分方程化成等价的一阶微分方程组这是一个含n个未知函数的一阶微分方程组 初值 例9 3试求二阶微分方程的数值解 状态变量初值和 起始时刻 终止时刻 该二阶微分方程可以通过下面的变换 得到等价的一阶方程组 编制如下函数m文件 functiondy weifen t x dy zeros 2 1 dy 1 x 2 x 2 zdy 2 sin t x 2 x 1 x 1 y 然后用函数ode45求解微分方程 t y ode45 weifen 0 20 0 6 plot t y 1 holdon 绘制y的图像plot t y 2 r 绘制y 的图像xlabel bft ylabel bfy title bf二阶微分方程的数值解 legend 实线 y的图像 虚线 Dy的图像 运算结果如图8 3所示 例9 4试求二阶微分方程的数值解 该二阶微分方程可以通过下面的变换 得到等价的一阶方程组 已知状态变量的初值起始时刻 终止时刻 1 计算二阶微分方程的解析解symsty dsolve D2y 3 cos 2 t 2 sin t t 3 5 y 0 0 Dy 0 0 t ezplot y 010 绘制解析解xlabel bft ylabel bfy title bf二阶微分方程的解析解和数值解 gtext bfy 3cos2t 4 2sint t 3 6 7t 2 4 2t 3 4 holdon 2 计算二阶微分方程的数值解fun inline x 2 3 cos 2 t 2 sin t t 3 5 t x tx1 ode45 fun 0 10 0 0 plot t x1 1 r 绘制数值解gridonlegend 曲线 解析解 星号 数值解 运算结果如图8 4所示 8 2 4刚性微分方程的数值解刚性方程是指其Jacobian矩阵的特征值相差十分悬殊 在解的性态上表现为 其中一些解变化缓慢 另一些变化快 且相差较悬殊 这类方程常常称为刚性方程 又称为Stiff方程 刚性方程和非刚性方程对解法中步长选择的要求不同 刚性方程一般不适合由ode45这类函数求解 而应该采用ode15s等函数 如果不能分辨是否是刚性方程 先试用函数ode45 再用函数ode15s 求解刚性微分方程数值解常用函数有 ode15s 它是基于Gears反向数值微分的多步算法 具有中等计算精度 ode23s 它是基于Rosebrock的单步算法 计算精度较低 函数ode15s和ode23s的调用格式与函数ode45相同 例9 5求刚性微分方程组的数值解已知 初始条件和 起始时刻和终止时刻 使用函数dsolve可求出解析解为 可见 当时 两个分量均趋于0 下降极快 而下降很慢 编制描述刚性微分方程组的M文件 刚性微分方程组的函数文件functionfun stiff t y fun 0 01 y 1 99 99 y 2 100 y 2 调用函数ode15s求解 t y ode15s stiff 0 400 2 1 tstep length t 总步数minh min diff t 最小步长maxh max diff t 最大步长plot t y xlabel bft ylabel bfy title bf刚性微分方程组的数值解 gtext bfy 1 gtext bfy 2 8 2 5微分代数方程的数值解微分代数方程 differentialalgebraionequation DAE 是指微分方程中某些变量受到某些代数方程条件的约束 它的一般形式是 其中 描述的方法与普通的微分方程相同 为奇异矩阵 由M文件的求解选项变量Mass说明 在无法严格证明是奇异矩阵时 可以先假设它是奇异矩阵 如果在MATLAB求解过程中没有出现是非奇异矩阵的警告信息 则求解是有效的 否则无效 例9 6求微分代数方程组的数值解 已知 初始条件和 起始时刻和终止时刻 该方程组的第3个是代数方程 给定3个状态变量之间的约束关系 用矩阵形式表示该微分代数方程 编写描述微分代数方程的函数文件 functiondx DAE t x dx 0 2 x 1 x 2 x 3 0 3 x 1 x 2 2 x 1 x 2 5 x 2 x 3 2 x 2 x 2 x 1 x 2 x 3 1 Options是用命令odeset设置的可选积分参数 将常系数矩阵M赋给质量矩阵Mass M 1 0 0 0 1 0 0 0 0 x0 0 8 0 1 0 1 options odeset options Mass M t x ode15s DAE 0 20 x0 options plot t x 1 r t x 2 b t x 3 g title 微分代数方程的数值解 gtext bfx 1 t gtext bfx 2 t gtext bfx 3 t xlabel bft ylabel bfx gridon