淺談反應擴散方程

淺談反應擴散方程

倪維明

從物理、 化學到生物, 反應擴散方程在近代科學中已被廣泛地用來描述各種現象。 擴散, 如 我們下面第一節所談的, 似乎是一種 「平凡化」 的過程: 任何初始狀態, 經過長時間擴散, 總是 達到處處常數的境界。 因此, 在 1952 年, Alan Turing 用反應擴散方程來描述自然界中的種種 pattern formation, 是個劃時代的創舉!

原始形態學中最早的一個例子, 應該是 A. Trembley 在 1744 年發現的水螅 (hydra) 的 再生現象: 當水螅的頭部被切掉 48 小時後, 又能長出一個新的頭。 為了瞭解這一奇特的現象, Turing 在 1952 年提出了一個關於 activator-inhibitor 的數學模型。 在這個理論中, Turing 假設 activator 擴散的很慢但 inhibitor 擴散的很快, 而他充分利用這兩者之間擴散速度的差 異, 給出一個 「直觀」 的論證, 來說明水螅頭部的再生現象。 他的觀點在 1972 年經過數值模擬, 終於被 formulated 成一個非線性方程組, 並引發了大量數學上的研究。

1998 年在 Notices of AMS 上筆者有一篇文章 [13], 介紹了在 「自律 (autonomous)」 系 統中, 反應與擴散的互動及其所產生的各種不同的機制, 包括了有趣的 「凝聚 (condensation)」

現象, 例如: 尖峰解 (peak solutions), 也簡要說明了上述 Alan Turning 原創性觀點的數學 證明。

在這篇短文中, 我們將由擴散的機制談起, 但我們的後續討論著重於 「非自律」 系統中, 反 應、 擴散與空間異質性 (spatial heterogeneity) 之間的交互作用。 我們將以生態學中古老的 logistic equation 為例, 來闡明這三者之間的互動關係。

要特別一提的是, 雖然本文以引發讀者的興趣為主, 但文中所選的題材, 都是基本、 簡單而 又重要的, 同時包含了當前研究前沿中較受關注的 open problems. 感興趣的讀者, 歡迎與筆 者聯繫, 或直接嘗試解決這些問題。

1. 首先我們由熱方程 (heat equation) 談起

假設在空間 Rⁿ 中我們有一有界區域 Ω, 它光滑的邊界是由良好的絕熱材料製成。 若將 Ω

17

內部在點 x, 時間 t > 0 的溫度, 記為 u(x, t), 則 u(x, t) 滿足下列方程











ut= ∆u in Ω × (0, ∞), u(x, 0) = u0(x) in Ω,

∂νu = 0 on ∂Ω × (0, ∞),

(1)

此處 ∆ = Pn i=1

∂²

∂x²i 是 Laplace 運算元, ν 代表 ∂Ω 上的單位外法線方向, 邊值條件 ∂νu = 0 即代表了這是個封閉系統: 無熱流入也無熱流出。 為了簡單起見, 在本文中, 我們將始終假設 u0 ≥ 0 且 u0 6≡ 0。

顯而易見的 (我們也已提過), 當我們讓時間 t → ∞ 時, 不論初始值 u0(x) 如何分佈, u(x, t) 必定趨於一個常數 u0 — 即 u0(x) 在 Ω 上的平均值。 數學上來說, 我們有

u(x, t) = u0+Z

Ω

u0ψ2

e^−µ2tψ2(x) + · · · (2)

其中 0 = µ₁ < µ2 ≤ µ3 ≤ . . . 是 ∆ 的一串固有值(eigenvalues), 它們對應的正規化固有函 數 (normalized eigenfunctions) 為 ψ₁ ≡ 常數, ψ₂, ψ₃· · ·







∆ψk+ µkµk = 0 in Ω,

∂νψk= 0 on ∂Ω. (3)

事實上, (2) 式不但告訴了我們解 u(x, t) 趨於一個常數 u₀ (初始值 u₀(x) 的平均值), 它 還指出趨於 u0 的速度 (rate) 是由 ∆ 的第一個非零的固有值 µ2 所決定。

