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

倭算数理研究所

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

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

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

ニュートンの冷却の法則

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

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

ここで

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

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

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

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

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 属性に指定している数値は、ファイル中の列番号です。

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


今回は冷却定数の値を r = 0.05 にしましたが、実験とはあまりあってませんね。

計算物理学入門

計算物理学入門

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

プログラミングGROOVY

*1:サポートしているファイルタイプは、現在の所 CSV, SSV, TSV の3つです。 区切り文字はそれぞれコンマ、スペース、タブです。