So-net無料ブログ作成
計算物理学入門 ブログトップ
前の3件 | 次の3件

2次元の写像 〜 『計算物理学入門』読み(その6) [計算物理学入門]

以下の2本の連立差分方程式で表される動的な系を描画してみる。 simple_map_x.png simple_map_y.png

プログラムは以下の通り。stream-takeの第2引数の値を変えると、プロットする点の数を簡単に変更できる。

(use graph)
(use util.stream)

(graph-init 800 600 :title "Simple Map" :display 2 :x-min -6 :x-max 10 :y-min -4 :y-max 8)
(graph-grid 16 12 'gray)

(define (simple-map-stream x y)
  (stream-cons (cons x y)
	       (simple-map-stream (+ 1 (- y) (abs x))
				  x)))

(graph-plot-point-stream
 (stream-take (simple-map-stream -0.1 0) 10000)
 'red)

プロットした結果は以下の通り。 SimpleMap.png
なんとも、不思議な絵柄が現れてくる。


nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:学問

車輪の外周の1点の動きのプロット 〜 『計算物理学入門』読み(その5) [計算物理学入門]

車輪の回転に伴って、車輪の上の1点がどのような動きをするのかをプロットしてみる。
車輪上の点の運動は、以下の方程式で表される。
whell_equation.png
ここで、
t: 時間
r: 車輪の半径
ω: 車輪の回転の角速度
vcm: 車輪の中心の移動速度(r×ω)
となる。

プロットするための、プログラムは以下の通り。前回同様にプロットする点のリストを作ってそれをプロットしてもよいのだが、リストの長さが5000と長いリストになってしまうので(このぐらいでメモリの使用量が問題になる事はないが)、ストリームを使った場合を示す。

(use srfi-1)
(use graph)
(use util.stream)
(use math.const)

(graph-init 800 600 :title "Wheel" :display 2 :x-min -1 :x-max 50 :y-min -4 :y-max 4)
(graph-grid 51 8 'gray)

(define (wheel r omega t)
  (let ((vcm (* r omega)))
    (cons (+ (- (* r (cos (* omega t))))
             (* vcm t))
          (* r (sin (* omega t))))))

(define (wheel-stream t)
  (let* ((point (wheel 1 1 t))
         (x (car point))
         (y (cdr point)))
    (if (< x 50)
        (stream-cons point (wheel-stream (+ t 0.01)))
        stream-null)))

(define plot-wheel-stream
  (graph-plot-point-stream (wheel-stream 0) 'red))

プロットした結果は以下の通り。 Wheel.png

ストリームを受け取ってプロットできるようにした、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 srfi-1)
  (use util.stream)
  (use rfb)
  (export graph-init graph-clear graph-close
          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 graph-plot-point-stream))

(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))

(define (graph-plot-point-stream p-str c)
  (stream-for-each (lambda (point)
                     (graph-set-point (car point) (cdr point) c))
                   p-str))

(provide "graph")


nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:学問

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))))

プロットした結果は以下の通り。画像をクリックすると原寸大表示されます。 CoolingOfCoffee.png

グラフに罫線を引けると便利なので、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")


nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:学問
前の3件 | 次の3件 計算物理学入門 ブログトップ

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。

×

この広告は1年以上新しい記事の更新がないブログに表示されております。