幸或不幸, 有趣的是, 到目前為止, 數學家們還沒有找到決定 µ2 大小的有效方法, 但是, 許多人甚至以為決定 Neumann 固有值 µ2 大小的方式與決定 Dirichlet 固有值大小的方式是 一樣或類似的。 ([15, pp.308-309], [10, p.239])

Dirichlet 固有值問題:







∆ψk+ λkψk= 0 in Ω,

ψk = 0 on ∂Ω , (4)

滿足一個簡單的性質:

「區域單調性 (Domain-Monotonicity Property)」: 若 Ω1 ⊆ Ω2, 則 λk(Ω1) ≥ λk(Ω2), 對所有 k = 1, 2, . . . 都成立。

在 2007 年筆者與王學鋒的文章 [14] 裡, 對這類問題有一些初步的討論。 當然, 上述教科 書 ([15], [10]) 在 [14] 一文發表後, 都做了適度的更正。

圖 1

在 [14] 中, 作者提出了幾個相關的幾何量; 例如:「內在直徑 (intrinsic diameter)」, . . . 等。 我們先來看個例子。 設

Ωa =







{x ǫ Rⁿ| a < |x| < 1} if 0 < a < 1 {x ǫ Rⁿ| 1 < |x| < a} if 1 < a.

在 [14] 中, 證明了 µ2 (Ωa) 對 a > 0 是嚴格遞減的。

由上圖我們看到, 當 a 增加時, (-1, 0) 與(1, 0) 之間 「在 Ωa 之內的距離」 是遞增的, 直 觀上來說, 這意味著, 在 Ωa 的內部 應該需要更長的時間, 來 「平均」(-1, 0) 與 (1, 0) 附近的 溫度或熱量之差異。

同樣的, 幸或不幸, 事情沒那麼簡單! 還有許多其它的幾何量在這裡也扮演了重要的角色。

在 [14] 一文中, 還給出了許多其它例子, 有些甚至 「似乎」 違反了我們的直觀 · · ·

如何決定 µ₂ 的大小, 這個問題至今仍未能解決, 有興趣的讀者不妨參閱 [14] 一文。

綜上所述, 我們知道, 對於 Dirichlet 邊值問題,











ut= ∆u in Ω × (0, ∞) , u(x, 0) = u0(x) in Ω ,

u = 0 on ∂Ω × (0, ∞) ,

(5)

不論初值 u₀ 為何, 當 t → ∞, 解 u(x, t) 一定趨於 0, 而且滿足一個簡單的規律: 區域 Ω 愈 大, u(x, t) → 0 的速度就愈慢!

對於 Neumann 邊值問題 (1), 情形就複雜多了, 而且最令人困惑的是, 到現在我們還不 知道解 u(x, t) → u0 的速度依賴於那些 Ω 的幾何量, 這仍有待大家的繼續努力。

那麼, 全空間的解又如何呢? 為簡便計, 我們只考慮有界的初值u0, 以及有界的解u(x, t):







ut= ∆u in Rⁿ× (0, ∞),

u(x, 0) = u0(x) on Rⁿ. (6)

自然, 首先我們要問: 當 t → ∞ 時, 解 u(x, t) 是否一定會收斂 ?

有趣的是, 與 Dirichlet 或 Neumann 邊值問題都不一樣, 在全空間, (6) 的解並不一定 會收斂! 這個問題也有相當長的歷史。 我們在下一結果中, 給出了解收斂的充分必要條件。

Proposition 1. 當 t → ∞ 時, (6) 的解 u(x, t) → u∞(x) 若且唯若初始值 u0 滿足 :

R→∞lim 1 BR(x)

Z

BR(x)

u0(y)dy = u∞(x), 其中 BR(x) 是以 x 為中心 R 為半徑的球。 而且 u∞ 必定得是常數。

對這個結果有興趣的讀者, 可以參閱 [11] —其中有更廣的敘述。 最後, 我們用一個例子 ([2]) 來結束本節關於擴散的討論。

例: 令 u0ǫ C^∞(R) 為一偶函數滿足 |u0|L^∞ = 1 且

u0(x) = (−1)^k for x ǫ [k! + 2^k, (k + 1)! − 2^k+1].

則 (6) 的唯一有界解 u(x, t) 滿足 lim inf

t→∞ u(0, t) = −1 and lim sup

