Julia言語を使って常微分方程式を解く（捕食者－非捕食者モデル）

こんにちは！Otepipiです。 今日はプログラミング言語Juliaを使用して微分方程式を解いてみます。 例題は下記「微分方程式で数学モデルを作ろう」から『捕食者ー非捕食者モデル』です。

微分方程式で数学モデルを作ろう

例題

$$捕食者－非捕食者の相互作用の単純なモデルは下記で与えられる。$$ $$\frac{dx}{dt}=x+y$$

$$\frac{dy}{dt}=-x+y$$

$$x………捕食者の個体数$$ $$y………非捕食者の個体数$$

$$初期値を下記のように仮定する。$$ $$x(0)=y(0)=1000$$

$$未来のxとyを決定せよ。$$

Juliaでのプログラミング

まず、常微分方程式を解くためのODEパッケージとグラフ描画のパッケージを使用することを宣言します。

using ODE
using Plots

次に今回求める微分方程式を記述します。 コーティングの方法については下記サイトを参考にしています。

github.com

function f(t, v)
     (x, y) = v
     dx_dt=x + y
     dy_dt=-x + y
    
    [dx_dt ; dy_dt]
    end;

初期値とタイムスパンを指定します。

v_in = [1000.0; 1000.0];

time = 0.0:0.1:1;

常微分方程式を解きます。

t, v = ode45(f, v_in, time);

プロットの準備をします。

x = map(v -> v[1], v);
y = map(v -> v[2], v);

結果をプロットします。

plot(t,y,seriestype=:scatter,label="y")
plot!(t,x,seriestype=:scatter,label="x")

解析解との比較

結果が正しいか確認するため解析解と比較します。 この微分方程式系の解析解は以下になります。 $$ x = 1000e^t(\cos t + \sin t) $$ $$ y = 1000e^t(\cos t - \sin t) $$

今回のシミュレーション結果と解析解を比較します。

x_th=map(x -> 1000*exp(x)*(cos(x)+sin(x)), t);
y_th=map(x -> 1000*exp(x)*(cos(x)-sin(x)), t);

plot!(t,y_th,label="y_th")
plot!(t,x_th,label="x_th")

結果は良く一致していました。