Elmでやってみるシリーズ10:ボールの衝突回数で円周率を計算する
id:wamanさんからひっそりと提案もらいましたので、やってみました。
今回は記事自身もElmで書きましたので、github-pages上の全画面(github上のソース)からどうぞ。
以下にも一応hatena blog記事中にiframeで表示しましたが、スクロールや文字の配置が読みにくい・もしくはJSの実行速度が遅いです。hatena blogではiframeにseamless属性つけられば良いのだが。iPhoneのSafariでもiframe中だとうまくいきませんが、全画面なら動くようです。
ソースはこちら。
PiByBall.elm
module PiByBall where import Window import Debug (log) import Graphics.Input (input, Input, button, dropDown) import Keyboard data Status = Pause | Running type State = { stat:Status, x1:Float, x2: Float, v1: Float, v2: Float, count: Int, ratio: Float } data Event = Start | Stop | TimeTick | ChangeN Int frameRate = 320 -- 画面更新頻度 initialState : State initialState = { stat=Pause, x1=-200, x2=-100, v1=1, v2=0, count=0, ratio=1 } inputSignal : Signal Event inputSignal = let f running = if | running -> Start | otherwise -> Stop in merges [(f <~ inpRunning.signal), (ChangeN <~ inpN.signal), (sampleOn (fps frameRate) (constant TimeTick))] colide v1 v2 r = (((r-1)*v1 + 2*v2)/(r+1), (2*r*v1 - (r-1)*v2)/(r+1), 1) nextState : Event -> State -> State nextState = \event ({stat, x1, x2, v1, v2, count, ratio} as state) -> case event of Start -> (log "Start" {initialState|stat<-Running, ratio<-ratio}) Stop -> (log "Stop" {state|stat<-Pause}) ChangeN n -> (log "ChangeN" {initialState|ratio<-100^toFloat n}) TimeTick -> if stat == Running then let (new_v1, new_v2, countIncl) = if | x1+v1>= x2+v2 -> colide v1 v2 ratio | x2+v2 >= 0 -> (v1, -v2, 1) | otherwise -> (v1, v2, 0) in (log "timetick" State Running (x1+new_v1) (x2+new_v2) new_v1 new_v2 (count+countIncl) ratio) else (log "Pause" {initialState|ratio<-ratio}) currentState : Signal State currentState = foldp nextState initialState inputSignal description1 = [markdown| ## ボールをぶつけるだけで円周率がわかる? ### シミュレーションの舞台 以下のように2つの質点M1,M2と壁を考えます。<br/> |] description2 = [markdown| 表示上の判り易さのために、質点の大きさに差を付けていま<br/> すが、表示されている大きさは質点の質量の比率には対応し<br/> ていません。 ### 質量について M1とM2の質量をそれぞれm1,m2としたとき、m1とm2の比率<br/> を以下とします。 m1:m2 = 100^N : 1 ここでNは0以上の整数値です。Nに応じて上記の比率は具体<br/> 的には以下のようになります。 <table border="1"> <tr><th>N</th><th>M1の質量(100^N) : M2の質量</th></tr> <tr><td>0</td><td>1 : 1</td></tr> <tr><td>1</td><td>100 : 1</td></tr> <tr><td>2</td><td>10000 : 1</td></tr> <tr><td>3</td><td>1000000 : 1</td></tr> <tr><td>:</td><td>:</td></tr> </table> ### シミュレーション 前提として、質点および壁は完全弾性衝突するとします。<br /> そして質点M1に右向きの適当な初速を与え、M2のM1および<br /> 壁に対する衝突回数をカウントします。</p> 実際にやってみましょう。まず以下でNは変更せずに(N=0の<br /> まま)「開始」ボタンを押してみて下さい。 |] description3 = [markdown| N=0のとき、衝突回数は最終的に3になったはずです。<br /> 「最終的」といっても計算の打ち切り処理はしてませんので、<br /> 永久に衝突しなくなるだろう時点を適当に判断してください。<br /> さらにNを1,2..と変えてみると、以下の結果になるでしょう。<br /> <table border="1"> <tr><th>N</th><th>衝突回数</th></tr> <tr><td>0</td><td>3</td></tr> <tr><td>1</td><td>31</td></tr> <tr><td>2</td><td>314</td></tr> <tr><td>3</td><td>3141</td></tr> <tr><td>4</td><td>31415</td></tr> </table> 注意深い読者は気づいたでしょうが、この回数が円周率に対応します。 <table border="1"> <tr><th>N</th><th>衝突回数c</th><th>c/10^N</th></tr> <tr><td>0</td><td>3</td><td>3.0</td></tr> <tr><td>1</td><td>31</td><td>3.1</td></tr> <tr><td>2</td><td>314</td><td>3.14</td></tr> <tr><td>3</td><td>3141</td><td>3.141</td></tr> <tr><td>4</td><td>31415</td><td>3.1415</td></tr> </table> Nを増やせば増やすほど、精度が上っていきます。 ### 留意点など - 初速は結果には関係ありません。 - 質点と壁の具体的な初期位置は結果には関係ありません(M1,M2,壁の順序で並んでい<br />て、M1の初速が右向きである必要はあります) - 衝突による速度の変化だけが結果を決めます。 - 正しい表示のためには、一定の離散時間でプロットするのではなく、時間精度を適<br />宜細かくとかしていく必要がありますが、このシミュレーションでは時間間隔一定<br />でプロットしています。動きが変なのは、そのせいです。 ### 参考その他 この記事はElmを使って書いています。この記事を紹介している記事は[こちら](http://uehaj.hatenablog.com/entry/2014/08/03/234120)。 以下を参考に(計算式は丸パクリ)させて頂きました。 - [「2つのボールをぶつけると円周率がわかる」らしいのでシミュレーションしてみた](http://wasan.hatenablog.com/entry/2014/04/10/073638) - [「2つのボールをぶつけると円周率がわかる」のをしつこく確かめてみた・・・解析的に](http://wasan.hatenablog.com/entry/2014/04/15/045611) |] bkcolor = rgb 200 200 256 inpN : Input Int inpN = input 0 -- Nを選択。 selectN : Element selectN = plainText "N=" `beside` dropDown inpN.handle [ ("0", 0) , ("1", 1) , ("2", 2) , ("3", 3) , ("4", 4) , ("5", 5)] inpRunning = input False startButton : Element startButton = button inpRunning.handle True "開始" stopButton : Element stopButton = button inpRunning.handle False "停止" -- シミュレーションを表示 simulation : Int -> Int -> State -> Element simulation w h state = layers [ collage w (h `div` 2) <| [move (-(toFloat w / 4), 0) <| filled bkcolor (rect ((toFloat w)/2) ((toFloat h)/2))] , flow down [ flow right [ collage w (h `div` 2) <| [ traced {defaultLine|width<-4} (segment (0, 200) (0, -200)) , move (min 0 state.x1, 0) (filled red <|circle 5) , move (min 0 state.x2, 0) (filled red <|circle 2) ] ] ] ] -- 画面を表示 main : Signal Element main=let disp w h state = spacer 10 10 `beside` flow down [ description1 , image 610 362 "fig1.png" , description2 , selectN , if state.stat == Running then stopButton else startButton , flow down [ "M1,M2の質量の比率(m1:m2)= 100^N:1 = "++show state.ratio++":1" |> plainText , "M1の位置="++show state.x1 |> plainText , "M2の位置="++show state.x2 |> plainText , "衝突回数:"++show state.count |> plainText ] , simulation w h state , description3 ] in disp <~ Window.width ~ constant 400 ~ currentState
補足
- Elmの次回リリースでは、Markdown interpolationというのができるようになるので、この手のはもっとみやすく書けるようになるでしょう。
- asパターン無いと思ってたらあった。キーワードasを使用します。上ではnextStateの引数「{stat, x1, x2, v1, v2, count, ratio} as state)」で使用。
- Signalを1つのEventストリームにマージして、Event->State->Stateという関数を作ってfoldpに与えて…という基本構造は、いろいろ考えても、おちつくところにおちつく。あんまりバリエーションが生じない気がする。そのための、ある種のフレームワークがいくつか提案されているようだが、今後調査してみたい。(→Playground, automaton)
- 古典的FRPでは、離散イベントを扱うEvent、連続的な変化を抽象化したBehaviorの2つで整理するようだが、Elmの採用するArrowizedFRPにおけるSignalは離散的であるという意味で古典的FRPのEventに対応している。SignalはEventのようにタイムスタンプを保持しているわけではないが、Time.timeStampでタイムスタンプを持ったSIgnalを生成することができる。Signal.sampleOnする先も離散的でよい(シグナルは最後の値1個を常に保持しているので値がとれないということはない)。Elmでは本質的に連続的に変化する値を扱うことはない。
関連エントリ
「Elmでやってみる」シリーズのまとめエントリ - uehaj's blog
- 作者: Miran Lipovača,田中英行,村主崇行
- 出版社/メーカー: オーム社
- 発売日: 2012/05/23
- メディア: 単行本(ソフトカバー)
- 購入: 25人 クリック: 580回
- この商品を含むブログ (67件) を見る
- 作者: Miran Lipovaca
- 出版社/メーカー: オーム社
- 発売日: 2012/09/21
- メディア: Kindle版
- 購入: 4人 クリック: 9回
- この商品を含むブログを見る
- 作者: Graham Hutton,山本和彦
- 出版社/メーカー: オーム社
- 発売日: 2009/11/11
- メディア: 単行本(ソフトカバー)
- 購入: 14人 クリック: 503回
- この商品を含むブログ (117件) を見る