tobiom’s diary

備忘録

ノートに万年筆で書き込みたい

高校生のころからSEIYUのクロチェットRとかいう名前のノートを良く使ってたけど、この前SEIYUに行ったら売ってなかった。

SEIYUで取り扱わなくなってしまったのか、ちょっと離れた店舗も探してみたがどこも売ってない。

このクロチェットR、万年筆で書き込んでも裏写りしない上に見た目がシンプルで好みだっただけに残念だった。

更にB5の30枚で35円/冊なので気兼ねなく多用出来たのが強かったのに...

仕方なく他のノートをいくつか試してみたところ、DAISOのCompleteとかいう3冊組のノートが一番クロチェットRに近い品質だったので今後はこれを使っていく事になりそう。

f:id:tobiom:20170404222512j:plain

 

実際に万年筆で試し書き↓

f:id:tobiom:20170404222200j:plain

その裏↓

f:id:tobiom:20170404222320j:plain

若干透けて見えるけど大して気にならない程度(DAISOのカートリッジ式はフローが多いからちょっと透けやすい)。

 

2重振り子の数値シミュレーション

2重振り子はカオス性を示す物理系として最も有名な例の1つ(だと思っている)。

この2重振り子をシミュレーションしてみた。

解析方針

この2重振り子、微分方程式の数値解法として主流なRunge-Kutta法を適用するとエネルギーが発散してシミュレーションが破綻し易いという事を訊き、FPUのシミュレーションを作る時に利用したSymplectic法を用いる事を考えた。

だが、実装した事のあるSymplectic法はHamiltonianが座標xだけの関数V(x)と運動量pだけの関数T(p)に分離できる前提の陽解法であった。

そこでSymplectic性をもつ陰解法であるGauss-Legendre法を採用。

これによってエネルギーを比較的高精度で保存するシミュレーションを実装する事が可能となった。

しかし、シミュレーションする系がカオスである以上、原理的に現実の系をほぼ寸分のずれも無く再現する事は不可能であるため、あくまでも初期値鋭敏性を分かり易く視覚化するのみのシミュレータになる。

Hamiltonianの導出について

以下、Hamiltonianは運動エネルギー項とポテンシャル項の和となるとする。

2重振り子の系のHamiltonianは形が複雑なため、一旦Lagrangianを導出してからLegendre変換でHamiltonianを求めることにした。

Lagrangian自体は比較的簡単に求められて

L=T-V
T=\displaystyle\frac{1}{2}(m_1+m_2)l_1^2\dot{\theta_1}^2+\displaystyle\frac{1}{2}m_2l_2^2\dot{\theta_2}^2+\displaystyle\frac{1}{2}m_2l_1l_2\dot{\theta_1}\dot{\theta_2}\cos{\left(\theta_1-\theta_2\right)}
V=-(m_1+m_2)gl_1\cos{\theta_1}-m_2gl_2\cos{\theta_2}

である。

ここからLegendre変換p=\displaystyle\frac{\partial L}{\partial \dot{q}}を用いてp_i=\displaystyle\frac{\partial L}{\partial \dot{\theta_i}}を計算する。
このときL=T+Vにおいて\dot{\theta_i}Tの中にのみ現れるのでT偏微分を考えればよろしい。

p_1=(m_1+m_2)l_1^2\dot{\theta_1}+\displaystyle\frac{1}{2}m_2l_1l_2\dot{\theta_2}\cos{\left(\theta_1-\theta_2\right)}

p_2=m_2l_2^2\dot{\theta_2}+\displaystyle\frac{1}{2}m_2l_1l_2\dot{\theta_1}\cos{\left(\theta_1-\theta_2\right)}

この式を用いてL\theta_ip_iの関数に書き換えるとHが得られる。

直接p_i\dot{\theta_i}について解いて代入する方法だと面倒なので、ある程度計算をやり易くするために次の対称行列Aを使う。

