Help us understand the problem. What is going on with this article?

Haskell+OpenGLでロトカ=ヴォルテラ方程式を描画してみた

More than 3 years have passed since last update.

先日「ロトカ=ヴォルテラ方程式知ってる?」と聞かれて人工生命などに熱中していた懐かしい記憶が蘇りました。そう言えば昔はPythonで相図を書いた記憶があったのですが今だったらとHaskellを使ってロトカ=ヴォルテラ方程式の時間発展を描画してみました。描画にはOpenGLを使いました。

OpenGLを使ってゴニョゴニョするのにGLUTを使っています。OpenGL Utility ToolkitのHaskellバインディングです。予めインストールしておきましょう。

$ cabal install GLUT

ロトカ=ヴォルテラ方程式はwikiのものをそのまま使用しました。あとHaskellでOpenGLを使うときはrunhaskellではなくghcでちゃんとコンパイルしてから実行しないので気をつけてください。これなんで何でしょう。以下のコードは完全なコードなのでコピペして試すことができます。

import Data.IORef
import Graphics.UI.GLUT

type Line = [(GLfloat, GLfloat)]

-- 計算用パラメータ
dt    = 0.01
alpha = 5.0
beta  = 3.0
gamma = 2.0
delta = 2.0

start = (1.0, 1.0)

main :: IO ()
main = do
    (_progName, _args) <- getArgsAndInitialize
    initialDisplayMode $= [DoubleBuffered]
    _window <- createWindow "Lotka-Volterra Equation"
    orbit <- newIORef [start]
    idleCallback $= Just (idle orbit)
    displayCallback $= display orbit
    mainLoop

idle :: IORef Line -> IdleCallback
idle orbit = do
    ps <- get orbit
    let (x, y) = head ps
    let nx = (\x' -> x + x' * dt) $ x * (alpha - beta * y)
    let ny = (\y' -> y + y' * dt) $ (negate y) * (gamma - delta * x)
    orbit $= ((nx, ny) : take 500 ps)
    postRedisplay Nothing

display :: IORef Line -> DisplayCallback
display orbit = do 
    clear [ColorBuffer]
    ps <- get orbit
    let scale = 0.5
    let (ox, oy) = (1.2, 2.1)
    renderPrimitive Lines $
        mapM_ (\(x, y) -> vertex $ Vertex3 ((x - ox) * scale) ((y - oy) * scale) 0) ps
    swapBuffers

簡単に解説していきます。

main :: IO ()
main = do
    (_progName, _args) <- getArgsAndInitialize
    initialDisplayMode $= [DoubleBuffered]
    _window <- createWindow "Lotka-Volterra Equation"
    orbit <- newIORef [start]
    idleCallback $= Just (idle orbit)
    displayCallback $= display orbit
    mainLoop

最初の方はほとんどお決まりの処理です。大事なのはidleCallbackdisplayCallbackでそれぞれ待機時、描画時に呼ばれる関数を指定します。まず待機時に呼ばれる関数を見て行きましょう。

idle :: IORef Line -> IdleCallback
idle orbit = do
    ps <- get orbit
    let (x, y) = head ps
    let nx = (\x' -> x + x' * dt) $ x * (alpha - beta * y)
    let ny = (\y' -> y + y' * dt) $ (negate y) * (gamma - delta * x)
    orbit $= ((nx, ny) : take 500 ps)
    postRedisplay Nothing

この関数は毎フレーム呼ばれるのでここでロトカ=ヴォルテラ方程式の計算をしています。orbitというIORefを用意してここに計算した軌道を貯めていっています。postRedisplayを呼ぶとdisplayCallbackが実行されます。

display :: IORef Line -> DisplayCallback
display orbit = do 
    clear [ColorBuffer]
    ps <- get orbit
    let scale = 0.5
    let (ox, oy) = (1.2, 2.1)
    renderPrimitive Lines $
        mapM_ (\(x, y) -> vertex $ Vertex3 ((x - ox) * scale) ((y - oy) * scale) 0) ps
    swapBuffers

ここではorbitの軌道を実際に描画しています。実はmainでDisplayModeにDoubleBufferedを指定していたのでダブルバッファリングを使うことができます。この場合swapBuffersを呼べばバッファの内容が描画されます。

だいたいこんな感じです。何も難しいことはしていません。実際に実行すると軌道がぐるぐる動くのが見れて楽しいと思います。発散するのは計算誤差からでしょうか・・・

lotz
実用関数型プログラミング言語 Haskell の情報を発信しています
http://lotz84.github.io/
folio-sec
誰もがかんたんに資産運用することができるサービス「フォリオ」を作っているFinTech系スタートアップ
https://corp.folio-sec.com/
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away