[作ってみた] [プログラミング] 反応拡散方程式で遊ぶ


我々には少し先までしか見えないが、そこになすべきことをたっぷりと見付けることができる。
―アラン・チューリング

タテジマキンチャクダイ

タテジマキンチャクダイ

以前,セル・オートマトンに関係するプログラミングで作った反応拡散方程式シミュレーションのお話です.
眠らせておくのももったいないので,ここで紹介しておきます.

反応拡散方程式というのは,ある空間上での何かの物質の濃度の変化を表現した微分方程式です.もともとは化学の分野で反応をモデル化する際に考えられたものですが,その解が複雑な振る舞いをすることから様々な分野で興味の対象となり,実世界の模様との関連性が指摘されたりもしています.

今回はこの反応拡散方程式のシミュレータの紹介です.まずは動画をどうぞ.

この模様はチューリングパターンと呼ばれ,動物の模様や植生の分布などに類似の模様が見られることがあります.

今回扱う反応拡散方程式は1980年代に提案されたGray-Scottモデルと呼ばれるもので,以下の形に書き下されます.

  \frac{\partial u}{\partial t} = D_{u}\nabla^2u - uv^2 + F(1 - u)\\   \frac{\partial v}{\partial t} = D_{v}\nabla^2v + uv^2 - (F + k)v

 

uとvという2種類の物質(?)についての式で,左辺が変化の度合い,右辺第1項が拡散,第2項が反応,第3項が生成と消滅に対応しています.Fとかkとかはパラメータです.

解析的に解けなさそうな式ですが,差分方程式にしてしまえばコンピュータで解くのはそれほど困難ではありません.得られるパターンが独特かつパラメータによって多様に変化するので,数値解析の例題としてはなかなか楽しい部類に入るのではないかと思います.

私がこの反応拡散方程式にはじめて出会ったのは,まともにプログラミングを学んだ少し後の頃でしたが,当時紹介されていた動画などを見てすっかり気に入ってしまった記憶があります.

1つは「ぬるぬる」している点です.特にはじめて触った頃のプログラミングというのはテキストや数式を扱うことが多いため,どうにも真四角でお固い世界だと錯覚しがちです.こんな「ぬるぬる」は当時の私にとってあまりコンピュータっぽくなかったというのがあります.
もう1つは,式が簡単であることです.すごく複雑なものが複雑な式から出てくるならまあわかりますが,この式自体はそれほど複雑ではありません.
各部分では単に隣接した点の情報のみから変化しているだけなのですが,それでも全体としてパターンが現れるというのはなかなか面白いと思います.

プログラムを置いておきます.手元の 64bit Windows 用(7で動作確認済み)なので,実行できない方は一番下のソースコードからコンパイルしてみて下さい.
マウスでドラッグした部分が消せます.キーボードの左2列(1,2,q,w,a,s,z,x)でパラメータの変更,iキーで分布を初期化できます.

rd_system_2 [追記: 録画用でウィンドウが小さいままだったので大きくしました.]

rd_system_32bit [さらに追記: 32bit版も用意しました]


私の大学の機械工学科には「シミュレーション & ビジュアライゼーション演習」という授業があります.この授業では何らかのシミュレーションの可視化に挑戦します.

実際に扱う対象は連続体というよりはどちらかというと剛体系の人が多いかと思いますが,何にしても可視化してみるというのは楽しいもので,あれこれ計算した結果が絵になって出てくる達成感がある気がします.
私も昨年度TAで参加したときに楽しませてもらいました.

プログラミングでちょっと遊んでみたいという人は,マイコンのような動くものにつなげてみるか,あるいはこういう可視化のできるもので何か作ってみると手ごたえを感じやすいと思います.


実装はC++/OpenGL+freeglut,方程式は1次の陽解法です.私はこういうものを作るときはたいていC# & 速度の必要な部分のみC++で実装してしまうのですが(実はこれもそうでした),今回は公開する前に移植しやすいようにC++で書き直してみました.

ソースはgithubにあげておきます.
githubを使うのははじめてなので,何か不具合があったら申し訳ございません.

ソース


[作ってみた] [プログラミング] 反応拡散方程式で遊ぶ」への6件のフィードバック

  1. ピンバック: [作ってみた] [プログラミング] 剛体3重振り子 | ぞうさんの何でもノート

  2. ピンバック: [フラクタル] マンデルブロー集合 10の200乗倍への挑戦 | ぞうさんの何でもノート

  3. ゼタ

    初めまして、ゼタと申します。
    forkさせていただき、現状はリファクタリングのみですが、
    今後機能追加等、いろいろ弄ってみようと思ってます。
    ところで、権利まわりは特に主張されていないという認識でよろしいでしょうか?
    コピーライトを書きたかったのでどうしようかと思いまして。
    (少なくとも、今後もGitHub上ではfork元としてoniku64さんの名前は残りますね。)
    なお、特に商用利用の予定はございません。
    ご査証のほど、よろしくお願いいたします。

    返信
    1. takamoto 投稿作成者

      ゼタ様
      ご連絡ありがとうございます。

      権利まわりは私が不勉強でした。GitHubに上げた時点で明記しておくべきでしたね。大変申し訳ございません。

      ソフトウェアは好きにして頂いて構いません。
      MITライセンスとしてGitHubに反映しておきました。ご確認ください。

      返信
      1. ゼタ

        ライセンスの件、確認いたしました。
        ご対応ありがとうございます!

        あと、実はわざわざforkしたのには理由がありまして、
        アプリケーションとしての出来はすごく良いのですが、
        プログラムソースとして会社だとぶん殴られるレベルだったので、
        見て頂きたい狙いがありました。
        具体的にいうと、mallocで確保している領域をdelete解放していたり、
        glutMainLoop()の後に処理が書いてある(デフォルトでは関数中でexit()が呼ばれる)などです。

        自分の気づいた点に関しては、forkリポジトリの方に反映させましたので、
        コミットログをご参考くださいませ。
        (ご無礼大変失礼しました。)

        返信
        1. takamoto 投稿作成者

          返事が遅くなってしまって申し訳ございませんでした。

          かなり昔に書いたコードを見直さずに引っ張ってきたので、いろいろと酷い箇所があったようです。ご指摘ありがとうございます。
          mallocの件は書いた時点で気づくべきでした。

          少しだけ手直しをしてGithubの方へ反映させておきました。

          ありがとうございました。

          返信

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です