2.6 コーヒーの冷却の問題：問題 2.1 コーヒーの冷却問題のためのプログラム

今回は、前回の「ニュートンの冷却の法則」に従う系のシミュレーションに関する補足問題です。

問題

a. 電卓を使ってオイラー法でニュートンの冷却の法則の数値解を計算せよ。
b. 実験データに合う r の近似値を探せ。
c. 設問 b で見つけた r の値を使って、温度の時間変化をグラフに示せ。
d. シミュレーションで使った時間幅の値が結果に影響しないほど小さいことを確かめよ。
e. r の値はブラック・コーヒーとクリーム入りコーヒーでは同じか。
f. ニュートンの冷却の法則はコーヒーの冷却に適用できるか。
g. 温度差が 1/2, 1/4, 1/8 になるまでコーヒーが冷めるにはどれほどの時間がかかるか。

a.

誰がやるか。

b, c.

前回のシミュレーションを、いろいろな r （冷却定数）の値について実行してみましょう。　シミュレーションタイプは並列シミュレーションです（こちらを参照）。

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

def rList = [0.01d, 0.02d, 0.03d, 0.04d, 0.05d]
def tmax = 46.0
def dt = 0.1

new SimulationBuilder().parallel{

    def pngout = png(fileName:'problem2_1b.png', overwrite:true){
        chart(title:'problem 2.1 b', domainLabel:'Time', rangeLabel:'Temperature', range:[30d, 85d]){
            // 各 r についての系列を定義
            for(r in rList){
                point(label:"r = $r", x:'time', y:'T')
            }
            marker(label:'Experiment', fileName:'data2_2.csv', x:0, y:1)
        }
    }

    // 各 r についてのシミュレーションを定義
    for(r in rList){
        composite{
            simulation{
                system(new CoolSystem(), r:r, T_room:17d, dt:dt)
                init(time:0.0, T:82.3d)

                count(total:10)
            }
            continueWhile { it.time <= tmax }
            dataOutputter(pngout["r = $r"])
        }
    }

}.simulate()

class CoolSystem extends AnnotationObservableSystem{

    @State BigDecimal time
    @State double T
    @Parameter BigDecimal dt
    @Parameter double r, T_room

    @Override
    void evolve(Map<String, Object> params) {
        double change = -r * (T - T_room) * dt.doubleValue()
        T += change
        time += dt
    }
}

problem2_1b.groovy
CoolSystem クラスは前回と全く同じです。
シミュレーションをしたい r の値を List にして、各値に対応する point ノードと simulation ノードを for 文で作成しています。

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

実験データは 0.02 ≦ r ≦ 0.03 の間にあるので、この間をもう少し細かくやってみると（シミュレーション結果は線にしてます）：

となって、時間全体にわたって実験データとシミュレーション結果が一致する r の値はなさそうです。　キチンと評価するには最小2乗法とかやればいいのかも知れませんが、これは第7章5節でやるようなのでここではしません。　そもそも、初期温度を固定してるので最小2乗法をやってうまく評価できるのかもよく分かんないし、r = 0.03 にするとある程度の温度まではそれなりに合っていて、60℃くらいからズレが大きくなっていく感じなので、ニュートンの冷却の法則によるモデリングがよくないんでしょう。

まぁ逆に言えば、50〜60℃くらいまでは r = 0.03 でそれなりに合っているとしていいでしょう。　問題2.2 で75℃くらいまででシミュレーションさせる問題があるので、著者の望んでる答えはそんなところでしょう（笑）　少し気になるのは、データがズレ出して以降の挙動は、r の値を変えればまたそれなりにニュートンの冷却の法則に従うのか？それとも関数形自体が（たとえば冪関数などに）変わってしまうのか？というところですが、シミュレーションする気が起こらないのでこの辺で・・・

d.

dt と t = 46.0 までシミュレーションを行った結果の温度 T とのグラフは以下のようになりました：

例えば dt = 0.01 （上記の r を求めるのに使った値）と dt = 0.005 や dt = 0.02 での値は、温度の変化量に比して1%もないので、ここで使った dt の値は結果に影響しないとしていいでしょう。　このあたりの話は、問題2.3でまた詳しくやるようなので、ここではこれ以上評価しません。

次回へ続く。