技術メモ

役に立てる技術的な何か、時々自分用の覚書。幅広く色々なことに興味があります。

カオス理論入門(pythonでアトラクタの軌跡を描く)

カオスについて勉強していた時に書いて見たかったもの。アトラクターの軌道を見てみたい。
目標としては3次元プロットして一定軌道をグルグル回っているところをプロットしてみる。有名なレスラーアトラクターとローレンツアトラクターを描いていきたい。

実際のコード

汎用性の高いものになるよう意識したので、色々カスタマイズできるようにした。
Attractorクラスを継承すれば、描画や実行は新しく書く必要がないようにした。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pdb


class Attractor:
    def __init__(self):
        self.dt = 0.01
        self.x = self.y = self.z = 0.1
        self.output = []

    def initialize(self,x,y,z):
        self.x = x
        self.y = y
        self.z = z

    def run(self,T):
        for t in range(T):
            self.step()
            self.output.append([self.x,self.y,self.z])
        return self.output

    def step(self):
        pass

    def draw(self,display=True,plotfig=None):
        if not self.output:
            print("not run...")
            return 0
        if display:
            fig = plt.figure()
            ax = fig.gca(projection='3d')
            ax.plot(*zip(*self.output))
            plt.show()
        else:
            try:
                plotfig.plot(*zip(*self.output))
            except:
                print("draw Error...")
                print("USAGE: draw(display=False,plotfig=<object plt.figure().gca>)")


class Lorenz(Attractor):
    def __init__(self,p,r,b):
        # p=10,r=28,b=8/3 in Lorenz(1963)        
        Attractor.__init__(self)
        self.p = float(p)
        self.r = float(r)
        self.b = float(b)

    def step(self):
        x,y,z = self.x,self.y,self.z
        self.x += self.dt*(-self.p*(x-y))
        self.y += self.dt*(-x*z + self.r*x - y)
        self.z += self.dt*(x*y - self.b*z)
        return self.x, self.y, self.z


class Rossler(Attractor):
    def __init__(self,a,b,c):
        # a=0.2,b=0.2,c=5.7 in Rossler(1976)
        Attractor.__init__(self)
        self.a = float(a)
        self.b = float(b)
        self.c = float(c)

    def step(self):
        x,y,z = self.x,self.y,self.z
        self.x += self.dt*(-y-z)
        self.y += self.dt*(x+self.a*y)
        self.z += self.dt*(self.b+x*z-self.c*z)
        return self.x, self.y, self.z


if __name__=="__main__":
    fig = plt.figure(figsize=plt.figaspect(0.5))
    ax = fig.add_subplot(1,2,1,projection="3d")
    la = Lorenz(10,28,8.0/3)
    la.run(10000)
    la.draw(display=False,plotfig=ax)
    ax2 = fig.add_subplot(1,2,2,projection="3d")
    rl = Rossler(0.2,0.2,5.7)
    rl.run(100000)
    rl.draw(display=False,plotfig=ax2)
    plt.show()

実行結果

f:id:swdrsker:20161208195826p:plain

参考にしたサイト:
Lorenz system - Wikipedia
Rössler attractor - Wikipedia