Euler法とRunge-Kutta法のカオス時系列データの比較

スポンサードリンク

モチベーション

先日、非線形力学系ローレンツ方程式のアトラクタ描画をルンゲクッタ法で計算しました。x,y,zに関する時系列データも簡単に計算できたので、描画してみたところ以下のようなグラフになりました。以下はx-tグラフです。

上記のグラフがカオス的な振る舞いを示しているというのは前回の記事に示した通りですが、これがルンゲクッタではなく、違う近似解法の場合どういった挙動を示すのか気になりました。

ネットを見る限りオイラー法とルンゲクッタ法の解の誤差について、示している記事もいくつかあったのですが、１変数の一階微分方程式や単純なバネのモデルのような１変数の二階微分方程式しか見つかりませんでした。

なので、今回は３変数関数のローレンツ方程式を用いて、数値解法が異なることによる、誤差の程度を検証したいと思います(定量評価はしてません)。

Eular法とは

微分方程式の近似解法の中でも比較的計算が簡単なものですが、その分、とても粗い近似解法です。

研究等ではまず使われないですが、数値解放の概念を理解する手助けにはなります。以下のような方程式を考えます。
\begin{eqnarray}
\frac{dy}{dt}=\dot{y}=f(t,y(t))
\end{eqnarray}

次に微分の定義を用いて、
\begin{eqnarray}
\dot{y}=\lim_{\Delta{t} \to 0} \frac{y(t+\Delta t)-y(t)}{\Delta t}
\end{eqnarray}
0<t<<1の時
\begin{eqnarray}
\dot{y}=\frac{y(t+\Delta t)-y(t)}{\Delta t}
\end{eqnarray}
と表せる。これを最初の式を用いて式変形していくと、
\begin{eqnarray}
\dot{y}\Delta t = y(t+\Delta t)-y(t) \\
y(t+\Delta t) = y(t)+\Delta tf(t,y(t))
\end{eqnarray}
となり、t秒の情報(y(t)、f(t,y(t))を用いて、Δt秒後のy(t+Δt)を求めるのがEular法です。

計算ではこのΔtを小さくするほど、微分に近づくことになるので精度はよくなりますが、その分、計算時間が長くなります。

Runge-Kutta法とは

ルンゲクッタ法は、計算に少々時間がかかりますが、その精度の高さから、微分方程式の数値解析には広く用いられています。

今回計算に用いているのは、４次のルンゲクッタ法です。4次のルンゲクッタ法には公式が存在するのですが、その公式の導出をしている人はどうやらいないみたいです。

公式の導出にはオイラー法とホイン法（２次のルンゲクッタ法）を用いなければいけないみたいですが、計算量がえげつないみたいで、そんなものを書こうとするマゾな人はいないらしいです。

簡単に言うとテイラー展開の４次までの項の係数と一致するように計算するらしいです。まあいろんな人が書いているので、ここで詳しく書く必要もないでしょう。

どうしても知りたい人はこのサイトを参考にしてください笑。→http://eman-physics.net/math/differential21.html

スポンサードリンク

計算結果比較

以下に、オイラー法とルンゲクッタ法で計算したx-tグラフを示します。

最初の方（t＝０から５００）くらいでは、二つの数値解法にはあまり差がありませんが、徐々に差は大きくなり、やがて全く異なる挙動になります。

試しにt＝０から５００の範囲のグラフを載せます

解に大きなズレはなくこの程度の範囲ならば、オイラー法でも良さそうです。

しかし実際の科学計算ではオイラー法はまず用いられません。研究で使用されている方程式はより高次で、高変数のものが多く、計算範囲ももっと広い場合がほとんどです。

なので、この程度の方程式かつこの程度の時間範囲しか解の正確性を保てないオイラー法はまず使われません。

オイラー法は比較的計算も楽で、簡単に理解できるので、数値解法の玄関的な役割として、教科書等に書かれているだけです。

ついでに、１変数関数であるならば、ある程度計算時間が長くてもオイラー法とルンゲクッタ法の解に差はあまり見られません。以下は一次元の減衰系振動モデルのx-tグラフです。

若干紫っぽくなってるのは、赤と青が重なっているためです。一次元の簡単な運動方程式程度であれば、二つの解放に誤差が出ないことが観察できます。

今回用いたプログラム

ローレンツ方程式

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

# lorenz equation
def lorenz(x, y, z, p=10, r=28, b=8/3):
    x_dot = -p*x +p*y
    y_dot = -x*z +r*x-y
    z_dot = x*y -b*z
    return np.array([x_dot, y_dot, z_dot])


t=0             ##### 初期時間######
dt = 0.01       ##### 微分間隔######
stepCnt =500 ######時間ステップ##




########ルンゲクッタ####################
# Need one more for the initial values
xs = np.empty((stepCnt + 1,))
ys = np.empty((stepCnt + 1,))
zs = np.empty((stepCnt + 1,))

# Setting initial values
xs[0], ys[0], zs[0] = (0.1, 0.5, 1.2)
k0=[0,0,0]
k1=[0,0,0]
k2=[0,0,0]
k3=[0,0,0]
listData=[0] 


for i in range(stepCnt):
    x,y,z=xs[i],ys[i],zs[i]

    k0 = dt * lorenz(x,y,z)
    k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)
    k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)
    k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

    dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0
    dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0
    dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0
    xs[i+1] = xs[i] + dx
    ys[i+1] = ys[i] + dy
    zs[i+1] = zs[i] + dz
      
    listData.append(i)  


