Haskell で振動子

Haskell Physics numerical analysis

Haskell を使って，振動子を数値的に解こう。

まずは振動子の基本的な道具立て。

-- Oscillator.hs

module Oscillator where


-- まず、1自由度系の状態を表す型 State を定義しよう。
-- これは 時刻 (t) 位置座標 (x) 速度 (v) からなる。

data State = State { t, x, v :: Double }

-- 例えば、時刻 0.0 位置 1.0 速度 0.0 の状態は、
-- 次のようにして作ることができる。
--   s = State 0.0 1.0 0.0
-- この状態の位置座標を取得するには次のようにする。
--   x s

-- 質量を表す型 Mass を定義しておこう。これは Double と同じで良い。

type Mass = Double

-- 力を表す型 Force を定義しよう。
-- 力は、一般に、状態の関数であるから、次のようにすると良いだろう。
-- これを力関数と呼ぶことにする。

type Force = State -> Double

-- 系は質量と力によって完全に特徴づけられるので、
-- これらを合わせて System という型を作ろう。

type System = (Mass, Force)

-- 系がある状態にあるときの、加速度を得る関数を用意する。
-- これは、ニュートンの運動方程式であり、
-- これによって、系の時間発展が実現される。

acceleration :: System -> State -> Double
acceleration (m,f) state = (f state) / m


-- 調和振動を引き起こす力関数を生成する関数を定義する。
-- これは、バネ定数を受け取って、力関数を返す。

type SpringConstant = Double

makeHarmonicForce :: SpringConstant -> Force
makeHarmonicForce k state = - k * (x state)

-- 例えば、バネ定数 1.0 のときの力関数を得るためには、次のようにする。
--   makeHarmonicForce 1.0

-- 摩擦力を追加する関数を定義しよう。
-- 摩擦係数と元になる力関数を受け取って，新しい力関数を返す。

addDamping :: Double -> Force -> Force
addDamping coefficient f0 state =
  (f0 state) - coefficient * (v state)

-- 強制力を追加する関数も定義しよう。
-- 振幅、角振動数、位相、元になる力関数を受け取って、
-- sin関数型の強制力を付け加えられた力関数を返す。

addDriving :: Double -> Double -> Double -> Force -> Force
addDriving amp angFreq phase f0 state =
  (f0 state) + amp * sin (angFreq * (t state) + phase)


-- 系の全エネルギーを得る関数も用意しておこう。
-- 運動エネルギーとポテンシャルエネルギーの和であり、
-- 状態の関数として表される。
-- また、パラメーターとして、バネ定数と質量も必要だ。

energy :: SpringConstant -> Mass -> State -> Double
energy k m state = (kinetic m state) + (potential k state)
  where
    kinetic :: Mass -> State -> Double
    kinetic m state = (m / 2.0) * ((v state) ** 2)
    potential :: SpringConstant -> State -> Double
    potential k state = (k / 2.0) * ((x state) ** 2)

次に、差分法を実装しよう。基本的な Euler法と Runge-Kutta法（4次）を実装する。

-- FiniteDifferenceMethods.hs

module FiniteDifferenceMethods where

import Oscillator


-- 差分法を特徴づける型 FiniteDifferenceMethod を定義したい。
-- これは、状態を受け取って、新しい状態を返す。
-- 新しい状態を作るには、時間刻み幅と、時間発展に関する系の情報が必要だ。
-- そこで、次のように定義できるだろう。

type TimeDiff = Double

type FiniteDifferenceMethod = TimeDiff -> System -> State -> State


-- Euler法は、次のように定義できる。

euler :: FiniteDifferenceMethod
euler dt sys (State t0 x0 v0) = State t1 x1 v1
  where
    t1 = t0 + dt
    x1 = x0 + dt * v0
    v1 = v0 + dt * a0
      where a0 = acceleration sys (State t0 x0 v0)


-- Runge-Kutta法（4次）は、次のように定義できる。

