[Python][Matplotlib] 変な形のアトラクターをつくる差分方程式系を適当に探してみる
差分方程式系を作る
二変数の差分方程式系を作れば良いので、つまりは
の形式であれば何でも良いということになる。しかし、大抵の写像は平衡点に収束するか爆発するかである。手を動かしていろいろな写像を作ってみると、「平衡点に収束しないが爆発もしない」という微妙なさじ加減がなかなかむずかしいとわかった。
そこで方針を変えて、既存の写像をちょっと変えて新たな写像を作ることにする。ここではロジスティック写像に変更を加えることにした。ロジスティック写像は次のような一変数の差分方程式。
ロジスティック写像は生物の個体数の変化を記述する数理モデルでもある。その場合、 は一個体が最大で生む子の数を表す。個体数が増えてくると、餌や住み場所がなくなるので、生物は増えにくくなる。この混み合いの効果は負の二次項によって表現されている。
ロジスティック写像をもとに、次のような差分方程式系を作ってみた。
も もロジスティック写像に似たような格好になっているが、二次項の部分に和 が入っているのがポイントである。後で示すパラメータでは となっている。ロジスティック写像がそうであったように、生物の個体数変化に即して考えると、 は増えるのは早いが他個体から負の影響を受けやすく、 は増えるのはゆっくりだが他個体からの負の影響を受けにくいという構造になっており、それでうまくバランスして2種の生物が共存してくれることを期待した。
作図する
最初に1000回計算した後で、 回計算して点を打つ方法で作図した。パラメータ と の値をいろいろいじってみたが、 のときが面白かった。
なんだか奇妙な形の範囲で点が動きつづけていることがわかる。大きく見て、上と下の二つ領域があるようにも見える。 のときの図を見ると、真ん中に大きな穴が二つあり、この領域には点が移動しないこともわかる(理屈はわからないが・・・)。
定量的な話をすると、xの値はときどき負の値をとっている。このことは、いま指定したパラメータ()では、この差分方程式を生物の個体数変化として解釈することはできないといういことを意味する。「うまくバランスして2種の生物が共存してくれることを期待した」のだが、少なくともいま指定したパラメータ()に関しては、この企みは失敗したわけである。ロジスティック写像では、 が不安定平衡点になっていて、 が負の値に突入すると は常に負になるので、はどんどん小さくなってに発散してしまうはずである。今回の差分方程式系でも同じようにに発散してしまうかと思いきや、なぜか 周辺に戻ってきている。
もう少し様子を見るために、次のような三次元の差分方程式に直して作図してみる。
式中のという値がポイントになると思ったので、これを各時刻で計算して、第三軸にプロットした。形が理解しやすいよう、z軸で回転したGIFアニメーションを作成した。
点の集合はひしゃげた長方形の形をしているようにも見える。z軸で上側の部分と下側の部分で大きく二つに分かれており、その間を行き来するルートが複数あるのが分かる。この図をz軸でつぶして上から見たのが最初の二次元の図であるから、xy平面で見たときに上と下の2領域があるという見立ては間違いではなかったことになる。
安定性解析
平衡点
がともに0でない平衡点なら
が成り立つ。ここで について整理すると
となる。この等式は特殊なパラメータのときでないと満たされない。当然、 では満たされない。したがって、ともに0でない平衡点 は存在しない。
適当に手計算をすると、平衡点は次の三つであることが分かる。
ソースコード
2次元
from tqdm import tqdm import matplotlib.pyplot as plt import numpy as np x0 = 0.5 y0 = 0.5 N = 10000 N_run = 1000 FILENAME = 'data.csv' a = 3.65 b = 3.6 def func(x, y): x_next = a * x * (1 - (x + y) / 2) y_next = b * y * (1 - (x + y) / 3) return (x_next, y_next) def calc(): x, y = (x0, y0) with open(FILENAME, mode='w') as f: for _ in range(N_run): x, y = func(x, y) for _ in tqdm(list(range(N))): f.write('{:f}, {:f}\n'.format(x, y)) x, y = func(x, y) return def plot_matplotlib(filename): data = np.loadtxt(FILENAME, delimiter=',') fig = plt.figure() ax = fig.add_subplot() ax.scatter(x=data[:, 0], y=data[:, 1], s=0.5, color='C0', antialiased=True) fig.savefig(filename, dpi=600) plt.close(fig) return def main(): calc() plot_matplotlib(filename='figure_a{:.2f}_b{:.2f}_N{:.1e}.jpg'.format(a, b, N)) return if __name__ == '__main__': main()
3次元
import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # NOQA import matplotlib.animation as animation import numpy as np from tqdm import tqdm x0 = 0.5 y0 = 0.5 N = 10000 N_run = 1000 FILENAME = 'data.csv' a = 3.65 b = 3.6 def func(x, y): x_next = a * x * (1 - (x + y) / 2) y_next = b * y * (1 - (x + y) / 3) z_next = x + y return (x_next, y_next, z_next) def calc(): x, y, z = (x0, y0, 0) with open(FILENAME, mode='w') as f: for _ in range(N_run): x, y, z = func(x, y) for _ in tqdm(list(range(N))): f.write('{:f}, {:f}, {:f}\n'.format(x, y, z)) x, y, z = func(x, y) return def plot_gif(): data = np.loadtxt(FILENAME, delimiter=',') fig = plt.figure() ax = fig.add_subplot(111, projection='3d') def plot(angle): ax.scatter(xs=data[:, 0], ys=data[:, 1], zs=data[:, 2], s=0.5, color='C0', antialiased=True) ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.view_init(10, angle) anim = animation.FuncAnimation( fig, plot, interval=300, repeat_delay=2000, frames=np.arange(0, 360, 10) ) anim.save('anim_a{:.2f}_b{:.2f}_N{:.1e}.gif'.format(a, b, N), writer="imagemagick") plt.close(fig) return def main(): calc() plot_gif() return if __name__ == '__main__': main()