2.6 コーヒーの冷却の問題 cool

今回はニュートンの冷却の法則に従う系のシミュレーション。　やっと物理っぽい系が出てきたよ?。

ニュートンの冷却の法則

ニュートンの冷却の法則とは「周囲と温度が高い物体は、周囲との温度差に比例して冷却していく」という法則です。　『計算物理学入門』では、部屋に置いた熱いコーヒーの冷却に適用しています。　この冷却の法則は微分方程式で書くと以下のようになります：

　　 $\frac{dT}{dt} = -r (T - T_s)$

ここで

状態変数
- ${ t }$ : 時間
- ${ T }$ : 温度
パラメータ
- ${ T_\textrm{room} }$ : 周囲の温度
- ${ r }$ : 冷却定数

です。　解析解に関してはこちらを参照。

Groops によるシミュレーション

では Groops でシミュレーションを行うコードを見ていきましょう。　特に必須ではありませんが、ここでは時間に関する数値を（double ではなく） BigDecimal として実装しています。　また、実験結果と比較できるようにデータが載せられてたので、そちらもプロットするようにします。　データ・ファイルはこちら：

data2_2.csv

1列目が時刻（分）、2列目がブラック・コーヒーの温度（℃）、3列目がクリーム入りコーヒーの温度（℃）のデータです。　ここではブラック・コーヒーのデータを使います。

@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.*

// 時間関連の数値は BigDecimal でいっちゃうヨ！
def tmax = 46.0
def dt = 0.1
double r = 0.05d

new SimulationBuilder().composite{
//new SimulationBuilder().simulation(type:'composite'){

    // 反復を10回繰り返すシミュレーションを・・・
    simulation{
        system(new CoolSystem(), r:r, T_room:17.0d, dt:dt)
        init(time:0.0, T:82.3d)

        count(total:10)
    }

    // ・・・この反復条件の間だけ繰り返す
    continueWhile { it.time <= tmax }

    console()
    png(fileName:'cool.png', overwrite:true){
        chart(title:'Cool', domainLabel:'Time', rangeLabel:'Temperature'){
            point(label:"Simulation (r = $r)", x:'time', y:'T')
            marker(label:'Experiment', fileName:'data2_2.csv', x:0, y:1)
                        // ファイルからデータを読み込んでプロット
        }
    }
}.simulate()

/** ニュートンの冷却の法則に従う物理系 */
class CoolSystem extends AnnotationObservableSystem{

    @State BigDecimal time
    @State double T

    @Parameter BigDecimal dt
    @Parameter double T_room
    @Parameter double r

    @Override
    void evolve(Map<String, Object> params) {
        // オイラー法の適用箇所
        double dT = -r * (T - T_room) * dt.doubleValue()    // dT = -r(T - T_room)dt
        T += dT
        time += dt
    }
}

cool.groovy
1回の反復毎にデータを出力しているとパフォーマンス劣化やデータ量の増大で効率が悪くなるので、ここでは10回に1回だけ必要なデータを出力するようにしています。　Groops では composite タイプの simulation ノード（もしくは単に composite ノードでも可）に simulation ノードをネストすることで、これを実現できます。
チャートに他ファイルからデータを読み込んでプロットしたい場合は、chart ノード化に marker ノード を書き、ファイル名などを指定します*1。　x, y 属性に指定している数値は、ファイル中の列番号です。

シミュレーション結果は以下のようになります：