数学 #13 数値計算 #6 Lotka-Volterra方程式

Lotka-Volterra(ロトカ=ヴォルテラ)方程式

Lotka-Volterra方程式とは，捕食者と被食者の増減関係をモデル化した非線形微分方程式です．

dx/dt=x(α-βy)
dy/dt=-y(γ-δx)

x:被食者
y:捕食者
t:時間の個体数
α, β, γ, δ:正の実数パラメータ

被食者の増殖速度
dx/dt=αx-βxy
αxは，被食者が多ければ多いほど，その増加速度は大きいことを意味します，つまりαは増加に関するパラメータ．
-βxyは，被食者，捕食者が多い程，被食者の個体数が減少することを意味します(被食者が食べられる量，捕食者が食べる量はそれぞれの個体数に比例)，つまりβは減少に関するパラメータ．

捕食者の増殖速度
dy/dt=δxy-γy
δxyは，被食者，捕食者が多い程，捕食者の個体数が増加することを意味します(↑の-βxyと対になっている)，つまりδは増加に関するパラメータ．
-γyは，捕食者の自然現象の速度を表します，つまりγは減少に関するパラメータ．

今回はルンゲ=クッタ法で数値的に解いてみましょう．
被食者-赤
捕食者-青


Fig.1
初期値:
x₀=1
y₀=1
パラメータ:
α=0.9
β=0.5
γ=0.75
δ=0.25


Fig.2
初期値:
x₀=1
y₀=1
パラメータ:
α=1.5
β=1.5
γ=0.1
δ=0.5


Fig.3
初期値:3
x₀=γ/δ=3
y₀=α/β=1.8
パラメータ:
α=0.9
β=0.5
γ=0.75
δ=0.25

Fig.1を見ると，周期解になっていることが分かります．また，被食者，捕食者の位相が若干ずれていることが分かります．これは，

魚が増える→サメが食べることの出来る量が増えるので，サメは増え，魚は急激に減る→サメが食べることの出来る魚の量が減ったので，サメは自然減少する他無い→サメが減ったので，魚は増える→…

というサイクルを表しているのです．
※ここで，
魚=被食者
サメ=捕食者

Fig.2では，周期がFig.1のおよそ2倍になっています．また，個体数のピーク時は，被食者よりも捕食者の方が多いことが分かります．

Fig.3では，未来永劫 個体数が変わりません．これは，初期値がx₀=γ/δ，y₀=α/βでは，

dx/dt=0, dy/dt=0

となるためです．

現実におこる現象は人智を超えた複雑怪奇なものですが，こんな簡単な式によって決まっているのかもしれません，しかしそこには無秩序と秩序の狭間にあるカオスが渦巻いているのでしょう．

数学 #12 数値解析 #5 van der Pol方程式

van der Pol方程式

van der Pol方程式とは，次の微分方程式です．

y''+ε(y²-1)y'+y=0

Wikipediaより，
「ファン・デル・ポール振動子とは、非線形の減衰を受けた非保存系の振動子である。」
だそうです．

今回はこれを解いてみましょう．

数値計算で求める場合，色々な手法がありますが，ここでは，
v=y'
とおき，

y'=v
v'=ε(1-y²)v-y

とすることで，2変数にしてみましょう．
さらに数値計算手法として，3次のテイラー級数を用いてみます．
即ち，

y_n+1=y_n+hy'_n+(h²/2)y''_n+(h³/6)y'''_n
v_n+1=v_n+hv'_n+(h²/2)v''_n+(h³/6)v'''_n

y'=vやy''=v'は上にありますし，v''=y'''やv'''は気合で求められます(少しのdと公式 そして私はあなたと微分を覚えた)．

初期条件の設定
基本的に，数値計算で微分方程式を解くときは初期条件が必要です．
具体的には最初に
y₀と，v₀=y'₀
を決めれば良いのです．

実際の計算例

以下のグラフでは，横軸にy，縦軸にv=y'をとっています．
Fig.1にベクトル場((y, v)での傾き)，Fig.2に数値計算で求めた解を示します．


Fig.1


Fig.2

パラメータ:
ε=0.3
h(刻み幅)=0.01
初期値
赤(y₀, v₀)=(2.00092238555422, 0)
青(y₀, v₀)=(0.5, 0)
緑(y₀, v₀)=(2.5, -2.5)

さて，計算結果から以下の事が分かります．
赤の軌道…周期解になっている(以後ずっと回り続ける)．
青の軌道…閉曲線(赤軌道)の“内側”を初期値としているが，赤軌道に落ち込んでいる．
緑の軌道…閉曲線(赤軌道)の“外側”を初期値としているが，赤軌道に落ち込んでいる．

このようにvan der Pol方程式の解は，十分な時間が経つとこの閉曲線上(赤軌道)を回っている点の運動とすることが出来ます．このようなattractorを“limit cycle”というようです．

今回で使った3次のテーラー級数では，刻み幅hを大きくし過ぎると，発散してしまうようなので注意．

数学 #11 数値解析 #4 台形則

