生態系のバランスや、バタフライエフェクトの問題をODE45で解く：無償ソフトで技術計算しよう【シミュレーション応用編】（1）（1/2 ページ）

「微分方程式を解く」といえば、無味乾燥なイメージ。そこで今回は、生態系のバランスや、気象など私たちの生活に関わっているような数式や問題を取り上げる。

[伊藤孝宏，MONOist]

	FreeMatでエンジニアリング計算：
⇒	入門編
⇒	グラフィックス編
⇒	プログラミング基礎編
⇒	プログラミング処理編
⇒	プログラミング応用編
⇒	シミュレーション基礎編

　今回から応用編として3回に渡って、ode45を使って微分方程式を解いてみます。「微分方程式を解く」といった場合、無味乾燥なイメージを持たれるかと思いますが、例題として、私たちの生活に関わっているような問題を取り上げてみます。

筆者注：FreeMatはコマンド入力後にエンターキーを押すとコマンドを実行します。本連載ではエンターキーの記述を省略しますが、操作の際にはコマンド入力後にエンターキーが必要です。

ロジスティック方程式

　生物の増え方に関する研究は、害虫や細菌などの増殖を予測することに応用できるため、古くから研究されてきました。例えば、固体数をN、世代当たりの増加率をrとすると、世代当たりrNだけ増加します。

　すると、

　となり、初期値がN₀として解くと、

　と爆発的に増加することになります。しかし現実には、食料や場所の制約から、ある程度以上には増加できないと予想されます。制約を入れた状態での個体数の増加を表したのがロジスティック方程式です。ロジスティック方程式は下記で表されます。

　ここで、rは増加率で、Kは環境収容力で前述の制約に相当するものです。

　ロジスティック方程式をFreeMatで解いてみましょう。増加率rを0.1、環境収容力Kを10⁶、初期個体数N₀を100として、0〜200世代の個体数の変化を計算してみます。下記のように、関数名をTNとして、ロジスティック方程式を無名関数で定義します。次にode45で0〜200の範囲で数値解法をしてみます。得られた結果をプロットすると、図1が得られます。ある世代（今回の場合は70世代当たり）を超えると爆発的に増加し、環境収容力に近づくにつれて、増加は緩やかになることが分かります。

--> TN=@(t,N) 0.1*(1e6-N)/1e6*N;
--> [t,N]=ode45(TN,[0,200],100);
--> plot(t,N);grid('on');

図1：個体数の増加曲線

ロトカ・ヴォルテラ方程式

　生態系の頂点に立つものを除けば、生物は弱肉強食の関係にあります。一般に肉食動物から食べられてしまう草食動物は、繁殖力が大きく、逆に、肉食動物は子どもの数が少ないというふうに、バランスを取っています。この関係を表したのが、下記のロトカ・ヴォルテラ方程式です。

　ここで、xは被食者の個体数、yは捕食者の個体数です。

　下記のように式を書き直します。草食動物xは植物を食べて自身が死ぬよりも多く子孫を残すため、a*xの割合で増加します。同時に、肉食動物yから食べられてしまい、これは草食動物xと肉食動物yの数に比例し、b*x*yの割合で減少します。一方、肉食動物yは、草食動物xを食べて子孫を増やします。これは、えさとなる草食動物xとも比例するため、d*x*yの割合で増加します。同時に、自身は寿命で死滅するため、c*yの割合で減少します。

　係数をa=0.5、b=0.02、c=0.5、d=0.01、初期値をx0=50、y0=1として、FreeMatで解いてみると、図2の結果が得られます。

図2：個体数の変化

　青線が草食動物xの個体数の変化で、緑線が肉食動物yの変化を示します。草食動物xが増加すると、えさが増えるため、肉食動物yも増えます。一方、草食動物xは肉食動物yから捕食されるため、肉食動物yの増加に伴い、草食動物xは減少します。えさとなる草食動物xが減少するため、肉食動物yも減少するといったサイクルを繰り返します。係数a,b,c,dを変えてみると、変化の様子が変わります。

　ex404.mがFreeMatのプログラムで、function dydt=rv(t,y)で係数a,b,c,dの設定と方程式の定義を行い、ode45に渡しています。計算結果は、配列yの1列目にx、2列目にyが格納されます。

function ex404
    clear;
    tspan=[0,100];
    t0=[50;1];
    [t,y]=ode45(@rv,tspan,t0);
    plot(t,y);
    grid('on');
    xlabel('Generation');ylabel('Number');
function dydt=rv(t,y)
    a=0.5;b=0.02;c=0.5;d=0.01;
    dydt=[a*y(1)-b*y(1)*y(2);-c*y(2)+d*y(1)*y(2)];

ex404.m

＞＞「ex404.m」ダウンロード