rk4 :: FiniteDifferenceMethod
rk4 dt sys (State t0 x0 v0) = State t1 x1 v1
  where
    t'  = t0 + dt / 2.0
    t1  = t0 + dt

    dx1 = dt * v0
    dv1 = dt * (acceleration sys (State t0 x0 v0))

    x'1 = x0 + dx1 / 2.0
    v'1 = v0 + dv1 / 2.0

    dx2 = dt * v'1
    dv2 = dt * (acceleration sys (State t' x'1 v'1))

    x'2 = x0 + dx2 / 2.0
    v'2 = v0 + dv2 / 2.0

    dx3 = dt * v'2
    dv3 = dt * (acceleration sys (State t' x'2 v'2))

    x'3 = x0 + dx3
    v'3 = v0 + dv3

    dx4 = dt * v'3
    dv4 = dt * (acceleration sys (State t' x'3 v'3))

    x1 = x0 + (dx1 + 2.0 * dx2 + 2.0 * dx3 + dx4) / 6.0
    v1 = v0 + (dv1 + 2.0 * dv2 + 2.0 * dv3 + dv4) / 6.0

そして、系の時間発展を得るためのソルバーを実装する。

-- Solver.hs

module Solver where

import Oscillator
import qualified FiniteDifferenceMethods as F


-- 系の時間発展をシミュレートする関数を定義しよう。
-- これは、初期状態を受け取って、状態のリストを返すものとする。

solve ::
  F.TimeDiff -> Double -> F.FiniteDifferenceMethod -> System ->
  State -> [State]

solve dt endtime fdm sys state0 = solve' nsteps state0
  where
    nsteps = ceiling (endtime / dt)

    solve' :: Int -> State -> [State]
    solve' n state
      | n <= 0 = [state]
      | otherwise = state : (solve' n' state')
           where n' = n - 1
                 state' = fdm dt sys state


-- 状態のリストから、特定の物理量のリストを得るための関数を用意しよう。
-- 時刻と物理量の組のリストを返すとする。

makeData :: (State -> Double) -> [State] -> [(Double,Double)]

makeData operator states = map f states
  where
    f :: State -> (Double, Double)
    f state = (t state, operator state)

-- 物理量のリストを出力する関数を用意しよう。
-- gnuplotで表示することを念頭において、
-- 時刻と物理量を空白で区切って出力する。

printData :: [(Double,Double)] -> IO ()

printData [] = putStrLn ""

printData ((t,q):rest) =
  do putStrLn $ show t ++ " " ++ show q
     printData rest

ちなみに、ここまでのファイルの行数をカウントすると次のようになる。

$ egrep -v '^\s*(--|$)' -c Oscillator.hs FiniteDifferenceMethods.hs Solver.hs
Oscillator.hs:23
FiniteDifferenceMethods.hs:32
Solver.hs:25

さて、それでは、適当にパラメーターを決めて、計算してみよう。

-- Main.hs

import Oscillator
import qualified FiniteDifferenceMethods as F
import qualified Solver as S

m       = 1.0
k       = 1.0
f0      = makeHarmonicForce k
sys     = (m, f0)
state0  = State { t = 0.0, x = 1.0, v = 0.0 }
endtime = 20.0
dt      = 0.01

例えば、次の1行を追加すると、Euler法の計算結果を得られる。

main = S.printData $ S.makeData x $ S.solve dt endtime F.euler (m, f0) state0

Euler法、Runge-Kutta法、解析解をプロットすると、次のようになる。

また、減衰振動の結果を得るためには、次のようにすれば良い。

main = S.printData $ S.makeData x $ S.solve dt endtime F.rk4 (m, addDamping 0.1 f0) state0

摩擦係数を変えたときの計算結果は次のようになる。

とりあえず、今日はここまで。

Haskell を使うと、概念が

  • 自然に
  • 簡潔に

表現できると思う。