気体分子運動論（Joule, Kronigの理論）

趣味で統計力学を勉強し始めました。アインシュタインのブラウン運動についての記事を読んでいるうちにもっと知りたいと思ったので。

今回は統計力学に入る前の気体分子運動論の取り扱いについて基本的な考えを知るため「Jouke,Kronigの理論」について書こうと思います。

参考文献は以下の記事で用いているものと同じです。

断面積S、高さhの容器（簡単のため立方体）を考えます。その中にN個の分子が入っているとします。そしてこれらの分子同士に相互作用はないものとします。

2. 分子はばらばらに動いていますが、平均をとるとばらつきはなく

N/3個はx軸方向に
N/3個はy軸方向に
N/3個はz軸方向に

動いているものとします。

3. 速さは皆$v = \sqrt{v_{x}^{2}+v_{y}^{2}+v_{z}^{2}}$であるとします。

4. 鉛直方向に速さvで移動している分子の数はこれまでの仮定からN/3個です。そのN/3個の分子は単位時間あたりに$\frac{v}{2h}$回、上下の壁に衝突します。

5. 単位時間あたりの一つの分子の運動量変化は$2mv \times \frac{v}{2h} $です

6. 壁が受ける力は単位時間あたりのN/3個の分子の運動量変化に等しいので、以下が成り立つはずです。なおpは圧力, Vは体積です。
\begin{eqnarray}
pS &=& \frac{N}{3} 2mv \times \frac{v}{2h} \\
\leftrightarrow p\frac{V}{h} &=& \frac{Nmv^2}{3h}\\
\leftrightarrow pV &=& \frac{Nmv}{3}
\end{eqnarray}

今回は一番簡易なモデルなため雑な計算方法ですが、異なる速度成分、異なる質量、異なる分子数として計算したステファンボルツマンのモデルを用いても厳密には意味が異なりますが
\begin{eqnarray}
pV =\frac{Nmv}{3}
\end{eqnarray}

の関係式は得られます。そしてステファンボルツマンの理論を用いることで「ボルツマン定数」を定義できたり、熱力学では得ることのできなかった関係式を得ることができます。

次回は「Jouke,Kronigの理論」を広げた「ステファンボルツマンの理論」について述べていきたいと思います。

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

「統計力学」勉強しようとしたけど、マクスウェルの速度分布が難しすぎる件ww

先に言っておきますが、マクスウェルの速度分布について綺麗にまとめているドキュメントではありません。

その過程というかあまりにもはやすぎる段階でつまづいたので、ネタとして記録しておこうと思ったまでです。誰得な内容ですが、寝る前や研究の休憩途中で「こんなこともわからない奴が東工大の大学院生か！」とでも私を罵ることで、読者様に優越感を与えられるなら幸いですよ。

参考図書は以下。すごく細かく書いてあるのでわかりやすいです。物理学科の学部生ならヨユーで理解できると思います。（自分は機械系なのでお許しをww）。

ついでにこの本は東北大で宇宙物理学(MHD)の研究をしているドクター（吾輩の兄）から譲りウケたもの。

　　>>>熱力学・統計力学
 　　

マクスウェルの速度分布の粗い導出

本書には

「分布が適当になると、分子同士の衝突にも関わらず分布が変わらなくなるであろうが、それはどのような分布であろうか。これは分子運動論で大切な問題である（p185から抜粋）」

この速度分布を雑に導くためにマクスウェルが最初においた仮定が

「分子の3方向の速度成分の分布は、違う方向について互いに独立である」というもの。

xとyとz方向の速度分布はそれぞれ独立だと。なんだか偏微分とかしなそうで楽そう！

で厳密にはこの仮定さえも証明が必要らしいのだが、これを証明なしで仮定することで簡単に速度分布が求まるみたい。

だがしかし…

N個の分子のうち、速度成分が

\begin{eqnarray}
(v_x, v_y , v_z)　\mathrm{と}　(v_x + dv_x, v_y + dv_y, v_z + dv_z)
\end{eqnarray}

の間にあるものの数を

\begin{eqnarray}
F(v_x, v_y , v_z)　dv_x dv_y dv_z
\end{eqnarray}

とする。

この仮定はぎり理解できました。

速度(v_x,v_y,v_z)時点での個数を$F(v_x, v_y , v_z)$とおけば速度空間で考えた場合、以下の図のように体積をかけることで、$(v_x, v_y , v_z)　\mathrm{と}　(v_x + dv_x, v_y + dv_y, v_z + dv_z)$間の数が計算できるのは想像できます。

問題はこの後の文章。

「この現象についてどの方向も全く同等であるから、$F(v_x, v_y , v_z)$は速度ベクトル$\vec{v}$の大きさだけ、したがってv²によるはずである。これを$F(v^2)$と書こう。(p185から引用)」

どの方向も同等→それは仮定からわかる
$F(v_x, v_y , v_z)$は速度ベクトル$\vec{v}$はv²による→????

