医学常见的微分方程的数值解法matlab(四阶龙格—库塔法).pptVIP

微分方程的数值解法 四阶龙格—库塔法 （The Fourth-Order Runge－Kutta Method） 常微分方程（Ordinary differential equations, ODE) 初值问题---给出初始值 边值问题---给出边界条件 一.解ODE的基本机理: 3. 根据式(2.2)编写计算导数的M函数文件-ODE文件 例题1：著名的Van der Pol方程 四阶Runge-Kutta公式 三. Runge-Kutta 法解Van der Pol 方程的Matlab 程序结构主程序：RK_vanderpol.m 子程序：RK_sub.m（函数文件） 解法2：采用Runge_Kutta法编程计算 主程序：RK_vanderpol.m t0=0; tN=20; y0=[0.25; 0]; h=0.001; t = t0 : h : tN; N = length (t); j = 1; for i = 1 : N t1 = 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; end subplot (121), plot (t, yy1, t, yy2); grid subplot (122), plot (yy2, yy1); grid 四. Matlab对应命令：ode23,ode45 五. 动力学系统的求解 1. 动力学方程 3. Matlab 程序(主程序：ZCX) t0;Y0;h;N;P0,w; %输入初始值、步长、迭代次数、初始激励力； for i = 1 : N t1 = t0 + h P=[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) next i 输出数据或图形 Matlab 程序(子程序：ZCX_sub.m) 例题2：三自由度质量弹簧系统 * * 与初值常微分方程解算有关的指令 ode23 ode45 ode113 ode23t ode15s ode23s ode23tb 2. 把高阶方程转换成一阶微分方程组 1. 列出微分方程 初始条件 令 (2.1) (2.2) (2.3) 例：著名的Van der Pol方程 令 降为一阶 初始条件 把t,Y作为输入宗量，把 作为输出宗量 %M function file name: dYdt.m  function Yd = f (t, Y) Yd = f (t,Y) 的展开式 例Van der Pol方程 %M function file name: dYdt.m  function Yd = f (t, Y) Yd=zeros(size(Y)); 4. 使编写好的ODE函数文件和初值 供微分方程解算指令（solver）调用 Solver解算指令的使用格式 [t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol); ode45 输出宗量形式 说明： t0：初始时刻；tN：终点时刻Y0：初值； tol：计算精度 %