Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

捕食者と被食者の方程式の求解

この例では、ode23ode45 の両方を使用して捕食者/被食者モデルを表す微分方程式を解く方法を示します。これらはルンゲ・クッタ積分手法の可変ステップ サイズを使用する常微分方程式の数値解に対する関数です。ode23 は、中規模な精度に対して単純に 2 次、3 次の公式の組み合わせを使用し、ode45 は、高い精度に対して 4 次、5 次の組み合わせを使用します。

Lotka-Volterra の方程式または捕食者と被食者のモデルとして知られる 1 次の常微分方程式の組み合わせについて考えます。

dxdt=x-αxydydt=-y+βxy.

変数 x および y は、被捕食者と捕食者のそれぞれの個体群のサイズを表します。2 次のクロス項は、種の間の相互関係を説明します。捕食者が存在しない場合に被捕食者の個体群は増加し、被捕食者が少なくなると捕食者の個体群は減少します。

方程式のコード作成

方程式系のシミュレーションを実行するために、状態と時間の値を与えて、状態導関数の列ベクトルを返す関数を作成します。2 つの変数 x および y は、MATLAB® でベクトル y の最初の 2 つの値として表すことができます。同様に、導関数はベクトル yp の最初の 2 つの値です。この関数は t および y の値を受け入れて、yp の方程式によって生成された値を返さなければなりません。

yp(1) = (1 - alpha*y(2))*y(1)

yp(2) = (-1 + beta*y(1))*y(2)

この例では、方程式は lotka.m というファイルに含まれています。このファイルでは、α=0.01 および β=0.02 のパラメーター値を使用します。

type lotka
function yp = lotka(t,y)
%LOTKA  Lotka-Volterra predator-prey model.

%   Copyright 1984-2014 The MathWorks, Inc.

yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y;

方程式系のシミュレーション

ode23 を使用して、lotka で定義された微分方程式を区間 0<t<15 について解きます。x(0)=y(0)=20 の初期条件を使用して、捕食者と被捕食者の個体群が等しくなるようにします。

t0 = 0;
tfinal = 15;
y0 = [20; 20];   
[t,y] = ode23(@lotka,[t0 tfinal],y0);

結果のプロット

時間に対して結果の個体群をプロットします。

plot(t,y)
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
legend('Prey','Predators','Location','North')

Figure contains an axes object. The axes object with title Predator/Prey Populations Over Time, xlabel t, ylabel Population contains 2 objects of type line. These objects represent Prey, Predators.

次に、相互に個体群をプロットします。結果の位相面プロットにより、個体群間の巡回的な関係が非常に明らかになります。

plot(y(:,1),y(:,2))
title('Phase Plane Plot')
xlabel('Prey Population')
ylabel('Predator Population')

Figure contains an axes object. The axes object with title Phase Plane Plot, xlabel Prey Population, ylabel Predator Population contains an object of type line.

さまざまなソルバーの結果の比較

2 度目は ode23 ではなく ode45 を使用して方程式系を解きます。ode45 ソルバーでは、各手順において時間がかかりますが、手順もより大きくなります。それにもかかわらず、既定の設定で、ソルバーは連続的な拡張式を利用して、使用した各ステップの期間において 4 つの等間隔な時間点で出力を生成するため、ode45 の出力は滑らかです ('Refine' オプションを使用して点の数を調整できます)。比較のために両方の解をプロットします。

[T,Y] = ode45(@lotka,[t0 tfinal],y0);

plot(y(:,1),y(:,2),'-',Y(:,1),Y(:,2),'-');
title('Phase Plane Plot')
legend('ode23','ode45')

Figure contains an axes object. The axes object with title Phase Plane Plot contains 2 objects of type line. These objects represent ode23, ode45.

結果は、異なる数値的手法を使用して微分方程式を解くと、若干異なる解が生成されることを示します。

参考

|

関連するトピック