########オイラー法####################

# Need one more for the initial values
xs2 = np.empty((stepCnt + 1,))
ys2 = np.empty((stepCnt + 1,))
zs2 = np.empty((stepCnt + 1,))
# Setting initial values
xs2[0], ys2[0] ,zs2[0]= (0.1, 0.5,1.2)
k02=[0,0,0]
k12=[0,0,0]
listData1=[0] 
for i in range(stepCnt):
    x,y,z=xs2[i],ys2[i],zs2[i]

    k02 = dt * lorenz(x,y,z)
    k12 = dt * lorenz(x+k02[0], y+k02[1] ,z+k02[2])
    dx = (k02[0]+k12[0])/2.0
    dy = (k02[1]+k12[1])/2.0
    dz = (k02[2]+k12[2])/2.0
    xs2[i+1] = xs2[i] + dx
    ys2[i+1] = ys2[i] + dy
    zs2[i+1] = zs2[i] + dz
    listData1.append(i)      
    

###########all-time-datasiries(graph-setting)##################     
plt.plot(listData, xs, c='b', label="Runge-Kutta")
plt.plot(listData1, xs2, c='r', label="Eular")
ax = plt.gca() 
ax.set_ylabel("X(t)")
ax.set_xlabel("time")
ax.set_title("x-t graph")
plt.legend() 
plt.show()

# -*- coding: utf-8 -*-

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

# lorenz equation

def lorenz(x, y, z, p=10, r=28, b=8/3):

x_dot = -p*x +p*y

y_dot = -x*z +r*x-y

z_dot = x*y -b*z

return np.array([x_dot, y_dot, z_dot])

t=0 ##### 初期時間######

dt = 0.01 ##### 微分間隔######

stepCnt =500 ######時間ステップ##

########ルンゲクッタ####################

# Need one more for the initial values

xs = np.empty((stepCnt + 1,))

ys = np.empty((stepCnt + 1,))

zs = np.empty((stepCnt + 1,))

# Setting initial values

xs[0], ys[0], zs[0] = (0.1, 0.5, 1.2)

k0=[0,0,0]

k1=[0,0,0]

k2=[0,0,0]

k3=[0,0,0]

listData=[0]

for i in range(stepCnt):

x,y,z=xs[i],ys[i],zs[i]

k0 = dt * lorenz(x,y,z)

k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)

k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)

k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0

dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0

dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0

xs[i+1] = xs[i] + dx

ys[i+1] = ys[i] + dy

zs[i+1] = zs[i] + dz

listData.append(i)

########オイラー法####################

# Need one more for the initial values

xs2 = np.empty((stepCnt + 1,))

ys2 = np.empty((stepCnt + 1,))

zs2 = np.empty((stepCnt + 1,))

# Setting initial values

xs2[0], ys2[0] ,zs2[0]= (0.1, 0.5,1.2)

k02=[0,0,0]

k12=[0,0,0]

listData1=[0]

for i in range(stepCnt):

x,y,z=xs2[i],ys2[i],zs2[i]

k02 = dt * lorenz(x,y,z)

k12 = dt * lorenz(x+k02[0], y+k02[1] ,z+k02[2])

dx = (k02[0]+k12[0])/2.0

dy = (k02[1]+k12[1])/2.0

dz = (k02[2]+k12[2])/2.0

xs2[i+1] = xs2[i] + dx

ys2[i+1] = ys2[i] + dy

zs2[i+1] = zs2[i] + dz

listData1.append(i)

###########all-time-datasiries(graph-setting)##################

plt.plot(listData, xs, c='b', label="Runge-Kutta")

plt.plot(listData1, xs2, c='r', label="Eular")

ax = plt.gca()

ax.set_ylabel("X(t)")

ax.set_xlabel("time")

ax.set_title("x-t graph")

plt.legend()

plt.show()

一次元減衰振動モデル

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


k = 1.0
m = 1.0
w = 0.01

tr=[]
xr=[]
te=[]
xe=[]
#############ルンゲクッタ###################
def g(t, x):
    return np.array([x[1], (-k/m*x[0])-2*w*x[1]])

x = np.array([0.01, 0.0])    # 初期条件 x(0) = 1, v(0) = 0
dt = 0.1

