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_0{e}^{\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)),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

Exercise

試描述下列方程其解的長時間行為

$$\frac{dX(t)}{dt}=X-X^3$$

---

2nd-order ordinary differential equation

Example - 彈簧震動

虎克(R. Hooke, 1635–1703)在研究彈簧的振動時, 寫下了有名的簡諧運動方程式

$$m\frac{d^{2}X(t)}{d{t}^2}+kX(t)=u(t),t>0,$$
其中

$m$ 是懸掛物質量,

$k$ 為彈簧彈性係數,

$X(t)$ 為彈簧在時間

$t$ 時的位移,

$u(t)$ 為外力.

給定初始位置跟速度,

$X(0)=X_0$,

$X^{\prime }(0)=Y_0$, 就可以解以上方程得到彈簧位移函數

$X(t)$.

若無外力,

$u(t)=0$, 此系統的解為

$$X(t)=C_1\cos \omega t+C_2\sin \omega t,\omega =\sqrt{k/m}.$$

---

一般而言, 二階常微分方程可以被改寫為兩個聯立的一階常微分方程, 例如以上簡諧運動, 我們引進

$Y(t)=X^{\prime }(t)$ 可得

$$\begin{gathered} \frac{dX(t)}{dt}=Y(t), \\ \frac{dY(t)}{dt}=-\frac{k}{m}X(t)+u(t), \end{gathered}$$

接著我們可以寫一簡單程式, 以 ode45 來解這微分方程.

我們先定義函數: (例如

$k=1$,

$m=2$,

$u(t)=\sin (t)$)

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$)

[t,y]=ode45(@shm,[0 100],[2 1]');   % 解上列 simple harmonic motion
plot(t,y(:,1));                    % 畫出解來

---

Autonomous system

若方程式無外力項, 或是獨立變數

$t$ 並沒有出現在方程式裡, 則稱為 Autonomous system. 這樣的二階方程可以寫為

$$\frac{d^{2}X}{d{t}^2}=f(X,X^{\prime })$$
例如無外力的簡諧運動其方程為

$$\frac{d^{2}X}{d{t}^2}=-{\omega }^{2}X,$$
而有阻尼的簡諧運動其方程為

$$\frac{d^{2}X}{d{t}^2}=-{\omega }^{2}X-\gamma X^{\prime }$$

二階 Autonomous sytem 可以畫出其相圖(phase portrait), 告訴我們這個系統其解運動的方向及行為, 以及決定系統的穩定性.

Matlab program for 2nd-order ODE

將 2nd-order autonomous system 改寫成 (引進

$Y=X^{\prime }$)

$$\begin{gathered} \frac{dX}{dt}=f_1(X,Y), \\ \frac{dY}{dt}=f_2(X,Y), \end{gathered}$$
例如無外力的簡諧運動其方程為

$$\begin{gathered} \frac{dX}{dt}=Y, \\ \frac{dY}{dt}=-{\omega }^{2}X, \end{gathered}$$
而有阻尼的簡諧運動其方程為

$$\begin{gathered} \frac{dX}{dt}=Y, \\ \frac{dY}{dt}=-{\omega }^{2}X-\gamma Y, \end{gathered}$$

• pplane

---

2nd-order boundary value problem

二階邊界值問題(boundary value problem)相對初始值問題(initial value problem)要困難許多. 一般會以 shooting method 或是將方程式離散化求根的方式來解.

• Lecture - ODE solver
• Lecture - fsolve

---

Higher order ordinary differential equation with initial condition

Example - 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)$.

---

Partial differential equation

Heat equation

$u=u(x,t)$ 描述在某個時間

$t$ 某個位置

$x$ 的溫度, 其方程為 (

$D$ 是常數)

$$u_t=D{u}_{xx}$$

wave equation

$u=u(x,t)$ 描述在某個時間

$t$ 某個位置

$x$ 波的高度, 其方程為 (

$c$ 是常數)

$$u_{tt}=c^2{u}_{xx}$$

• Lecture - 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+ϵY(t)$, 他必須滿足原方程式:

$\frac{d(X_p+ϵY(t))}{dt}=f(X_p+ϵY(t))$.

而我們有

$$\frac{d(Xp+ϵY)}{dt}=ϵ\frac{dY}{dt}$$

$$f(X_p+ϵY)\approx f(X_p)+ϵf^{\prime }(X_p)Y=ϵf^{\prime }(X_p)Y$$
因此可以得到

$$\frac{dY}{dt}=f^{\prime }(X_p)Y.$$
因此若

$f^{\prime }(X_p)>0$, 那一點點小擾動就會指數成長, 所以這不動點就是不穩定的.

反之, 若

$f^{\prime }(X_p)<0$, 那小擾動都會指數遞減至零, 所以這不動點就是穩定的.

Example

$$\frac{dX(t)}{dt}=X-X^3$$
這個系統

$f(X)=X-X^3$, 所以

$f^{\prime }(X)=1-3X^2$. 將三個不動點套入得到

$$f(0)=1,f(1)=-2,f(-1)=-2.$$
所以有兩個穩定的不動點,

$X=\pm 1$, 一個不穩定的不動點,

$X=0$.

---

以上這些推導也適用於 system of 1st-order ODE, 只是這時

$f^{\prime }(X_p)$ 會是個矩陣. 而判斷系統的穩定與否就要看

$e^{f^{\prime }(X_p)}$ 的行為. 所以就要看

$f^{\prime }(X_p)$ 的特徵值長怎樣. 如果有個特徵值大於零, 那

$e^{f^{\prime }(X_p)}$ 就會長大, 也就是一點點擾動就會變大, 所以不動點就是不穩定的.

所以通常這裡我們需要判斷一個矩陣,

$f^{\prime }(X_p)$, 他特徵值會不會有大於零的, 以及大於零的這些裡最大的是誰. 這時需要用到 shift-inverse power method 來求最大正數的特徵值.

• Lecture - Power method

---

Example - Allen-Cahn equation

$$\frac{\partial u(x,t)}{\partial t}=\delta u_{xx}+u-u^3,\delta \ll 1$$
這個系統有三個不動解,

$u(x,t)\equiv 0$,

$1$ 以及

$-1$.