從物理、化學到生物, 反應擴散方程在近代科學中已被廣泛地用來描述各種現象。 擴散, 如我們下面第一節所談的, 似乎是一種「平凡化」的過程: 任何初始狀態, 經過長時間擴散, 總是達到處處常數的境界。因此, 在1952年, Alan Turing 用反應擴散方程來描述自然界中的種種 pattern formation, 是個劃時代的創舉!
原始形態學中最早的一個例子, 應該是 A. Trembley 在 1744 年發現的水螅 (hydra) 的再生現象:當水螅的頭部被切掉 48 小時後, 又能長出一個新的頭。 為了瞭解這一奇特的現象, Turing 在 1952 年提出了一個關於 activator-inhibitor 的數學模型。 在這個理論中, Turing 假設 activator 擴散的很慢但 inhibitor 擴散的很快, 而他充分利用這兩者之間擴散速度的差異, 給出一個「直觀」的論證, 來說明水螅頭部的再生現象。 他的觀點在 1972 年經過數值模擬, 終於被 formulated 成一個非線性方程組, 並引發了大量數學上的研究。
1998 年在 Notices of AMS 上筆者有一篇文章
在這篇短文中, 我們將由擴散的機制談起, 但我們的後續討論著重於「非自律」系統中, 反應、擴散與空間異質性 (spatial heterogeneity) 之間的交互作用。我們將以生態學中古老的 logistic equation 為例, 來闡明這三者之間的互動關係。
要特別一提的是, 雖然本文以引發讀者的興趣為主, 但文中所選的題材, 都是基本、 簡單而又重要的, 同時包含了當前研究前沿中較受關注的 open problems. 感興趣的讀者, 歡迎與筆者聯繫, 或直接嘗試解決這些問題。
1. 首先我們由熱方程 (heat equation) 談起
假設在空間 ${\Bbb R}^{n}$ 中我們有一有界區域 $\Omega$, 它光滑的邊界是由良好的絕熱材料製成。 若將 $\Omega$ 內部在點 $x$, 時間 $t\gt 0$ 的溫度, 記為 $u(x,t)$, 則 $u(x,t)$ 滿足下列方程 \begin{equation} \begin{cases} u_{t} =\Delta u \ \ &\hbox{in} \ \ \Omega \times (0,\infty), \\ u(x,0)= u_{0}(x) \ \ &\hbox{in} \ \ \Omega,\\ \partial_{u} u= 0 \ \ &\hbox{on} \ \ \partial \Omega \times (0,\infty), \end{cases} \end{equation} 此處 $\Delta = \sum^{n}_{i=1}\frac{\partial^2 }{\partial x_i^2}$ 是 Laplace 運算元, $u$ 代表 $\partial\Omega$ 上的單位外法線方向, 邊值條件 $ \partial_{u} u=0$ 即代表了這是個封閉系統: 無熱流入也無熱流出。 為了簡單起見, 在本文中, 我們將始終假設 $u_{0}\geq 0$ 且 $u_{0} \not\equiv 0$。
顯而易見的 (我們也已提過), 當我們讓時間 $t\rightarrow \infty$ 時, 不論初始值 $u_{0}(x)$ 如何分佈, $u(x,t)$ 必定趨於一個常數 $\overline{u_{0}}$ --- 即 $u_{0}(x)$ 在 $\Omega$ 上的平均值。 數學上來說, 我們有 \begin{equation} u(x,t)= \overline{u_{0}}+\Big(\int_{\Omega }u_0\psi_2 \Big)e^{-\mu_2 t}\psi _2(x)+\cdots \label{2} \end{equation} 其中 $ 0=\mu_1\lt \mu_2 \leq \mu_3 \leq \dots$ 是 $\Delta $ 的一串固有值(eigenvalues), 它們對應的正規化固有函數 (normalized eigenfunctions) 為 $ \psi_1 \equiv$ 常數, $\psi_2, \psi_3 \cdots$ \begin{equation} \begin{cases} \Delta \psi _{k}+\mu _{k} \mu _{k}=0 \quad &\text{ in $\Omega$,}\\ \partial _{u}\psi _{k}=0 \quad &\text{ on $\partial \Omega$.} \end{cases} \end{equation}
事實上, \eqref{2} 式不但告訴了我們解 $u(x,t) $ 趨於一個常數 $ \overline{u_{0}} $ (初始值 $u_{0}(x)$ 的平均值), 它還指出趨於 $ \overline{u_{0}}$ 的速度 (rate) 是由 $\Delta$ 的第一個非零的固有值 $\mu_2$ 所決定。
幸或不幸, 有趣的是, 到目前為止, 數學家們還沒有找到決定
$\mu_2$ 大小的有效方法, 但是, 許多人甚至以為決定 Neumann 固有值 $\mu_2$ 大小的方式與決定 Dirichlet 固有值大小的方式是一樣或類似的。
([pp.308-309]
Dirichlet固有值問題: \begin{equation} \begin{cases} \Delta \psi _{k}+\lambda _{k}\:\psi _{k}=0 \quad &\text{ in $\Omega$,}\\ \psi _{k}=0 \quad &\text{ on $\partial \Omega$ ,} \end{cases} \end{equation} 滿足一個簡單的性質:
「區域單調性 (Domain-Monotonicity Property)」: 若 $\Omega_{1} \subseteq \Omega_{2}$, 則 $\lambda_{k}(\Omega_{1}) \geq \lambda_{k}( \Omega_{2})$, 對所有 $k=1, 2,\ldots $ 都成立。
在 2007 年筆者與王學鋒的文章
在

由上圖我們看到, 當 $a$ 增加時, (-1, 0) 與(1, 0) 之間「在 $\Omega_a$ 之內的距離」是遞增的, 直觀上來說, 這意味著, 在 $\Omega_a$ 的內部 應該需要更長的時間, 來「平均」(-1, 0) 與 (1, 0) 附近的溫度或熱量之差異。
同樣的, 幸或不幸, 事情沒那麼簡單! 還有許多其它的幾何量在這裡也扮演了重要的角色。
在
如何決定 $\mu_2$ 的大小, 這個問題至今仍未能解決, 有興趣的讀者不妨參閱
綜上所述, 我們知道, 對於 Dirichlet 邊值問題, \begin{equation} \begin{cases} u_{t}=\Delta u \quad &\text{ in $\Omega \times (0, \infty)$ ,}\\ u(x, 0)\:=\:u_{0}(x) \quad &\text{ in $\Omega$ ,}\\ u=0 \quad &\text{ on $\partial\Omega\times (0, \infty)$ ,} \end{cases} \end{equation} 不論初值 $u_{0}$ 為何, 當 $t\rightarrow \infty$, 解 $u(x,t)$ 一定趨於 0, 而且滿足一個簡單的規律: 區域 $\Omega$ 愈大, $u(x,t)\rightarrow 0$ 的速度就愈慢!
對於 Neumann 邊值問題 (1), 情形就複雜多了, 而且最令人困惑的是, 到現在我們還不知道解 $u(x,t)\rightarrow \overline{u_0}$ 的速度依賴於那些 $\Omega$ 的幾何量, 這仍有待大家的繼續努力。
那麼, 全空間的解又如何呢? 為簡便計, 我們只考慮有界的初值$u_{0}$, 以及有界的解$u(x,t)$: \begin{equation} \begin{cases} u_{t}=\Delta u \quad &\text{ in $\mathbb{R}^{n} \times (0,\infty)$,} u(x, 0)=u_{0}(x) \quad &\text{ on $\mathbb{R}^{n}$.} \end{cases}\label{6} \end{equation}
自然, 首先我們要問: 當 $t\rightarrow \infty$ 時, 解 $u(x,t)$ 是否一定會收斂 ?
有趣的是, 與 Dirichlet 或 Neumann 邊值問題都不一樣, 在全空間, \eqref{6} 的解並不一定會收斂! 這個問題也有相當長的歷史。 我們在下一結果中, 給出了解收斂的充分必要條件。
proposition 1. 當 $t \rightarrow \infty$ 時, \eqref{6} 的解 $u(x,t)\rightarrow{u_{\infty}(x)}$ 若且唯若初始值 $u_{0}$ 滿足 $:$ $$ \lim_{R\rightarrow\infty} \frac{1}{B_{R}(x)} \int_{B_{R} (x)} u_{0}(y)dy \quad = \quad u_{\infty}(x), $$ 其中 $B_{R}(x)$ 是以 $x$ 為中心 $R$ 為半徑的球。而且 $u_{\infty}$ 必定得是常數。
對這個結果有興趣的讀者, 可以參閱
例: 令 $u_{0}\:\epsilon \: C^{\infty}(\mathbb{R})$ 為一偶函數滿足 $|u_{0}|_{L^{\infty}}=1 \quad $ 且 $$u_{0}(x) = (-1)^{k} \quad \hbox{for}\ x \: \epsilon \: [k!+2^{k}, (k+1)!-2^{k+1}].$$ 則 \eqref{6} 的唯一有界解 $u(x, t)$ 滿足 $$\liminf_{t\rightarrow\infty}\quad u(0, t) = -1 \quad \hbox{and} \quad \limsup_{t\rightarrow\infty} u(0,t) =1.$$
2. 當有反應項時, 情形就大不相同了
首先, 我們從最簡單的反應項談起: 不依賴於時間或空間的非線性反應項; 也就是說, 一個「自律(autonomous)」的方程 \begin{equation} \begin{cases} u_{t}=d\Delta u + f(u) \quad &\text{ in $\Omega \times (0, T)$ },\\ u(x, 0)\:=\:u_{0}(x) \quad &\text{ in $\Omega$ },\\ \partial_{u}u =0 \quad &\text{ on $\partial\Omega\times (0, T)$ },\\ \end{cases}\label{7} \end{equation}
此處 $d \gt 0$ 代表擴散係數, $0\lt T\leq\infty$ 是解存在的最大時間段, 對這最簡單的情形, \eqref{7} 的解即使是全域存在而且有界, 當 $t\rightarrow \infty$ 時, 解也不知道是否一定會收斂! (更多細節請參看 [N;Chapter 2])。
但是, 如果 \eqref{7} 的解是收斂的, 那麼它收斂的極限一定是 \eqref{7} 的平衡態 (equilibrium, steady state)。
在 1978 年, Casten 與 Holland
proposition 2. 當 $\Omega$ 是凸區域時, \eqref{7} 沒有局部穩定的非常數的平衡態。
換言之, 只要 $\Omega$ 是凸的, \eqref{7} 的任何穩定的平衡態, 必定是個常數, 且必須是反應項 $f(s)$ 的零點! 也就是說, 此時, 自律系統 \eqref{7} 是難以用來描述或刻劃 pattern formation 的。
非自律系統可以描述或刻畫的現象, 就非常豐富了; 比方說, 當反應項依賴於空間變數 $x \: \epsilon\: \Omega$ 時, 情形就已經大不相同。 在本節剩下的篇幅裡, 我們用一個古老而簡單的方程來闡明這個情形。
在大學二年級「常微方程」的課程裡, 我們都學過 logistic equation \begin{equation} u_{t} = ru(1-\frac{u}{K}),\label{8} \end{equation}
此處 $r\gt 0$, $K\gt 0$ 都是常數, 這個方程是生物數學家 Verhulst 在 1837 年提出來, 用來描述一個種群的人口總量。 $r$ 代表種群人口的增長率, $K$ 是環境資源所允許的最大人口總量。 此方程中的 $u=u(t)$ 代表在時間 $t\geq0$ 時的人口總量, 這個方程比 C. Darwin 的「進化論」還要早 22 年!
如果我們將 \eqref{8} 式寫成 $$u_{t} =\frac{r}{K}u(K-u),$$
再重新調整 $t$ 的大小, 吸收掉常數 ($ r/K$),就得到了
\begin{equation} u_{t} = u(K-u).\label{9} \end{equation}加入空間變量 $x\:\epsilon\:\Omega$, 並假設物種遷移是隨機的, 我們即有
\begin{equation} \begin{cases} u_{t}=d\Delta u + u(K-u) \quad &\text{in $\Omega \times (0,T)$ },\\ \partial_{u}u =0 \quad &\text{on $\partial\Omega\times (0,T)$ },\\ \end{cases}\label{10} \end{equation}此處$d \gt 0$ 是常數, 表示種群中個體遷移的強度。
數學上我們容易證明, 當 $K$ 是常數時, $V(x)\equiv\:K$ 是個全域穩定的平衡態; 也就是說, 不論 \eqref{10} 中$u$ 的初始值$u(x,0)$為何(只要不全為0), 最終當 $t\rightarrow \infty$ 時, $u(t,x)$ 一定趨於 $V(x)\equiv\:K$! 同樣的, 當 $K\not\equiv$ 常數時, 用上、 下解方法, 我們可以得到
proposition 3. $K \not\equiv$ 常數, 則 \eqref{10} 有唯一的正平衡態 $u_{d};$ 同時, $u_{d}$ 是全域穩定的。
2006年, 樓元
由 $K \not\equiv$ 常數, 馬上得到 $u_{d} \not\equiv$ 常數, 所以
\begin{equation} \int_{\Omega}K\lt \int_{\Omega}u_{d}.\label{11} \end{equation}這個簡單的不等式 \eqref{11} 卻有不尋常的生物意義, 它告訴我們: 總人口一定嚴格大於環境資源承載總容量 (total carrying capacity)! 這似乎違反了我們的直觀! 請不要忘記, \eqref{11} 式成立的兩個重要因素是: 種群要有遷移(dispersal), 不論大小; 再就是資源分佈必須是不均勻的(反映在「$ K \not\equiv$常數」上)!
不難證明(直觀上也非常清楚), 當$d\rightarrow\:0$ 時, $u_d$ 一致收斂到$K$, 而當 $d\rightarrow\: \infty$ 時, $u_d$ 一致收斂到 $\overline{K} \equiv\frac{1}{|\Omega|}\int_{\Omega}K(x)dx. $ 由此推得, 當 $d\rightarrow\:0$ 或 $\infty$ 時, \begin{equation} \int_{\Omega}u_{d}\rightarrow\int_{\Omega}K.\label{12} \end{equation}
\eqref{11} 式及 \eqref{12} 式表明了, 若將 $\int_{\Omega}u_{d}$ 看成 $d$ 的函數, 在某一個 $d_{0}$ 它會達到它的最大值。 自然, 我們要問:
(i) 如何決定 $d_{0}$?
(ii) $\int_{\Omega}u_{d_{0}}$ 到底能有多大 (與 $\int_{\Omega}K$ 相比)?
這是兩個簡單而基本的問題, 卻極具挑戰性, 到目前為止, 我們仍然不知道答案。更進一步, 令$$E(K)\equiv\max_{0\: \lt \: d \:\lt \:\infty }\frac{\int_{\Omega}u_{d}}{\int_{\Omega}K}.$$ 我們知道 $E(K)\gt 1$, 接下來考慮 $$\rho \: \equiv\: sup \left \{ \: E(K)|K\geq 0, \: \not\equiv 0\right \}.$$
問題: $\rho$ 是否有限? 如果是, 它的值是多少?
幾年前, 我有個猜測: 當 $n=1$ 時, $\rho=3$。
這個猜測最近被白學利、何小清和李芳
我們提到的這幾個簡單的數學問題有一定的生態學中的意義。 如問題 (i) 與 (ii), 它們的答案可以告訴我們, 當資源在空間分佈不均勻時, 種群應如何選擇它的擴散率, 以使得它的人口總量達到最大值 $E(K)$, 以及 $E(K)$ 的值; $\rho$ 的意義在於: 指出對於給定的資源總量, 如何選擇資源在空間的分佈, 使得 $E(K)$ 達到最大, 以及比較這個最大值是環境資源承載總容量 $\int_{\Omega}K$ 的若干倍。
關於 $u_{d}$ 與 $K$ 的關係, 還有許多有趣的性質與問題, 在此我們就不一一敘述了。
3. 異質環境中的物種競爭
\eqref{11} 這個簡單的不等式在異質環境中, 對兩個物種的競爭有深遠的影響。 假設兩個物種 $U$, $V$ 之間的競爭滿足經典的 Lotka-Volterra 方程組 \begin{equation} \begin{cases} U_{t}=d_{1}\Delta U + U(K(x)-U-c\: V) &\text{ in $\Omega \times (0,\infty)$ },\\ V_{t}=d_{2}\Delta V + V(K(x)-b\: U-V) &\text{ in $\Omega \times (0,\infty)$ },\\ U(x,0) = U_{0}(x),\: V(x, 0) = V_{0}(x)\quad &\text{ in $\Omega$ },\\ \partial_{u}U =\partial_{u}V =0 \quad &\text{ on $\partial\Omega\times (0,\infty)$ }. \end{cases}\label{13} \end{equation}
在 \eqref{13} 中我們假設了這兩種群 $U$, $V$ 需要完全相同的資源 $K(x)$, 而他們彼此的競爭能力以 $b$, $c$ 來表示。 當 $b$, $c$ 都介於 0 與 1 之間時, 我們說 $U$, $V$ 之間的競爭是「弱競爭」。
在弱競爭時, 如果資源在空間中是均勻分佈的 (即 ${K(x)}\equiv\: \text{常數}$), 一個眾所周知的事實是:
不論初始值 $U_{0}$, $V_{0}$ 的大小, $U$, $V$ 必定共存! 可以證明 (
但是當資源分佈不再均勻時, 樓元
其實, 早在 1998 年,
在最近的一篇文章
在何小清與筆者的一系列文章中
4.
現在讓我們回到原先的 logistic equation \eqref{8} 及其簡化後的 \eqref{9}。 從常微分方程的角度來看, \eqref{8} 與 \eqref{9} 基本上是沒有差異的。 但是在空間異質 (heterogeneous) 時, \eqref{8} 所對應的偏微分方程就應該是 \begin{equation} \begin{cases} u_{t}= d\:\Delta \: u +r(x)u(1-\dfrac{u}{K(x)})\quad &\text{ in $ \Omega \times (0, T)$ },\\ u(x, 0)=u_{0}(x) \quad &\text{ in $\Omega$ },\\ \partial _{u}u=0 \quad &\text{ on $\partial \Omega \times (0, T)$ }.\\ \end{cases}\label{15} \end{equation} 那麼, \eqref{15} 與 \eqref{10} 有差異嗎? 差異何在? 特別要關心的是, \eqref{15} 的平衡態是否也滿足 \eqref{11}呢?
首先, 對許多生態學家來說, 研究 \eqref{15} 是更為必要的, 因為他們相信內在成長速率(intrinsic growth rate) $r(x)$ 與 環境資源承載容量 (carrying capacity) $K(x)$ 是兩個不同的概念, 不同的量, 對一個種群是兩種不同的特質, 應該被區分開來, 但這並不表示二者之間沒有關係。 顯而易見的, 當 $r(x)\equiv\: cK(x)$ 時 ($c$ 是個正常數), \eqref{15} 與 \eqref{10} 幾乎是一樣的, 接下來, 我們來討論一般的情形。
標準的上下解方法可以用來證明: 對每一個 $d\gt 0$, \eqref{15} 具有唯一的一個全域穩定的平衡態 $u_{d}$, 在
下一個問題是, $L$ 與 $\overline{K}$ 孰大? 如果我們更進一步, 將 $r$ 設成 $K$ 的函數,
(i) 如果 $r(K)$ 是一個 凸函數, 則 $L \gt \overline{K}$;
(ii) 如果 $r(K)$ 是一個 凹函數, 則 $L \lt \overline{K}$.
一般來說, $r$ 依賴於許多變數, 而不是只依賴於 $K$, 但是如果將其它變數固定, 只讓 $K$ 變動, 從直觀上來看, 許多時候 $r$ 應該是 $K$ 的凹函數。 所以當 $d$ 很大時, \eqref{15} 的解通常不具備 \eqref{11} 式的性質。
在 § 3 中談到 Lotka-Voltera 競爭---擴散 (competition-diffusion) 的種群有趣現象大多依賴於 \eqref{11} 式, 因此, 如果我們考慮對應於 \eqref{15} 式的 Lotka-Voltera 競爭---擴散 (competition-diffusion), § 3 中談到的各種現象是否依然會發生? 在什麼條件下可以發生? 這些問題都是當下我與何小清正在探索的方向。
5.
我們已經談過 \eqref{11} 式對種群競爭的影響, 但是也不要忽略了 \eqref{11} 式本身在生態學中的意義: 在資源分佈不均的情況下, 對有隨機擴散的種群, 它能支援的人口總量超過它的環境承載容量!} 直觀上, 這似乎令人難以置信!
近年來, 有部分資深的生態學家注意到這個「純數學」預測的理論結果, 開始設計一系列的實驗來證實它。
在
DeAngelis 與他的團隊也已開始用實驗來探討 $r$ 與 $K$ 的關係, 我們十分期待瞭解 $r$ 與 $K$ 之間關係的實驗資料。
6.
反應擴散方程一直是極為豐富而且活躍的領域, 許多對於理解自然十分根本的問題, 仍懸而未解。 在這短短的篇幅中, 竟能包羅這麼多有趣甚或「出人意表」的結果, 似乎顯示出, 我們對於反應和擴散間交互作用的機制, 真正所知甚少。 這是個廣闊而吸引人的領域, 有待我們進一步深入探索。
後記
本文是基於作者於 2015 年 6 月初在臺灣清華大學對高年級 (大三、大四) 以上的學生 (包括研究生) 作的一次通俗演講的講稿整理而成。
感謝陳兆年教授的邀請, 宋宥松同學的筆記, 以及理論中心的資助。 特別感謝中研院的李宣北博士, 也是作者大學時代的老同學的邀稿、 建議與幫助, 以及助理編輯黃馨霈小姐的耐心, 本文得以在「數學傳播」上刊登。
最後, 還要感謝孫毓孜小姐的周到, 使新竹之行圓滿愉快!
參考資料
---本文作者任教中國華東師範大學及美國明尼蘇達大學---