数値計算ネタ．

積分の近似計算です．今回も対象として，フレネル積分を用います．
しかし，計算方法が違います．今回は台形則という方法を使います．


参考までにフレネル積分(sin版)．


台形則

式はこんな感じです．



式だけだと意味不明なので，図をFig.1に示します．


Fig.1 trapezoid

赤線:y=sin(x²)
青い台形一つの面積:h/2(f(xⁱ)+f(xⁱ⁺¹)

台形のパラメタは以下の様に対応．Fig.1ではh=0.25です．
上底:f(xⁱ)
下底:f(xⁱ⁺¹
高さ:h=x軸の刻み幅
ですので，青い台形の面積の総和をとる式を構成すると，始めに示した式が出来上がります．

また，h(刻み幅)が小さければ小さいほど，計算の精度が上がる事も分かると思います(台形がどんどん細くなり誤差が減っていく)．例えば，h=1/16の時はFig.2のようになります．青の面積を全て足し合わせたものが近似値なので，精度が高くなる事が納得できると思います．


Fig.2

実際にFig.1のパラメタで計算したのがFig.3．


Fig.3
0≦x≦8
h=1/4
刻み幅が大きいので，かなりがたついています(ただ，前回の様に発散してはいません)．

hを小さくすると(Fig.2と同じ)Fig.4のような感じ．


Fig.4
0≦x≦8
h=1/16

非常に滑らかになり，精度が大幅に向上しています．

数学 #10 数値解析 #3 ルンゲ・クッタ法

先日紹介したオイラー法．今回は，それよりも精度の高いルンゲ・クッタ法を紹介したいと思います．上には上がいるものです．

ルンゲ・クッタ法(Runge-Kutta method)

初期値を，



とすると，(hは刻み幅)







となります．

一般化したりできますが，省略します．

オイラー法の様にしてこれを計算すると，非常に精度良く数値解を求める事が出来ます．

数学 #7 数値解析 #2 積分の近似計算

今回は，積分の近似計算について．

フレネル積分

フレネル積分は以下のように定義されています．





残念なことに，この積分は初等関数で表すことは出来ないようです．
というわけで近似計算．





↑では，
1.三角関数をマクローリン展開
2.項毎に積分
を行い，無限級数の形で解を出しています．

では，この級数のグラフを示しましょう．


Fig.1 Fresnel(cos)


Fig.2 Fresnel(sin)

第n項までの和のグラフを表示しています(nは1から40近くまで取ってあります)．
この級数展開の問題は，tを大きくすると，より多くの項を計算に入れないといけないという点にあります．
実際，nを大きく取っても，tの値を大きくすると直ぐに発散してしまいます．


Fig.3 Fresnel

x(t)とy(t)をxy平面にプロットすると，上のような螺旋模様が描かれます．


Fig.4

発散した時の図，アヴァンギャルドな感じが楽しめます．．

先に述べた問題が発生しない，他の近似計算法も沢山あるのですが，それは次回以降…

数学 #2 数値解析 #1 オイラー法

微分方程式が趣味の方は居るでしょうか．
日々閉じた解を追い求める，これ以上のロマンはありません．

しかし，行き詰ってしまうことも少なくないはず．そこはそれ，人生は旅です．気楽に行きましょう．

そんなわけでオイラー法です．

例えば，
y'=x²+y²
この微分方程式は，解を初等関数で表すことはできないようです．

ではどうするか．数値解を求めれば良いのです．

まず，初期値x₀，y₀，刻み幅hを決めます．

さて，
x₁=x₀+h
y₁=y(x₁)
としましょう．y(x₁)はy()が分からないので求められませんが，テイラー展開することが出来ます．

y₁=y(x₀+h)
=y₀+hy'₀+(h²/2!)y''_n…

始めの2項にのみ注目しましょう．(hが十分に小さいとすれば，それ以降を無視しても誤差は小さいです)
つまり，
x₁=x₀+h
y₁=y₀+hy'₀
と書くことが出来ます．

ちなみに，ここで，y'₀は
y'₀=y'(x₀, y₀)
ですのでご注意を．

最終的に，これを再帰的に記述すれば，
x_n+1=x_n+h
y_n+1=y_n+hf(x_n, y_n)
(但し，f(x, y)=y'(x, y))

となります．
これらから，初期値を決めてしまえば，数値的振る舞いが分かるのです．

下に例を示しましょう．


Fig.1
先の微分方程式について，それぞれの座標での傾きを表したものです(-2≦x≦2, -2≦y≦2)．これで，解の外形をイメージすることが出来ます．


Fig.2 x₀=-1.6, y₀=-1.8, h=1/4, 1/16
初期値を設定して"解いて"みました(-2≦x≦2, -2≦y≦2)．
hが小さいほど，誤差は小さくなります(hが1/2になると，誤差も1/2になる)．

ちなみにこれ，精度はそれ程良くないです．
もっと精度を良くしたいというのであらば，テーラー級数の項を増やすか，あの方法を使いましょう…と，次回へ引き延ばしてみたり．