t→∞

u(0, t) = 1.

2. 當有反應項時, 情形就大不相同了

首先, 我們從最簡單的反應項談起: 不依賴於時間或空間的非線性反應項; 也就是說, 一個

「自律(autonomous)」 的方程











ut= d∆u + f (u) in Ω × (0, T ) , u(x, 0) = u0(x) in Ω ,

∂νu = 0 on ∂Ω × (0, T ) ,

(7)

此處 d > 0 代表擴散係數, 0 < T ≤ ∞ 是解存在的最大時間段, 對這最簡單的情形, (7) 的 解即使是全域存在而且有界, 當 t → ∞ 時, 解也不知道是否一定會收斂! (更多細節請參看 [N;Chapter 2])。

但是, 如果 (7) 的解是收斂的, 那麼它收斂的極限一定是 (7) 的平衡態 (equilibrium, steady state)。 在 1978 年, Casten 與 Holland [3] 證明了

Proposition 2. 當 Ω 是凸區域時, (7) 沒有局部穩定的非常數的平衡態。

換言之, 只要 Ω 是凸的, (7) 的任何穩定的平衡態, 必定是個常數, 且必須是反應項 f (s) 的零 點! 也就是說, 此時, 自律系統 (7) 是難以用來描述或刻劃 pattern formation 的。

非自律系統可以描述或刻畫的現象, 就非常豐富了; 比方說, 當反應項依賴於空間變數 xǫΩ 時, 情形就已經大不相同。 在本節剩下的篇幅裡, 我們用一個古老而簡單的方程來闡明這個情 形。

在大學二年級 「常微方程」 的課程裡, 我們都學過 logistic equation ut= ru(1 − u

K), (8)

此處 r > 0, K > 0 都是常數, 這個方程是生物數學家 Verhulst 在 1837 年提出來, 用來描 述一個種群的人口總量。 r 代表種群人口的增長率, K 是環境資源所允許的最大人口總量。 此 方程中的 u = u(t) 代表在時間 t ≥ 0 時的人口總量, 這個方程比 C. Darwin 的 「進化論」

還要早 22 年!

如果我們將 (8) 式寫成

ut= r

Ku(K − u), 再重新調整 t 的大小, 吸收掉常數 (r/K), 就得到了

ut = u(K − u). (9)

加入空間變量 x ǫ Ω, 並假設物種遷移是隨機的, 我們即有







ut= d∆u + u(K − u) in Ω × (0, T ) ,

∂νu = 0 on ∂Ω × (0, T ) , (10) 此處d > 0 是常數, 表示種群中個體遷移的強度。

數學上我們容易證明, 當 K 是常數時, V (x) ≡ K 是個全域穩定的平衡態; 也就是說, 不論 (10) 中u 的初始值u(x, 0)為何(只要不全為 0), 最終當 t → ∞ 時, u(t, x) 一定趨於 V (x) ≡ K! 同樣的, 當 K 6≡ 常數時, 用上、 下解方法, 我們可以得到

Proposition 3. K 6≡ 常數, 則 (10) 有唯一的正平衡態 ud; 同時, ud 是全域穩定的。

2006 年, 樓元 [12] 發現了一個簡單而出色的結果: 將 ud 的方程通除以 ud 再積分, 用一 次散度 (Divergence) 定理, 得到

d Z

Ω

| ▽ ud| u²_d

2

+ Z

Ω

(K − ud) = 0.

由 K 6≡ 常數, 馬上得到 ud6≡ 常數, 所以 Z

Ω

K <

Z

Ω

ud. (11)

這個簡單的不等式 (11) 卻有不尋常的生物意義, 它告訴我們: 總人口一定嚴格大於環境 資源承載總容量 (total carrying capacity)! 這似乎違反了我們的直觀! 請不要忘記, (11) 式 成立的兩個重要因素是: 種群要有遷移(dispersal), 不論大小; 再就是資源分佈必須是不均勻 的 (反映在 「K 6≡常數」 上)!

不難證明 (直觀上也非常清楚), 當d → 0 時, ud 一致收斂到K, 而當 d → ∞ 時, ud 一 致收斂到 K ≡ _|Ω|¹ R

