HOME > natural science Laboratory > コンピュータ・シミュレーション講座 > 計算物理学

計算物理学 第1章序章
1-2.Duffing振動子の位相空間上でのカオス的運動のアニメーション

文責:遠藤 理平 (2011年4月 1日) カテゴリ:計算物理学(126)

本稿は、凝縮系物理を中心として広範囲にわたる計算手法に関して、そのアルゴリズムの導出から実装までを丁寧にまとめられていることで知られている計算物理学 J.M.ティッセン (著), (訳) 松田 和典, 道廣 嘉隆, 谷村 吉隆, 高須 昌子, 吉江 友照を元に、 本書で取り扱われているテーマについて、具体的にC++言語にてプログラミングし、計算結果を gnuplot や OpenGL を用いて図や動画にすることを目的とします。

1-1.Duffing振動子のカオス的運動のシミュレーション では、横軸:時間、縦軸:位置として、プロットしました。 一定空間を運動するようなカオス的運動の特徴を調べる場合、 位相空間での軌道を調べるほうが顕著であることがあります。 位相空間とは、横軸:位置、縦軸:運動量で張られる空間をいいます。

本稿では、 gnuplot のみを使用してDuffing振子の位相空間上でのカオス的運動のアニメーションを作成します。

計算結果

初期値、パラメータは、1-1.Duffing振動子のカオス的運動のシミュレーション の一つ目と同じです。

位相空間上での軌道もカオス的運動であることが見られます。 はじめは軌道が重なっていますが、時間が経つにつれ軌道は離れていきます。

gnuplot プログラム

gnuplot のテンプレートを利用して、 1.gnuplot 上でルンゲクッタ法(4次)を用いて数値計算し、
2.各時刻ごとに結果をファイルに出力し、
3.出力したファイルを読み込み、各時刻ごとに2次元プロットし、
4.各時刻ごとにファイル名が連番になるように jpgファイルを出力します。
そして、連番番号のjpgファイルを結合して、アニメーションを作成します。

参考

【gnuplot でアニメーション】gnuplotでルンゲクッタ法
【gnuplot でアニメーション】サイクロイド曲線

gnuplot のテンプレートファイルの作り方

1.次の2つのファイルをgnuplot の実行ファイル(wgnuplot.exe)と同じフォルダに作ります。

ファイル名:Duffing.plt

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

set xr[-2:2]
set yr[-3:4]

DATA = "Duffing.data"
set print DATA

fx(x, vx)  = vx
fvx(x, vx) = -gamma/mass*vx + 2.0*a/mass*x - 4.0*b/mass*x**3 + F0/mass*cos(omega*t)

#パラメータ
F0    = 2.0
omega = 2.4
gamma = 0.1
mass = 1.0
a = 1.0/5.0
b = 1.0/2.0
x1 = 0.5
vx1 = 0
x2 = 0.5001
vx2 = 0
t=0

dt = 0.01    #時間刻み
n_max = 100   #静止画数
m_max = 100   #間引き数
n=0;
m=0
call "rk.plt"

ファイル名:rk.plt

 
k1x  = dt* fx(x1, vx1)
k1vx = dt*fvx(x1, vx1)
k2x  = dt* fx(x1+k1x/2,vx1+k1vx/2)
k2vx = dt*fvx(x1+k1x/2,vx1+k1vx/2)
k3x  = dt* fx(x1+k2x/2,vx1+k2vx/2)
k3vx = dt*fvx(x1+k2x/2,vx1+k2vx/2)
k4x  = dt* fx(x1+k3x,  vx1+k3vx)
k4vx = dt*fvx(x1+k3x,  vx1+k3vx)
x1  =x1  +(k1x  + 2*k2x  + 2*k3x  + k4x )/6
vx1 =vx1 +(k1vx + 2*k2vx + 2*k3vx + k4vx)/6

k1x  = dt* fx(x2, vx2)
k1vx = dt*fvx(x2, vx2)
k2x  = dt* fx(x2+k1x/2,vx2+k1vx/2)
k2vx = dt*fvx(x2+k1x/2,vx2+k1vx/2)
k3x  = dt* fx(x2+k2x/2,vx2+k2vx/2)
k3vx = dt*fvx(x2+k2x/2,vx2+k2vx/2)
k4x  = dt* fx(x2+k3x,  vx2+k3vx)
k4vx = dt*fvx(x2+k3x,  vx2+k3vx)
x2  =x2  +(k1x  + 2*k2x  + 2*k3x  + k4x )/6
vx2 =vx2 +(k1vx + 2*k2vx + 2*k3vx + k4vx)/6

print x1, mass*vx1, x2, mass*vx2
t = t + dt

m=m+1
if(m < m_max) reread
m=0

outfile(n) = sprintf("f/%d.jpg",n+10000)  #出力ファイル名
title(n) = sprintf("t = %4.2f",t )  #タイトル名
unset label 
set label title(n)  font 'Times,18'  at -1.9 , 3.6 
set output outfile(n)
plot DATA using 1:2 with line linecolor rgb "red", DATA using 3:4 with line linecolor rgb "green"

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

2.アニメーション作成用のjpgファイルを格納するフォルダ「f」を上記と同じフォルダに作成します。
3.gnuplot を起動し、次のコマンドを打ちます。

 
load "Duffing.plt"

これで、フォルダ名「f」の中に、連番jpg ファイルが作成されます。 そして、連番番号のjpgファイルをフリーソフト giam を利用して、結合することでアニメーションを作成します。

課題メモ

1.連番 jpg からのアニメーションではなく、直接gif アニメーションを作成する

参考

【gnuplotでアニメーション】直接 gif アニメーションを出力する

※ミスの指摘、質問、コメントなどがございましたら、HTMLフォーム、もしくはinfo:natural-science.or.jpから、お待ちしております。

目次

第1章 序章

第2章 対称ポテンシャルによる量子散乱

ゴール:水素原子のクリプトン原子による散乱に対する散乱断面積の計算

  • ・ヌメロフ法による常微分方程式の解法
  • ・球ベッセル関数
  • ・ルジャンドル多項式
  • ・位相シフト→断面積を計算するアルゴリズム

その他



▲このページのトップNPO法人 natural science トップ

関連記事

計算物理学







▲このページのトップNPO法人 natural science トップ