pythonで計算(一次元振動モデル、ロトカヴォルテラ方程式)



スポンサードリンク

今回は、pythonの勉強の一環でやってきた簡単な数値計算のソースコードを紹介していきます。

一次元減衰強制振動

x-t(振幅ー時間)グラフは以下のようになります。

両解法(オイラー法とルンゲクッタ法)において、以前にも記事に書いた通り、一次元で方程式が単純ということもあり、そんなに解に差は見られません。

>>>以前の記事はこちら

python初心者の方は、設定パラメータを動かして、どんなグラフになるか、遊んでみてください。

以下に、一次元減衰強制振動のモデルを示します。

 

ばね定数:k、減衰係数:c、質量:m、強制外力:m/u


スポンサードリンク

ロトカヴォルテラ方程式

自然界における捕食者と被食者の個体数変動の数理モデルです。

以下の数式で定義されています。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= ax-bxy\\
\frac{dy}{dt}= cxy-dy
\end{array}
\right.
\]

xが被食者の個体数、yは捕食者の個体数、tは時間、4つの係数a、b、c、dは正の実数パラメータです。

以下は、a=8, b=3, c=4 ,d=18の時の計算結果。

周期関数になることがわかりますが、直感的に何が起きているか不明なので、定性的に述べていきます。ライオンとシマウマの例で考えると、

  1. 捕食者であるライオンが増加
  2. それにより被食者であるシマウマが減少
  3. エサ(シマウマ)が減少するので、ライオンは減少
  4. 天敵(ライオン)が減るので、シマウマは増加
  5. 1以降を繰り返す

といったことを繰り返しています。よくできたモデルですが、簡易化されているため、実際はもっと複雑です。以下は上で述べたことを図にしたものです。赤が捕食者(ライオン)で、青が被食者(シマウマ)です。

 

ルンゲクッタ法とオイラー法の比較

オイラー法とルンゲクッタ法の誤差の確認として、被食者数(上で言えばシマウマ)をグラフにしてみました。

ぱっと見では、解に差が内容に思われますが、縮尺を変更すると

時刻t=122では,0.5くらいの差があります。

割合でいうと、0.5/12=4%くらいの差があります。この程度の計算回数で誤差4パーは論外ですね笑

まあルンゲクッタが正解かというとそうでもないですが。

ついでに、グラフが直線的になっているのは、計算する際の微分間隔が雑なためです。

以下にコードを載せますが、そこの「微分間隔」と書いてあるところを小さくしていけば、理論的な微分に近づき、滑らかな解(グラフ)が得られます。

計算に用いたプログラム

 

以上です。基礎的な微分方程式や運動方程式の数値解法の勉強する際には役立てみてください。


スポンサードリンク

記事が役に立ったらシェア!

投稿者:

中村 俊

中村 俊

1993/09/04生まれ。機械系大学院を休学し、ベンチャーでインターンしている最中。直近では、デカルトの「方法序説」に感銘を受けた。 趣味:読書、web開発の勉強、異分野の論文読んだり、記事書いたり。 最終的には経営者か研究者になりたい。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

CAPTCHA