A=\begin{pmatrix} (m_1+m_2)l_1^2 & \displaystyle\frac{1}{2}m_2l_1l_2\cos{\left(\theta_1-\theta_2\right)} \\ \displaystyle\frac{1}{2}m_2l_1l_2\cos{\left(\theta_1-\theta_2\right)} & m_2l_2^2 \end{pmatrix}

するとこの行列A\theta_iのベクトル\vec{\bf\theta}=\begin{pmatrix}\theta_1\\\theta_2\end{pmatrix}を用いて\vec{\bf p}=\begin{pmatrix}p_1\\p_2\end{pmatrix}
\vec{\bf p}=A\vec{\bf \dot{\theta}}、即ち\vec{\bf \dot{\theta}}=A^{-1}\vec{\bf p}と表せる。

ここで\theta_iの関数TT=\frac{1}{2}\dot{\vec{\bf\theta}}^tA\vec{\bf\dot{\theta}}であったから、ここに\vec{\bf \dot{\theta}}=A^{-1}\vec{\bf p}を代入すると

T=\frac{1}{2}\left(A^{-1}\vec{\bf p}\right)^tA\left(A^{-1}\vec{\bf p}\right)=\frac{1}{2}\vec{\bf p}^t A^{-1} \vec{\bf p}

こうして得られたp_iの関数としてのTにポテンシャル項Vを足すとHamiltonianが得られる。

運動方程式

以上の結果を正準方程式に適用すると、以下のようになる。

\dot{\vec{\bf p}}=-\displaystyle\frac{\partial (T+V)}{\partial\vec{\bf\theta}}=-\frac{1}{2}\vec{\bf p}^t \displaystyle\frac{\partial A^{-1}}{\partial\vec{\bf\theta}} \vec{\bf p} - \begin{pmatrix}(m_1+m_2)gl_1\sin{\theta_1}\\m_2gl_2\sin{\theta_2}\end{pmatrix}

\dot{\vec{\bf\theta}}=\displaystyle\frac{\partial T}{\partial\vec{\bf p}}=A^{-1} \vec{\bf p}

動画


2重振り子
何故かカクついてる

gnuplotを使ったCoulomb散乱の図示化

gnuplotでCoulomb散乱のモデルを計算してコマ毎に描画したのをWindowsムービーメーカーで繋げて動画にしてみました。
計算には2次のSymplectic法を用いています。

コードを書いて

Coulomb_Scattering.plt
#Coulomb_Scattering

set terminal jpeg enhanced font "Times" 20 size 600, 480
set tics font 'Times,12'
set nokey
set ticslevel 0

#画面
wx=10.0
wy=8.0
set xr[-wx*0.5:wx*0.5]
set yr[-wy*0.5:wy*0.5]

ASIZE=2*15

#配列にデータを格納する関数を定義するための文字列の生成
str_qx = "qx_set(i,x) = "
str_qy = "qy_set(i,x) = "
str_px = "px_set(i,x) = "
str_py = "py_set(i,x) = "
do for [k=0:ASIZE-1] {str_qx = str_qx . sprintf("(int(i)==%d ? _qx%d=x :", k, k)}
do for [k=0:ASIZE-1] {str_qy = str_qy . sprintf("(int(i)==%d ? _qy%d=x :", k, k)}
do for [k=0:ASIZE-1] {str_px = str_px . sprintf("(int(i)==%d ? _px%d=x :", k, k)}
do for [k=0:ASIZE-1] {str_py = str_py . sprintf("(int(i)==%d ? _py%d=x :", k, k)}
str_qx = str_qx. "1/0"
str_qy = str_qy. "1/0"
str_px = str_px. "1/0"
str_py = str_py. "1/0"
do for [k=0:ASIZE-1] {str_qx = str_qx . ")"}
do for [k=0:ASIZE-1] {str_qy = str_qy . ")"}
do for [k=0:ASIZE-1] {str_px = str_px . ")"}
do for [k=0:ASIZE-1] {str_py = str_py . ")"}
#文字列実行
eval str_qx
eval str_qy
eval str_px
eval str_py

