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

倭算数理研究所

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

2.8 簡単なグラフ表示 : 問題 2.4 冷却時間

今回はグラフ表示と緩和時間の計算です。 グラフ表示は今まで既に行っていますが、まぁ一応やっておきます。 緩和時間は放射性崩壊などの半減期と似た量です。

問題

a. プログラム cool を修正して、コーヒーの温度についての数値解を、数値として表示するのではなくグラフで示せ。

b. コーヒーの温度が、初期の温度と室温との差の { 1/e \sim 0.37 } に達するまでの時間を求めよ。 そのような時間を緩和時間という。 緩和時間は初期温度や室温に依存するか。 異なる { r } の値について計算を試み、緩和時間の { r } に対する定性的な依存性を求めよ。

a.

プログラム cool のグラフ表示は既にこちらで行っていますが、ここでももう一度やっておきましょう。 グラフを PNG ファイルとして生成するには png ノードとして書きます:

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

double r = 0.03d, T_room = 17.0d, T0 = 82.3d

new SimulationBuilder().simulation{

    system(new CoolSystem(), r:r, T_room:T_room, dt:1.0)
    init(time:0.0, T:T0)
    continueWhile { it.time <= 46.0 }

    console()
    png(fileName:'problem2_4a.png', overwrite:true){
        chart(title:'Problem 2.4 a', domainLabel:'Time', rangeLabel:'Temperature'){
            point(label:"Simulation", x:'time', y:'T')
            marker(label:'Analytic Solution', function:{ t -> T_room - (T_room - T0)*Math.exp(-r*t) })
        }
    }
}.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) {
        T += -r * (T - T_room) * dt.doubleValue()
        time += dt
    }
}

生成されるグラフは以下のようになります:

b.

次は緩和時間 (relaxation time) の計算。 まずは初期温度、室温、冷却定数を固定して行うシミュレーション・スクリプト

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

double r = 0.03d, T_room = 17.0d, T0 = 82.3d

double T_end = T_room + (T0 - T_room) / Math.E
println "r = $r, T_end = $T_end"

new SimulationBuilder().simulation{

    system(new CoolSystem(), r:r, T_room:T_room, dt:0.1)
    init(time:0.0, T:T0)
    continueWhile { it.T > T_end }

    console(dataEntries:['time', 'T'], outputTiming:OBSERVE_FINAL)
}.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) {
        T += -r * (T - T_room) * dt.doubleValue()
        time += dt
    }
}

緩和時間の精度は小数第1位までです(dt = 0.1 のため)。 このスクリプトでの実行結果より、緩和時間は 33.2 (min) でした。 ニュートンの冷却の法則を表す微分方程式の解析解は

  { \displaystyle
\begin{align*}
    T(t) = T_s - (T_s - T_0) e^{-rt}
\end{align*}
}

となるので(こちら参照)、緩和時間は冷却定数の逆数 { 1/r } になるはずです。 今のシミュレーションでは { r = 0.03 } としていたので、{ 1/r \sim 33.3 } となり、シミュレーション結果とだいたい一致しますね。

いろいろな室温 T_room, 初期温度 T0 で緩和時間を求めるシミュレーション・スクリプトは以下のようになります:

@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.*
import org.waman.groops.simulation.system.*

double r = 0.03d

new SimulationBuilder().parallel{
    for(T_room in (15..20)){
        for(T0 in (80..85)){
            simulation{
                system(new CoolSystem(), r:r, T_room:T_room, dt:0.1,
                       T_end:T_room + (T0 - T_room) / Math.E)
                init(time:0.0, T:T0)
                parameterObservableSet()
                console(format:"T_room : %s, T0 : $T0, Relaxation Time : %s", dataEntries:['T_room', 'time'],
                        outputTiming:FINAL, outputHeader:false)
            }
        }
    }
}.simulate()

class CoolSystem extends AnnotationControllableObservableSystem{

    @State BigDecimal time
    @State double T

    @Parameter BigDecimal dt
    @Parameter double T_room, r, T_end

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

    @Override
    boolean continueIteration(DataStore store){
        return T > T_end
    }
}
  • problem2_4b2.groovy 直
  • 物理系を少し修正して、反復を停止する温度 T_end ( = T_room + (T0 - T_room) / e ) もシミュレーション・パラメータとして付け加えています。 また、反復の停止を物理系に判断させるために、物理系を AnnotationControllableObservableSystem のサブクラスとして作成しています。
  • 標準出力にデータを出力をする console() ノードに format 属性を指定して、出力する文字列のフォーマットを設定しています。 「%s」 には dataEntries 属性に指定された観測量の観測値が順に渡されます。

シミュレーション結果から、室温、初期温度を変えても緩和時間は変わらないことが分かります。 また、冷却定数を変えて同様のシミュレーションを行うと、冷却定数を { r } として緩和時間が { 1/r } であることも分かります。

計算物理学入門

計算物理学入門

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

プログラミングGROOVY