常微分方程式に空間変数を導入して偏微分方程式へ拡張¶
Volterra 問題を例として¶
以前扱った Volterra 問題は、自由変数が時間だけの常微分方程式で、これは、生物の数がどこでも一定ということを前提としていることを意味していた． しかし、実際の生物は空間分布が一定ではなく、そしておそらくこの分布は生物の数量変化に大きな影響を与える結構重要なファクターなんじゃないかな．
そこで、空間変数を自由変数として導入し、問題を偏微分方程式にして取り扱ってみよう． しかも、こういう問題は空間次元が 1のときと 2以上のときとで根本的に挙動が異なるので、思い切って空間次元を 2以上にとって、取り扱ってみよう．
最初は, おとなしめなモデルで¶
空間変数を新たに導入するということは、問題のモデルに新しいファクターを導入する必要が生ずる． もちろんいろいろなファクターがあり得るので、本来は慎重に考えないといけないところだ．
しかし今回は数値計算への熟達が主眼であるので、モデルとしては少し安直に、「常微分方程式からのズレがあまり無い」モデルを導入し、数値計算結果の解釈に苦労しなくて良いようにしよう．
そこで、単なる「拡散」だけを導入することにして、被食者数を $u(x, t)$ で、捕食者数を $v(x, t)$ で表すとして、 $x \in [0, L] \times [0, L]$, $t > 0$ とし、以下のようなモデルを考えよう．
$\left\{\begin{array}{rcl} \displaystyle \frac{\partial u}{\partial t} & = & \epsilon_1 \Delta u - C_1 uv + D_1 u, \\ \displaystyle \frac{\partial v}{\partial t} & = & \epsilon_2 \Delta v + C_2 uv - D_2 v \end{array}\right.$
ただし、係数 $\epsilon_1$, $\epsilon_2$, $C_1$, $C_2$, $D_1$, $D_2$ は全て定数で正としておこう．
そして、空間領域に関する境界条件としては周期的境界条件を課しておくことにしよう．
method of line を適用する¶
ここでは、空間領域 $\Omega = [0, L] \times [0, L]$ に対し、その分割数を $N \times N$ とし、境界条件として周期的境界条件を離散化しよう．
よって、$\Delta x = \Delta y = L/N$ として $u(x,t)$ の離散化に相当する従属変数を $ u(t) = \{ u_{i,j}(t)\}_{i, j = 1}^N$ とし、境界条件は 外側にはみ出す点を仮想的に考えたとして
$$ \left\{ \begin{array}{rcl} \displaystyle u_{0, j} & := & u_{N, j} \mbox{ for } \forall j ,\\ \displaystyle u_{N+1, j} & := & u_{1, j} \mbox{ 同上 } , \\ \displaystyle u_{i, 0} & := & u_{i, N} \mbox{ for } \forall i, \\ \displaystyle u_{i, N+1} & := & u_{i, 1} \mbox{ 同上 } \end{array} \right. $$と解釈することになる．
こうしておいて、空間微分(ラプラシアン) $\Delta$ をたとえば差分近似で置き換える． たとえば $U$ に対する二階中心差分ならば
$$ (\Delta u)_{i,j} := \frac{ u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j+1} - 4 u_{i,j} }{\Delta x^2} $$となる． これらを使った場合は、全体は $i = 1,2,\cdots, N$, $j = 1,2,\cdots, N$ に対して
$$ \left\{ \begin{array}{rcl} \displaystyle \frac{d u_{i,j}(t) }{dt} & = & \epsilon_1 (\Delta u)_{i,j} - C_1 u_{i,j} v_{i,j} + D_1 u_{i,j}, \\ \displaystyle \frac{d v_{i,j}(t) }{dt} & = & \epsilon_2 (\Delta v)_{i,j} + C_2 u_{i,j} v_{i,j} - D_2 v_{i,j}, \\ \end{array} \right. $$という $N \times N$ 元連立常微分方程式で近似できることになる．
あとは少しずつプログラムを作っていくだけだ． まずは、問題のパラメータを設定してしまおう．
ϵ1, ϵ2 = 0.001, 0.001
C1 = C2 = D1 = D2 = 1.0
L = 1.0
N = 40
Δx = L/N
Δt = 0.1 * Δx^2
ラプラシアンに使う係数行列を下記のように定義して…
d = [ -2.0 for i in 1:N ]
v = [ 1 for i in 1:N-1 ]
d2 = diagm(d, 0) + diagm(v, 1) + diagm(v, -1)
d2[1,N] = d2[N, 1] = 1.0
d2 = (1/(Δx)^2) * d2
下にあるように、近似データそのものである行列 $u$ に対して、上で定義した行列 $d2$ を用いて $d2 * u + u * d2$ を計算すると、これがラプラシアンを適用した結果相当になっている(線形代数が苦手な人は、for 文を使ってプログラムしたほうが良いかも)．
function f(U)
u, v = U[1:N,:], U[N+1:end, :]
uv = u .* v
fu = ϵ1 * ( d2 * u + u * d2 ) - C1 * uv + D1 * u
fv = ϵ2 * ( d2 * v + v * d2 ) + C2 * uv - D2 * v
res = vcat( fu, fv )
return res
end
上を見て分かるように、 $u(x)$ の近似行列 $u_{i,j}$ と $v(x)$ の近似行列 $v_{i,j}$ を受け渡しする時は、上下につなげて、一つにして扱っている．
そこで、これで連立常微分方程式の右辺(の近似)が計算できたことになるので、Runge-Kutta 法を定義する．
function RK(u)
f1 = f(u)
f2 = f(u + Δt/2 * f1)
f3 = f(u + Δt/2 * f2)
f4 = f(u + Δt * f3)
return u + Δt * (f1 + 2*f2 + 2*f3 + f4)/6
end
これで計算手順の準備は整った．
次に、初期値を適当に作ろう． 今回は、ガウス型に穏やかに分布しているものを考えよう．
u0 = [ 0.3 + 0.4 * exp( - 4.6 * 4 * ((i-N/2)^2 + (j-N/2)^2) / (N^2) ) for i in 1:N, j in 1:N ]
v0 = [ 0.4 + 2.0 * exp( - 4.6 * 16 * ((i-N/4)^2 + (j-N/4)^2) / (N^2) ) for i in 1:N, j in 1:N ]
U0 = vcat( u0, v0 )
どんな初期値なのか、目で見ておこう．
using Plots
gr()
X = Y = Δx:Δx:L
plot(X, Y, u0 , st = :surface)
plot(X, Y, v0, st = :surface)
ふむ、こんな感じか．
さて、これで準備は整ったので、早速、計算してみよう
using ProgressMeter
U = copy(U0)
sq_U = copy(U0)
@showprogress for i in 1:200000
U = RK(U)
if i % 2000 == 0 # 一定数ステップおきにデータを格納するとして
sq_U = cat(3, sq_U, U) # 3次元目に追加していく
end
end
試しに結果の数字をちょろっと見てみる．
sq_U
あまり変なことにはなっていなさそうなので、プロットしてみよう．
プロットには、fill=true をつけて、色を塗って平面で表してもらおう．
minimum(sq_U), maximum(sq_U)
range = (0.19, 2.87)
contour(X, Y, sq_U[1:N,:,15], clims = range, fill=true, fmt = :png) # なんとなくこいつをプロットしてみた．
ふむ、なんとなくうまく計算できているようだな． しかし、いちいち plot コマンドを打って確認するのも疲れるので、一気に全データを画像にしてしまうことを考えよう．
まず、画像を入れておくディレクトリをコマンドで作らせよう．
run(`mkdir U`)
run(`mkdir V`)
次に、データを画像化するコマンドを作ろう． 基本的には、plot の直後に savefig コマンドを使うだけだ．
function figureU(num)
contour(X, Y, sq_U[1:N,:,num], clims = range, fill = true)
savefig( "U/" * @sprintf("%04d", num) * ".png" )
contour(X, Y, sq_U[N+1:end,:,num], clims = range, fill = true)
savefig( "V/" * @sprintf("%04d", num) * ".png" )
end
では、全データを画像にしよう．やや時間がかかりそうだな．
@showprogress for n in 1:size(sq_U)[3]
figureU(n)
end
これで画像ファイルがたくさんできたはずだ． julia/jyupiter のファイル操作画面から、画像を確認しよう．
ffmpeg などの、動画作成ツールが有るならば、動画にしてみよう¶
フリーのツールで ffmpeg というものがあり、これを使うと連番の画像ファイルを無理なくきれいに動画にできる． 古い juliabox にはインストールされていないが、新しい juliabox にはインストールされているようだ.
例えばこれを julia から直接呼び出して使うならば、今回は以下のようにすればよいだろう．
(詳しい人向け: terminal からコマンドを打っても良い)
run(`ffmpeg -r 5 -i "U/%04d.png" -pix_fmt yuv420p -s 1200x800 U/out.mp4`) # U ディレクトリのデータ
run(`ffmpeg -r 5 -i "V/%04d.png" -pix_fmt yuv420p -s 1200x800 V/out.mp4`) # V ディレクトリのデータ
これでディレクトリ U とディレクトリ V の下に、それぞれ out.mp4 という動画ファイルができたはずだ． julia/jyupiter のファイル操作画面で、その動画ファイルをクリックして確認しよう．
レポート問題¶
今回は空間方向だけ差分化して、あとは常微分方程式として扱うという「線の方法」をとった． 時間方向も差分化したり、有限要素法を使ったり(かな-り面倒だが)他の方法で計算してみよう．
レポート問題(チャレンジング)¶
今回は空間変数の導入にあたり「拡散効果」だけを導入したが、おそらく、「敵から逃げる」「餌を追う」という効果を導入したほうが生物モデルとしては妥当だろう． それにはどうしたらよいか．いろいろ考えて、試してみよう．