今回はfortranとpythonを使ってシミュレーションをしてみましょう。
何回かに分け、基本のところからやりたいと思います。
初回は「相互作用のない二つの粒子」のモデルを作ってどのような運動をするのか見てみます。
手順としては
- fortran で運動の計算を行う
- python(vpythonというライブラリを使う) でシミュレーション結果を見る
という方法を取りたいと思います。
ではさっそくやっていていきましょう。
相互作用のない2粒子のシミュレーション
相互作用のない粒子を考えます。
つまり、反発もしないので互いにすり抜けてしまう場合を考えています。
少し奇妙なモデルですが、これを基礎としていろいろな場合についてシミュレーションしていきます。
fortranで運動の計算を行う
ではまず、fortranで運動の計算を行います。
状況としては
重さ1の粒子2つが速さ1と-1でぶつかる場合を考えます。
このとき、外力、相互作用ともに0です。
このことを踏まえ、コードを書くと次のようになります。
program main
implicit none
integer,parameter :: N=2
real :: x(N),vx(N)
real :: Fx(N),Fx_next(N)
real :: dt = 0.01
integer :: i
integer :: i_max = 2000
integer :: i_rec = 10
integer :: fout = 11
x(1)= 1.0
vx(1)=1.0
Fx(1) = 0.0
x(2)= 7.0
vx(2)=-1.0
Fx(2) = 0.0
open(fout,file = "two-zero-force.txt")
write(fout,*)0.0,x(1),x(2),vx(1),vx(2)
do i=1,i_max
x(1)=x(1) + vx(1) * dt + 0.5 * Fx(1) * dt * dt
Fx_next(1) = 0.0
vx(1) = vx(1) +0.5 * (Fx(1) + Fx_next(1)) *dt
Fx(1) = Fx_next(1)
x(2)=x(2) + vx(2) * dt + 0.5 * Fx(2) * dt * dt
Fx_next(2) = 0.0
vx(2) = vx(2) +0.5 * (Fx(2) + Fx_next(2)) *dt
Fx(2) = Fx_next(2)
if(mod(i,i_rec) == 0) then
write(fout,*) i*dt,x(1),x(2),vx(1),vx(2)
endif
enddo
close(fout)
end program main
少しわかりにくいと思うので説明を入れます。
ファイルを作りデータを記入
計算結果を
open(fout,file = "two-zero-force.txt") write(fout,*)0.0,x(1),x(2),vx(1),vx(2)
で「two-zero-force.txt」という名前のファイルを作成させます。
そこに左から「時間、粒子1の座標x(1)、粒子2の座標x(2)、粒子1の速度vx(1)、粒子2の速度vx(2)」という順番で得られたこれらのデータを入れていきます。
実際につくられたファイルを見ると
のようにデータが書き込まれていることが分かります。
速度ベルレ法
肝心の計算方法ですが、今回は速度ベルレ法を使いました。
以下の部分です。
do i=1,i_max x(1)=x(1) + vx(1) * dt + 0.5 * Fx(1) * dt * dt Fx_next(1) = 0.0 vx(1) = vx(1) +0.5 * (Fx(1) + Fx_next(1)) *dt Fx(1) = Fx_next(1) x(2)=x(2) + vx(2) * dt + 0.5 * Fx(2) * dt * dt Fx_next(2) = 0.0 vx(2) = vx(2) +0.5 * (Fx(2) + Fx_next(2)) *dt Fx(2) = Fx_next(2) enddo
単位時間(dt)当たりの位置x(1),x(2)、速度vx(1),vx(2)、力Fx(1),Fx(2)etc…から次の時間における位置、速度を決定する方法です。
詳しくは
をどうぞ。
pythonでシミュレーションできているか確認
ではこのデータがほんとうに粒子の動きを再現できているか確認してみましょう。
以下のpythonコードで確認できます。
import numpy as np
from vpython import *
data = np.loadtxt('two-zero-force.txt')
t = data[:, 0]
x1 = data[:, 1]
x2 = data[:, 2]
# シミュレーションの描画スペースの確保と設定
mydisplay = canvas(autoscale=False) # 描画スペースを作る。
mydisplay.x = mydisplay.y = 0 # パソコン画面上での描画スペースの位置
mydisplay.width = mydisplay.height = 640 # 描画スペースの幅と高さ
mydisplay.background = color.white# 背景色
mydisplay.foreground = [0, 0, 0] # デフォルトの前景色。
mydisplay.center = vector (5, 0, 0)# 描画の原点を移動
mydisplay.range = 20 # 描画する座標の範囲
particle1 = sphere(pos=vector (x1[0], 0, 0), radius=1, color=vector(0, 1, 0))# 粒子1を初期位置に描く
particle2 = sphere(pos=vector(x2[0], 0, 0), radius=1, color=vector(1, 0, 0))# 粒子2を初期位置に描く
myclock = label(pos=vector(0, -10, 0), text='time is %5.1f' % t[0]) # 時計を初期時刻で描く
for k in range(1, t.size): # ループを用いてパラパラ漫画を作る
rate(25) # パラパラ漫画をめくる速さ
particle1.pos = vector(x1[k], 0, 0) # 粒子1の位置を変更して描く
particle2.pos = vector(x2[k], 0, 0) # 粒子2の位置を変更して描く
myclock.text = ('time is %5.1f' % t[k])
ちなみに今回のシミュレーションに使ったライブラリはvpython7で、vpython6を使っている方は文法が若干違うので注意してください。
これを実行すると
まず、互いに近づいていき…↓
相互作用がないので貫通し…↓
互いに逆方向に離れていきます↓
うまくシミュレーションできました！！
次回は相互作用がある場合を考えたいと思います。
コメント
[…] […]