---
title: 數學軟體實作 - Differential equations
tags: 2020 Fall - 數學軟體實作
GA: G-77TT93X4N1
---
# Lecture - Differential equations
## Ordinary differential equation
### 1st-order ordinary differential equation
#### Example - 人口成長模型 1
> 假設人口的成長率與總人口數成正比
假設在時間 $t$ 的人口數為 $P(t)$, 人口成長率與人口數的比值為 $\lambda$, 則有
$$
\frac{dP(t)}{dt} = \lambda P(t)
$$
給定某地區的初始人口數 $P(t_0)=P_0$, 可以很容易得出
$$
P(t) = P_0e^{\lambda(t-t_0)}.
$$
所以根據這個模型人口會指數成長無上限.
#### Example - 人口成長模型 2
由於人口不可能會無上限成長, 因此
> 假設人口的成長不能超過某最大容量 $M$
並且有以下關係
$$
\frac{dP(t)}{dt} = \lambda P(t)(M-P(t))
$$
給定某地區的初始人口數 $P(t_0)=P_0$, 需要解以上微分方程才能得到人口成長曲線.
---
一般而言, 一階常微分方程可以寫成
$$
\frac{dX(t)}{dt} = f(X(t)), \quad X(t_0)=X_0,
$$
其中要解的函數為 $X(t)$, 而 $f$ 是一個函數.
例如以上兩例子
* Model 1: $f(X)=\lambda X$
* Model 2: $f(X)=\lambda X(M-X)$
#### Matlab program for 1st-order ODE
* [dfield](https://www.mathworks.com/matlabcentral/fileexchange/71408-dfield8)
#### Exercise
試描述下列方程其解的長時間行為
$$
\frac{dX(t)}{dt} = X - X^3
$$
---
### 2nd-order ordinary differential equation
#### Example - 彈簧震動
虎克(R. Hooke, 1635–1703)在研究彈簧的振動時, 寫下了有名的簡諧運動方程式
$$
m\frac{d^2X(t)}{dt^2} + kX(t) = u(t), \quad t>0,
$$
其中 $m$ 是懸掛物質量, $k$ 為彈簧彈性係數, $X(t)$ 為彈簧在時間 $t$ 時的位移, $u(t)$ 為外力.
給定初始位置跟速度, $X(0)=X_0$, $X'(0)=Y_0$, 就可以解以上方程得到彈簧位移函數 $X(t)$.
若無外力, $u(t)=0$, 此系統的解為
$$
X(t) = C_1\cos\omega t + C_2\sin\omega t, \quad
\omega=\sqrt{k/m}.
$$
---
一般而言, 二階常微分方程可以被改寫為兩個聯立的一階常微分方程, 例如以上簡諧運動, 我們引進 $Y(t) = X'(t)$ 可得
$$
\frac{dX(t)}{dt} = Y(t),\\
\frac{dY(t)}{dt} = -\frac{k}{m}X(t) + u(t),
$$
接著我們可以寫一簡單程式, 以 `ode45` 來解這微分方程.
我們先定義函數: (例如 $k=1$, $m=2$, $u(t)=\sin(t)$)
```matlab=
function dydt = shm(t,y)
% shm: simple harmonic motion
% 存成 shm.m
k=1;
m=2;
dydt = [y(2); -(k/m)*y(1)+u(t)];
function value = u(t)
value = sin(t);
end
end
```
接著以 `ode45` 來解: (假設要解的區間為 $t\in[0,100]$, 初始值 $X_0=2$, $Y_0=1$)
```matlab=
[t,y]=ode45(@shm,[0 100],[2 1]'); % 解上列 simple harmonic motion
plot(t,y(:,1)); % 畫出解來
```
---
#### Autonomous system
若方程式無外力項, 或是獨立變數 $t$ 並沒有出現在方程式裡, 則稱為 Autonomous system. 這樣的二階方程可以寫為
$$
\frac{d^2X}{dt^2} = f(X, X')
$$
例如無外力的簡諧運動其方程為
$$
\frac{d^2X}{dt^2} = -\omega^2X,
$$
而有阻尼的簡諧運動其方程為
$$
\frac{d^2X}{dt^2} = -\omega^2X-\gamma X'
$$
二階 Autonomous sytem 可以畫出其相圖([phase portrait](https://en.wikipedia.org/wiki/Phase_portrait)), 告訴我們這個系統其解運動的方向及行為, 以及決定系統的穩定性.
#### Matlab program for 2nd-order ODE
將 2nd-order autonomous system 改寫成 (引進 $Y=X'$)
$$
\frac{dX}{dt} = f_1(X, Y),\\
\frac{dY}{dt} = f_2(X, Y),
$$
例如無外力的簡諧運動其方程為
$$
\frac{dX}{dt} = Y,\\
\frac{dY}{dt} = -\omega^2X,
$$
而有阻尼的簡諧運動其方程為
$$
\frac{dX}{dt} = Y,\\
\frac{dY}{dt} = -\omega^2X-\gamma Y,
$$
* [pplane](https://www.mathworks.com/matlabcentral/fileexchange/61636-pplane)
---
### 2nd-order boundary value problem
二階邊界值問題(boundary value problem)相對初始值問題(initial value problem)要困難許多. 一般會以 shooting method 或是將方程式離散化求根的方式來解.
* [Lecture - ODE solver](https://hackmd.io/@teshenglin/ms_ode)
* [Lecture - fsolve](https://hackmd.io/@teshenglin/ms_fsolve)
---
### Higher order ordinary differential equation with initial condition
#### Example - 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)$.
---
## Partial differential equation
#### Heat equation
$u = u(x,t)$ 描述在某個時間 $t$ 某個位置 $x$ 的溫度, 其方程為 ($D$ 是常數)
$$
u_t = Du_{xx}
$$
#### wave equation
$u = u(x,t)$ 描述在某個時間 $t$ 某個位置 $x$ 波的高度, 其方程為 ($c$ 是常數)
$$
u_{tt} = c^2u_{xx}
$$
* [Lecture - Chebfun](https://hackmd.io/@teshenglin/ms_chebfun)
---
## Analysis
### Stationary solution for 1st-order system
#### Example - 1st-order ODE
$$
\frac{dX(t)}{dt} = f(X(t))
$$
如果有個數字 $X_p$ 滿足 $f(X_p)=0$, 那 $X(t)\equiv X_p$ 會是這系統的一個解, 而且是個不動的常數解.
要找出所有這些不動的常數解, 就是要解 $f(X_p)=0$, 也就是**求根問題**.
#### Example
$$
\frac{dX(t)}{dt} = X - X^3
$$
這個系統有三個不動點, $X=0$, $X=1$ 以及 $X=-1$.
### Linear stability of the stationary solution
找出這些不動的解之後, 我們會想知道他們的穩定性, 也就是若對他做一點點擾動會不會跑掉. 最常見的手法就是在這個點附近做線性化, 也就是做泰勒展開式.
假設對 $X_p$ 擾動一點點的解長這樣: $X(t) = X_p+\epsilon Y(t)$, 他必須滿足原方程式: $\frac{d(X_p+\epsilon Y(t))}{dt} = f(X_p+\epsilon Y(t))$.
而我們有
$$
\frac{d(Xp+\epsilon Y)}{dt} = \epsilon\frac{d Y}{dt}
$$
$$
f(X_p+\epsilon Y) \approx f(X_p) + \epsilon f'(X_p)Y = \epsilon f'(X_p)Y
$$
因此可以得到
$$
\frac{dY}{dt} = f'(X_p)Y.
$$
因此若 $f'(X_p)>0$, 那一點點小擾動就會指數成長, 所以這不動點就是不穩定的.
反之, 若 $f'(X_p)<0$, 那小擾動都會指數遞減至零, 所以這不動點就是穩定的.
#### Example
$$
\frac{dX(t)}{dt} = X - X^3
$$
這個系統 $f(X)=X-X^3$, 所以 $f'(X) = 1-3X^2$. 將三個不動點套入得到
$$
f(0)=1, \quad f(1)=-2, \quad f(-1)=-2.
$$
所以有兩個穩定的不動點, $X=\pm 1$, 一個不穩定的不動點, $X=0$.
---
以上這些推導也適用於 system of 1st-order ODE, 只是這時 $f'(X_p)$ 會是個矩陣. 而判斷系統的穩定與否就要看 $e^{f'(X_p)}$ 的行為. 所以就要看 $f'(X_p)$ 的特徵值長怎樣. 如果有個特徵值大於零, 那 $e^{f'(X_p)}$ 就會長大, 也就是一點點擾動就會變大, 所以不動點就是不穩定的.
所以通常這裡我們需要判斷一個矩陣, $f'(X_p)$, 他特徵值會不會有大於零的, 以及大於零的這些裡最大的是誰. 這時需要用到 shift-inverse power method 來求**最大正數**的特徵值.
* [Lecture - Power method](https://hackmd.io/@teshenglin/ms_power)
---
#### Example - Allen-Cahn equation
$$
\frac{\partial u(x,t)}{\partial t} = \delta u_{xx} + u - u^3, \quad \delta\ll 1
$$
這個系統有三個不動解, $u(x,t)\equiv 0$, $1$ 以及 $-1$.
Sign in with Wallet
Connect another wallet