ロジスティック方程式を使った適当な数値シミュレーショングラフを作ってみる - 日報 #117
最近めっきりと技術的な内容がなくなってきております。良くない傾向です。
作業的なコーディングをしなければならない時はもちろんありますが、
かといって新しい知識のインプットや、
試行錯誤のような脳の運動を絶やしてはならないわけです。
数値シミュレーションをすることがちょっと増えたので、
本日はその辺りの話を。
数値計画とそのシミュレーションをやっていて、
すっっっっごい雑に曲線を作りたくなる時ってあると思います。
例えばじゃぁ、、、
1月に1だったものを、
12月までに10にしまぁす！
みたいな「最初の値と最後の値と、それを適当な数のプロットで結んで、何となく曲線にする」ような話。
意味のある関数で描ければそりゃぁ最高なんですが、
とにかくまぁ一旦それっぽく繋げたくなる時ってほら、あるじゃないですか。
あったらあんまよくないんですけどね！
と、とにかく、そういう時に、
皆さんどうやって曲線を描いてますかね？
前段の例で考えていきましょう。
通る２点は決定しています。
横軸を月、縦軸を値としたら、 (1, 10) と (12, 10) ですね。
これを直線で引くと、
y = 9/11 * x + 2/11
が直線成長するモデルになります。
グラフはこんな感じ。
では、これを二次関数にしてみます。
y = ax^2 + bx + c
なので、3点が決まらないとa,b,cは確定しませんね。
ということで、a = 1として計算してみますと、
y = x^2 - (134/11)*x + 134/11
もちろんこれでは大変おかしなグラフになります。
変曲点が二点の間にあって下に凸ってことは、そりゃそうなっちゃいます。
だから下に凸にするなら、変曲点を左より、例えば(0,0)にもってくるようにします。
するとb=0なので、これで計算すると、
y = (9/143) x^2 + 134/143
だいたいこんな感じです。
「2次曲線的な成長」とかいう、訳分からん成長ならこんな感じですね。
しかし工夫がありませんね。
もう少し便利な関数はないものでしょうか。
俗にいう、「指数関数的な成長」ならどうでしょうか。
丁度2点を通るためには2元1次の方程式程度になっている必要があるので、
y = x^a * e^b
とか適当に考えてみます。
計算しますと、
y = x^(ln10/ln12)
こんな感じです。
うーん、イマイチですね。。。
もっとパラメータをいじって自由に形を変える関数のほうがよさそうです。
「ある値からある値へのデジタルな変化を平滑化する」
といえば電子回路でお馴染みシグモイド関数です。
ということで、ロジスティック曲線を使ってみます。
ロジスティック曲線の式は、
dN/dt = r ((K-N)/K ) N
です。生物群の個体数Nの変化量は、その時のNの値に「何らかの値」を積算したもの、という発想から生まれた方程式です。
その「何らかの値」を増殖率rとかで表現しきるのもいいのですが、
個体には必ず「ブレーキ」が存在します。
「ブレーキ」は、その個体が存在しうる資源、収容能力が沢山余っているほどゆるく、
満杯に近づくほど強く、かかればいいわけなので、
増殖率rに対して、「収容能力の余力」を積算すればいいわけです。
「収容能力の余力」は「限界収容能力に対する空きの割合」で単純にOKなわけなので、
限界収容能力をKとしたら、(K-N)/K となるわけですね。
なので、それに増殖率rをかけて、冒頭の「何らかの値」とすればいいので、結局、
dN/dt = r ((K-N)/K) N
となるわけです。
んで、これを解いていくわけですが、
一般にこの微分方程式を解くとこんな感じになります。
N(t) = K / (1 + exp(rK(t0 - t))
今回のケースで言うと、
増殖率rは、適当な値でOK。なんとなーく、月ごとの成長率くらいに考えて、0.2とかにしておきましょう。
限界収容能力Kは、じゃぁ最終的に到達する値10にしましょうか。
初期値t0が厄介ですが、今回は1月で1というのを初期値としたいので、N(1)=1を解いて、
t0 = 1 + ln(K-1)/rK
⇔ t0 = 1 + ln(9)/2
なので書き下すと、
y = 10 / (1 + exp(2(1+ln(9)/2 - x))
なんとなくそれっぽいですね。
ただ、結構勢いが急ですね。。。
なので、適当にrをいじっていきましょう。
r=0.05とかにしてみると、一気にこんな感じ。
ところでrは元々下記の式において定性的意味を付与されています。
dN/dt = r * N
要は、その時の個体数Nに対して、
次の瞬間「どの程度の割合増えるのか？」
ということです。
なにもブレーキがなかったらどの程度の成長率でしょうか。
N = exp(rt+c)
を(1,1), (12,10)で解くと、
r = ln(10)/11 = 0.21
くらいなんですね。
ちなみにその時の式はこんな感じです。
y = exp((ln(10)/11)*x - ln(10)/11)
これだと、急激に成長してしまいます。
イメージこんな感じ。
これをもっとゆるやかに成長させて、かつ、ゆるやかにブレーキさせることを考えると、
だいたいその半分から1/4程度にしてみるといい感じになります。
ここまでやっといてなんですけど、
この数値計算、マジで、意味は無いです。
とにかく、それっぽく見せたい！って時は、
y = K / (1 + exp(Kr(1+ln(K-1)/(Kr) - x))
という式を作って、Kを目標値にして、rを成長率にで適当にやれば、
まぁなんとなくいい感じのグラフ、かけまっせ！！
以上！