for i in range(10000):
    t = i*dt

    k1 = g(t, x)
    k2 = g(t + dt/2, x + k1*dt/2)
    k3 = g(t + dt/2, x + k2*dt/2)
    k4 = g(t + dt, x + k3*dt)
    x += (k1 + 2*k2 + 2*k3 + k4)*dt/6
    xr.append(x[0])
    tr.append(t)   

############オイラー法#########################
x1 = np.array([0.01, 0.0])    # 初期条件 x(0) = 1, v(0) = 0 
for i in range(10000):
    t = i*dt

    k1 = g(t, x1)
    k2 = g(t + dt, x1 + k1*dt)
    x1 += (k1 + k2)*dt/2
    xe.append(x1[0])
    te.append(t)
           
fig = plt.figure()
ax = plt.gca() 
plt.plot(tr, xr, c='b', label="Runge-Kutta")
plt.plot(te, xe, c='r', label="Eular")
ax = plt.gca() 
ax.set_ylabel("X(t)")
ax.set_xlabel("time")
ax.set_title("x-t graph")
plt.legend() 
plt.show()

# -*- coding: utf-8 -*-

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

import math

k = 1.0

m = 1.0

w = 0.01

tr=[]

xr=[]

te=[]

xe=[]

#############ルンゲクッタ###################

def g(t, x):

return np.array([x[1], (-k/m*x[0])-2*w*x[1]])

x = np.array([0.01, 0.0]) # 初期条件 x(0) = 1, v(0) = 0

dt = 0.1

for i in range(10000):

t = i*dt

k1 = g(t, x)

k2 = g(t + dt/2, x + k1*dt/2)

k3 = g(t + dt/2, x + k2*dt/2)

k4 = g(t + dt, x + k3*dt)

x += (k1 + 2*k2 + 2*k3 + k4)*dt/6

xr.append(x[0])

tr.append(t)

############オイラー法#########################

x1 = np.array([0.01, 0.0]) # 初期条件 x(0) = 1, v(0) = 0

for i in range(10000):

t = i*dt

k1 = g(t, x1)

k2 = g(t + dt, x1 + k1*dt)

x1 += (k1 + k2)*dt/2

xe.append(x1[0])

te.append(t)

fig = plt.figure()

ax = plt.gca()

plt.plot(tr, xr, c='b', label="Runge-Kutta")

plt.plot(te, xe, c='r', label="Eular")

ax = plt.gca()

ax.set_ylabel("X(t)")

ax.set_xlabel("time")

ax.set_title("x-t graph")

plt.legend()

plt.show()

ルンゲクッタ法とオイラー法で別々にリストを作って、別々にfor文を回しているので、すごい効率の悪いプログラムです。

if文等用いて変数を減らし、一つのfor文ですむように書き直したいと思います。

スポンサードリンク

最後に

今回は３変数関数におけるオイラー法とルンゲクッタ法の数値計算の誤差について述べました。

数値解法は他にもホイン法、ベルレ法、かえる跳び法、修正オイラー法等あるので、それらを用いて計算し、解の比較をしていきたいと思います。

ついでにルンゲクッタ法よりも僕は、計算回数と長期正確性を加味して、速度ベルレ法が一番おすすめです笑。

記事が役に立ったらシェア！

ローレンツ方程式をRunge-Kutta法（python）で計算

スポンサードリンク

前回は、アトラクターの勉強ということで、非線形力学系のレスラー方程式のアトラクター再構成の数値実験を行いました。

>>>前回の記事「Rossler方程式を Runge-Kutta法を用いて計算してみた（言語はpython）。」

>>>前々回の記事「多変量時系列解析に関する論文」

今回はローレンツ方程式について実施した数値実験の結果と、時間遅れの大きさについて論じて行きます。

時間遅れ座標系を用いた、アトラクタの再構成

まず、初めにローレンツ方程式を示します。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= -px+py\\
\frac{dy}{dt}= -xz+rx-y\\
\frac{dz}{dt}= xy-bz
\end{array}
\right.
\]
x,y,z 全ての変数の時系列データが観測できた場合のアトラクタを以下に示します。

アトラクタ再構成のやり方には

