【計算物理】サンプルコード#1：固い壁に囲まれた粒子【FORTRAN】

【計算物理】サンプルコード#1：固い壁に囲まれた粒子【FORTRAN】

■概要

数値計算のためのプログラミングを解説していこうと思うが、コードの抽象的な説明から入るのは嫌なので、単純なサンプルコードを紹介しながら、各部分を説明していく形式にしようと思う。言語は個人的に最も慣れ親しんでいるという理由でとりあえずFORTRAN90を用いる。

具体的なビジョンもないし、特に深い理由はないが、最初の例としてなんとなく視覚的にわかりやすそうであったため、固い壁に囲まれた粒子の2次元の運動を追う非常に単純な計算を選んでみた。

【数字に強いビジネスパーソンになる】外資系投資銀行流のビジネス・シミュレーション術講座

上のアニメーションを描くのに使った計算コードは以下の通りだ。以下で各部分の意味を解説する。

program main

 implicit none
 integer :: i 
 integer :: fstp
 real(8) :: x,y
 real(8) :: vx, vy
 real(8) :: dt
 real(8) :: a

 open(10,file='sample.dat')

 a    = 1.0d0
 fstp = 1000
 dt   = 0.01d0

 call random_number(vx)
 call random_number(vy)
 call random_number(x)
 call random_number(y)

 write(10,*) '# x, y'
 write(10,*) x, y

 write(10,*)
 write(10,*) 

!----------------------------
! loop start (!で始まる部分は処理されません)
!----------------------------
do i=1, fstp

  x = x + vx*dt
  y = y + vy*dt

  if (x >= a) then
   x = a
   vx = -vx
  end if

  if(x <= 0 ) then
   x=0.0d0
   vx =-vx
  end if

  if (y >= a) then
   y = a
   vy = -vy
  end if

  if(y <= 0) then
   y=0.0d0
   vy =-vy
  end if

  write(10,'(4e19.9)') x, y

  write(10,*)
  write(10,*)

end do

end program

■コードの説明

物理的な内容は単純で、扱っているのは二次元平面に拘束された一つの粒子の初期条件$(x,y, v_x, v_y)$が乱数で与えられ、$x \ (\text{or} \ y )=0,a$で壁に当たり弾性衝突するというだけの系である。

まず、最初の部分

program main

 implicit none
 integer :: i 
 integer :: fstp
 real(8) :: x,y
 real(8) :: vx, vy
 real(8) :: dt
 real(8) :: a

は、プログラミングの開始を宣言した後、どんな記号にどんな種類の値を対応付けるか、ということを宣言している。implicit noneはここではとりあえずおまじないと思ってもらえればいい。

integerというのが整数型、realというのが実数型であることを示している。real(8)の(8)は変数の精度を表しているのだが、これも最初は深く考えなくていい。

次に

open(10,file='sample.dat')

の部分で、計算を行った結果を書き出すファイルを呼び出している。数字の10はファイルの識別番号として適当に与えており、sample.datをファイル名としている。

続いて、先ほど使いますよ、と宣言していた変数に、実際の値を入れている。

 a    = 1.0d0
 fstp = 1000
 dt   = 0.01d0

 call random_number(vx)
 call random_number(vy)
 call random_number(x)
 call random_number(y)

ここでは箱の一片の長さ$a$として1.0という値を入れている。fstpはステップ数で、dtは時間幅の刻みある。

その後、初期条件x, y, vx, vyは乱数で決定している。call random_number(...)というのがそれ用に組み込まれている手続きで、こうだけ書いておけば(...)内の変数に0から1の間の乱数を入れてくれる。

その次の

 
 write(10,*) '# x, y'
 write(10,*) x, y

 write(10,*)
 write(10,*) 

では、先ほど呼び出した10番のファイルに書き込みを行っている。「'」で挟んだ '# x, y' はxとyの値ではなく、そのまま文として書き込まれる。#を先頭に入れておくことで、値を読み込むときにスキップしてくれる。そして次の段の write(10,*) x, yで実際に値を書き込み、 write(10,*)だけの指示は、は空の一列を書き込むことを意味する。二列空けることで、それをはさむ二つの文字列が二つの独立なブロックとなる。

そしていよいよ

 
do i=1, fstp

   x = x + vx*dt
   y = y + vy*dt

...

end do

の部分で具体的な計算を行っている。ここでの指示は、do ... end doで挟まれた部分を、i=1からfstp回繰り返すということを意味している。

そして、各ステップごとに、運動方程式$dx/dt=v_x$を積分している（$y$成分も同様）。掛け算には*という記号を用いる。

x=x+vx*dtは数学的には違和感のある式だが、右辺にあるxがその処理に至る前の値を持たされたxで、vx*dtを加えたものを新たなxの値として、左辺のxに担わせる、ということを意味している。

途中に挟まれている

 
   if (x >= a) then
   x = a
   vx = -vx
   end if

   if(x <= 0 ) then
   x=0.0d0
   vx =-vx
   end if

の部分が壁にぶつかったときの処理だ。if文と呼ばれるこの指示は、if (条件) thenとし、その条件が満たされた場合、次に書かれた処理を行う、というものだ。そして最後にend ifとして処理を閉じる。

ここでは、xがa以上（あるいは0以下）になった場合、x=a（あるいは0）として、x方向の速度vxの符号を反転させよ、という指示になっている。y成分も同様である。

これらの指示によって最初のアニメーションのような計算結果が得られる。このコードだけでは描画はできず、別の処理が必要となるので、それについては別の記事で説明する。

以降、粒子数を増やす、粒子間の衝突を含める、粒子の分布の変化を追う、などこのコードを発展させながら遊んでみようと思う。

内定率94％！完全無料のITスクール【ITCE Academy】