《常微分方程式培训教材》由会员分享,可在线阅读,更多相关《常微分方程式培训教材(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