微分座標系(V(t)=($x(t),\dot x(t),\ddot y(t)$)を用いる
時間遅れ座標系(V(t)=($x(t),x(t+\delta),x(t+2\delta)$))を用いる

方法があるが、ノイズがある場合には、微分座標系はノイズを拡大させる原因にもなるので、時間遅れ座標系を用いる場合が多いです。

以下に観測できるデータがx(t)、時間遅れの大きさ$\delta=1,10,100,500$の場合での再構成したアトラクタを示します。

そして力学系のアトラクタは平衡点、リミットサイクル、kトーラス、ストレンジアトラクタ、ランダムに分類されます。

状態空間（位相構造）を観察する限りでは、ローレンツ方程式のアトラクタはストレンジアトラクタというものに分類されます。

このアトラクターが「ストレンジアトラクターになっている」というのが大切で、これが観察できれば、この非線形力学系は「カオス的な振る舞いをする」と言えます。

ターケンスの埋め込み定理

しかし、上記のような時間遅れ座標系を用いてアトラクタを再構成するには、どうやら以下のような定理が成り立たないといけないみたいです。

「d 次元のコンパクトな多様体 M と C ¹ 級微分同相写像 f : M → M，C ¹ 級関数 g : M → R¹ が与えられたとき， m = 2 d + 1 であれば，次式の写像 V : M → R ^m は，埋め込みである」

だめだ、わけワカンねぇ。本当に同じ言語なのかと疑うほどであるが、m=2d+1が成り立てば良いらしい。

ここでmとはアトラクタの次元のことで、ここではm=3になります。

問題はdです。この「ターケンスの埋め込み定理」が成り立つには、d＝１であれば良いのですが、僕はこのdが何を意味するのかわかりません。可観測な状態変数の数であるならば、今回は可観測なデータがx(t) のみを仮定しているので、成り立ちますがどうなのでしょう。

ここら辺、追求したいですが、体力が切れてきたので、ここはまた今度にします。

しかし前回読んだ論文では、ターケンスの埋め込み定理はあくまで十分条件らしいので、必ず成り立っていないといけないみたいではないようです。

スポンサードリンク

時間遅れの大きさの最適化

時間遅れの大きさの最適化には以下の二つの方法があります。

線形自己相関関数が最初に０になる点を求める
平均相互情報量が最初に０なる点を求める

以下では平均相互情報量の方法について述べて行きます。まずは定義から

\begin{eqnarray}
\hat{I}_{AB}(T)=\sum_{a_{i},b_{k}} P_{AB}(a_{i},b_{k})I_{AB}(a_{i},b_{k})
\end{eqnarray}
A={x(n)},B={x(n+T)}とすれば、
\begin{eqnarray}
\hat{I(T)}=\sum_{n=1}^N P(x(n),x(n+T))\log2(\frac{P(x(n),x(n+T))}{P(x(n))P(x(n+T))})
\end{eqnarray}
となります。そしてこの$\hat{I(T)}$ の最初（シグマなので）の最小値をとる時のTを「最適な遅れ時間」として選択します。

そして、この遅れ時間の最適化の計算はプログラムを書いて計算し、また次回以降載せます。

計算するとδ=18くらいになるみたいです。上の再構成アトラクタ（図２〜図５）もδ=10が一番元の位相構造に近いので、適切な手法と言えそうです。

計算に用いたプログラム

以下に計算に用いたプログラムを載せます。

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


# rorenz equation
def lorenz(x, y, z, p=10, r=28, b=8/3):
    x_dot = -p*x +p*y
    y_dot = -x*z +r*x-y
    z_dot = x*y -b*z
    return np.array([x_dot, y_dot, z_dot])


t=0             ##### 初期時間######
dt = 0.01       ##### 微分間隔######
stepCnt = 10000 ######時間ステップ##
delta = 500      ######遅延時間＃####


# Need one more for the initial values
xs = np.empty((stepCnt + 1,))
ys = np.empty((stepCnt + 1,))
zs = np.empty((stepCnt + 1,))

# Setting initial values
xs[0], ys[0], zs[0] = (0.1, 0.5, 1.2)
k0=[0,0,0]
k1=[0,0,0]
k2=[0,0,0]
k3=[0,0,0]

########x(t)のデータ作成####################
f = open('xdata.csv', 'w') 
csvWriter = csv.writer(f)
# Stepping through "time".


for i in range(stepCnt-2*delta):
    x,y,z=xs[i],ys[i],zs[i]

    k0 = dt * lorenz(x,y,z)
    k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)
    k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)
    k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

    dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0
    dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0
    dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0
    xs[i+1] = xs[i] + dx
    ys[i+1] = ys[i] + dy
    zs[i+1] = zs[i] + dz
    listData = []   
    listData.append(x)  
    csvWriter.writerow(listData)  
    
    
########x(t-δ)のデータ作成####################
f = open('xdata2.csv', 'w') 
csvWriter = csv.writer(f)
# Stepping through "time".
    
for i in range(stepCnt-delta):
    x,y,z=xs[i],ys[i],zs[i]

    k0 = dt * lorenz(x,y,z)
    k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)
    k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)
    k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

    dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0
    dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0
    dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0
    
    xs[i+1] = xs[i] + dx
    ys[i+1] = ys[i] + dy
    zs[i+1] = zs[i] + dz
    listData2 = [] 
    if i >= 500: 
        listData2.append(x)  
        csvWriter.writerow(listData2)  

    
