数値解析法及び演習 第九回

常微分方程式
  1. Δtの二次の精度
  2. ルンゲクッタ法
  3. 連立1階微分方程式
  4. gnuplotでアニメ
  5. 今日の宿題

[1]常微分方程式その2 「Δtの二次の精度」

オイラー法における誤差はΔtに比例する(1ステップの誤差はΔtの二乗に比例するため)．すなわち誤差を10分の1にしたければΔtを1/10にする必要がある．これは計算にかかる時間が10倍になることを意味する．これでは時間がかかってしょうがない．そこでΔtの2乗，4乗に誤差が比例するようなアルゴリズムが存在する．4乗に比例するものはルンゲクッタ法と呼ばれ広く用いられている．

ルンゲクッタ法に入る前に誤差がΔtの2乗に比例するアルゴリズムについて考えてみよう．


\begin{displaymath} {dy\over dt}=f(t,y) \end{displaymath} (1)


という微分方程式が与えられたものとする．タイムステップがΔtのとき，n+1ステップ目のyの値yn+1は，ynを用いて

\begin{displaymath} y_{n+1}=y_n+\Delta t y_n'+{1\over 2}\Delta t^2 y_n'' \end{displaymath} (2)

と書くことができる．ここでy''は二次の微係数で


\begin{displaymath} y''={\partial f\over \partial t}+{\partial f\over \partial y}{dy\over dt} =f_t+f_yf \end{displaymath} (3)

となる．2項出てくることに注意．これを用いると先ほどの式は


\begin{displaymath} y_{n+1}=y_n+\Delta t f++{1\over 2}\Delta t^2 (f_t+f_yf) \end{displaymath} (4)

のようにかきかえられる．ここでポイントは，我々はft,fyを知らないと言うことである．なので別の式でさらに書き替えないといけない．そのため式(8)が，以下のように書けるものと仮定する．


\begin{displaymath} y_{n+1}=y_n+\lambda_1f_1+\lambda_2f_2 \end{displaymath} (5)

ここでf1，f2は


$\displaystyle f_1$ $\textstyle =$ $\displaystyle \Delta t f(t_n,y_n)$ (6)
$\displaystyle f_2$ $\textstyle =$ $\displaystyle \Delta t f(t_n+\alpha\Delta t,y_n+\beta f_1)$ (7)

である．f2の中身に注意．f1，f2はそれぞれtnでの傾き，および傾きk1でΔt秒進んだとしたときのt+Δt秒後の傾きを表す．この二つの傾きを適当に重ね合わせることによりΔtの二乗の精度にすることができる．(11)式をテーラー展開すると，λ1，λ2，α，βの間の関係式が


$\displaystyle \lambda_1+\lambda_2$ $\textstyle =$ $\displaystyle 1$ (8)
$\displaystyle \lambda_1\alpha$ $\textstyle =$ $\displaystyle {1\over 2}$ (9)
$\displaystyle \lambda_2\beta$ $\textstyle =$ $\displaystyle {1\over 2}$ (10)

と求めることができる．求めるべき数が4つあるのに式が３つしか存在しない．なのでいろいろな値の組み合わせが


\begin{displaymath} \lambda_1=\lambda_2={1\over 2}, \alpha=1, \beta=1 \end{displaymath} (11)


\begin{displaymath} \lambda_1={2\over 3}, \lambda_2={1\over 3}, \alpha=\beta={3\over 2} \end{displaymath} (12)

といったようにありうる．いずれの場合も二次の精度まで正しく計算できる．つまり一ステップの誤差はΔtの3乗に比例する. 今回は(11)式を採用し，以下の(13)-(15)式で積分を行う．


\begin{displaymath} y_{n+1}=y_n+{f_1+f_2\over 2} \end{displaymath} (13)


$\displaystyle f_1$ $\textstyle =$ $\displaystyle \Delta t f(t_n,y_n)$ (14)
$\displaystyle f_2$ $\textstyle =$ $\displaystyle \Delta t f(t_n+\Delta t,y_n+f_1)$ (15)

練習問題1

