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 に対しても +, * などの演算子が使用できるので書くのは簡単です。