COMSOL Multiphysics® ソフトウェアのバージョン 6.3 のリリースにより, 時間依存スタディの過去の解を現在の時間ステップで使用できるようになりました. つまり, 以前の時刻の解を現在の時刻の解に反映させることができるのです. これは, 遅延微分方程式のモデリングなど, 多くの用途に利用できます. この新機能を, 少し変わった例で見てみましょう.
山の草原の思索
少し前に, 山でハイキングをする機会があり, マーモット を見つけました. 数枚の写真を撮った後, この愛らしい草食動物のライフサイクルについて考え始めました. するとすぐに, 山の草原で草を食むマーモットの個体数の変化を記述する方程式を思いつくというアイデアが浮かびました.
このブログのインスピレーションの源となった, ふわふわのマーモット.
大学の数値解析の教科書を少し読み返してみると, Lotka–Volterra 方程式やロジスティック微分方程式など, 役に立つかもしれない様々な方程式が浮かび上がってきました. これらの方程式を参考に, 2つの常微分方程式 (ODE) を書こうとしましょう. 1つは牧草地のマーモットの数の変化率 M(t), もう1つは可食バイオマス量の変化率 B(t) です.
マーモットの個体数の変化率を表す最初の式には, 繁殖に関する項 R(t) と死亡に関する項 D(t) が含まれており, これらについては後ほど詳しく説明します. 2番目の式には3つの定数があります: G_0 は成長率, B_{max} は牧草地で生育できる植物の最大質量 E_0 はマーモットの1日あたりの消費量です.
さて, マーモットの個体数を表す式に移りましょう. まず, 現在の個体数に線形依存する繁殖に関する項から始めます:
一方, マーモットの死亡率は利用可能な資源と関連しているはずなので, 中間ステップとして, 1日あたりの餌面積を定義します:
ここで, A_0 は牧草地の面積です. 餌場面積が大きいほど, 捕食動物 に遭遇する可能性が高くなると仮定します. また, 面積あたりの捕食動物の数 P(t) = P_0 は一定であると仮定します. したがって, マーモットの死亡率は, 現在のマーモットの数 M(t), 各マーモットの1日の餌場面積, そして捕食動物の密度に比例します. このことから, 死亡率の式は, 以下の2つの未知数に依存します:
さらに, マーモットのコロニー全体が一斉に冬眠から目覚めると仮定し, これらの方程式を解いて夏季の個体数変化を予測します. その後, コロニーは冬季に備えて巣穴に戻ります. マーモットの数または植生の量が臨界値を下回り, コロニーまたは環境の崩壊を示唆した場合にソルバーを停止する停止条件も設定します.
スクリーンショットは, 一連の常微分方程式をモデル化する方法を示しています. 常微分方程式ごとに単位が異なる場合があります. 時間依存ソルバー 機能にも 停止条件 機能が追加されています.
これらの方程式は, 上記のスクリーンショットに示すように, COMSOL Multiphysics® の グローバル常微分方程式および微分代数方程式 インターフェースを使用して解くことができます. このインターフェースは, 解くべき2つの変数 M(t) と B(t), そしてそれぞれの時間微分と初期値 M_\text{init}=M(t=0) と B_\text{init}=B(t=0) に関する方程式を定義します. これらの方程式を解くことで, 夏季の個体数変化を予測できます.
植生バイオマス (緑) とマーモット個体群 (茶色) の動態を表す方程式の解. 夏の終わりの個体群の注釈付き. マーモットの個体数は当初急激に増加しますが, その後, 捕食の増加により減少します.
結果は, マーモットの個体数が当初急増し, 利用可能な食物を急速に消費することを示しています. これは, 捕食の増加による個体数の減少という修正をもたらします. さらにしばらくすると, マーモットの個体数は均衡状態に達します. (興味のある読者への補足として, これが 安定した 均衡状態であるかどうかを判断するのは非常に困難です.)
遅延微分方程式による追加要因の考慮
この時点では, 妥当と思われる結果が得られましたが, いくつかの懸念事項もあります. 考慮すべき重要な要因が少なくとも2つあります:
- マーモットは冬眠から目覚めると交尾を始めますが, すぐに繁殖するわけではありません. 妊娠期間があるからです.
- マーモットの個体数が十分に長い期間にわたって増加すると, より多くの捕食者が草原に引き寄せられます.
では, これら2つの要因を数学的にどのように記述し, モデルにどのように実装するかを見てみましょう. まず, 再生産率を受胎時の個体群の関数とし, 再生産率がゼロとなる妊娠期間 T_g を考慮します. したがって, 再生産率の式は次のようになります:
ここで H(t-T_g) はヘヴィサイドの階段関数です. この割合は, 妊娠期間中の個体群減少率の割合 \frac{M\left( T_g \right)}{M_\text{init}} によっても減少します. 若いマーモットが餌を探し始めると, 年長のマーモットによる捕食は停止すると仮定します. したがって, この割合は期間 t=T_g 後も一定のままです.
これで, if() 文と at(time,variable) 演算子を組み合わせてソフトウェアに入力できる遅延微分方程式ができました. 繁殖項は次のように記述できます:
R0*if(t>Tg,at(t-Tg,M)*at(Tg,M)/Minit,0)
at() 演算子の最初の引数は, 2番目の引数を評価する前の時刻です. 求解中に過去の時刻を参照するこの機能を有効にするには, 以下のスクリーンショットに示すようにソルバー設定を変更する必要があります. 出力時間ステップが遅延時間を超えない 厳密な時間ステップ を使用することで, 出産期の開始を確実に捉えることができます. それ以外の場合は, 以前と同じ方法でモデルを解きます.
求解中に時間演算子の使用を有効にする方法を示すスクリーンショット.
この一定でない繁殖率を組み込んだ場合の個体群を見てみましょう. まず, マーモットが生まれていない期間に個体群が減少し, その後, 個体群が増加して横ばいになり始めることがわかります.
出生率に遅延期間を導入した場合のマーモットの個体数 (茶色) と植生の量 (緑).
最後に, マーモットの平均個体数が十分に長い期間にわたって増加した場合に発生する捕食者の増加を考慮して, 死亡率を修正しましょう. つまり, 過去数日間におけるマーモットの個体数の時間平均を追跡する必要があります.
この積分を期間 T_p にわたるすべての過去の時間ステップにわたって評価するのではなく, 微積分の基本定理を適用して, 過去数日間の平均人口の変化率を表す式を書くことができます:
この方程式は, モデルに別の全体方程式を追加するだけで解くことができますが, t=0 から評価されるため, 別の if() 文を使用して M(t<0) = M_\text{init} を定義します. 移動平均を使用して捕食者密度 P(t) = \left( P_0 /M_\text{init}\right) M_\text{ave}(t) をスケールアップし, 死亡率に影響を与えます. もう一度解くと, 夏の終わりの個体群への影響がわかります.
朗報です! 個体群が増加し, マーモットは長い冬眠のために穴を掘ることができます.
結果はマーモットの個体群 (茶色) とバイオマス量 (緑) を示しており, 繁殖率と死亡率はどちらも以前の解に依存しています.
実験的観察により, マーモットの個体群の増加が確認されています.
個体群バランスモデリングに関する簡単な説明
ここで強調しておきたいのは, ここで紹介した例は, 分かりやすいモデルを用いて, いくつかの画期的な新機能を実証することを目的としているということです. 実際の個体群動態モデルは, これよりもはるかに複雑で, 多くの場合, コンパートメント化アプローチや区間ベースアプローチを用いており, 異なる変数を用いて全個体群の異なる割合を追跡します. このアプローチの入門例としては, 疫学で用いられる SEIR モデル (感受性, 曝露, 感染, 回復と免疫) が挙げられます. これについては, ブログ “COMSOL Multiphysics® による COVID-19 の拡散モデル化” で解説しています.
この区間ベースアプローチは, 特に化学工学の分野で注目されており, 液滴や沈殿物のサイズを追跡するためによく用いられています. 以下のチュートリアルモデルは, このタイプのモデリングに焦点を当てています:
注意事項とまとめ
COMSOL Multiphysics® バージョン 6.3 における遅延微分方程式の実装と解法の基礎を説明しました. このシンプルなモデルは非常に興味深いダイナミクスを生み出す可能性がありますが, この例は学習目的であることをご留意ください. このシンプルな例を理解すれば, より高度なモデリングタスクに取り組む準備が整います. 例えば, 時間に関する遅延微分方程式と空間に関する微分方程式を組み合わせることも可能です. ソフトウェアのすべてのインターフェースは, これらの新しい演算子をサポートしています. マルチフィジックスモデリングで使用するインターフェースの組み合わせに関わらず, 以前の解を参照する方法は, ここで示したものと全く同じです.
上記のモデルを試してみませんか? アプリケーションギャラリから関連する MPH ファイルをダウンロードしてください:
コメント (0)