読者です 読者をやめる 読者になる 読者になる

技術メモ

役に立てる技術的な何か、時々自分用の覚書。統計学、機械学習、神経科学、暗号理論...と幅広く興味があります。

カオス理論入門(プログラム編)

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

実際のコード

汎用性の高いものになるよう意識したので、色々カスタマイズできるようにした。

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