matlab 常微分方程数值解法

科学计算与MATLAB 主讲 唐建国中南大学材料科学与工程学院2012 第七讲常微分方程数值解法 内容提要 引言欧拉近似方法龙格 库塔 R K 方法MATLAB的常微分方程函数小结 1 引言 物理学所研究的各种物质运动中 有许多物质运动的过程是用常微分方程来描述的 例如 质点的加速运动 简谐振动等 简单问题可以求得解析解 多数实际问题靠数值求解 一阶常微分方程 ODE 初值问题 数值解法就是求y x 在某些分立的节点xn上的近似值yn 用以近似y xn ODE OrdinaryDifferentialEquation 2 1简单欧拉 L Euler 1707 1783 方法 欧拉数值算法就是由初值通过递推求解 递推求解就是从初值开始 后一个函数值由前一个函数值得到 关键是构造递推公式 2 欧拉近似方法 欧拉数值算法递推公式构造 2 1 1差分法 差分法就是用差商近似代替微商 即 代入微分方程得到 对于等间隔节点 可以得到 在xn节点上 微分方程可以写为 作如下近似 则得到欧拉解法递推公式的一般形式 具体求解过程为 简单欧拉方法程序function outx outy MyEuler fun x0 xt y0 PointNum MyEuler用前向差分的欧拉方法解微分方程 fun表示f x y x0 xt表示自变量的初值和终值 y0表示函数在x0处的值 其可以为向量形式 PointNum表示自变量在 x0 xt 上取的点数ifnargin 5 PointNum 0 如果函数仅输入4个参数值 则PointNum默认值为100PointNum 100 endifnargin 4 y0默认值为0y0 0 endh xt x0 PointNum 计算步长hx x0 0 PointNum h 自变量数组y 1 y0 将输入存为行向量 输入为列向量形式fork 1 PointNumf feval fun x k y k 计算f x y 在每个迭代点的值f f y k 1 y k h f 对于所取的点x迭代计算y值endouty y outx x plot x y 画出方程解的函数图 例题 2 1 2折线法 几何意义 从 x0 y0 作曲线y x 的切线 得切线方程 此切线与x x1交点纵坐标为 y1 y y t 从 x1 y1 作曲线y x 在 x1 y x1 的切线的平行线 得直线方程 与x x2交点纵坐标为 y1 y y t y2 x1x2x3x4x5x6x7 按照相似的方法 从 xn yn 作曲线y x 在 xn y xn 的切线的平行线 得直线方程 与x xn 1交点纵坐标为 Q1 折线近似曲线 n增大 误差变大 x1x2x3x4x5x6x7 y y t 2 1 3数值积分 对微分方程作从x0到x积分得 在 x0 x1 积分采用矩形近似 得 同样 在 x0 x2 积分采用矩形近似 得 同样 在 x0 xn 1 积分采用矩形近似 得 作如下近似 得 2 1 4欧拉法误差 利用泰勒级数得 作如下近似 局部截断误差由泰勒展式可以看出 在欧拉方法中 假定yn是精确的 由yn计算yn 1所产生的误差的数量级为O h2 整体截断误差如果再考虑到yn本身的误差 计及误差的累积效果和局部截断误差 这种误差称为欧拉方法整体截断误差 数量级为O h 欧拉法是数值方法解微分方程最简单的一种 精度不高 实际应用不多 但它可以反映数值求解微分方程的特征 2 2改进的欧拉近似方法 节点 步长 对于一般的常微分方程 积分法求欧拉公式时 矩形法采用前一点函数值求积 若利用后一点 即为向后的欧拉方法 向后的欧拉方法递推公式为 即 向后的欧拉方法 隐式方法 预报 校正法 用欧拉方法预报用向后的欧拉方法校正 用欧拉方法预报 用向后的欧拉方法校正 迭代 积分法求欧拉公式时 积分采用梯形近似 即得到改进的欧拉方法 2 2 2改进的欧拉方法 改进的欧拉方法递推公式为 隐式方法 预报 校正法 用欧拉方法预报用改进的欧拉方法校正 计算公式为 例题 3 龙格 库塔 R K 方法 1欧拉方法 对于一般的常微分方程 几种方法的比较 取前一点的斜率作为平均斜率 即 整体截断误差 O h 2向后的欧拉方法 取后一点的斜率作为平均斜率 即 整体截断误差 O h 3改进的欧拉方法 取两点斜率平均值 即 整体截断误差 O h2 欧拉方法 前 后 和改进的欧拉方法优点是算法简单 但计算精度较低 不能满足实际问题求解对精度的要求 所以使用较少 龙格 库塔 R K 方法就是一种精度较高的较为实用计算常微分方程的方法 从以上三种方法的比较 可以推测 若取多点处斜率的加权平均作为平均斜率 精度会更高 这就是龙格 库塔法的基本思想 四阶龙格 库塔法计算公式为 整体截断误差 O h4 例题 time Euler 0 2500time EulerPro 0 5000time RK4 1 0150 4 MATLAB的常微分方程函数 ode45ode23ode113ode15sode23sode23tode23tb 格式 x y ode45 fun x0 xn y0 option 说明 适用于求解一阶常微分方程组fun定义微分方程组的函数文件名 x0 xn 求解区域y0初始条件向量option可选参数 由ODESET函数设置 比较复杂x输出自变量向量 y输出 y y y 没有一种算法可以有效地解决所有的ODE问题 因此MATLAB提供了多种ODE函数 例题 二 三阶龙格库塔法 例题 四 五阶龙格库塔法 t1 0 0310t2 0 0160 例题 用MATLAB的符号解法 求解常微分方程 dsolve Dy 3 x y x exp x 2 ans 1 3 exp x x 3 t C1 exp 3 x t dsolve Dy 3 x y x exp x 2 x ans exp x 2 exp 3 2 x 2 C1 例题 用MATLAB的符号解法 求解常微分方程的特解 dsolve x Dy 2 y exp x 0 y 1 2 exp 1 x ans exp x x exp x 2 exp 1 x 2 例题 采用ODE45求解描述某刚性问题的微分方程 functiondy rigid t y dy zeros 3 1 行向量dy 1 y 2 y 3 dy 2 y 1 y 3 dy 3 0 51 y 1 y 2 options odeset RelTol 1e 4 AbsTol 1e 41e 41e 5 T Y ode45 rigid 012 011 options plot T Y 1 T Y 2 T Y 3 legend y1 y2 y3 例题 用MATLAB函数ode23 ode45 ode113 求解多阶常微分方程 functiondy myfun03 x y dy zeros 3 1 初始化变量dydy 1 y 2 dy 1 表示y的一阶导数 其等于y的第二列值dy 2 y 3 dy 2 表示y的二阶导数dy 3 2 y 3 x 3 3 y 2 x 3 3 exp 2 x x 3 dy 3 表示y的三阶导数 用ode23ode45ode113解多阶微分方程clear clc x23 y23 ode23 myfun03 1 10 11030 x45 y45 ode45 myfun03 1 10 11030 x113 y113 ode113 myfun03 1 10 11030 figure 1 第一幅图plot x23 y23 1 r x45 y45 1 ob x113 y113 1 g 作出各种函数所得结果legend ode23解 ode45解 ode113解 title ODE函数求解结果 figure 2 plot x45 y45 以ode45为例作出函数以及其各阶导数图legend y y一阶导数 y两阶导数 title y y一阶导数 y二阶导数函数图 例题 MATLAB在导热计算 传热过程涉及面广 数学模型复杂 计算过程中涉及到许多运算方法 导热虽是容易做数学处理的一种热量传递方式 但其过程往往涉及常微分 偏微分方程 线性 非线性 方程组的求解 有一外径为4cm 内径为1 5cm 载有电流密度I为5000A cm2的内冷钢质导体 导体单位时间发出的热量等于流体同时带走的热量 导体内壁面的温度为70 假定外壁面完全绝热 试确定 导体内部的温度分布 已知钢的导热系数k 0 38Kw m K 电导率 2 1011k m 解 这是圆柱坐标中常物性一维稳态导热问题 结合本题具体条件导热微分方程式可简化为 结合边界条件 可得这一导热问题的数学描述为 此常微分方程的分析解 可调用MATLAB符号工具箱中的dsolve函数来实现 在命令窗口执行下面的代码 clear d equat D2t Dt r 131579 0 condition t0 0075 70 D t0 02 0 边界条件 t dsolve d equat condition r 程序执行结果 程序执行结果 即求出温度分布方程为 工程上遇到的导热问题 往往由于物体的几何形状复杂或边界条件难以描述 无法求出分析解 此时可采用数值方法进行求解 常微分方程 ODE 包括初值问题和边值问题两种 初值问题ODE的数值解法常调用ode45 和ode23 函数实现 在MATLAB编辑器中编写函数BZ m 存盘 functionBZclearallclca 0 0075 b 0 02 r值的范围solinit bvpini tlinspace a b 10 7080 用初始值对解初始化sol bvp4c ODEfun BCfun solinit 解