[理学]11-常微分方程式.ppt

MATLAB 程序设计进阶篇常微分方程式 張智星 jang@.tw .tw/~jang 清大资工系 多媒体检索实验室 11-1 ODE 指令列表 MATLAB用于求解常微分方程式的指令:  指令 方法 適用ODE類別 ode45 Explicit Runge-Kutta (4, 5) pair of Dormand-Prince Nonstiff ODE ode23 Explicit Runge-Kutta (2, 3) pair of Bogacki and Shampine Nonstiff ODE ode113 Variable order Adams-Bashforth-Moulton PECE solver Nonstiff ODE ode15s Numerical differentiation formulas (NDFS) Stiff ODE ode23s Modified Rosenbrock formula of order 2 Stiff ODE ode23t Trapezoidal rule with a “free” interpolant Stiff ODE ode23tb Implicit Runge-Kutta formula with a backward differentiation formula of order two Stiff ODE ODE 指令列表 指令项目繁多，最主要可分两大类 适用于Nonstiff系统 一般的常微分方程式都是Nonstiff系统 直接选用ode45、ode23或ode113来求解 适用Stiff系统 速率（即微分值）差异相常大 使用一般的ode45、ode23或ode113来求解，可能会使得积分的步长（Step Sizes）变得很小，以便降低积分误差至可容忍范围以内，会导致计算时间过长 专门对付Stiff系统的指令，例如ode15s、ode23s、ode23t及ode 23tb 提示 使用 Simulink 來求解常微分方程式 Simulink是和MATLAB共同使用的一套软件 可使用拖拉的方式来建立动态系统 可直接产生C程序代码或进行动画显示 功能非常强大 11-2 ODE 指令基本用法 使用ODE指令时，必须先将要求解的ODE表示成一个函式 输入为t（时间）及y（状态变量，State Variables） 输出则为dy（状态变量的微分值） ODE函式的档名为odeFile.m，则呼叫ODE指令的格式如下： [t, y] = solver (odeFile, [t0, t1], y0); [t0, t1]是积分的时间区间 y0代表起始条件（Initial Conditions） solver是前表所列的各种ODE指令 t是输出的时间向量 y则是对应的状态变量向量 ODE 指令基本用法 以van der Pol微分方程为例，其方程式为： 化成标准格式 可用向量来表示成一般化的数学式 为一向量，代表状态变量 ODE 指令基本用法 假设=1，ODE档案（vdp1.m）可显示如下： type vdp1.m function dy = vdp1(t, y) mu = 1; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; 有了ODE档案，即可选用前述之ODE指令来求解 ODE 指令基本用法範例-1 (I) 在 =1 時，van der Pol 方程式并非Stiff系统，所以使用ode45来画出积分的结果 範例11-1：odeBasic01.m ode45(vdp1, [0 25], [3 3]); ODE 指令基本用法範例-1 (II) [0, 25] 代表积分的时间区间，[3 3]’ 则代表起始条件（必须以行向量来表示） 因为没有输出变量，所以上述程序执行结束后，MATLAB只会画出状态变量对时间的图形 ODE 指令基本用法範例-2 (I) 要取得积分所得的状态变量及对应的时间，可以加上输出变量以取得这些数据 范例11-2：odeGetData01.m [t, y] = ode45(vdp1, [0 25], [3 3]); plot(t, y(:,1), t, y(:,2), :); xlabel(Time t); ylabel(Solution y(t) and y(t)); legend(y(t), y(t)); ODE 指令基本用法範例-2 (II) ODE 指令基本用法範例-3 (I) 也可以画出 及 在 相位平面（Phase Plane）的运动情况 範例11-3：odePhastPlot01.m [t, y] = ode45(vdp1, [0