《常微分方程式培训教材》由会员分享，可在线阅读，更多相关《常微分方程式培训教材（46页珍藏版）》请在金锄头文库上搜索。
1、MATLAB 程序设计进阶篇常微分方程式,張智星 jangcs.nthu.edu.tw http:/www.cs.nthu.edu.tw/jang 清大资工系 多媒体检索实验室,11-1 ODE 指令列表,MATLAB用于求解常微分方程式的指令:,ODE 指令列表,指令项目繁多，最主要可分两大类 适用于Nonstiff系统 一般的常微分方程式都是Nonstiff系统 直接选用ode45、ode23或ode113来求解 适用Stiff系统 速率（即微分值）差异相常大 使用一般的ode45、ode23或ode113来求解，可能会使得积分的步长（Step Sizes）变得很小，以便降低积分误差至可容
2、忍范围以内，会导致计算时间过长 专门对付Stiff系统的指令，例如ode15s、ode23s、ode23t及ode 23tb,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是输出的时间向
3、量 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
4、-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); ylabe
5、l(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 25, 3 3); plot(y(:,1), y(:,2), -o); xlabel(y(t); ylabel(y(t);,ODE 指令基本用法範例-3 (II),当 值越来越大时，van der Pol方程式就渐渐变成一个Stiff的系统，此时若要解此方程式，就必须改用专门
6、对付Stiff系统的指令,ODE 指令基本用法範例-4 (I),将 值改成1000，ODE档案改写成（vdp2.m）： type vdp2.m function dy = vdp2(t, y) mu = 1000; dy = y(2); mu*(1-y(1)2)*y(2)-y(1);,ODE 指令基本用法範例-4 (II),选用专门对付Stiff系统的ODE指令，例如ode15s，来求解此系统并作图显示 範例11-4：ode15s01.m,t, y= ode15s(vdp2, 0 3000, 2 1); subplot(2,1,1); plot(t, y(:,1), -o); xlabel(T
7、ime t); ylabel(y(t); subplot(2,1,2); plot(t, y(:,2), -o); xlabel(Time t); ylabel(y(t);% 注意單引號的重覆使用,ODE 指令基本用法範例-4 (III),由上图可知， 的变化相当剧烈（超过 ），这就是Stiff系统的特色,ODE 指令基本用法範例-5 (I),若要画出二维平面相位图，可以使用下列范例： 范例11-5：ode15s02.m 若要产生在某些特定时间点的状态变量值，则呼叫ODE指令的格式可改成： t, y = solver(odeFile, t0, t1, tn, y0); 其中t0, t1, tn
8、即是特定时间点所形成的向量,t, y= ode15s(vdp2, 0 3000, 2 1); subplot(1,1,1); plot(y(:, 1), y(:, 2), -o); xlabel(y(t); ylabel(y(t),ODE 指令基本用法範例-5 (II),11-3 ODE 指令的選項,ODE指令可以接受第四个输入变量，代表积分过程用到的各种选项（Options），此种ODE指令的格式为： t, y = solver(odeFile, t0, tn, y0, options); 其中options是由odeset指令来控制的结构变量 结构变量即包含了积分过程用到的各种选项 ode
9、set的一般格式如下： options = odeset(name1, value1, name2, value2,) 其字段name1的值为value1，字段name2的值为value2，依此类推 未被设定的字段，其域值即为默认值,ODE 指令的選項,也可以只改变一个现存的options结构变量中，某个字段的值，其格式如下： newOptions = odeset(options, name, value); 若要读出某个字段的值，可用odeget，其格式如下： value = odeget(otpions, name); 其中name为域名，value即为对应的域值 当使用odeset指令
10、时，若无任何输入变量，则odeset会显示所有的域名及域值，并以大括号代表默认值,ODE 指令的選項, odeset AbsTol: positive scalar or vector 1e-6 RelTol: positive scalar 1e-3 NormControl: on | off NonNegative: vector of integers OutputFcn: function_handle OutputSel: vector of integers Refine: positive integer Stats: on | off InitialStep: positive
11、 scalar MaxStep: positive scalar BDF: on | off MaxOrder: 1 | 2 | 3 | 4 | 5 Jacobian: matrix | function_handle JPattern: sparse matrix Vectorized: on | off Mass: matrix | function_handle MStateDependence: none | weak | strong MvPattern: sparse matrix MassSingular: yes | no | maybe InitialSlope: vecto
12、r Events: function_handle ,由 odeset 產生的 ODE 選項,由 odeset 產生的 ODE 選項,若F(t, y, Jacobian) 傳回,若 F( , , JPattern) 傳回,，且,由 odeset 產生的 ODE 選項,常用到的欄位來進行說明,f在积分误差容忍度方面，每一次积分所产生的局部误差e(i)，必须满足下列方程式： max(RelTol* , AbsTol(i) 其中i代表第i个状态变量 降低RelTol及AbsTol来求得更精确的积分结果,範例-1 (I),範例11-6：odeRelTol01.m,subplot(2,1,1); ode
13、45(vdp1, 0 25, 3 3); title(RelTol=0.01); options = odeset(RelTol, 0.00001); subplot(2,1,2); ode45(vdp1, 0 25, 3 3, options); title(RelTol=0.0001);,範例-1 (II),第一個圖所使用的相對誤差值是0.01（預設值），第二個圖所使用的相對誤差值是0.00001，因此我們得到較細密的點，但所花的計算時間也會比較長,積分輸出方面說明,积分输出的相关处理方面 选用一个OutputFcn 当ODE指令没有输出变量时，此输出函式OutputFcn会被MATLAB
14、呼叫 OutputFcn的默认值是”odeplot”，其功能为画出所有的状态变量 其它可用的函式 odephas2：画出2-D的相位平面（Phase Plane） odephas3：画出3-D的相位平面 odeprint：随时在指令窗口印出计算结果,積分輸出方面說明,以 Lorenz 渾沌方程式（Lorenz Chaotic Equation）為例 type lorenzOde.m function dy = lorenzOde(t, y) % LORENZODE: ODE file for Lorenz chaotic equation IGMA = 10.; RHO = 28; BETA
15、= 8./3.; A = -BETA 0 y(2) 0 -SIGMA SIGMA -y(2) RHO -1 ; dy = A*y;,範例-2,使用 ode45 對Lorenz 渾沌方程式進行數值積分 範例11-7：odeLorenz01.m 上列圖中，共有三條曲線，代表三個狀態變數隨時間變化的圖形,ode45(lorenzOde, 0 10, 20 5 -5);,範例-3,上述範例畫三度空間之相位圖形 範例11-8：odeLorenz02.m 圖形中只出現一條曲線，此曲線代表以三個狀態變數為座標、以時間為參數的一條三度空間中的曲線,options = odeset(OutputFcn, ode
16、phas3); % 使用 odephas3 進行繪圖 ode45(lorenzOde, 0 10, 20 5 -5, options);,提示,要觀看 Lorenz 渾沌方程式隨時間而變的動畫，可在 MATLAB 指令視窗下直接輸入lorenz 指令,範例-4 (I),假設 OutputFcn 設成“myfunc”： options = odeset(OutputFcn, myfunc) ODE 指令會呼叫 myfunc(tspan, y0, init) 讓 myfunc 進行各種初始化動作 積分步驟中，ODE 指令會持續呼叫status=myfunc(t, y) 若 status=1，則停止積分 積分結束時，ODE 指令會呼叫 myfunc( , , done)，讓 myfunc 進行收尾動作 Output
