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

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] として結果を考察せよ。

        《レポート》

1. 演習４のプログラムと自分なりに考察したことをメールで以下のアドレスに送ること。
   
   提出先：　ryoichi@r02.mbox.media.kyoto-u.ac.jp
   締切り：　１１月９日（火）
   
   氏名、学籍番号を忘れずに。授業に関する要望・感想があれば遠慮せずに書いてください。

レポート提出者（出した人は番号があるか確認してください）：

提出者：

0500-16-3529

0500-15-0173

0500-15-0389

0500-15-0404

0500-15-0647

0500-15-0656

0500-15-0638

0500-15-0807

0500-15-0852

0500-15-1279

0500-15-1528

0500-15-1537

0500-15-1804

0500-15-1921

0500-15-1958

0500-15-1976

0500-15-2033

0500-15-2230

0500-15-2697

0500-15-2991

0500-15-3077

0500-15-3246

0500-15-3255

0500-15-3308

0500-15-3596

0500-15-3676

0500-15-3845

0500-15-3999

0500-15-4191

0500-15-4252

0500-15-4271

0500-15-4539

0500-15-4557

0500-15-4978

0500-15-5142

0500-15-5385

0500-15-6211

0500-15-6220

0500-15-6337

0500-15-6892

0500-15-6991

0500-15-7236

0500-15-7281

0500-15-7451

0500-15-7782

0500-15-8091

0500-15-8233

0500-15-8289

0500-15-8313

0500-15-8725

0500-15-8805

0500-15-8823

0500-15-9179

0500-15-9419

0500-15-9464

0500-15-9473

0500-15-9571

0500-15-9581

0500-15-9689

0500-15-9894

0500-14-0721

0500-14-1030

0500-14-1746

0500-14-1245

0500-14-2350

0500-14-2476

0500-14-2529

0500-14-2645

0500-14-2734

0500-14-2841

0500-13-0420

0500-13-0466

0500-13-1697

0500-13-2326

0500-13-2952

0500-11-2154

0500-11-3071

6100-15-2893

6100-14-1218

[2004.11.10, 22:30]

＊受け取りの確認の意味で番号（なければ名前）を載せています。掲載を希望しない人があれば対応しますので、レポート毎にメールに希望を記載してください。