式(13),(14),(15)を用いて先週の宿題(dy/dt=-y, t=0でy=1）をt=2まで解け．オイラー法の場合と精度を比較せよ．なお，プログラムを作成するときは関数副プログラムを用い，メインプログラムの中ではf(t,x)と書き，f(t,x)の具体的な形を関数副プログラムで定義すること．

[2] ルンゲクッタ法

先ほどのアルゴリズムの誤差はΔtの二乗である．これを4乗に高めたのがルンゲクッタ法である．

さきほどと同じ考え方を使うと，ルンゲクッタ法のアルゴリズムは以下のようになる．


\begin{displaymath} y_{n+1}=y_n+{1\over 6}({f_1+2f_2+2f_3+f_4}) \end{displaymath} (16)


$\displaystyle f_1$ $\textstyle =$ $\displaystyle \Delta t f(t_n,y_n)$  
$\displaystyle f_2$ $\textstyle =$ $\displaystyle \Delta t f(t_n+{1\over 2}\Delta t,y_n+{1\over 2}f_1)$  
$\displaystyle f_3$ $\textstyle =$ $\displaystyle \Delta t f(t_n+{1\over 2}\Delta t,y_n+{1\over 2}f_2)$  
$\displaystyle f_4$ $\textstyle =$ $\displaystyle \Delta t f(t_n+\Delta t,y_n+f_3)$ (17)

練習問題2

式(16),(17)を用いて，ルンゲクッタ法を用いて　dy/dt=-y, t=0でy=1 をt=2まで解け．３つの場合の誤差を比較せよ．

[3] 連立微分方程式

先ほどのルンゲクッタ法は多変数に拡張できる．例えば，x=x(t)とy=y(t)がそれぞれ


$\displaystyle {dy\over dt}$ $\textstyle =$ $\displaystyle f(t,y,z)$  
$\displaystyle {dz\over dt}$ $\textstyle =$ $\displaystyle g(t,y,z)$ (18)

にしたがって進化する場合は，以下のようにアルゴリズムを組めばよい．


$\displaystyle y_{n+1}$ $\textstyle =$ $\displaystyle y_n+{1\over 6}(f_1+2f_2+2f_3+f_4)$  
$\displaystyle z_{n+1}$ $\textstyle =$ $\displaystyle z_n+{1\over 6}(g_1+2g_2+2g_3+g_4)$ (19)


$\displaystyle f_1$ $\textstyle =$ $\displaystyle \Delta t f(t_n,y_n,z_n)$  
$\displaystyle g_1$ $\textstyle =$ $\displaystyle \Delta t g(t_n,y_n,z_n)$ (20)


$\displaystyle f_2$ $\textstyle =$ $\displaystyle \Delta t f(t_n+{1\over 2}\Delta t,y_n +{1\over 2}f_1,z_n+{1\over 2}g_1)$  
$\displaystyle g_2$ $\textstyle =$ $\displaystyle \Delta t g(t_n+{1\over 2}\Delta t,y_n +{1\over 2}f_1,z_n+{1\over 2}g_1)$ (21)


$\displaystyle f_3$ $\textstyle =$ $\displaystyle \Delta t f(t_n+{1\over 2}\Delta t,y_n +{1\over 2}f_2,z_n+{1\over 2}g_2)$  
$\displaystyle g_3$ $\textstyle =$ $\displaystyle \Delta t g(t_n+{1\over 2}\Delta t,y_n +{1\over 2}f_2,z_n+{1\over 2}g_2)$ (22)


$\displaystyle f_4$ $\textstyle =$ $\displaystyle \Delta t f(t_n+\Delta t,y_n+f_3,z_n+g_3)$  
$\displaystyle g_4$ $\textstyle =$ $\displaystyle \Delta t g(t_n+\Delta t,y_n+f_3,z_n+g_3)$ (23)

練習問題3

質点の投げ上げの数値シミュレーションを行う．今回解く方程式は以下の連立微分方程式となる．


\begin{displaymath} {d^2x\over dt^2}=-g \end{displaymath} (24)


$\displaystyle {dx\over dt}$ $\textstyle =$ $\displaystyle y$  
$\displaystyle {dy\over dt}$ $\textstyle =$ $\displaystyle -g$ (25)

式(25)で与えられる連立微分方程式を解け．t=0でx=0.D0, y=1.D0 とすること．Δt=1.D-2でt=2となるまで積分せよ．xが初期位置，yが初速度に対応する．yを変化させて結果がどう変わるか以下のアニメで観察せよ．

[4] gnuplotでアニメ

gnuplotを用いて簡単なアニメーションを作ることができる．先ほどの練習問題3の結果を可視化しよう．各タイムステップにおいて 

write(1,*) 0.D0 , x
と2つの数値を並べてファイルに書き出す（0.D0はxy平面上でx座標を固定するため）．さらに，以下の内容でanim.txtというファイルを作成する． 

  set xrange [-1:1]
  set yrange [0:1]

  if (n>200) exit

  plot "ファイル名" every 1:1:n:0:n:0 with points pt 7 ps 2

  pause 0.05

  n=n+1

  reread

gnuplotを起動し，以下のコマンドによってアニメーションを見ることができる．

gnuplot> n=1
gnuplot> load 'anim.txt'

プログラム中で，y=1.D0　としている箇所が初速に相当する． y=0.5D0, y=2.D0として計算をやり直し，アニメを再生してみること．yrangeの値を変化させると縦方向の範囲を変えられる．

[5] 今日の宿題 (1.5点：12/9日まで)

練習問題3のプログラム，および初速を3通り変えた場合のそれぞれのt=2.D0におけるx座標を送付すること．

[6] 追加の宿題 (0.5点：12/9日まで)

式(26)で与えられる微分方程式を解け．t=0でx=1.D0, dx/dt=0.D0　とすること．Δt=1.D-2でt=10となるまで積分せよ．

write(1,*) t, x などとしてファイルにデータを書き出し，gnuplotを用いて横軸t, 縦軸xでグラフをかく．Aを変化させると，グラフの形状が変化する．プログラムおよびグラフの形状が変化する境界のAの値，それぞれのAの領域でのグラフの形状を記述した文章（もしくは画像）を送付すること．


\begin{displaymath} {d^2x\over dt^2}=-Ax-2{dx\over dt} \end{displaymath} (26)




日程表へ戻る <<