２−６．常微分方程式のまとめと演習

Hamiltonianが

H = (p²+q²)/2        (a1)

で与えられる１次元調和振動子を数値的に解いてみよう。

運動方程式は

dq/dt - p = 0           (a2)

dp/dt + q = 0           (a3)

        と２つの連立１次常微分方程式で表される。

        それぞれの数値計算アルゴリズムを以下にまとめ、初期条件 [p₀, q₀] = [1, 0]、 dt=0.1 としたサンプルプログラムを与える。

• Euler法  [fortran, C]
  
  q_n+1 = q_n + dt p_n
  p_n+1 = p_n - dt q_n
  
  　
• Leap-Frog法  [fortran, C]
  
  q_n+1 = q_n-1 + 2 dt p_n
  p_n+1 = p_n-1 - 2 dt q_n
  
  　
• Runge-Kutta法
  
  (2次)  [fortran, C]
  q_n+1/2 = q_n + 1/2 dt p_n
  p_n+1/2 = p_n - 1/2 dt q_n
  q_n+1    = q_n + dt p_n+1/2
  p_n+1    = p_n - dt q_n+1/2
  
  (4次)  [fortran, C]
  q' = q_n + 1/2 dt p_n
  p' = p_n - 1/2 dt q_n
  q'' = q_n + 1/2 dt p'
  p'' = p_n - 1/2 dt q'
  q''' = q_n + dt p''
  p''' = p_n - dt q''
  q_n+1    = q_n + 1/6 dt (p_n + 2p' + 2p'' + p''')
  p_n+1    = p_n - 1/6 dt (q_n + 2q' + 2q'' + q''')
  
  　
• Predictor-Corrector法  [fortran, C]
  
  q'_n+1 = q_n + dt p_n
  p'_n+1 = p_n - dt q_n
  q_n+1    = q_n + 1/2 dt [ p'_n+1+ p_n ]
  p_n+1    = p_n - 1/2 dt [ p'_n+1+ p_n ]
  
  　
• Symplectic法
  
  （１次）  [fortran, C]
  q_n+1 = q_n + dt p_n
  p_n+1 = p_n - dt q_n+1

        《演習》

1. Euler法のサンプルプログラムを参考にして、Leap-Frog法、Runge-Kutta法、Predictor-Corrector法、Symplectic法 で調和振動子の時間発展を求めよ。それぞれの方法のサンプルプログラムを参考にしてもよい。
   　
2. 各種方法による結果をgnuplotを使ってグラフ化し、比較考察せよ。
   
   gnuplot
   .....
   gnuplot> set data style lines
   gnuplot> plot 'fort.10' u 1:2           t - p(t) のプロット
   gnuplot> plot 'fort.10' u 1:3           t - q(t) のプロット
   gnuplot> plot 'fort.10' u 1:4           t - e(t) のプロット（eは全エネルギー＝保存すべき量であることに注意）
   gnuplot> plot 'fort.10' u 2:3           p(t) - q(t) のパラメトリックプロット（厳密解は(p²+q²)/2=e(0)の円であることに注意）
   　
3. Symplectic法による結果の特徴を考察せよ。影のハミルトニアン
   
   H' = H + dt/2 pq
   
   が（十分正確に）保存していることを確認せよ。
   
   　
4. ４次のRunge-Kutta法、及び１次のSymplectic法を用いて以下の非調和振動子の時間発展を求めよ。
   
   H = p²/2 + q⁴/4
   
   初期条件 [p₀, q₀] = [2ⁿ, 0], n=0,1,2,3,... として結果を考察せよ。
   　
5. （アドバンストな自習問題）電場(E)磁場(B)中での電子（質量m、電荷-q）の運動を考える。適当な摩擦（係数g）を取り入れた運動方程式は
   
   dx/dt = v
   m dv/dt = qv x B - gv + qE
   
   で与えられる。
   B = [0, 0, 1], E = [cos(t), 0, 0], m = q = g = 1
   初期条件 x₀ = [0,0,0], v₀ = [v₀, 0, 0], [0, v₀, 0] として結果を考察せよ。