[Python][Matplotlib] 変な形のアトラクターをつくる差分方程式系を適当に探してみる

はじめに

以前、グモウスキー・ミラの写像を作図したときに、差分方程式系のアトラクターはなかなか綺麗に描けるものだとわかった。

cyanatlas.hatenablog.com

そこで、自分でも適当に写像を作って、変な形のアトラクターを描いてみよう、と思ったのがきっかけ。

差分方程式系を作る

二変数の差分方程式系を作れば良いので、つまりは
x_{n+1}=f(x,y)
y_{n+1}=g(x,y)
の形式であれば何でも良いということになる。しかし、大抵の写像は平衡点に収束するか爆発するかである。手を動かしていろいろな写像を作ってみると、「平衡点に収束しないが爆発もしない」という微妙なさじ加減がなかなかむずかしいとわかった。

そこで方針を変えて、既存の写像をちょっと変えて新たな写像を作ることにする。ここではロジスティック写像に変更を加えることにした。ロジスティック写像は次のような一変数の差分方程式。

x_{n+1}=a x_n (1 - x_n)

ロジスティック写像は生物の個体数の変化を記述する数理モデルでもある。その場合、a は一個体が最大で生む子の数を表す。個体数が増えてくると、餌や住み場所がなくなるので、生物は増えにくくなる。この混み合いの効果は負の二次項によって表現されている。

ロジスティック写像をもとに、次のような差分方程式系を作ってみた。

x_{n+1}=a x_n (1 - (x_n + y_n) / 2)
y_{n+1}=b y_n (1 - (x_n + y_n) / 3)

xy もロジスティック写像に似たような格好になっているが、二次項の部分に和 x+y が入っているのがポイントである。後で示すパラメータでは a>b となっている。ロジスティック写像がそうであったように、生物の個体数変化に即して考えると、x は増えるのは早いが他個体から負の影響を受けやすく、y は増えるのはゆっくりだが他個体からの負の影響を受けにくいという構造になっており、それでうまくバランスして2種の生物が共存してくれることを期待した。

作図する

最初に1000回計算した後で、N 回計算して点を打つ方法で作図した。パラメータ ab の値をいろいろいじってみたが、a=3.65, b=3.6 のときが面白かった。

f:id:cyanatlas:20200718005128j:plain
N = 10^3のとき
f:id:cyanatlas:20200718005218j:plain
N = 10^4のとき
f:id:cyanatlas:20200718005246j:plain
N = 10^5のとき
f:id:cyanatlas:20200718005313j:plain
N = 10^6のとき

なんだか奇妙な形の範囲で点が動きつづけていることがわかる。大きく見て、上と下の二つ領域があるようにも見える。N=10^6 のときの図を見ると、真ん中に大きな穴が二つあり、この領域には点が移動しないこともわかる(理屈はわからないが・・・)。

定量的な話をすると、xの値はときどき負の値をとっている。このことは、いま指定したパラメータ(a=3.65, b=3.6)では、この差分方程式を生物の個体数変化として解釈することはできないといういことを意味する。「うまくバランスして2種の生物が共存してくれることを期待した」のだが、少なくともいま指定したパラメータ(a=3.65, b=3.6)に関しては、この企みは失敗したわけである。ロジスティック写像では、x=0 が不安定平衡点になっていて、x が負の値に突入すると  a x(1-x)は常に負になるので、xはどんどん小さくなって-\inftyに発散してしまうはずである。今回の差分方程式系でも同じように-\inftyに発散してしまうかと思いきや、なぜか0 周辺に戻ってきている。

もう少し様子を見るために、次のような三次元の差分方程式に直して作図してみる。

x_{n+1}=a x_n (1 - (x_n + y_n) / 2)
y_{n+1}=b y_n (1 - (x_n + y_n) / 3)
z_{n+1}=x_n + y_n

式中のx_n + y_nという値がポイントになると思ったので、これを各時刻で計算して、第三軸にプロットした。形が理解しやすいよう、z軸で回転したGIFアニメーションを作成した。

f:id:cyanatlas:20200718014654g:plain
N = 10^4のとき
f:id:cyanatlas:20200718025602g:plain
N = 10^5のとき

点の集合はひしゃげた長方形の形をしているようにも見える。z軸で上側の部分と下側の部分で大きく二つに分かれており、その間を行き来するルートが複数あるのが分かる。この図をz軸でつぶして上から見たのが最初の二次元の図であるから、xy平面で見たときに上と下の2領域があるという見立ては間違いではなかったことになる。

安定性解析

平衡点

 (x,y) がともに0でない平衡点なら
x=a x (1 - (x + y) / 2)
y=b y (1 - (x + y) / 3)
が成り立つ。ここで x+y について整理すると
2(1-\frac{1}{a}) = x + y = 3 (1 - \frac{1}{b})
となる。この等式は特殊なパラメータのときでないと満たされない。当然、a=3.65, b=3.6 では満たされない。したがって、ともに0でない平衡点 (x,y) は存在しない。

適当に手計算をすると、平衡点は次の三つであることが分かる。

 x=0, y=3
 x=2, y=0
 x=0, y=0

安定性

ヤコビ行列の固有値を計算すると、

 x=0, y=3 (-3.6, -1.825)
 x=2, y=0 (-3.65, 1.2)
 x=0, y=0 (3.6, 3.65)

となって、いずれの平衡点でも固有値の絶対値が1より大きく、不安定である。

f:id:cyanatlas:20200718023700j:plain
解析に使ったMathematicaのコード

ソースコード

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()