なんで方向依存性がないと速度分布が大きさに依存するのー？

マジでこんなことつまづく奴おるんかいなって感じです。

さらに

「$v_xとv_x + dv_x$の間に数は$v_x$だけの関数$f(v_x)$に比例していて…」

またここにきて新たな関数$f(v_x)$。もうなんかね。さすが「物理」だよって感じですね。

これを職業としてる人は頭どうかしてるんじゃww

最近東大が熱力学第二法則を量子力学使って証明したみたいな感じの記事がありましたね。

東大の同い年の理論物理の研究している大学院生と話してみたいなww

頭の回転早すぎて、会話にならなさそうww

趣味で「物理学」はこれだから楽しいですね。アウトプットが求められてるわけではないからわからんくてもストレスかからないですし。

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

地球のトンネル -反対側まで何分?-

スポンサードリンク

地球の重力

皆さん常日頃から感じている重力。この重力のおかげで、我々は大地の上に立つことができています。

また飛行機からスカイダイビングしたり、ジェットコースターで落ちるスリリングを楽しんでいるわけですね。

ではこの重力の向き、我々日本の反対側にいる人たちはどのように感じているのでしょうか。

答えは簡単、というよりも常識とでもいいましょうか。

地球の反対側にいる人だって普通に立っているのです。答えはやはり地面に向かって重力を感じている、です。

要するに、我々地球の上に住む生命は皆、地球の中心に向かって重力で引っ張られているのです。

スポンサードリンク

地球の中心を貫くトンネルを作ろう!

図１を見ればわかる通り、地球の反対側では重力の向きは反対ですね。

ということは、地球の中心までは重力で加速されてトンネルを落ちていき、中心を通り過ぎた後は重力で減速しながらトンネルを登っていく、そんなトンネルが作れるはずです。

作れるかどうかという現実問題はさておき、この夢のようなトンネルができれば、なんのエンジンも必要なくおよそ45分で反対側に到着できます。

必要なのは母なる大地、地球の重力エネルギーだけです。

図の上半分では、地球の重力 (青矢印) で加速され地球の中心へ向かう。中心を通り過ぎた後は重力で減速されながら反対側の地表に向かって上昇していく。

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

pythonで計算（一次元振動モデル、ロトカヴォルテラ方程式）

スポンサードリンク

今回は、pythonの勉強の一環でやってきた簡単な数値計算のソースコードを紹介していきます。

一次元減衰強制振動

# -*- 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.1
u = 5
###################################################
tr=[]
xr=[]
te=[]
xe=[]
def g(t, x):
    return np.array([x[1], (-k/m*x[0])-2*w*x[1]+u])

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


######################ルンゲクッタ#########################
for i in range(800):
    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(800):
    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() 
ax.plot(tr, xr, c='b', label="Runge-Kutta")
ax.plot(te, xe, c='r', label="Eular")
ax.set_xlabel("time")
ax.set_ylabel("x(t)")
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.1

u = 5

###################################################

tr=[]

xr=[]

te=[]

xe=[]

def g(t, x):

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

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

dt = 0.1

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

for i in range(800):

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

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

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

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

ax.set_xlabel("time")

ax.set_ylabel("x(t)")

ax.set_title("x-t graph")

plt.legend()

plt.show()

x-t(振幅ー時間)グラフは以下のようになります。

両解法（オイラー法とルンゲクッタ法）において、以前にも記事に書いた通り、一次元で方程式が単純ということもあり、そんなに解に差は見られません。

＞＞＞以前の記事はこちら

python初心者の方は、設定パラメータを動かして、どんなグラフになるか、遊んでみてください。

以下に、一次元減衰強制振動のモデルを示します。

ばね定数：k、減衰係数：c、質量：m、強制外力：m/u

スポンサードリンク

ロトカヴォルテラ方程式

自然界における捕食者と被食者の個体数変動の数理モデルです。

以下の数式で定義されています。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= ax-bxy\\
\frac{dy}{dt}= cxy-dy
\end{array}
\right.
\]

xが被食者の個体数、yは捕食者の個体数、tは時間、4つの係数a、b、c、dは正の実数パラメータです。

以下は、a=8, b=3, c=4 ,d=18の時の計算結果。

周期関数になることがわかりますが、直感的に何が起きているか不明なので、定性的に述べていきます。ライオンとシマウマの例で考えると、

捕食者であるライオンが増加
それにより被食者であるシマウマが減少
エサ（シマウマ）が減少するので、ライオンは減少
天敵（ライオン）が減るので、シマウマは増加
１以降を繰り返す

といったことを繰り返しています。よくできたモデルですが、簡易化されているため、実際はもっと複雑です。以下は上で述べたことを図にしたものです。赤が捕食者（ライオン）で、青が被食者（シマウマ）です。

ルンゲクッタ法とオイラー法の比較

オイラー法とルンゲクッタ法の誤差の確認として、被食者数（上で言えばシマウマ）をグラフにしてみました。