ΩK(x)dx. 由此推得, 當 d → 0 或 ∞ 時, Z

Ω

ud→ Z

Ω

K. (12)

(11) 式及 (12) 式表明了, 若將 R

Ωud看成 d 的函數, 在某一個 d0 它會達到它的最大值。

自然, 我們要問:

(i) 如何決定 d₀? (ii) R

Ωud0 到底能有多大 (與 R

ΩK 相比)?

這是兩個簡單而基本的問題, 卻極具挑戰性, 到目前為止, 我們仍然不知道答案。

更進一步, 令

E(K) ≡ max

0 < d < ∞

R

Ωud

R

ΩK. 我們知道 E(K) > 1, 接下來考慮

ρ ≡ sup { E(K)|K ≥ 0, 6≡ 0} . 問題: ρ 是否有限? 如果是, 它的值是多少?

幾年前, 我有個猜測: 當 n = 1 時, ρ = 3。 這個猜測最近被白學利、 何小清和李芳 [1] 證 明了, 但是 n > 1 的情形則完全沒有頭緒, 但是我們相信 ρ 會遠大於 3 !

我們提到的這幾個簡單的數學問題有一定的生態學中的意義。 如問題 (i) 與 (ii), 它們的 答案可以告訴我們, 當資源在空間分佈不均勻時, 種群應如何選擇它的擴散率, 以使得它的人口 總量達到最大值 E(K), 以及 E(K) 的值; ρ 的意義在於: 指出對於給定的資源總量, 如何選擇 資源在空間的分佈, 使得 E(K) 達到最大, 以及比較這個最大值是環境資源承載總容量 R

ΩK 的若干倍。

關於 ud 與 K 的關係, 還有許多有趣的性質與問題, 在此我們就不一一敘述了。

3. 異質環境中的物種競爭

(11) 這個簡單的不等式在異質環境中, 對兩個物種的競爭有深遠的影響。 假設兩個物種 U, V 之間的競爭滿足經典的 Lotka-Volterra 方程組



















Ut= d1∆U + U(K(x) − U − c V ) in Ω × (0, ∞) , Vt= d2∆V + V (K(x) − b U − V ) in Ω × (0, ∞) , U(x, 0) = U₀(x), V (x, 0) = V₀(x) in Ω ,

∂νU = ∂νV = 0 on ∂Ω × (0, ∞) .

(13)

在 (13) 中我們假設了這兩種群 U, V 需要完全相同的資源 K(x), 而他們彼此的競爭能 力以 b, c 來表示。 當 b, c 都介於 0 與 1 之間時, 我們說 U, V 之間的競爭是 「弱競爭」。

在弱競爭時, 如果資源在空間中是均勻分佈的 (即 K(x) ≡ 常數), 一個眾所周知的事實 是: 不論初始值 U0, V0 的大小, U, V 必定共存! 可以證明 ([9]): 當 t → ∞ 時,

(U, V )(x, t) → ( 1 − c

1 − bc, 1 − b 1 − bc)K.

而且這個結論對任意 d₁, d2 都成立。

但是當資源分佈不再均勻時, 樓元 [12] 在 2006 年證明了一個有趣而驚人的結果: 存在某 些 d1 < d₂, 不論初始值 U₀, V₀ 的大小, U 最終總是將 V 完全消滅掉! 也就是說, 即使在 「弱 競爭」 的情況下, 只要資源分佈是不均勻的, 對於適當的遷移擴散速率 d₁, d₂, 共存不再是可能 的了。

其實, 早在 1998 年, [5] 就曾經考慮過一個 (13) 的特殊情形



















Ut= d1∆U + U(K(x) − U − V ) in Ω × (0, ∞) , Vt= d2∆V + V (K(x) − U − V ) in Ω × (0, ∞) , U(x, 0) = U₀(x), V (x, 0) = V₀(x) in Ω ,

∂νU = ∂νV = 0 on ∂Ω × (0, ∞) ;

(14)

也就是說, 如果除了擴散速度不同, U, V 這兩個種群在任何其他方面都是完全相同的, 那麼 [5]

也得到了一個驚人的結論: 只要 d₁ < d₂, 不論初始值 U₀, V₀ 的大小, 最終 U 總是將 V 完全 消滅掉! 這就是著名的 「Slower diffuser always prevails!」