########x(t-2δ)のデータ作成####################
f = open('xdata3.csv', 'w') 
csvWriter = csv.writer(f)
# Stepping through "time".
    
for i in range(stepCnt):
    x,y,z=xs[i],ys[i],zs[i]
    
    k0 = dt * lorenz(x,y,z)
    k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)
    k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)
    k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

    dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0
    dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0
    dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0
    
    xs[i+1] = xs[i] + dx
    ys[i+1] = ys[i] + dy
    zs[i+1] = zs[i] + dz
    listData3 = []   
    if i >= 1000:
        listData3.append(x)  
        csvWriter.writerow(listData3)  

  
###########all-time-datasiries(graph-setting)##################     
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(xs, ys, zs,c='m')
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor")
f.close()
plt.show()
plt.savefig("Lorenz-alltimesireies.png")    

f = open('xdata.csv', 'rb')
f = open('xdata2.csv', 'rb')
f = open('xdata3.csv', 'rb')

##########embedding-theory-datasireis(graph-setting)#############
x = np.loadtxt("xdata.csv", unpack=True)
y = np.loadtxt("xdata2.csv", unpack=True)
z = np.loadtxt("xdata3.csv", unpack=True)
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(x, y, z,c='m')
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Takens embedding '$\delta=500$'")
f.close()
plt.show()
plt.savefig("Lorenz-embedding-500.png")

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

# -*- coding: utf-8 -*-

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

import csv

import matplotlib.animation as animation

# rorenz equation

def lorenz(x, y, z, p=10, r=28, b=8/3):

x_dot = -p*x +p*y

y_dot = -x*z +r*x-y

z_dot = x*y -b*z

return np.array([x_dot, y_dot, z_dot])

t=0 ##### 初期時間######

dt = 0.01 ##### 微分間隔######

stepCnt = 10000 ######時間ステップ##

delta = 500 ######遅延時間＃####

# Need one more for the initial values

xs = np.empty((stepCnt + 1,))

ys = np.empty((stepCnt + 1,))

zs = np.empty((stepCnt + 1,))

# Setting initial values

xs[0], ys[0], zs[0] = (0.1, 0.5, 1.2)

k0=[0,0,0]

k1=[0,0,0]

k2=[0,0,0]

k3=[0,0,0]

########x(t)のデータ作成####################

f = open('xdata.csv', 'w')

csvWriter = csv.writer(f)

# Stepping through "time".

for i in range(stepCnt-2*delta):

x,y,z=xs[i],ys[i],zs[i]

k0 = dt * lorenz(x,y,z)

k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)

k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)

k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0

dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0

dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0

xs[i+1] = xs[i] + dx

ys[i+1] = ys[i] + dy

zs[i+1] = zs[i] + dz

listData = []

listData.append(x)

csvWriter.writerow(listData)

########x(t-δ)のデータ作成####################

f = open('xdata2.csv', 'w')

csvWriter = csv.writer(f)

# Stepping through "time".

for i in range(stepCnt-delta):

x,y,z=xs[i],ys[i],zs[i]

k0 = dt * lorenz(x,y,z)

k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)

k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)

k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0

dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0

dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0

xs[i+1] = xs[i] + dx

ys[i+1] = ys[i] + dy

zs[i+1] = zs[i] + dz

listData2 = []

if i >= 500:

listData2.append(x)

csvWriter.writerow(listData2)

########x(t-2δ)のデータ作成####################

f = open('xdata3.csv', 'w')

csvWriter = csv.writer(f)

# Stepping through "time".

for i in range(stepCnt):

x,y,z=xs[i],ys[i],zs[i]

k0 = dt * lorenz(x,y,z)

k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)

k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)

k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0

dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0

dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0

xs[i+1] = xs[i] + dx

ys[i+1] = ys[i] + dy

zs[i+1] = zs[i] + dz

listData3 = []

if i >= 1000:

listData3.append(x)

csvWriter.writerow(listData3)

###########all-time-datasiries(graph-setting)##################

fig = plt.figure()

ax = fig.gca(projection='3d')

ax.plot(xs, ys, zs,c='m')

ax.set_xlabel("X Axis")

ax.set_ylabel("Y Axis")

ax.set_zlabel("Z Axis")

ax.set_title("Lorenz Attractor")

f.close()

plt.show()

plt.savefig("Lorenz-alltimesireies.png")

f = open('xdata.csv', 'rb')

f = open('xdata2.csv', 'rb')

f = open('xdata3.csv', 'rb')

##########embedding-theory-datasireis(graph-setting)#############

x = np.loadtxt("xdata.csv", unpack=True)

y = np.loadtxt("xdata2.csv", unpack=True)

z = np.loadtxt("xdata3.csv", unpack=True)

fig = plt.figure()

ax = fig.gca(projection='3d')

ax.plot(x, y, z,c='m')

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.set_zlabel("Z")

