Newtonの冷却のシミュレーションのグラフ化 〜 『計算物理学入門』読み(その4)
Newtonの冷却のシミュレーションについてグラフ化してみる。(graphモジュールは、後述のように拡張している)
普通の手続き型の言語を使う場合は、1つのループの内部で、値の計算とプロットの両方が行われるように、コーディングされる事が多い。それに対して、Lisp系の言語の場合、結果の外部への表示のようなプログラムの本質に関係ない事は、最後におこなうようにプログラミングする事が多い。
下のプログラムでは、fold関数の部分でシミュレーション結果のリストを作って、その結果をfor-eachに渡してプロットしている。このようにすると、for-eachの第1引数を変えるだけで、色々な事ができる。(例えば、結果をファイルに出力するなど) ただし、このやり方だと、メモリを多量に消費してしまう。プロットしたい点の対のリストを受け取ってグラフにプロットする関数があると便利なので、それをgraphモジュールに追加する。
(use srfi-1)
(use graph)
(graph-init 800 600 :title "Cooling of coffee" :display 2 :x-min 0 :x-max 60 :y-min 0 :y-max 100)
(graph-grid 6 10 'gray);; 冷却のシミュレーションの結果をグラフ化
(for-each (lambda (point)
(let ((x (car point))
(y (cdr point)))
(graph-set-point x y 'green)))
(reverse (fold (lambda (time result-lst)
(let ((time-ini (caar result-lst))
(Temp-ini (cdar result-lst)))
(cons (cons time
(euler time-ini Temp-ini time
(lambda (time Temp) (- (* 0.026 (- Temp 17))))
0.001))
result-lst)))
'((0 . 82.3))
(iota 23 2 2)))
'green);; 実測値のグラフ化
(for-each (lambda (x y)
(graph-set-point x y 'cyan))
(iota 24 0 2)
'(82.3 78.5 74.3 70.7 67.6
65.0 62.5 60.1 58.1 56.1
54.3 52.8 51.2 49.9 48.6
47.2 46.1 45.0 43.9 43.0
41.9 41.0 40.1 39.5));; graphモジュールの機能を使用して、シミュレーション結果をプロット
(graph-plot-points
(reverse (fold (lambda (time result-lst)
(let ((time-ini (caar result-lst))
(Temp-ini (cdar result-lst)))
(cons (cons time
(euler time-ini Temp-ini time
(lambda (time Temp) (- (* 0.026 (- Temp 17))))
0.001))
result-lst)))
'((0 . 82.3))
(iota 23 2 2))))
プロットした結果は以下の通り。画像をクリックすると原寸大表示されます。
グラフに罫線を引けると便利なので、graphモジュールに機能を追加した。
;; -*- coding: utf-8; mode: scheme -*-
;;
;; graph.scm - Drawing graph module with Gauche-rfb
;;
;; Copyright (c) 2008 Ettem
;; All rights reserved.
;;(define-module graph
(use rfb)
(use srfi-1)
(export graph-init graph-clear
graph-set-scale!
graph-set-x-min! graph-set-x-max! graph-set-y-min! graph-set-y-max!
graph-set-point graph-line graph-box graph-grid graph-circle
graph-plot-fx graph-plot-points))(select-module graph)
;; グラフのインスタンス
(define *graph* #f);; グラフの表示範囲
(define-class()
((x-min :init-keyword :x-min
:getter get-x-min
:setter set-x-min!)
(x-max :init-keyword :x-max
:getter get-x-max
:setter set-x-max!)
(y-min :init-keyword :y-min
:getter get-y-min
:setter set-y-min!)
(y-max :init-keyword :y-max
:getter get-y-max
:setter set-y-max!)
(window-x :init-keyword :window-x
:getter get-window-x
:setter set-window-x!)
(window-y :init-keyword :window-y
:getter get-window-y
:setter set-window-y!)
(inset-ratio :init-keyword :inset-ratio
:init-value 0.05
:getter get-inset-ratio
:setter set-inset-ratio)))(define-method get-x-len ((graph
))
(- (get-x-max graph) (get-x-min graph)))(define-method get-y-len ((graph
))
(- (get-y-max graph) (get-y-min graph)))(define-method get-x-pixel ((graph
))
(round->exact
(* (get-window-x graph)
(- 1 (* 2 (get-inset-ratio graph))))))(define-method get-y-pixel ((graph
))
(round->exact
(* (get-window-y graph)
(- 1 (* 2 (get-inset-ratio graph))))))(define-method get-inset-left ((graph
))
(round->exact (* (get-window-x graph) (get-inset-ratio graph))))(define-method get-inset-top ((graph
))
(round->exact (* (get-window-y graph) (get-inset-ratio graph))));; グラフの値から、画面上のピクセルへの変換
(define (x->x-pixel x)
(+ (get-inset-left *graph*)
(round->exact
(* (/ (abs (- x (get-x-min *graph*)))
(get-x-len *graph*))
(get-x-pixel *graph*)))))(define (y->y-pixel y)
(- (get-window-y *graph*)
(get-inset-top *graph*)
(round->exact
(* (/ (abs (- y (get-y-min *graph*)))
(get-y-len *graph*))
(get-y-pixel *graph*)))));;; API
;; グラフの表示
;; win-x: ウィンドウの幅
(define (graph-init win-x win-y . restarg)
(let-keywords restarg ((title #f) (display 0) (port #f)
(x-min 0) (x-max win-x)
(y-min 0) (y-max win-y))
(set! *graph* (make
:x-min x-min :x-max x-max
:y-min y-min :y-max y-max
:window-x win-x :window-y win-y))
(let ((rfb-arg '()))
(if title
(set! rfb-arg (append `(:title ,title) rfb-arg)))
(if port
(set! rfb-arg (append `(:port ,port) rfb-arg)))
(apply rfb-init win-x win-y :display display rfb-arg))));; グラフウィンドウを閉じる
(define graph-close
rfb-close);; グラフの消去
(define (graph-clear c)
(rfb-clear c));; グラフの表示範囲
(define (graph-set-scale! x-min x-max y-min y-max)
(set-x-min! *graph* x-min)
(set-x-max! *graph* x-max)
(set-y-min! *graph* y-min)
(set-y-max! *graph* y-max))(define (graph-set-x-min! x)
(set-x-min! *graph* x))(define (graph-set-x-max! x)
(set-x-max! *graph* x))(define (graph-set-y-min! y)
(set-y-min! *graph* y))(define (graph-set-y-max! y)
(set-y-max! *graph* y))
;; 点の描画
(define (graph-set-point x y c)
(rfb-set-pixel (x->x-pixel x)
(y->y-pixel y)
c));; 線の描画
(define (graph-line x1 y1 x2 y2 c)
(rfb-line (x->x-pixel x1)
(y->y-pixel y1)
(x->x-pixel x2)
(y->y-pixel y2)
c));; 矩形描画
(define (graph-box x1 y1 x2 y2 c . restarg)
(let-keywords restarg ((filled? #f))
(rfb-box (x->x-pixel x1)
(y->y-pixel y1)
(x->x-pixel x2)
(y->y-pixel y2)
c
:filled? filled?)));; グリッド線の描画
(define (graph-grid n-x n-y c)
(for-each (lambda (x)
(rfb-line (x->x-pixel x)
(y->y-pixel (get-y-min *graph*))
(x->x-pixel x)
(y->y-pixel (get-y-max *graph*))
c))
(iota (+ n-x 1) (get-x-min *graph*) (/ (get-x-len *graph*) n-x)))
(for-each (lambda (y)
(rfb-line (x->x-pixel (get-x-min *graph*))
(y->y-pixel y)
(x->x-pixel (get-x-max *graph*))
(y->y-pixel y)
c))
(iota (+ n-y 1) (get-y-min *graph*) (/ (get-y-len *graph*) n-y))));; 円の描画
;; 目盛にあった、円は描けない場合がある...
;; Windowのアスペクト比と、目盛のアスペクト比は一致させる必要がある。
(define (graph-circle x y r c)
(rfb-circle (x->x-pixel x)
(y->y-pixel y)
(x->x-pixel r)
c));; y=f(x)の描画
(define (graph-plot-fx f c)
(let* ((x-min (get-x-min *graph*))
(x-max (get-x-max *graph*))
(x-pixel (* (get-window-x *graph*)
(- 1 (* 2 (get-inset-ratio *graph*)))))
(delta-x (/ (- x-max x-min) x-pixel)))
(for-each (lambda (x)
(graph-set-point x (f x) c))
(iota (+ x-pixel 1) x-min delta-x))));; 点のリストのプロット
(define (graph-plot-points p-lst c)
(for-each (lambda (point)
(graph-set-point (car point) (cdr point) c))
p-lst))(provide "graph")