在最近的一篇文章 [6] 裡, 何小清與筆者利用了 §2 中的 E(K), 在 bc ≤ 1 時, 將 (13) 的全域動力學性態完全刻劃清楚了。 這個結果不但同時包含了 「弱競爭」 及 「Slower diffuser always prevails」, 而且將這兩個重要而又似乎無關的定理統一起來, 使得它們之間的關聯也清 楚了。

在何小清與筆者的一系列文章中 [6]−[8], 我們依照解的全域動力學性態, 對 d₁, d2 做了 分類刻劃; 同時更進一步, 也考慮了更為一般的情形。

4.

現在讓我們回到原先的 logistic equation (8) 及其簡化後的 (9)。 從常微分方程的角度 來看, (8) 與 (9) 基本上是沒有差異的。 但是在空間異質 (heterogeneous) 時, (8) 所對應的偏 微分方程就應該是











ut= d ∆ u + r(x)u(1 − u

K(x)) in Ω × (0, T ) , u(x, 0) = u0(x) in Ω ,

∂νu = 0 on ∂Ω × (0, T ) .

(15)

那麼, (15) 與 (10) 有差異嗎? 差異何在? 特別要關心的是, (15) 的平衡態是否也滿足 (11) 呢?

首先, 對許多生態學家來說, 研究 (15) 是更為必要的, 因為他們相信內在成長速率(in- trinsic growth rate) r(x) 與 環境資源承載容量 (carrying capacity) K(x) 是兩個不同的 概念, 不同的量, 對一個種群是兩種不同的特質, 應該被區分開來, 但這並不表示二者之間沒有 關係。 顯而易見的, 當 r(x) ≡ cK(x) 時 (c 是個正常數), (15) 與 (10) 幾乎是一樣的, 接下 來, 我們來討論一般的情形。

標準的上下解方法可以用來證明: 對每一個 d > 0, (15) 具有唯一的一個全域穩定的平衡 態 ud, 在 [4] 一文中, 當 d 足夠小的時候, 推導出

Z

Ω

ud = Z

Ω

K − dhZ

Ω

▽K · ▽(1

r) + o(1)i

. (16)

所以, 當 r 與 K 是 「正相關 (positively correlated)」(也就是說, 當 K 增加時, r 不會減少) 的時候, 我們有

▽K · ▽(1

r) = −1

r² ▽ K · ▽r ≤ 0, 代回 (16) 式, 即得, 對足夠小的 d,

Z

Ω

ud>

Z

Ω

K;

也就是說, 如果 r 與 K 是 「正相關」 的話, 只要 d 足夠小, (11) 式依然成立! 那麼對於足夠大 的 d 又如何呢? 在 [4] 一文中, 也推出了如下的結論:

ud→ L ≡ R

Ωr(x) R

Ω r(x) K(x)

as d → ∞. (17)

下一個問題是, L 與 K 孰大? 如果我們更進一步, 將 r 設成 K 的函數, [4] 就給出了下 面的解答:

(i) 如果 r(K) 是一個 凸函數, 則 L > K;

(ii) 如果 r(K) 是一個 凹函數, 則 L < K.

一般來說, r 依賴於許多變數, 而不是只依賴於 K, 但是如果將其它變數固定, 只讓 K 變 動, 從直觀上來看, 許多時候 r 應該是 K 的凹函數。 所以當 d 很大時, (15) 的解通常不具備 (11) 式的性質。

在 §3 中談到 Lotka-Voltera 競爭—擴散 (competition-diffusion) 的種群有趣現象大多 依賴於 (11) 式, 因此, 如果我們考慮對應於 (15) 式的 Lotka-Voltera 競爭—擴散 (competition- diffusion), §3 中談到的各種現象是否依然會發生? 在什麼條件下可以發生? 這些問題都是當 下我與何小清正在探索的方向。

5.

我們已經談過 (11) 式對種群競爭的影響, 但是也不要忽略了 (11) 式本身在生態學中的 意義: 在資源分佈不均的情況下, 對有隨機擴散的種群, 它能支援的人口總量超過它的環境承載 容量! 直觀上, 這似乎令人難以置信!

