scipyで常微分方程式を解く(ローレンツ方程式)

お久しぶりです

しばらくブログを放置ぎみだったのでたまには更新します．

はじめに

常微分方程式をscipyのodeintを用いて解くことが多いのですが，使い方を忘れてしまうことが多いためメモとして残しておきます． 本業でのシミュレータをお見せするわけにはいかないので，今回はカオス理論で有名なローレンツ方程式を解くようなコードを書いてみました．というのも，先日，金田邦彦さんの「カオスの紡ぐ夢の中で」を久々に読み返したときに表紙にローレンツ方程式が描かれており，綺麗な図形だなと思って自分でも書いてみたくなったからです．

ローレンツ方程式

wikipediaを参考にして．．． ローレンツ方程式は以下の3つの変数を解くことで表現できます．

dx/dt=-px+py

dy/dt=-xz+rx-y

dz/dt=xy-bz

このシンプルな3つの微分方程式はもともと天気の気象モデルを現わすためのものですが，コンピュータシミュレーション上でカオス的な振舞いをすることが判明してから，ローレンツカオスという呼ばれかたをされています．

ローレンツが方程式を論文に掲載した際には

p=10

r=28

b=8/3

というパラメータにしたそうです． 今回のコードもこのパラメータにしたがってコードを書いていきます． なお時間と3変数の初期値についてはdt=0.01,t=100,x=y=z=0.1にしました．

オイラー法によるシミュレーション

まずは，僕がよく常微分方程式を大雑把に解きたいときに使うオイラー法での解法を載せておきます．

#coding:utf-8
import matplotlib.pyplot as plt
import numpy
from mpl_toolkits.mplot3d import Axes3D

class LorenzAtractor(object):
    def __init__(self,p,r,b,t,dt,x,y,z):
        self._p=p
        self._r=r
        self._b=b
        self._t=t
        self._dt=dt
        self._x=x
        self._y=y
        self._z=z
    def calc(self):
        resultx=[]
        resulty=[]
        resultz=[]
        F_dx=lambda value: -self._p*value[0]+self._p*value[1]
        F_dy=lambda value: -value[0]*value[2]+self._r*value[0]-value[1]
        F_dz=lambda value: value[0]*value[1]-self._b*value[2]
        for i in range(0,len(self._t)):
            xx=F_dx([self._x,self._y,self._z])*self._dt
            yy=F_dy([self._x,self._y,self._z])*self._dt
            zz=F_dz([self._x,self._y,self._z])*self._dt
            self._x+=xx
            self._y+=yy
            self._z+=zz
            resultx.append(self._x)
            resulty.append(self._y)
            resultz.append(self._z)
        self.show(resultx,resulty,resultz)
    def show(self,x,y,z):
        fig=plt.figure()
        ax=fig.add_subplot(111,projection='3d')
        ax.plot(x,y,z,'o')
        plt.show()

p=10.0
r=28.0
b=8/3.0
x=0.1
y=0.1
z=0.1
dt=0.01
t=numpy.arange(0,100,dt)
LorenzAtractor(p,r,b,t,dt,x,y,z).calc()

結果

scipyのodeintによるシミュレーション

次にscipyのodeintを使ったコードです．

#coding:utf-8
import scipy.integrate
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

class LorenzAtractor(object):
    def __init__(self,p,r,b,t,x,y,z):
        self._p=p
        self._r=r
        self._b=b
        self._t=t
        self._x=x
        self._y=y
        self._z=z

    def Func(self,result,t0):
        F_dx=lambda result: -self._p*result[0]+self._p*result[1]
        F_dy=lambda result: -result[0]*result[2]+self._r*result[0]-result[1]
        F_dz=lambda result: result[0]*result[1]-self._b*result[2]
        return [F_dx(result),F_dy(result),F_dz(result)]
       
    def calc(self):
        initval=[self._x,self._y,self._z]
        result=scipy.integrate.odeint(self.Func,initval,self._t)
        self.show(result[:].T[0],result[:].T[1],result[:].T[2])

    def show(self,x,y,z):
        fig=plt.figure()
        ax=fig.add_subplot(111,projection='3d')
        ax.plot(x,y,z,'o')
        plt.show()
p=10.0
r=28.0
b=8/3.0
x=0.1
y=0.1
z=0.1
t=numpy.arange(0,100,0.01)
LorenzAtractor(p,r,b,t,x,y,z).calc()

結果

少し考察してみる

見比べると同じパラメータでも各円にある中央の穴のサイズが違うものなんですね． オイラー法ではやっぱり曲線を描くのが苦手な側面が現われているのでしょうか．

配列の定義をしてみたり，オイラー法で解くところではfor文を利用したりと，もっとpythonicな書きかたもできるような気もします．それに比べてodeintのほうでは一行で済んでますね． このあたりの差が長い時間計算するときに効いてきそうです．

最後に

このくらいの計算量なら数秒で描画されますし，Pythonで書くとずいぶん楽に常微分方程式の近似解を得ることができるものだと関心しました． カオスの面白いところはこのようなシンプルな式に様々なパラメータを代入するだけで予想もつかないような何か意味がありそうな図形が出てくるところですね．もっとも美しい式として知られる式に「オイラーの公式」がありますが，このようなカオス・フラクタル図形は誰が見ても美しいと感じますね．気が向いたら他の図形も描写してみようと思います．

それから今回はmatplotlibを使ったのですが，簡単にグラデーションが効いたような綺麗なグラフは作りたいです．Processingあたりを見ていると点が輝いていているので，そんなグラフも作ってみたいなと思いました．