ここでは、個体群動態の古典モデルの一つであるRicker Mapを用いて、差分方程式（時間離散モデル）の数値解法を解説する．目標は以下の分岐図を描くこと．

サンプルコードは"test_numpy.py"．

• モデルの概略
• 時系列データの生成
• 分岐図の作成

モデルの概略

Ricker modelは密度効果を考慮した1種類の生物に対する個体群動態モデルである．観測時刻間での個体群サイズの変動は、たった一つの生態パラメータrで特徴付けられる．下図の方程式で、Ntは観測時刻（タイムステップ）tにおける個体群サイズを表す．タイプステップ間の個体群サイズ変化率N_t+1/N_tはexp(r(1-N_t))ということで，個体数が大きいほど小さくなる負の密度効果を表している．グラフ中の赤い曲線は時刻が１進むときにN_tとN_t+1の間の関係を表している．緑色の線はN_t+1 = N_tとなる１：１の直線なので、赤線と緑線の交点はN_t+1 = N_tとなる「平衡点」となる．図はr = 1.5のときを表している．

ためしに最初の１０ステップの個体群動態をグラフにすると以下の通り．平衡点N*=1.0に近づいていくのが分かります．

時系列データの生成

このモデルの面白いところは、増殖率を表すrの値が多くなると、平衡点に近づいていくとは限らないこと．

では，実際にプログラムコートを書いていきます．まず3つのクラスをインポートします．

import numpy as np   #for numerical calculations

import math    #for mathematical functions

import matplotlib.pyplot as plt   #for graphics

実はnumpyを読み込んだだけでは，数学関数（たとえば、指数関数）を呼び出すのにいちいちメッソドの形を使う必要があります．

np.exp(1.0)

しかし、以下のようなスクリプトを追加するとダイレクトにexp()が使えます．

from numpy import exp

exp(1.0)

上のように単純に時系列を発生させるには，結果を保存する配列(v3)，モデルパラメータ(r)，個体群サイズ(x)の初期値を用意し，ループ関数の中で方程式をまわしていけばいい．各タイムステップでの個体群サイズ(x)は常に配列v3に保存するので、過去のデータを保存する必要は無く，x = x*exp(r*(1.0 - x))のように上書きすればよい。

#Ricker map

v3 = np.zeros([100,2])  #note that the index starts from zero, arrey to save results

r = 1.5  #growth parameter

x = 0.1 #initial population size

#save the initial conditions

v3[0,0] = 0

v3[0,1] = x

for i in range(1, 100):  #note that this is the loop for 1 <=i < 100

    x = x*exp(r*(1.0 - x))

    #print(i,"\t", x, "\n")

    v3[i,0] = i

    v3[i,1] = x

上の例ではv3は100行2列の二次元配列で、一列目にタイプステップの値(i)を、二列目に個体群サイズの値(x)を保存していく．for loopのインデックスiの始まりと終わりに注意します（i=100は含まれません）．ちなみにpythonでデフォルトの浮動点小数のサイズは64ビット（C言語のdouble相当）で，これは変数の一覧でfloat64として確認できる．

単純な二次元散布図、二次元線グラフはそれぞれ以下のメソッドで描写可能．これらは０～９までの値をプロットする場合．

#plot

plt.scatter(v3[0:10,0], v3[0:10,1], color="red")

plt.show() #rstudio上でプロットを表示するには必要らしい

plt.plot(v3[0:10,0], v3[0:10,1],color="blue", marker="o")

plt.show() #rstudio上でプロットを表示するには必要らしい

０～９９ステップまでのすべてのデータを使い，グラフのタイトルなどを追加する場合は以下のスクリプトを実行する．

plt.plot(v3[:,0], v3[:,1],color="blue", marker="o")

plt.title('population dynamics with r = 1.5')

plt.xlabel('time step t')

得られるグラフは以下の通り．

これらの数値計算をr=2.2, 2,6, 3.0の値にそれぞれ変更してやり直すと，以下のように2周期，4周期，カオス的変動へと変わっていくのが分かります．生態学的な示唆としては，増殖速度が大きくなるにしたがって個体数変動が大きくなるという不安定化現象を示しています．

分岐図(Bifurcation Map)の作成

それでは，rの値を小刻みに変えていき，どのように個体群動態の「最終状態」が変遷するかを見ていくアルゴリズムを下図を使ってまず説明する．

rの値を決めるごとに異なる個体群動態が生じる．力学系においては，初期値から出発してしばらくのあいだはトランジエントなフェーズで変動するが、時刻が十分大きくなったとき，何らかの「最終状態」(omega limit set)に近づく．これは平衡点、周期解、カオスアトラクター等のうちのどれかである．数値的にこの「最終状態」をみつけるにはtransient phaseの長さ (0 <= t < initial_record)を十分長く設定し，この期間の情報は捨て（配列などに保存せず）、final phase（initial_record <=t < end_time ）の情報だけを記録すればよい．このfinal phaseの長さも短すぎず（短すぎるとカオスアトラクターを捉えきれない）、長すぎず（長すぎるとアトラクター内の構造が見えなくなってしまう）ちょうどいい期間を取る必要がある．分岐図を描く際には，最終状態だけに注目すればよく時刻情報には興味が無いので、保存すべきは(r, x)の値のみである．

以下のスクリプトではまず数値計算する全期間(end_time)，アトラクタの情報保存を開始する時刻(initial_record)を指定すると共に、探索するrの値の範囲(min_r<=r<=max_r)も指定し、rを変更するステップ幅(step_r)も合わせて決める．きれいな分岐図が描けるようにこれらのパラメータ値の決定には若干の試行錯誤が必要である．上図の左右を比較すれば，たとえばr=2.6の時にはfinal phaseを十分長くとっても実質的に4つの値しか取らない．なぜなら4周期の解に収束するからである．

#Ricker map & bifurcation

end_time = 1000

initial_record = 800

min_r = 1.0

max_r = 4.5

step_r = 0.005

ここではそれぞれのrの値について最終(end_time - initial_record)ステップの値を保存するのことにするので、データ保存用の配列のサイズは以下のように決めることができる．

size_result = int((end_time-initial_record)*(max_r-min_r)/step_r)

配列サイズは整数で定義する必要があるのでint()関数を使っている．以下のコードでデータ保存用配列を定義する．

v4 = np.zeros([size_result,2])  #note that the index starts from zero

あとは二つのfor loopを組み合わせればOK. 外側のiに対するループはパラメータrを変更するためのもので，内側のk（=時刻）に対するループはそれぞれ指定されたrの値において，Ricker modelを逐次的に解いて最終状態だけを保存していく．

j = 0

r = min_r

for i in range(0, int((max_r-min_r)/step_r)):

    x = 0.1  #initial population size

    for k in range(1, end_time):

        x = x*exp(r*(1.0 - x))

        if k >= initial_record:

            v4[j, 0] = r

            v4[j, 1] = x

            j += 1 #renew j

    r += step_r #renew r

すべての結果は二次元配列v4に保存されており、一点一点のサイズに気をつけながら単純な散布図を描けばbifurcation mapのできあがり．

plt.scatter(v4[:,0], v4[:,1], s=0.01, color="blue", marker='x')

plt.title('Bifurcation map of Ricker model')

plt.xlabel('r value')

差分モデルの解説は以上で終わり．常微分方程式を用いた連続時間モデルについても、基本的には離散時間へと近似するので、基本的な数値計算のアルゴリズムはこれでOKということになる．