數學軟體實作 - ODE solver

--- title: 數學軟體實作 - ODE solver tags: 2020 Fall - 數學軟體實作 GA: G-77TT93X4N1 --- Recall: [Function handle and pseudo code](https://hackmd.io/@teshenglin/ms_pseudo) # 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' + y = 0, \, y(0)=2, \, y'(0)=0, $$ 其中 $y=y(t)$, $\mu$ 是個常數. 我們引進 $y_1=y$, $y_2=y'$, 則可以將上式改寫成兩個一階的常微分方程 $$ y_1' =y_2, \\ y_2' = \mu(1-y_1^2)y_2 - y_1, $$ 其初始值為 $y_1(0)=2$, $y_2(0)=0$. 上面這個微分方程寫成函數就是以下這個樣子. (這是 `matlab` 內建的程式, 不需要另外存可以直接用, 其中 $\mu=1$) ```matlab= 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 輸入 ```matlab= [t,y]=ode45(@vdp1,[0 20],[2 0]'); % 解上列 van der pol ODEs plot(t,y(:,1)); % 畫出解來 ``` --- #### Remark: 1. 若我們直接在 command window 輸入 ```matlab= ode45(@vdp1,[0 20],[2 0]') ``` 會直接將 $y_1$ 以及 $y_2$ 畫出來, 不過沒有 output. 3. 若要畫在 phase plane (y(t), y'(t)) 上, 可以這樣做 ```matlab= [t,y]=ode45(@vdp1,[0 20],[2 0]'); plot(y(:,1), y(:,2)) ``` --- #### Exercise 解以下的 Lorenz system $$ \frac{dy_1}{dt} = \sigma(y_2 - y_1)\\ \frac{dy_2}{dt} = y_1(\rho-y_3) - y_2\\ \frac{dy_3}{dt} = y_1y_2 - \beta y_3, $$ 其中 $\rho = 28$, $\sigma = 10$, $\beta = 8/3$. 系統初始條件為 $(y_1(0), y_2(0), y_3(0)) = (0, 1, 0)$. 並以如下指令畫出其在三維的 phase plane ```matlab= plot3(y(:,1), y(:,2), y(:,3)) ``` --- ## Two point boundary value problems `ode45` 只能解常微分方程初始值問題, 因此, 若我們要解的問題是邊界值問題, 那就需要用其他方式來做. 比如說我們要解以下 Fisher-KPP equation $$ y'' + y - y^2=0, \quad y(0)=1, \quad y(8)=0; $$ 跟之前一樣我們引進 $y_1=y$, $y_2=y'$, 則可以將上式改寫成兩個一階的常微分方程 $$ y_1' =y_2, \\ y_2' = -y_1 + y_1^2, $$ 然後我們知道 $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, \quad y(0)=1, \quad y(8)=0; $$