オイラー法でローレンツアトラクタが見たい on Haskell

こんにちは! やきつか(@sy_t_)です。
最近やらなきゃいけないことが多いのですが、全然関係無いことをしてしまう自分に絶望してます。

はじめに
準備
オイラー法を実装する
ローレンツアトラクタを見る
- 実行結果
まとめ

はじめに

Haskellの勉強中に、「オイラー法の実装ならいい練習になるやん!」と思って実装しました。
プログラミングの一番の応用って数値計算だと思うんですよね(過激派)
オイラー法の実装だけでは面白くないので、適当な常微分方程式解いて終わります。
解析解のない微分方程式にはいろいろありますが、個人的に好きなローレンツ方程式にしました。

準備

それでは、対戦よろしくお願いします。
最近ハマっているパイプラインっぽい演算子を使わせて下さい。

(|>) :: (a -> b) -> (b -> c) -> (a -> c)
f |> g = g.f

タプルを使いたくないので、時空座標を表すCoordinate型を作ってShowのinstanceにします。

data Coordinate = Crd {getT :: Double, getX :: Double, getY :: Double, getZ :: Double}

instance Show Coordinate where
    show (Crd t x y z) =  map show |> foldl1 (\s0 s1 -> s0 ++ "," ++ s1) $ [t, x, y, z]

ローレンツ方程式もついでに定義しておきます。

lorenz :: Coordinate -> Coordinate
lorenz (Crd t x y z) = Crd t x' y' z'
    where
        x' = (-p)*x + p*y
        y' = (-x)*z + r*x - y
        z' = x*y - b*z

        p = 10
        r = 28
        b = 8/3

この段階で
f:id:yakituka:20190826182615p:plain
こんな感じのことができます。

オイラー法を実装する

準備がととのったのでオイラー法を実装します。

eulerMethod :: Double -> Double -> (Coordinate -> Coordinate) -> Maybe Coordinate -> Maybe Coordinate
eulerMethod _ _ _ Nothing = Nothing
eulerMethod h tmax f (Just vector@(Crd t x y z))
    | t < tmax = Just vector'
    | otherwise = Nothing
    where
        vector' = Crd t' x' y' z'

        t' = t + h
        x' = x + h * getX (f vector)
        y' = y + h * getY (f vector)
        z' = z + h * getZ (f vector)

刻み幅h、数値解の見たい上限時間tmax
微分方程式が $\frac{d\vec{x}}{dt}=f(t,\vec{x})$ と表される時のf
前ステップでの座標vectorを引数にとります。

Maybe Coordinate -> Maybe Coordinate

となっているのはiterate関数を使いやすくしたかっただけです。
t < tmaxの時はオイラー法で次のステップの座標を計算してJustで包んで返します。
tmax以降の時間はいらないので、その時はNothingを返します。

全体をNumericalCalcモジュールとします。

module NumericalCalc where

f |> g = g.f

data Coordinate = Crd {getT :: Double, getX :: Double, getY :: Double, getZ :: Double}

instance Show Coordinate where
    show (Crd t x y z) =  map show |> foldl1 (\s0 s1 -> s0 ++ "," ++ s1) $ [t, x, y, z]

eulerMethod :: Double -> Double -> (Coordinate -> Coordinate) -> Maybe Coordinate -> Maybe Coordinate
eulerMethod _ _ _ Nothing = Nothing
eulerMethod h tmax f (Just vector@(Crd t x y z))
    | t < tmax = Just vector'
    | otherwise = Nothing
    where
        vector' = Crd t' x' y' z'

        t' = t + h
        x' = x + h * getX (f vector)
        y' = y + h * getY (f vector)
        z' = z + h * getZ (f vector)

lorenz :: Coordinate -> Coordinate
lorenz (Crd t x y z) = Crd t x' y' z'
    where
        x' = (-p)*x + p*y
        y' = (-x)*z + r*x - y
        z' = x*y - b*z

        p = 10
        r = 28
        b = 8/3

ローレンツアトラクタを見る

全てが整ったので、Main.hsを実装します。先に全体を記載します。

import NumericalCalc

h = 0.01
tmax = 10000
initialCoordinates = Crd 0 1 1 1

main = do
    let trace = iterate (eulerMethod h tmax lorenz) (Just initialCoordinates)
    let n = 10000
    let path = "result.dat"

    writeFile path ""
    sequentiallyWrite path (take n trace)

sequentiallyWrite filePath [] = return ()
sequentiallyWrite filePath (Nothing:xs) = return ()
sequentiallyWrite filePath ((Just x):xs) = appendFile filePath (show x ++ "\n") >> sequentiallyWrite filePath xs

import NumericalCalc

で先程作ったモジュールを使えるようにします。

h = 0.01
tmax = 10000
initialCoordinates = Crd 0 1 1 1

は刻み幅と上限と初期条件です。
sequentiallyWrite関数はファイルパスとMaybeなリストを引数に、失敗(Nothing)するか終端にたどり着くまでひたすらファイルに印字する関数です。(同機能のものが既にありそうな気がする、誰か教えて下さい)
main部はオイラー法の反復合成リストtraceとステップ数n、パスを宣言してsequentiallyWrite関数に渡してます。

Main.hsを実行して、これによって得られたresult.datをgnuplotで適当に表示します。
f:id:yakituka:20190827231432p:plain