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

微分方程式が解けない [計算物理学入門]

原子核の崩壊のシミュレーションで核種が2種類ある場合、解析的に値を求めるには連立微分方程式を解く必要がある。でも、その微分方程式が解けなかった。微分法定式なんて、かれこれ10年以上解いていないから解き方をすっかり忘れてしまっていた。数学は、高校の途中から苦手意識を持ち始めていたから、忘れるのも無理はないかも。それでも、やっぱりショックだった。

出先で本を読みながら、「どうやって解くんだっけ?」となり、数学の教科書は家に置いてあるので、解き方をインターネットで調べた。(便利な世の中になった)

たどりついた先が、化学反応の計算のサイトだった。大学大学院は、化学専攻だったのに。大学の1年の時習っていたはずだった。さらにショック。


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

原子核の崩壊のシミュレーション 〜 『計算物理学入門』読み(その8) [計算物理学入門]

放射性原子核の崩壊の基本法則は以下のように書ける。

AtomDecayEquation.png

ここで、λは崩壊定数である。
この方程式の形は、ニュートンの冷却の法則のシミュレーションとほぼ同じ形をしている。 それ故、放射性崩壊のシミュレーションのプログラムは以下の通りの短いものになる。
(define (decay Lamda time N dt)
  (euler 0.0 N time (lambda (t N) (- (* Lambda N))) dt))

λを直接指定するよりも、半減期を指定できた方が便利である。半減期は以下で与えられる。

HalfLifeEquation.png
半減期を指定できるようにしたプログラムは以下の通り。単にラッパーの関数を作っただけ。
(define (decay-half-reduce T_1/2 time N dt)
  (decay (/ 2 (exp T_1/2)) time N dt))

オイラー法は1次の精度しかないので、2次の精度で計算してみる。2次の方程式は以下の通り。

SecondPrecisionEquation.png

2次の精度で求めるプログラムは以下の通り。

(define (decay-2nd-precision Lambda time N dt)
  (let loop ((x-i 0)
             (y-i N))
    (if (< x-i time)
        (loop (+ x-i dt) (- y-i (* (/ Lambda 2)
                                   (- (* 2 y-i) (* Lambda y-i dt))
                                   dt)))
        y-i)))

様々なdtで、解析解、オイラー法と比較してみると

gosh> (* 100 (exp (- (* 1 1))))
36.787944117144235
gosh> (map (lambda (dt)
	     (decay 1 1 100 dt))
	   '(0.1 0.01 0.001 0.0001 0.00001 0.000001 0.0000001))
(31.381059608999998
 36.603234127322985
 36.76954247709651
 36.78242603282879
 36.78739229905556
 36.78792572316515
 36.787938598960544)
gosh> (map (lambda (dt)
	     (decay-2nd-precision 1 1 100 dt))
	   '(0.1 0.01 0.001 0.0001 0.00001 0.000001 0.0000001))
(33.352959127436435
 36.78856187161925
 36.78795025306916
 36.78426556798378
 36.7875762401565
 36.78794411715042
 36.78794043835247)
あれ? なんか値が解析解に収束していかないぞ。誤差が溜まってきてしまうのかなあ。まあ、途中までの結果をみると、2次の精度で求められているようだ。


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

リサージュ図形のプロット 〜 『計算物理学入門』読み(その7) [計算物理学入門]

以下の方程式で表される図形をプロットしてみる。

LissajousCurveEquation.png
プログラムは以下の通り。
(use graph)
(use util.stream)
(use math.const)

(graph-init 400 400 :title "Lissajous curve" :display 2 :x-min -2 :x-max 2 :y-min -2 :y-max 2)
(graph-grid 4 4 'gray)

(define (lissajous A omega-x phi-x
                   B omega-y phi-y
                   t)
  (cons (* A (sin (+ (* omega-x t) phi-x)))
        (* B (sin (+ (* omega-y t) phi-y)))))

(define (lissajous-curve t dt)
  (stream-cons (lissajous 1 2 (/ pi 6)
                          1 3 (/ pi 4)
                          t)
               (lissajous-curve (+ t dt))))

(graph-plot-point-stream
 (stream-take (lissajous-curve 0 0.01) 350)
 'red)

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

LissajousCurve.png
綺麗だなあ。


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

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

×

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