Recall: Function handle and pseudo code

Lecture - ODE solver

Basic Matlab ODE solvers

Matlab 解常微分方程的指令中, 最常用的有以下兩個:

• ode45: 解 non-stiff 微分方程
• ode15s: 解 stiff 微分方程

其基本指令就是

[TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0)

其中

• ODEFUN 是個 function handle, 定義了要解的 ODE
• TSPAN = [T_0 T_FINAL], 定義另要解的時間範圍
• Y0 是初始值

也就是說, 這是在解一階常微分方程的初始值問題.

---

ODEFUN

要解一階常微分方程的初始值問題, 我們將方程寫成

$$\frac{dY}{dt}=F(t,Y),$$
其中

$Y=Y(t)$. 若

$Y$ 是個向量那就是個一階常微分方程系統.

接著我們將它寫成 matlab 函數:

dydt = myODE(t, Y)
% Input:  t - time, scalar
%         Y - nx1 vector
% Output: dydt: nx1 vector

dydt = F(t, Y)

定義出這個函數之後, 就可以把它以 @ 做成 function handle 輸入至 ode45 來解了.

---

van der Pol equation:

舉例來說, 我們想要解這個二階常微分方程的初始值問題:

$$y^{″}-\mu (1-y^2)y^{\prime }+y=0,y(0)=2,y^{\prime }(0)=0,$$
其中

$y=y(t)$,

$\mu$ 是個常數.

我們引進

$y_1=y$,

$y_2=y^{\prime }$, 則可以將上式改寫成兩個一階的常微分方程

$$\begin{gathered} y_1^{\prime }=y_2, \\ {y}_2^{\prime }=\mu (1-y_1^2)y_2-y_1, \end{gathered}$$
其初始值為

$y_1(0)=2$,

$y_2(0)=0$.

上面這個微分方程寫成函數就是以下這個樣子. (這是 matlab 內建的程式, 不需要另外存可以直接用, 其中

$\mu =1$)

function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%

dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

好了之後我想要以 ode45 來解,

$t\in [0,20]$. 這樣我需要在 matlab 的 command window 輸入

[t,y]=ode45(@vdp1,[0 20],[2 0]');   % 解上列 van der pol ODEs
plot(t,y(:,1));                    % 畫出解來

---

Remark:

1. 若我們直接在 command window 輸入
   ​​​​ode45(@vdp1,[0 20],[2 0]')
   

   會直接將
   $y_1$ 以及
   $y_2$ 畫出來, 不過沒有 output.
2. 若要畫在 phase plane (y(t), y'(t)) 上, 可以這樣做
   ​​​​[t,y]=ode45(@vdp1,[0 20],[2 0]');
   ​​​​plot(y(:,1), y(:,2))
   

---

Exercise

解以下的 Lorenz system

$$\begin{gathered} \frac{d{y}_1}{dt}=\sigma (y_2-y_1) \\ \frac{d{y}_2}{dt}=y_1(\rho -y_3)-y_2 \\ \frac{d{y}_3}{dt}=y_1{y}_2-\beta y_3, \end{gathered}$$
其中

$\rho =28$,

$\sigma =10$,

$\beta =8/3$. 系統初始條件為

$(y_1(0),y_2(0),y_3(0))=(0,1,0)$.

並以如下指令畫出其在三維的 phase plane

plot3(y(:,1), y(:,2), y(:,3))

---

Two point boundary value problems

ode45 只能解常微分方程初始值問題, 因此, 若我們要解的問題是邊界值問題, 那就需要用其他方式來做.

比如說我們要解以下 Fisher-KPP equation

$$y^{″}+y-y^2=0,y(0)=1,y(8)=0;$$

跟之前一樣我們引進

$y_1=y$,

$y_2=y^{\prime }$, 則可以將上式改寫成兩個一階的常微分方程

$$\begin{gathered} y_1^{\prime }=y_2, \\ {y}_2^{\prime }=-y_1+y_1^2, \end{gathered}$$
然後我們知道

$y_1(0)=1$.

Shooting method

假設我們令

$y_2(0)=\alpha$, 比如說

$\alpha =1$, 那這樣我們在

$t=0$ 就有兩個初始值, 就可以以 ode45 來解這問題, 得到

$y_1(t)$,

$y_2(t)$. 然後我們可以算

$y_1(8)$ 是多少, 假設

$y_1(8)=\beta$.

也就是說, 給定

$\alpha$, 我們就可以算出一個

$\beta =y_1(8)$. 因此我們可以將這想成是一個函數,

$f(\alpha )=\beta$.

而原問題是要找

$y(8)=0$, 也就是

$y_1(8)=0$, 也就是

$\beta =0$.
也就是說, 我們要找

$\alpha$, 使得

$\beta =0$.
也就是說, 這是個求根問題, 我們要找

$\alpha$, 使得

$f(\alpha )=0$.

Assignment

以 shooting method 解以下的 two-point boundary value problem:

$$y^{″}+y-y^2=0,y(0)=1,y(8)=0;$$