ax.set_title("Takens embedding '$\delta=500$'")

f.close()

plt.show()

plt.savefig("Lorenz-embedding-500.png")

今回は、遅れ時間の時系列データもプログラムできました。しかし、遅れ時間のデータをfor文で別々に書いているので、ここを一つにまとめられたらよかったのです。

それと、遅れ時間も毎回、δの部分を書き換えて計算しているので、δ の配列を作り、配列の要素ごとに計算するというのも、やってみたいです。

アトラクタの軌跡のアニメーションも作ろうと思って計算したのですが、PCが低スペックでフリーズしたこともあり、今回はやめておきました。

そこらへんも、次回以降書いて行きたいと思います。

スポンサードリンク

記事が役に立ったらシェア！

Rossler方程式を Runge-Kutta法を用いて計算してみた（言語はpython）。

スポンサードリンク

多変量時系列解析の勉強その２
前回の記事で「アトラクター」という単語が出てきたので、今回はその用語を中心に書いていこうと思う。

＞＞＞前回の記事「多変量時系列解析に関する論文」。

とりあえずググると、アトラクターとは、「ある力学系がそこに向かって時間発展をする集合のことである。カオスな力学系に対してアトラクターを描写することは、現在においてもカオス理論における一つの研究課題である。」とか「実世界における力学系は散逸的であることが多いであろう。すなわち、もし力学系に運動の駆動力が無ければ、運動は停止するものと考えられる（そのような散逸は、様々な原因による内部摩擦や熱力学的損失、物質の損失などにより生じうる）。散逸と駆動力が組み合わさることにより、初期の摂動を鎮め、その力学系の振る舞いを典型的なものへと落ち着かせる傾向にある。そのような典型的な振る舞いに対応している力学系からなる位相空間の一部分はattracting section または attractee と呼ばれる。」とかウィキペディアには書いてある。

しかしこれだけではもちろん意味不なので、論文を読んでみた。

＞＞＞読んだ論文はこちら。

今回はこちらの内容を元に、勉強したことを書いていこうと思います。

埋め込みとアトラクタ表示

例えばある３変数関数f(x,y,z)があるとする。しかし得られたデータは状態変数xの時系列データ(x(t))であるとする。そしてここから遅れ時間$/delta$なる変数を新たに導入し、新たにn次元ベクトル${\bf V}= (x(t),x(t+\delta),…,x(t+(n-1)\delta))(t = 0, 1, · · · , N − 1 − (n − 1)\delta)$を生成する。

以下では、例として、連続力学系のRossler方程式を取り上げる。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= -y-z\\
\frac{dy}{dt}= x+ay\\
\frac{dz}{dt}= b+z(x-c)
\end{array}
\right.
\]
なお今回は、a = 0.398, b = 2, c = 4 とし，ルンゲクッタ法を用いて数値計算を行なった.

まず、全ての状態変数（x、y、z）を用いて軌道を描いたアトラクターが以下である。

次に、３つの状態変数のうちxのみが観測可能であると仮定し、遅れ時間$\delta$を用いて。アトラクタを再構成したものが以下である。なお遅れ時間$\delta=1,10,500$を用いて時系列データを作成した。

図２　遅れ時間[latex]\delta=1[\mathrm{s}][/latex]の再構成アトラクタ

図３　遅れ時間[latex]\delta=10[\mathrm{s}][/latex]の再構成アトラクタ

図４　遅れ時間[latex]\delta=500[\mathrm{s}][/latex]の再構成アトラクタ

図からわかる通り、図３のアトラクタの位相は全ての状態変数が観測できた時に描けるアトラクタ（図１）と位相構造が似ている。

そして今回の計算ではTakensの埋め込み定理は満たされていないが、それでも元の力学系の位相的構造が保存されていることがわかる。

本来、アトラクタの再構成を適切におこなうには、時間遅れ$\delta$は時系列データの「自己相関関数」や「パワースペクトラム」を用いるらしいが、今回のような状態変数が少なく、データを視覚的に捉えやすい場合は、必ずしもそうする必要がないことがわかる。

スポンサードリンク

Runge-Kutta法によるR$\ddot{\mathrm{o}}$ssler方程式の計算

以下に、今回の計算に用いたプログラムを載せる。

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


# rossler equationの関数作成
def lorenz(x, y, z, a=0.398, b=2, c=4):
    x_dot = -y -z
    y_dot = x + a*y
    z_dot = b + z*(x-c)
    return np.array([x_dot, y_dot, z_dot])

#初期値指定
t=0
dt = 0.05
stepCnt = 7000
xs = np.empty((stepCnt + 1,))
ys = np.empty((stepCnt + 1,))
zs = np.empty((stepCnt + 1,))

xs[0], ys[0], zs[0] = (0.1, 0.5, 1.2)
k0=[0,0,0]
k1=[0,0,0]
k2=[0,0,0]
k3=[0,0,0]

