2−6.常微分方程式のまとめと演習
Hamiltonianが
H = (p2+q2)/2 (a1)
で与えられる1次元調和振動子を数値的に解いてみよう。
運動方程式は
dq/dt - p = 0 (a2)
dp/dt + q = 0 (a3)
と2つの連立1次常微分方程式で表される。
それぞれの数値計算アルゴリズムを以下にまとめ、初期条件 [p0, q0] = [1, 0]、 dt=0.1 としたサンプルプログラムを与える。
《演習》
- Euler法のサンプルプログラムをコピーし、コンパイル、実行、結果のグラフ表示を行う。具体的には、Linux端末で以下のコマンド を実行する。
gfortran -o eu1 eu1.f
./eu1
(このFrotranのサンプルプログラムでは自動的に"fort.10"という名前のファイルを作成するようになっている。Cのサンプルプログラム では同様のデータを端末上に出力するので、それをコピー&ペーストしてemacsなどで"fort.10"なるデータを別途作成する。)
gnuplot
gnuplotのプロンプト"gnuplot>"が表示されるので、以下を入力。
.....
gnuplot> plot 'fort.10' u 1:2 w line t - p(t) のプロット
gnuplot> plot 'fort.10' u 1:3 w line t - q(t) のプロット
gnuplot> plot 'fort.10' u 1:4 w line t - e(t) のプロット(eは全エネルギー=保存すべき量であることに注意)
gnuplot> plot 'fort.10' u 2:3 w line p(t) - q(t) のパラメトリックプロット(厳密解は(p2+q2)/2=e(0)の円である ことに注意)
- Euler法のサンプルプログラムを参考にして、Leap-Frog法で調和振動子の時間発展(p,qの時間変化(t依存性))を求めよ。それぞ れの方 法のサンプルプログラムを参考にしてもよいが、必ずプログラムソースを自分で読んで理解してから使うこと。 結果をgnuplotを使ってグラフ化し、オイラー法の結果と比較考察せよ。
- Symplectic法による結果の特徴を考察せよ。影のハミルトニアン
H' = H + dt/2 pq
が(十分正確に)保存していることを確認せよ。
- 4次のRunge-Kutta法[解法のヒン ト]、
及び1次のSymplectic法を用いて以下の非調和振動子の時間発展を求めよ。
H = p2/2 + q4/4
初期条件 [p0, q0] = [2n, 0], n=0,1,2,3,... として結果を考察せよ。
- (アドバンストな自習問題)電場(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
初期条件 x0 = [0,0,0], v0 = [v0, 0, 0], [0, v0, 0] として結果を考察せよ。