近年來, 有部分資深的生態學家注意到這個 「純數學」 預測的理論結果, 開始設計一系列的 實驗來證實它。 在 [16] 中, 作者們用水萍 (duckweed) 檢驗了對應於 (10) 的 「patch model」

的解 — 他們的實驗證明了 (11) 在 「patch model 」 中確實會發生!

DeAngelis 與他的團隊也已開始用實驗來探討 r 與 K 的關係, 我們十分期待瞭解 r 與 K 之間關係的實驗資料。

6.

^{反應擴散方}程一直是極為豐富而且活躍的領域, 許多對於理解自然十分根本的問題, 仍懸 而未解。 在這短短的篇幅中, 竟能包羅這麼多有趣甚或 「出人意表」 的結果, 似乎顯示出, 我們 對於反應和擴散間交互作用的機制, 真正所知甚少。 這是個廣闊而吸引人的領域, 有待我們進一 步深入探索。

後記

本文是基於作者於 2015 年 6 月初在臺灣清華大學對高年級 (大三、 大四) 以上的學生 (包括研究生) 作的一次通俗演講的講稿整理而成。

感謝陳兆年教授的邀請, 宋宥松同學的筆記, 以及理論中心的資助。 特別感謝中研院的李 宣北博士, 也是作者大學時代的老同學的邀稿、 建議與幫助, 以及助理編輯黃馨霈小姐的耐心, 本文得以在 「數學傳播」 上刊登。

最後, 還要感謝孫毓孜小姐的周到, 使新竹之行圓滿愉快!

參考資料

1. X. Bai, X. He and F. Li, An optimization problem and its application in population dynamics, Proc. AMS, to appear.

2. P. Collet and J.-P. Eckmann, Space-time behavior in problems of hydrodynamic type:

A case study, Nonlinearity 5(1992), 1265-1302.

3. R. G. Casten and C. J. Holland, Instability results for reaction-diffusion equations with Neumann boundary conditions, J. Diff. Equns. 27(1978), 266-273.

4. D. L. DeAngelis, Wei-Ming Ni and B. Zhang, Dispensal and heterogeneity: Single species, J. Math. Biol. to appear.

5. J. Dockery, V. Hutson, K. Mischaikow and M. Pernarowski, The evolution of slow dispersal rates: A reaction-diffusion model, J. Math. Biol. 37(1998), 61-83.

6. X. He. and Wei-Ming Ni, Global dynamics of the Lotka-Voltea competition-diffusion system: Diffusion and spatial heterogeneity, I, Comm. Pur Appli. Math., to appear.

7. X. He. and Wei-Ming Ni, Global dynamics of the Lotka-Voltea competition-diffusion system with equal amount of total resources, II, preprint.

8. X. He. and Wei-Ming Ni, Global dynamics of the Lotka-Voltea competition-diffusion system with equal amount of total resources, III, preprint.

9. S.-B. Hsu, Limiting behavior for competing species, SIAM J. Appl. Math. 34(1978), 760-763.

10. J. Jost, Partial Differentail Equations, Springer, 2002.

11. S. Kamin, On the stabilisation of solutions of the Cauchy problem for parabolic equa- tions, Proc. Roy. Soc. Edinburgh (Sect. A) 76(1976/77), 43-53.

12. Y. Lou, On the effects of migration and spatial heterogeneity on single and multiple species, J. Diff. Eqns. 223(2006),400-426.

13. Wei-Ming Ni, The Mathematics of Diffusion, CBMS-NSF Regimal Conf. Ser. Appl.

Math. 82, SIAM, Philadelphia, 2011.

14. Wei-Ming Ni and X. Wang, On the first Neumann eigenvalue, Discrete Contin. Dyn.

Syst. 17(2007), 1-19.

15. W. A. Strauss, Partial Differential Equations, John Wiley and Sons, Inc. New York, 1992.

16. B. Zhang, X. Li, D. DeAngelis, Wei-Ming Ni and G. Wang, Effects of dispersal on total biomass in a patchy, heterogeneous system: Analysis and Experiment, Math. Biosci., 264(2015), 54-62.

—本文作者任教中國華東師範大學及美國明尼蘇達大學—