# 書き込みファイルの指定
f = open('xdata.csv', 'w') 
csvWriter = csv.writer(f)

#ルンゲ-クッタ法
for i in range(stepCnt):
    x,y,z=xs[i],ys[i],zs[i]

    k0 = dt * lorenz(x,y,z)
    k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)
    k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)
    k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

    dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0
    dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0
    dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0

    xs[i+1] = xs[i] + dx
    ys[i+1] = ys[i] + dy
    zs[i+1] = zs[i] + dz
    listData = []   #listの初期化
    listData.append(x)  
    csvWriter.writerow(listData)  #1行書き込み   

#グラフ描画設定               
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(xs, ys, zs)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor")
f.close()
plt.show()

# -*- coding: utf-8 -*-

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

import csv

# rossler equationの関数作成

def lorenz(x, y, z, a=0.398, b=2, c=4):

x_dot = -y -z

y_dot = x + a*y

z_dot = b + z*(x-c)

return np.array([x_dot, y_dot, z_dot])

#初期値指定

t=0

dt = 0.05

stepCnt = 7000

xs = np.empty((stepCnt + 1,))

ys = np.empty((stepCnt + 1,))

zs = np.empty((stepCnt + 1,))

xs[0], ys[0], zs[0] = (0.1, 0.5, 1.2)

k0=[0,0,0]

k1=[0,0,0]

k2=[0,0,0]

k3=[0,0,0]

# 書き込みファイルの指定

f = open('xdata.csv', 'w')

csvWriter = csv.writer(f)

#ルンゲ-クッタ法

for i in range(stepCnt):

x,y,z=xs[i],ys[i],zs[i]

k0 = dt * lorenz(x,y,z)

k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)

k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)

k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0

dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0

dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0

xs[i+1] = xs[i] + dx

ys[i+1] = ys[i] + dy

zs[i+1] = zs[i] + dz

listData = [] #listの初期化

listData.append(x)

csvWriter.writerow(listData) #1行書き込み

#グラフ描画設定

fig = plt.figure()

ax = fig.gca(projection='3d')

ax.plot(xs, ys, zs)

ax.set_xlabel("X Axis")

ax.set_ylabel("Y Axis")

ax.set_zlabel("Z Axis")

ax.set_title("Lorenz Attractor")

f.close()

plt.show()

なお、上のプログラムは全ての状態変数が観測できた時の、アトラクタを描画する際のプログラムである。

本来なら、このプログラムに時間遅れを加味した時系列データをいい感じに元のデータにappendでもしたり、二次元配列のデータの任意の一部のデータを除去したりできればよかったです。ですがそこまでの実力がなかったので（一応数時間奮闘はしました笑）、時間遅れの再構成のアトラクタは別にプログラムを書いた（とは言ってもデータ処理の部分のみ）。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import csv
import pandas as pd

f = open('xdata.csv', 'rb')
f = open('xdata-2.csv', 'rb')
f = open('xdata-3.csv', 'rb')

fig = plt.figure()
ax = fig.gca(projection='3d')
x = np.loadtxt("xdata.csv", unpack=True)
y = np.loadtxt("xdata-2.csv", unpack=True)
z = np.loadtxt("xdata-3.csv", unpack=True)
ax.plot(x, y, z)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Takens embedding '$\delta=1$'")
f.close()
plt.show()

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

import csv

import pandas as pd

f = open('xdata.csv', 'rb')

f = open('xdata-2.csv', 'rb')

f = open('xdata-3.csv', 'rb')

fig = plt.figure()

ax = fig.gca(projection='3d')

x = np.loadtxt("xdata.csv", unpack=True)

y = np.loadtxt("xdata-2.csv", unpack=True)

z = np.loadtxt("xdata-3.csv", unpack=True)

ax.plot(x, y, z)

ax.set_xlabel("X")

ax.set_ylabel("Y")

ax.set_zlabel("Z")

ax.set_title("Takens embedding '$\delta=1$'")

f.close()

plt.show()

xdata.csvは観測できる時系列データです。そして、xdata-2とxdata-3は時間遅れの新規時系列データx(t+$\delta$)とx(t+2$\delta$)ということになります。

これらのデータはとても原始的に用意しました笑。xdataの時系列データの時間遅れの分だけ上の行から取り除いたデータを、xdata2にコピペし、さらに、そのxdata２の時系列データから、時間遅れの分だけ上の行から取り除いたデータを、xdata3にコピペする。

最後に時間遅れの分、xdata2,xdata3は、要素数（ファイルの行数）が少ないので、一番データ行数の少ない、xdata3と行数が合うようにデータを消して、３つのファイルの行数を揃えます。

こんな原始的にやらなくても、pythonならこれくらいのこと瞬殺できそうですが、ここはまた今度にします（pythonの無駄遣い笑）。

