微分方程的数值解法matlab(四阶龙格—库塔法)ppt课件

微分方程的数值解法 四阶龙格 库塔法 TheFourth OrderRunge KuttaMethod 常微分方程 Ordinarydifferentialequations ODE 初值问题 给出初始值边值问题 给出边界条件 与初值常微分方程解算有关的指令ode23ode45ode113ode23tode15sode23sode23tb 一 解ODE的基本机理 2 把高阶方程转换成一阶微分方程组 1 列出微分方程 初始条件 令 2 1 2 2 2 3 例 著名的VanderPol方程 令 降为一阶 初始条件 3 根据式 2 2 编写计算导数的M函数文件 ODE文件 把t Y作为输入宗量 把作为输出宗量 Mfunction dYdt mfunctionYd f t Y Yd f t Y 的展开式 例VanderPol方程 Mfunction dYdt mfunctionYd f t Y Yd zeros size Y 4 使编写好的ODE函数文件和初值供微分方程解算指令 solver 调用 Solver解算指令的使用格式 输出宗量形式 说明 t0 初始时刻 tN 终点时刻Y0 初值 tol 计算精度 例题1 著名的VanderPol方程 主程序 程序名 VanderPol ex1 m t0 0 tN 20 tol 1e 6 Y0 0 25 0 0 t Y ode45 dYdt t0 tN Y0 tol subplot 121 plot t Y subplot 122 plot Y 1 Y 2 解法1 采用ODE命令 VanderPol方程 子程序 程序名 dYdt m functionYdot dYdt t Y Ydot Y 2 Y 2 Y 1 2 1 Y 1 或写为 functionYdot dYdt t Y Ydot zeros size Y Ydot 1 Y 2 Ydot 2 Y 2 Y 1 2 1 Y 1 各种solver解算指令的特点 二 四阶Runge Kutta法 对I a b 作分割 步长 单步法 Runge Kutta方法多步法 Admas方法 计算的近似值时只用到 是自开始方法 Runge Kutta法是常微分方程的一种经典解法MATLAB对应命令 ode45 四阶Runge Kutta公式 四阶Runge Kutta法计算流程图 开始 Plot 初始条件 积分步长 迭代次数 输出结果 子程序计算 End 三 Runge Kutta法解VanderPol方程的Matlab程序结构主程序 RK vanderpol m子程序 RK sub m 函数文件 解法2 采用Runge Kutta法编程计算 主程序 RK vanderpol mt0 0 tN 20 y0 0 25 0 h 0 001 t t0 h tN N length t j 1 fori 1 Nt1 t0 h K1 RK sub t0 y0 K2 RK sub t0 h 2 y0 h K1 2 K3 RK sub t0 h 2 y0 h K2 2 K4 RK sub t0 h y0 h K3 y1 y0 h 6 K1 2 K2 2 K3 K4 yy1 j y1 1 yy2 j y1 2 t0 t1 y0 y1 j j 1 endsubplot 121 plot t yy1 t yy2 gridsubplot 122 plot yy2 yy1 grid 子程序 RK sub mfunctionydot vdpol t y ydot zeros size y ydot 1 y 2 ydot 2 y 2 y 1 2 1 y 1 或写为 ydot y 1 y 2 y 1 2 1 y 1 四 Matlab对应命令 ode23 ode45 调用格式 t y ode23 函数文件名 t0 tN y0 tol t y ode45 函数文件名 t0 tN y0 tol 默认精度 ode23 1e 3ode45 1e 6 说明 t0 初始时刻 tN 终点时刻y0 初值 tol 计算精度 3月15日作业 1 VanderPol方程的两种解法 1 采用ode45命令2 Runge Kutta方法2 Duffing方程的求解 Runge Kutta方法 计算步长h 0 005 计算时间t0 0 0 tN 100 要求 写出程序体 打印所绘图形 图形标题用个人的名字 Duffing方程 五 动力学系统的求解 1 动力学方程 2 二阶方程转成一阶方程 1 令 2 其中 即 2 3 Matlab程序 主程序 ZCX t0 Y0 h N P0 w 输入初始值 步长 迭代次数 初始激励力 fori 1 Nt1 t0 hP P0 sin w t0 0 0 0 0 输入t0时刻的外部激励力K1 ZCX sub t0 Y0 P P 输入 t0 h 2 时刻的外部激励力K2 ZCX sub t0 h 2 Y0 hK1 2 P K3 ZCX sub t0 h 2 Y0 hK2 2 P P 输入 t0 h 时刻的外部激励力K4 ZCX sub t0 h y0 hK3 P Y1 y0 h 6 K1 2K2 2K3 K4 t1 Y1 输出t1 y1 nexti输出数据或图形 Matlab程序 子程序 ZCX sub m functionydot f t Y P M K C 输入结构参数P1 zeros 3 1 inv M P A zeros 0 0 eye n n M 1K M 1C ydot AY P1 例题2 三自由度质量弹簧系统 矩阵表示 其中 动力学方程 解析解 已知参数 m1 m2 m3 1 k1 2 k2 2 K3 1 K4 2 P0 1 要求 采用四阶龙格 库塔法编程计算三个质量的响应时程 计算时间0 50 例如 4阶龙格 库塔法的结果 ode45的结果 第一个质量的位移响应时程 结果完全一致 MATLAB程序 1 4阶RK方法 2 采用ode45 m chap2 ex2 1 m m chap2 ex2 1 sub m 例题3 蹦极跳系统的动态仿真 蹦极者系着一根弹性绳从高处的桥梁 或山崖等 向下跳 在下落的过程中 蹦极者几乎处于失重状态 按照牛顿运动规律 自由下落的物体由下式确定 其中 m为人体的质量 g为重力加速度 x为物体的位置 第二项和第三项表示空气的阻力 其中位置x的基准为桥梁的基准面 即选择桥梁作为位置的起点x 0 低于桥梁的位置为正值 高于桥梁的位置为负值 如果人体系在一个弹性常数为k的弹性绳索上 定义绳索下端的初始位置为0 则其对落体位置的影响为 空气的阻力 整个蹦极系统的数学模型为 设桥梁距离地面为50m 即h2 50 蹦极者的起始位置为绳索的长度30m 即h1 30 蹦极者起始速度为0 其余的参数分别为k 20 a2 a1 1 m 70kg g 10m s2 初始条件 已知参数 初始条件变为 y0 30 0 初始位移和初始速度 t y ode45 bengji sub 0 0 01 100 y0 x1 50 y 1 x1代表蹦极者与地面之间的距离plot t x1 gridplot t y 1 grid y 1 代表位移 主程序 程序名 bengji m Matlab程序 functionydot f t y m 70 k 20 a1 1 a2 1 g 10 x y 1 x代表蹦极者的位移x dot y 2 x dot代表x的速度ifx 0ydot 0 1 k m a1 m a2 m abs x dot y 0 g elseydot 0 1 0 a1 m a2 m abs x dot y 0 g end 子程序 程序名 bengji sub m y 1 x1 结果分析 右上图为蹦极者与地面之间的距离 从结果可看出 对于体重为70kg