ぱっと見では、解に差が内容に思われますが、縮尺を変更すると

時刻t=122では、0.05くらいの差があります。

割合でいうと、0.05/12=0.004（0.4%）くらいの差があります。

まあルンゲクッタが正解かというとそうでもないですが。

ついでに、グラフが直線的になっているのは、計算する際の微分間隔が雑なためです。

以下にコードを載せますが、そこの「微分間隔」と書いてあるところを小さくしていけば、理論的な微分に近づき、滑らかな解（グラフ）が得られます。

計算に用いたプログラム

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


# rorenz equation
def lotka(x, y, a=8, b=3, c=4 ,d=18):
    x_dot = a*x -b*x*y
    y_dot = c*x*y -d*y
    return np.array([x_dot, y_dot])


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


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

# Setting initial values
xs[0], ys[0] = (10, 7)
k0=[0,0,0]
k1=[0,0,0]
k2=[0,0,0]
k3=[0,0,0]
listData=[0] 
listData1=[0] 
########ルンゲクッタ####################
# Stepping through "time".


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

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

    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
    xs[i+1] = xs[i] + dx
    ys[i+1] = ys[i] + dy
    
      
    listData.append(i)  
  
       
# Need one more for the initial values
xs2 = np.empty((stepCnt + 1,))
ys2 = np.empty((stepCnt + 1,))

# Setting initial values
xs2[0], ys2[0] = (10, 7)
k02=[0,0,0]
k12=[0,0,0]
listData1=[0] 

########オイラー法####################
# Stepping through "time".


for i in range(stepCnt):
    x,y=xs2[i],ys2[i]

    k02 = dt * lotka(x,y)
    k12 = dt * lotka(x+k02[0], y+k12[1])
    dx = (k02[0]+k12[0])/2.0
    dy = (k02[1]+k12[1])/2.0
    xs2[i+1] = xs2[i] + dx
    ys2[i+1] = ys2[i] + dy
    listData1.append(i)  
                    
###########all-time-datasiries(graph-setting)##################     
fig = plt.figure()
ax = plt.gca() 
ax.plot(listData, xs, c='b', label="Runge-Kutta")
ax.plot(listData1, xs2, c='r', label="Eular")
ax.set_xlabel("time")
ax.set_ylabel("Populazione")
ax.set_title("Lotka-Volterra-model")
plt.legend() 
plt.show()

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

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

# rorenz equation

def lotka(x, y, a=8, b=3, c=4 ,d=18):

x_dot = a*x -b*x*y

y_dot = c*x*y -d*y

return np.array([x_dot, y_dot])

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

dt = 0.001 ##### 微分間隔######

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

# Need one more for the initial values

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

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

# Setting initial values

xs[0], ys[0] = (10, 7)

k0=[0,0,0]

k1=[0,0,0]

k2=[0,0,0]

k3=[0,0,0]

listData=[0]

listData1=[0]

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

# Stepping through "time".

for i in range(stepCnt):

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

k0 = dt * lotka(x,y)

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

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

k3 = dt * lotka(x+k2[0], y+k2[1])

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

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

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

listData.append(i)

# Need one more for the initial values

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

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

# Setting initial values

xs2[0], ys2[0] = (10, 7)

k02=[0,0,0]

k12=[0,0,0]

listData1=[0]

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

# Stepping through "time".

for i in range(stepCnt):

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

k02 = dt * lotka(x,y)

k12 = dt * lotka(x+k02[0], y+k12[1])

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

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

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

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

listData1.append(i)

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

fig = plt.figure()

ax = plt.gca()

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

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

ax.set_xlabel("time")

ax.set_ylabel("Populazione")

ax.set_title("Lotka-Volterra-model")

plt.legend()

plt.show()

以上です。基礎的な微分方程式や運動方程式の数値解法の勉強する際には役立てみてください。

スポンサードリンク

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

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のデータ処理と、ローレンツ方程式ではどんな計算結果が得られるのか勉強したいと思います。

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

Tag: 物理

気体分子運動論（Joule, Kronigの理論）

「統計力学」勉強しようとしたけど、マクスウェルの速度分布が難しすぎる件ww

マクスウェルの速度分布の粗い導出

だがしかし…

地球のトンネル -反対側まで何分?-

地球の重力

地球の中心を貫くトンネルを作ろう!

pythonで計算（一次元振動モデル、ロトカヴォルテラ方程式）

一次元減衰強制振動

ロトカヴォルテラ方程式

ルンゲクッタ法とオイラー法の比較

計算に用いたプログラム

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

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

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

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

最後に

マクスウェルの速度分布の粗い導出

だがしかし…

地球の重力

地球の中心を貫くトンネルを作ろう!

一次元減衰強制振動

ロトカヴォルテラ方程式

ルンゲクッタ法とオイラー法の比較

計算に用いたプログラム

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

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

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

最後に

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