ニュートンの冷却の法則をSchemeで実装 〜 『計算物理学入門』読み(その2)

ニュートンの冷却の法則のシミュレーションのプログラムをSchemeで実装してみた。(この法則の詳細はWikidediaの記事を参照) 本では、コーヒーの冷却を例として挙げている。

解析的に解ける問題なので、数値計算で求める必要はないのだが、練習と言う事で、本に載せているみたい。(計算のアルゴリズムの精度や、安定性の評価は、解析解と数値解を比較してみて、評価する事が多い。)

プログラムは以下の通り。

;; Euler

;; オイラー法を使って、初期条件と導関数dy/dx=f(x,y)からxの時のyの値を求める。

;; x-ini: xの初期値

;; y-ini: yの初期値

;; x:

;; fxy: f(x,y)

;; dx: x_nからx_n+1までの区間の幅

(define (euler x-ini y-ini x fxy dx)

    (let loop ((x-i x-ini)

               (y-i y-ini))

      (if (< x-i x)

          (loop (+ x-i dx) (+ y-i (* (fxy x-i y-i) dx)))

          y-i)))
;; コーヒーの冷却問題

;; r: 冷却定数

;; time: 時間

;; T-init: コーヒーの温度

;; T-room: 室温

;; dt: 時間の区間幅

;; この場合rは、0.026ぐらいがちょうどよい。

(define (cool r time T-init T-room dt)

  (euler 0.0 T-init time (lambda (t T) (- (* r (- T T-room)))) dt))

ちなみに、解析的に求める場合は以下の通り。

;;; 解析解

(define (cool-analytical r time T-init T-room)

  (- T-room (* (- T-room T-init) (exp (- (* r time))))))