このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
捕食者と被食者の方程式の求解
この例では、ode23
と ode45
の両方を使用して捕食者/被食者モデルを表す微分方程式を解く方法を示します。これらはルンゲ・クッタ積分手法の可変ステップ サイズを使用する常微分方程式の数値解に対する関数です。ode23
は、中規模な精度に対して単純に 2 次、3 次の公式の組み合わせを使用し、ode45
は、高い精度に対して 4 次、5 次の組み合わせを使用します。
Lotka-Volterra の方程式または捕食者と被食者のモデルとして知られる 1 次の常微分方程式の組み合わせについて考えます。
変数 および は、被捕食者と捕食者のそれぞれの個体群のサイズを表します。2 次のクロス項は、種の間の相互関係を説明します。捕食者が存在しない場合に被捕食者の個体群は増加し、被捕食者が少なくなると捕食者の個体群は減少します。
方程式のコード作成
方程式系のシミュレーションを実行するために、状態と時間の値を与えて、状態導関数の列ベクトルを返す関数を作成します。2 つの変数 および は、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
というファイルに含まれています。このファイルでは、 および のパラメーター値を使用します。
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
で定義された微分方程式を区間 について解きます。 の初期条件を使用して、捕食者と被捕食者の個体群が等しくなるようにします。
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')
次に、相互に個体群をプロットします。結果の位相面プロットにより、個体群間の巡回的な関係が非常に明らかになります。
plot(y(:,1),y(:,2)) title('Phase Plane Plot') xlabel('Prey Population') ylabel('Predator Population')
さまざまなソルバーの結果の比較
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')
結果は、異なる数値的手法を使用して微分方程式を解くと、若干異なる解が生成されることを示します。