最後に

今回は「多変量時系列解析」に関してしっかりと理解するため、論文のような計算結果が得られるか確認しました。ルンゲクッタとかはググったら瞬殺できたのでよかったのですが、pythonに慣れていなく、データ処理の所に時間がかかってしまいました。

今後はpythonのデータ処理と、ローレンツ方程式ではどんな計算結果が得られるのか勉強したいと思います。

記事が役に立ったらシェア！

多変量時系列解析に関する論文

最近、人工知能とかデータサイエンスとか、Pythonに興味があるので、勉強を開始しました。とりまテキトーに論文を漁る所から。

辿り着いた論文は”Analyzing multiple nonlinear time series with extended Granger causality”　PDFはこちら。

スポンサードリンク

論文の内容の前に

その前に日本語で書かれても全く理解出来ないです。何せ専門用語しらないのですから。

今回は私がこの論文を読むにあたり、調べた用語などのアウトプットとアブストラクトの要約を述べようと思います。

多変量解析

多数の変量を同時に解析すること。多変量解析があるならば、一変量解析、二変量解析も当然あります。

データ解析は数値計算をするだけでなく、データの可視化をすることで、データに対する

理解をする手助けをすることも含めます。

例えば、一変量解析ならば、ヒストグラムなど。二変量解析ならば、散布図、三変量ならば、3次元の散布図などです。

人間は4次元以上のものは、グラフのようなデータの可視化ができないので、多変量解析はまずは、次元を落としてやるところから始めます。または、方程式を可視化する方法もあるんだとか。
グレンジャー因果

ある変数に他の変数を説明変数として加えたときに、平均二乗誤差に変化が出る場合、グレンジャー因果があるという。

この論文ではこの手法を用いています。

このグレンジャ因果は線形確率系に用いることで、多変量の関係を見つける手法らしいです。
ターケンスの埋め込み定理

１つの状態変数から多数の状態変数を復元する際に用いられる。この定理によると、復元したデータの位相空間は復元する前の位相空間（アトラクター）と微分的に同相となります。
相空間

一般的な力学系では、運動を記述する際、位置と時間だけのグラフでは、完全に特徴を記述できず、N次元の空間を想定する。この空間のことを相空間あるいは、状態空間という。

アブストラクトの内容

多変量解析において、多変数の相関を特定することは重要である。そりゃ、そうだ！運動量が位置（例えばx）によって変化するのに、xの関数で表そうとしないなんて馬鹿な話はありませんね。この論文では、線形確率系において、使われたグレンジャー因果を非線形系に拡張し、その結果を述べるそうです。

特に３変量やそれ以上の時系列データに対して、変数の相関が見られるかもしくは、他の工程により分解されるかどうかを決定することが可能となる、方法を提案するらしい。

なるほど、非線形系に対してグレンジャー因果を応用したときに、変数の相関を確認することができるかどうかを明らかにしたいそうです。

まずは、グレンジャー因果と多変量解析の勉強をしていくところですかね！

スポンサードリンク

記事が役に立ったらシェア！

Category: 論文関連

Euler法とRunge-Kutta法のカオス時系列データの比較

モチベーション

Eular法とは

Runge-Kutta法とは

計算結果比較

今回用いたプログラム

ローレンツ方程式

一次元減衰振動モデル

最後に

ローレンツ方程式をRunge-Kutta法（python）で計算

時間遅れ座標系を用いた、アトラクタの再構成

ターケンスの埋め込み定理

時間遅れの大きさの最適化

計算に用いたプログラム

Rossler方程式を Runge-Kutta法を用いて計算してみた（言語はpython）。

多変量時系列解析の勉強その２
前回の記事で「アトラクター」という単語が出てきたので、今回はその用語を中心に書いていこうと思う。

埋め込みとアトラクタ表示

Runge-Kutta法によるR\(\ddot{\mathrm{o}}\)ssler方程式の計算

最後に

多変量時系列解析に関する論文

論文の内容の前に

アブストラクトの内容

モチベーション

Eular法とは

Runge-Kutta法とは

計算結果比較

今回用いたプログラム

ローレンツ方程式

一次元減衰振動モデル

最後に

時間遅れ座標系を用いた、アトラクタの再構成

ターケンスの埋め込み定理

時間遅れの大きさの最適化

計算に用いたプログラム

多変量時系列解析の勉強その２ 前回の記事で「アトラクター」という単語が出てきたので、今回はその用語を中心に書いていこうと思う。

埋め込みとアトラクタ表示

Runge-Kutta法によるR\(\ddot{\mathrm{o}}\)ssler方程式の計算

最後に

論文の内容の前に

アブストラクトの内容

多変量時系列解析の勉強その２
前回の記事で「アトラクター」という単語が出てきたので、今回はその用語を中心に書いていこうと思う。