【Julia入門】Juliaで分岐図を描画してカオスを感じる

今回はJuliaで分岐図を作成してみたいと思います。

分岐図についてですが例として以下のようなものになります。 これはロジスティック関数の分岐図でパラメーターrを増やしていくと値が発散してカオス減少が確認できるものになります。

ja.wikipedia.org

今回、同じくロジスティック関数を題材にJuliaで分岐図を作成してみます。

ロジスティック関数

題材を以下書籍から抜粋します。この本は計算モデルの作り方やその解析方法について初歩から詳しく解説されていて大変勉強になります。

A Biologist's Guide to Mathematical Modeling in Ecology and Evolution

和書での類書は以下になるかと思いますが、この洋書のほうが中身が圧倒的に濃いです。

数理モデリング入門 ―ファイブ・ステップ法― 原著第4版

以下の再帰式を使用します。この式は繁殖による生物の個体数増加をモデル化した式になります。

$$ n(t+1) = n(t) + r n(t) \left(1-\frac{n(t)}{K} \right) $$

n(t)はt時間経過後の個体数、rは内的自然増加率、Kは環境収容力になります。

この式をシミュレートして、分岐図を作成します。

分岐図の作成

今回の分岐図の作成は以下のように実施しました。

1. 時間t=1,2,3,……199,200（タイムステップ1で時間200まで計算）とし、 K=1000, 初期個体数N(0)=10とする。

2. rを0.01, 0.02, ……, 2.99, 3.00 と振り、その時の上記ロジスティック関数の個体数n(t)を各々計算する。

3. 各rに対し、n(t=181)~n(t=200)をそれぞれプロットする。

結果、以下のような分岐図が得られました。 無事カオス現象を確認することが出来ましたね。

簡素な内容ですが今回はここまでにします。 最後にシミュレーションに使用したコードを載せておきます。

using Plots
using LinearAlgebra

M = 200; #ステップ数
r = [0.01*x for x in 1:300]
K = 1000;
N0 = 10;#初期値

global Population = zeros(M,300)

function recursion(r,K,M,N0)
     A = zeros(M)
     A[1] = N0
    for i = 1:M-1
        A[i+1] = A[i] + r*A[i]*(1-A[i]/K)
    end
    return A
end

for i = 1:300
global Population[:,i] = recursion(r[i],K,M,N0)
end

for i = 1:300
global p = scatter!(r[i].*ones(20),Population[M-20+1:M,i],legend = false,color = "black",marker=:x, markersize = 0.3)
end

yaxis!("Population sizes")
xaxis!("Intrinsic growth rate (r)")
display(p)