読者です 読者をやめる 読者になる 読者になる

倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

2.5 オイラー法のプログラム example

計算物理学入門 シミュレーション

さて、プログラミング演習がソコソコこなせてきたので、だんだん物理シミュレーションの方に入っていきましょう。 微分方程式に従う物理系のシミュレーションを行うオイラー以前の記事で説明しました。 今回はそのオイラー法を簡単な微分方程式に適用してみます。

微分方程式

今回扱う微分方程式とその初期条件は

  { \displaystyle
\begin{align*}
    \frac{dy}{dx} = 2x, \qquad y(1) = 1
\end{align*}
}

です。 ちなみに解析解は

  { \displaystyle
\begin{align*}
    y = x^2
\end{align*}
}

ですね。

物理系を表すクラスを実装する場合

まずは、物理系を表すクラスを定義した方法。 オイラー法のアルゴリズムを適用するのは PhysicalSystem のサブクラスに定義する evolve() メソッドです。

@GrabResolver('https://github.com/waman/groops-core/tree/master/repo')
@Grab('org.waman.groops:groops-core:0.1-beta')

import org.waman.groops.builder.SimulationBuilder
import org.waman.groops.simulation.system.*

new SimulationBuilder().simulation{

    system(new ExampleSystem(), dx:0.1)
    systemObservableSet()    // @State 変数を物理量として登録
    init(x:1.0, y:1d)    // 初期条件

    continueWhile{ it.x <= 2.0 }    // x が2以下の間、反復
    console(dataEntries:['x', 'y'])    // x, y を出力
}.simulate()

class ExampleSystem extends AnnotationPhysicalSystem{

    @State BigDecimal x
    @State double y
    @Parameter BigDecimal dx

    @Override
    void evolve(Map params){
        // オイラー法による状態更新
        double slope = 2d*x.doubleValue()    // dy/dx = 2x
        double change = slope*dx
        y += change
        x += dx
    }
}
  • example1.groovy 直
  • slope 変数に代入している部分が、微分方程式の具体的な形から計算する部分です。
  • change 変数に代入している部分は、微分方程式の詳細によらない、オイラー法の部分です(もちろん1階の微分方程式という前提で)。
  • 最後に @State 変数(状態を指定する変数) x, y の値を更新してます。

といった感じです。 特に難しくないですね。

ノードだけで構築した場合

今回の物理系もあまり「物理系」といった感じのものではないので、物理系を表すクラスを定義せずに、ノードだけで構築してみましょう。 ついでに結果をチャートに出力してみます。

@GrabResolver('https://github.com/waman/groops-core/tree/master/repo')
@Grab('org.waman.groops:groops-core:0.1-beta')

import org.waman.groops.builder.SimulationBuilder
import java.awt.Color

new SimulationBuilder().simulation{

    observableSystem{
        states = [x:1.0, y:1d]    // x は BigDecimal のインスタンスにしてます。
        parameters = [dx:0.1]          // dx も BigDecimal。
        evolve = { system, Map params ->
            system.with{
                // オイラー法による状態更新
                y += 2d*x.doubleValue()*dx.doubleValue()    // x, dx が BigDecimal に注意。
                x += dx
            }
        }
    }

    continueWhile{ it.x <= 2.0 }    // x が BigDecimal なので 2.0 も。
    console(dataEntries:['x', 'y'])
    // 結果をチャートに出力
    png(fileName:'example.png', overwrite:true){
        chart(title:'Example', domain:[1d, 2d]){    // 定義域の指定
            point(label:'Simulation', x:'x', y:'y')
            marker(label:'Analytic', function:{ double x -> x*x }, color:Color.BLUE)
        }
    }
}.simulate()
  • example2.groovy 直
  • オイラー法のアルゴリズムは observableSystem ノード下の evolve プロパティにセットしているクロージャに書かれています。
  • Groovy が Object クラスに追加している with() メソッドは、クロージャ内の(解決できない)プロパティやメソッドの処理を、with() を呼び出されているオブジェクトに委譲します。 具体的には、このクロージャ内でアクセスされている x, y, dx プロパティは system オブジェクトのものになります。
  • 1つ目のサンプルコードでは、x が double (java.lang.Double) の精度なので、dx をどんどん加えていくと少しずつ厳密な値からズレていきます。 もしこれが気になるなら、double (java.lang.Double) ではなく java.math.BigDecimal を使うことで解決できます。 Groovy では BigDecimal に対しても +, * などの演算子が使用できるので書くのは簡単です。

シミュレーション結果をチャートにすると以下のようになります:




計算物理学入門

計算物理学入門

  • 作者: ハーベイゴールド,ジャントボチニク,Harvey Gould,Jan Tobochnik,鈴木増雄,石川正勝,溜渕継博,宮島佐介
  • 出版社/メーカー: ピアソンエデュケーション
  • 発売日: 2000/12
  • メディア: 単行本
  • 購入: 1人 クリック: 28回
  • この商品を含むブログ (45件) を見る