用MATLAB解常微分方程

实验四实验四 求微分方程的解求微分方程的解 一 问题背景与实验目的一 问题背景与实验目的 实际应用问题通过数学建模所归纳而得到的方程 绝大多数都是微分方程 真正能得到代数方程的机会很少 另一方面 能够求解的微分方程也是十分有 限的 特别是高阶方程和偏微分方程 组 这就要求我们必须研究微分方程 组 的解法 既要研究微分方程 组 的解析解法 精确解 更要研究微分 方程 组 的数值解法 近似解 对微分方程 组 的解析解法 精确解 Matlab 有专门的函数可以用 本 实验将作一定的介绍 本实验将主要研究微分方程 组 的数值解法 近似解 重点介绍 Euler 折 线法 二 相关函数 命令 及简介二 相关函数 命令 及简介 1 dsolve equ1 equ2 Matlab 求微分方程的解析 解 equ1 equ2 为方程 或条件 写方程 或条件 时用 Dy 表示表示 y 关关 于自变量的一阶导数于自变量的一阶导数 用用用用 D2y 表示表示 y 关于自变量的二阶导数 依此类推 关于自变量的二阶导数 依此类推 2 simplify s 对表达式 s 使用 maple 的化简规则进行化简 例如 syms x simplify sin x 2 cos x 2 ans 1 3 r how simple s 由于 Matlab 提供了多种化简规则 simple 命令就 是对表达式 s 用各种规则进行化简 然后用 r 返回最简形式 how 返回形成 这种形式所用的规则 例如 syms x r how simple cos x 2 sin x 2 r cos 2 x how combine 4 T Y solver odefun tspan y0 求微分方程的数值解 说明 1 其中的 solver 为命令为命令 ode45 ode23 ode113 ode15s ode23s ode23t ode23tb 之一 之一 2 odefun 是显式常微分方程 00 yty ytf dt dy 3 在积分区间 tspan 上 从到 用初始条件求解 0f tt 0 t f t 0 y 4 要获得问题在其他指定时间点上的解 则令 tspan 210 ttt 要求是单调的 210f tttt 5 因为没有一种算法可以有效地解决所有的 ODE 问题 为此 Matlab 提供了多种求解器 Solver 对于不同的 ODE 问题 采用不同的 Solver 求解器求解器 Solver ODE 类型类型特点特点说明说明 ode45非刚性 单步算法 4 5 阶 Runge Kutta 方程 累计截断误差达 3 x 大部分场合的首选算 法 ode23非刚性单步算法 2 3 阶 Runge Kutta 方程 累计截断误差达 3 x 使用于精度较低的情 形 ode113非刚性多步法 Adams 算法 高低精 度均可到 63 10 10 计算时间比 ode45 短 ode23t适度刚性采用梯形算法适度刚性情形 ode15s刚性多步法 Gear s 反向数值微分 精度中等 若 ode45 失效时 可 尝试使用 ode23s刚性单步法 2 阶 Rosebrock 算法 低精度 当精度较低时 计算 时间比 ode15s 短 ode23tb刚性梯形算法 低精度当精度较低时 计算 时间比 ode15s 短 6 要特别的是 ode23 ode45 是极其常用的用来求解非刚性的标准形式 的一阶常微分方程 组 的初值问题的解的 Matlab 的常用程序 其中 ode23 采用龙格 库塔 2 阶算法 用 3 阶公式作误差估计来调节步长 具有 低等的精度 ode45 则采用龙格 库塔 4 阶算法 用 5 阶公式作误差估计来调节步长 具 有中等的精度 5 ezplot x y tmin tmax 符号函数的作图命令 x y 为关于参数 t 的符 号函数 tmin tmax 为 t 的取值范围 6 inline 建立一个内联函数 建立一个内联函数 格式 inline expr var1 var2 注意 括号里的表达式要加引号 例 Q dblquad inline y sin x pi 2 pi 0 pi 三 实验内容三 实验内容 1 几个可以直接用 Matlab 求微分方程精确解的例子 例 1 求解微分方程 并加以验证 2 2 x xexy dx dy 求解本问题的 Matlab 程序为 syms x y line1 y dsolve Dy 2 x y x exp x 2 x line2 diff y x 2 x y x exp x 2 line3 simplify diff y x 2 x y x exp x 2 line4 说明 1 行 line1 是用命令定义定义 x y 为符号变量为符号变量 这里可以不写 但为确保正确 性 建议写上 2 行 line2 是用命令求出的微分方程的解 1 2 exp x 2 x 2 exp x 2 C1 3 行 line3 使用所求得的解 这里是将解代入原微分方程将解代入原微分方程 结果应该为 0 但这里给出 x 3 exp x 2 2 x exp x 2 C1 2 x 1 2 exp x 2 x 2 exp x 2 C1 4 行 line4 用 simplify 函数对上式进行化简 结果为 0 表明 的确是微分方程的解 xyy 例 2 求微分方程在初始条件下的特解 并画出解函0 x eyxyey2 1 数的图形 求解本问题的 Matlab 程序为 syms x y y dsolve x Dy y exp x 0 y 1 2 exp 1 x ezplot y 微分方程的特解为 y 1 x exp x 1 x exp 1 Matlab 格式 即 x ee y x 解函数的图形如图 1 6 4 20246 30 20 10 0 10 20 30 40 50 x 1 x exp x 1 x exp 1 图 1 例 3 求微分方程组在初始条件下的特解 03 5 yx dt dy eyx dt dx t 0 1 00 tt yx 并画出解函数的图形 求解本问题的 Matlab 程序为 syms x y t x y dsolve Dx 5 x y exp t Dy x 3 y 0 x 0 1 y 0 0 t simple x simple y ezplot x y 0 1 3 axis auto 微分方程的特解 式子特别长 以及解函数的图形均略 2 用 ode23 ode45 等求解非刚性的标准形式的一阶常微分方程 组 的初值 问题的数值解 近似解 例 4 求解微分方程初值问题的数值解 求解范围为 1 0 222 2 y xxy dx dy 区间 0 0 5 fun inline 2 y 2 x 2 2 x x y x y ode23 fun 0 0 5 1 x y plot x y o x ans 0 0000 0 0400 0 0900 0 1400 0 1900 0 2400 0 2900 0 3400 0 3900 0 4400 0 4900 0 5000 y ans 1 0000 0 9247 0 8434 0 7754 0 7199 0 6764 0 6440 0 6222 0 6105 0 6084 0 6154 0 6179 图形结果为图 2 00 050 10 150 20 250 30 350 40 450 5 0 6 0 65 0 7 0 75 0 8 0 85 0 9 0 95 1 图 2 例 5 求解描述振荡器的经典的 Ver der Pol 微分方程 7 0 0 1 0 0 1 2 2 2 yyy dt dy y dt yd 分析 令则 1 21 dt dx xyx 1 12 2 1 2 2 1 xxx dt dx x dt dx 先编写函数文件 verderpol m function xprime verderpol t x global mu xprime x 2 mu 1 x 1 2 x 2 x 1 再编写命令文件 vdp1 m global mu mu 7 y0 1 0 t x ode45 verderpol 0 40 y0 x1 x 1 x2 x 2 plot t x1 图形结果为图 3 0510152025303540 2 5 2 1 5 1 0 5 0 0 5 1 1 5 2 2 5 图 3 3 用 Euler 折线法求解 前面讲到过 能够求解的微分方程也是十分有限的 下面介绍用 Euler 折 线法求微分方程的数值解 近似解 的方法 Euler 折线法求解的基本思想是将微分方程初值问题 00 yxy yxf dx dy 化成一个代数方程 即差分方程 主要步骤是用差商替代微商 h xyhxy 于是 dx dy 00 xyy xyxf h xyhxy kk kk 记 从而 则有 1kkkk xyyhxx 1 hxyy kk 1 2 1 0 1 1 00 nk yxhfyy hxx xyy kkkk kk 例 6 用 Euler 折线法求解微分方程初值问题 1 0 2 2 y y x y dx dy 的数值解 步长 h 取 0 4 求解范围为区间 0 2 解 本问题的差分方程为 1 2 1 0 2 4 0 1 0 2 1 1 00 nk y x yyxfyxhfyy hxx hyx kkkk kk 其中 相应的 Matlab 程序见附录 1 数据结果为 0 1 0000 0 4000 1 4000 0 8000 2 1233 1 2000 3 1145 1 6000 4 4593 2 0000 6 3074 图形结果见图 4 00 20 40 60 811 21 41 61 82 1 2 3 4 5 6 7 图 4 特别说明 本问题可进一步利用四阶 Runge Kutta 法求解 读者可将两个 结果在一个图中显示 并和精确值比较 看看哪个更 精确 相应的 Matlab 程序参见附录 2 四 自己动手四 自己动手 1 求微分方程的通解 0sin2 1 2 xxyyx 2 求微分方程的通解 xeyyy x sin5 2 3 求微分方程组 0 0 yx dt dy yx dt dx 在初始条件下的特解 并画出解函数的图形 0 1 00 tt yx yf x 4 分别用 ode23 ode45 求上述第 3 题中的微分方程初值问题的数值解 近似解 求解区间为 利用画图来比较两种求解器之间的差异 0 2 t 5 用 Euler 折线法求解微分方程初值问题 1 0 12 3 2 y y x yy 的数值解 步长 h 取 0 1 求解范围为区间 0 2 6 用四阶 Runge Kutta 法求解微分方程初值问题 1 0 cos y xeyy x 的数值解 步长 h 取 0 1 求解范围为区间 0 3 四阶 Runge Kutta 法的迭代公式为 Euler 折线法实为一阶 Runge Kutta 法 1 2 1 0 2 2 2 2 22 6 34 23 12 1 43211 1 00 nk hLyhxfL L h y h xfL L h y h xfL yxfL LLLL h yy hxx xyy kk kk kk kk kk kk 相应的 Matlab 程序参见附录 2 试用该方法求解第 5 题中的初值问题 7 用 ode45 方法求上述第 6 题的常微分方程初值问题的数值解 近似解 从而利用画图来比较两者间的差异 五 附录五 附录 附录 1 fulu1 m clear f sym y 2 x y 2 a 0 b 2 h 0 4 n b a h 1 x 0 y 1 szj x y for i 1 n 1 y y h subs f x y x y x x h szj szj x y end szj plot szj 1 szj