#配列のデータを参照する関数を定義するための文字列の生成
str_qx = "qx_get(i) = "
str_qy = "qy_get(i) = "
str_px = "px_get(i) = "
str_py = "py_get(i) = "
do for [k=0:ASIZE-1] {str_qx = str_qx . sprintf("(int(i)==%d ? _qx%d :", k, k)}
do for [k=0:ASIZE-1] {str_qy = str_qy . sprintf("(int(i)==%d ? _qy%d :", k, k)}
do for [k=0:ASIZE-1] {str_px = str_px . sprintf("(int(i)==%d ? _px%d :", k, k)}
do for [k=0:ASIZE-1] {str_py = str_py . sprintf("(int(i)==%d ? _py%d :", k, k)}
str_qx = str_qx. "1/0"
str_qy = str_qy. "1/0"
str_px = str_px. "1/0"
str_py = str_py. "1/0"
do for [k=0:ASIZE-1] {str_qx = str_qx . ")"}
do for [k=0:ASIZE-1] {str_qy = str_qy . ")"}
do for [k=0:ASIZE-1] {str_px = str_px . ")"}
do for [k=0:ASIZE-1] {str_py = str_py . ")"}
#文字列実行
eval str_qx
eval str_qy
eval str_px
eval str_py

DATA = "Coulomb_Scattering.data"
set print DATA

#初期値・係数
mass = 0.5
alpha = 1.5
R=wx*0.4
v0 = 3.0
do for [k=0:ASIZE-1] {dummy=qx_set(k,-wx*0.5)}
do for [k=0:ASIZE-1] {dummy=qy_set(k,(k-(ASIZE-1)*0.5)*0.2)}
do for [k=0:ASIZE-1] {dummy=px_set(k,mass*v0)}
do for [k=0:ASIZE-1] {dummy=py_set(k,0)}

#微分方程式
partial_p(p) = (p)/(mass)
partial_qx(qx,qy) = -alpha*(qx-R)/(((qx-R)**2+qy**2)**(3/2))
partial_qy(qx,qy) = -alpha*qy/(((qx-R)**2+qy**2)**(3/2))
dt = 0.001    #時間刻み
n_max = 500  #静止画数
m_max = 20   #間引き数
n = 0;
m = 0
call "si2.plt"
si2.plt
do for [k=0:ASIZE-1] {
	cpx=px_get(k)
	cpy=py_get(k)
	cqx=qx_get(k)
	cqy=qy_get(k)
	dummy=px_set(k,cpx-dt*0.5*partial_qx(cqx,cqy))
	dummy=py_set(k,cpy-dt*0.5*partial_qy(cqx,cqy))
}
do for [k=0:ASIZE-1] {
	cpx=px_get(k)
	cpy=py_get(k)
	cqx=qx_get(k)
	cqy=qy_get(k)
	dummy=qx_set(k,cqx+dt*partial_p(cpx))
	dummy=qy_set(k,cqy+dt*partial_p(cpy))
}
do for [k=0:ASIZE-1] {
	cpx=px_get(k)
	cpy=py_get(k)
	cqx=qx_get(k)
	cqy=qy_get(k)
	dummy=px_set(k,cpx-dt*0.5*partial_qx(cqx,cqy))
	dummy=py_set(k,cpy-dt*0.5*partial_qy(cqx,cqy))
}

do for [k=0:ASIZE-1] {print qx_get(k),qy_get(k)}
m=m+1
if(m<m_max) reread
m=0

outfile(n) = sprintf("f%d.jpg",n+10000) #出力ファイル名
title(n) = sprintf("t = %4.2f",n*dt) #タイトル名
unset label
set label title(n) font 'Times,12' at -wx/2, wy/2
set output outfile(n)
plot DATA with points pt 1 ps 0.1

n = n+1
if(n < n_max) reread
n=0

実行して画像を生成

f:id:tobiom:20161216234524p:plain
f:id:tobiom:20161216234445j:plain

参考文献

gnuplotでの配列(もどき)の使用 - 米澤進吾 ホームページ
gnuplotには配列が標準装備で無かったのでこのサイトのコードをコピペして改悪参考にさせて頂きました。