数値計算法
オイラー法とルンゲクッタ法の習得および比較
工学部・知能機械工学科 2 年 後期 テキスト No.2，Ver.6
担当教員：綴木1

馴

講義の目的
演習の目的は以下の通りである．

目的 1．オイラー法を導出する．
目的 2．オイラー法を理解する．
目的 3．オイラー法を用いて楕円関数を数値的に解く．
目的 4．オイラー法の欠点を知る．
目的 5．ルンゲクッタ法を理解する．
目的 6．ルンゲクッタ法を用いて楕円関数を数値的に解く．
目的 7．ルンゲクッタ法の特性を理解する．
目的 8．オイラー法とルンゲクッタ法を比較する．
目的 9．手際よい良い，報告書作成能力を身につける．
以上を以て，この演習に真摯に取り組み，研究・開発能力の向上を目指す．

2

種々の数値計算法
数値計算法の手法には大きく分けて● 決定論的手法（オイラー法，ルンゲクッタ法など）● 確率論的手法（モンテカルロ法など）がある．本講義では，まず，決定論的手法を取り扱う．時間があれば確率的手法としてモンテカルロ法も取り扱う．

2.1

参考：決定論的手法とは決定論的手法とは，系を記述する方程式（時間発展方程式や定常状態の状態方程式）に確率的要素（統計力学における熱ゆらぎのようなもの）が含まれておらず，この方程式を解くことによって，系の状態を数値的に計算する手法である．
このようなシミュレーションの例としては，(1) 古典力学のニュートン方程式における，運動方程式の数値積分による，粒子系のミクロな動力学の計算．
(2) 波動方程式：拡散方程式や流体方程式のような偏微分方程式の数値積分による，系の時間発展の計算．などがある．
ニュートンの運動方程式は多数の粒子（原子・分子）の位置及び運動量の連立常微分方程式であり，この数値積分によるシミュレーションは分子動力学シミュレーション法と呼ばれる．この方法は，原子・分子だけでなく相互作用の定められた古典力学に従う他，自由度系，例えば天体の多体問題にも適用できる．
これに対して，偏微分方程式で記述される系の時間発展の計算は，流体力学の分野では差分化法や有限要素法として知られている．また，この種の偏微分方程式の問題は，しばしば，その方程式の固有値問題に置き換えられ，行列の対角化の問題に帰着されることがある．
これらのシミュレーションの特徴は，系の初期の状態が与えられれば，計算によりその系の振る舞いが一意に決まると言う点である．量子力学では，波動方程式の解である波動関数には確率的解釈がなされるが，その時間発展を決める波動方程式自身は決定論的方程式である．

3

オイラー (Euler) 法の導出

ニュートン方程式の数値積分の方法はどのようにすればよいか．簡単な例として，１次元の運動をする１粒子の時間発展方程式

dx(t)
dt
dy(t)
dtdvy (t)dtdvx (t)dt= f (x(t), vx (t), t) ≡ vx (t)= h(y(t), vy (t), t) ≡ vy (t)= g(x(t), y(t), vx (t), t)= k(x(t), y(t), vy (t), t)(1)を例にとって考える．ここに f (x(t), vx (t), t) と g(x(t), y(t), vx (t), t) は x, y, vx と t の関数である．関数 f (x(t), vx (t), t) と g(x(t), y(t), vx (t), t) が時間 t に依存しないときは「自律系」と呼ばれる．
時刻 t における変数の値 x(t) と vx (t) が分かっているものとして，微分方程式 (1) に従って時間発展するときに，微少時間 ∆t 後の時刻 t + ∆t における変数の値 x(t + ∆t) および vx (t + ∆t)を計算する処方箋（アルゴリズムとかスキームと呼ばれる）を導出する．x(t + ∆t) と vx (t + ∆t)
を時間について t のまわりでテイラー (Tahlor) 展開すると，

1 d2 x(t)dx(t)∆t +(∆t)2 + · · ·dt
2 dt2dvx (t)1 d2 vx (t)vx (t + ∆t) = vx (t) +∆t +(∆t)2 + · · ·2dt
2 dtx(t + ∆t) = x(t) +(2)となる．したがって，微少時間 ∆t 後の変数の値を計算したければ，この右辺の展開式を適当な次数で打ち切った表現を用いればよい．展開の O(∆t) の項で打ち切る方法を オイラー法，より高次の展開項まで取り組む方法にルンゲクッタ（Runge-Kutta）法がある．本解説書では，まずオイラー法について解説する．
オイラー法は，常微分方程式を数値的に解く最も簡単な解法である．式 (2) の右辺の展開をO(∆t) の項で打ち切り，展開係数の dx(t)/dt と dvx (t)/dt に式 (1) の表現を代入すれば，x(t + ∆t) = x(t) + (∆t)f (x(t), vx (t), t)vx (t + ∆t) = vx (t) + (∆t)g(x(t), y(t), vx (t), t)(3)となる．この式は，x(t) と vx (t) が与えられたときに微少時間 ∆t 後の x と vx の値を与える計算式になっており，無限小時間間隔の時の微分を有限時間間隔 ∆t での差分に置き換えたことになるので，このような計算式は「差分スキーム」と呼ばれる．また，差分スキームに用いる微少時間 ∆t のことを「時間刻み幅」と呼ぶことにする．
式 (3) の差分スキームは Taylor 級数展開の打ち切りの為に O((∆t)2 ) の誤差を含んでおり，（この誤差を「打ち切り誤差」と呼ぶ），O(∆t) まで正しいスキームである．このような差分スキームを Euler 差分スキームと呼ぶ．式 (3) と同様にすれば，y(t + ∆t) = y(t) + (∆t)h(y(t), vy (t), t)vy (t + ∆t) = vy (t) + (∆t)k(x(t), y(t), vy (t), t)が得られる．

(4)

4

オイラー (Euler) 法の実装

■ここからが重要■与えられた初期値 x(0) と vx (0) からスタートして，式 (3) にしたがって１ステップずつ時間発展を繰り返し追って行けば，微分方程式 (1) の近似解が得られる．
式 (3) で無視した項の存在を考えれば，∆t の値が大きすぎると，近似によって得られた解が微分方程式 (1) の真の解を正しく近似しないことは容易に理解できる．しかしながら，逆に ∆tを小さくとればとるほど，計算された解は微分方程式 (1) の真の解に近づいて行くかと言えばそうでもなく，次のような困難が生じる．
● ∆t を小さくすればするほど，一定の時間 t の時間発展を追跡する為に必要なステップ数が多くなり，計算量の増大を招く．
● 計算機のなかでは実数の演算は有限の桁で打ち切られた表現（「数値を丸める」と表現される）を用いて行われる為，演算をする度にこのような「丸め」による誤差が生じる．
従って ∆t を小さくとったときのステップ数の増大は丸めの誤差の堆積をまねき，制度の点で不都合が生じる．
これらの問題がある為に，あまりにも小さな時間刻み幅 ∆t の値を用いることは不適切である．

5オイラー法の課題

オイラー差分法を用いて，次の中心重力場（重力場あるいはクーロン場に相当）の中で運動する２次元１粒子の系の時間発展方程式

dx(t) = vx (t)dtdh(y(t), vy (t), t) =y(t) = vy (t)dtdx(t)
g(x(t), y(t), vx (t), t) =vx (t) = −32dt(x(t) + y(t)2 ) 2dy(t)k(x(t), y(t), vy (t), t) =vy (t) = −3dt(x(t)2 + y(t)2 ) 2
f (x(t), vx (t), t) =(5)を数値的に解き，その解の時間変化の軌道を gnuplot を使ってグラフィックス表示してみよ．
初期値は x = 1.8, y = 0, vx = 0, vy = 1 としてみよ（自分なりに決めても良い）．
■注意■ c 言語では，xy を pow(x, y) と記述して計算する．また計算機においては， 32 と言う表現と 1.5 という表現は意味が異なるので注意すること．前者は整数演算の場合その値は 1 になってしまう．
■挑戦■さらに，式 (5) を解析的に解き（数学的に微分方程式を解く），数値解と解析解で得られた真の値との比較を行え．ただし，式 (5) では，簡単のため粒子の質量などの物理定数は1 にとってある．具体的には c 言語により，”dat”ファイルへデータを出力し，gnuplot で出力してみよ．

6ルンゲクッタ（Runge-Kutta）法

オイラー法が式 (2) の展開の O(∆t) まで正しいスキームであるのに対して，ルンゲクッタ法
mは ∆t より高次まで正しいスキームである．式 (2) の高階の展開係数 d dtx(t)などはたとえ解析にn書き下せたとしても，一般的には複雑な形をしており，数値計算には向かない．しかしルンゲクッタ法では，これら高階の展開係数は数値的に計算できる．
一般に広く用いられているルンゲクッタ法は，４次のルンゲクッタ法と言い，展開の O((∆t)4 )まで正しいスキームである．そのスキームを以下に示す．

x(t + ∆t) = x(t) + (α1 + 2α2 + 2α3 + α4 )∆t/6
y(t + ∆t) = y(t) + (γ1 + 2γ2 + 2γ3 + γ4 )∆t/6vx (t + ∆t) = vx (t) + (β1 + 2β2 + 2β3 + β4 )∆t/6vy (t + ∆t) = vy (t) + (δ1 + 2δ2 + 2δ3 + δ4 )∆t/6(6)
α1 = f (x(t), vx (t), t) γ = h(y(t), v (t), t)1yβ1 = g(x(t), y(t), vx (t), t)





δ1 = k(x(t), y(t), vy (t), t)
α2 = f (x(t) + α1 ∆t/2, vx (t) + β1 ∆t/2, t + ∆t/2) γ = h(y(t) + γ ∆t/2, v (t) + δ ∆t/2, t + ∆t/2)21y1
β2 = g(x(t) + α1 ∆t/2, vx (t) + β1 ∆t/2, t + ∆t/2)






δ2 = k(y(t) + γ1 ∆t/2, vy (t) + δ1 ∆t/2, t + ∆t/2)
α3 = f (x(t) + α2 ∆t/2, vx (t) + β2 ∆t/2, t + ∆t/2) γ = h(y(t) + γ ∆t/2, v (t) + δ ∆t/2, t + ∆t/2)
32y2
 β3 = g(x(t) + α2 ∆t/2, vx (t) + β2 ∆t/2, t + ∆t/2)





δ3 = k(y(t) + γ2 ∆t/2, vy (t) + δ2 ∆t/2, t + ∆t/2)
α4 = f (x(t) + α3 ∆t, vx (t) + β3 ∆t, t + ∆t) γ = h(y(t) + γ ∆t, v (t) + δ ∆t, t + ∆t)43y3

β4 = g(x(t) + α3 ∆t, vx (t) + β3 ∆t, t + ∆t)





(7)

δ4 = k(y(t) + γ3 ∆t, vy (t) + δ3 ∆t, t + ∆t)

ただし，ここでは楕円関数を用いて，dx(t) = vx (t)dtdx(t)
g(x(t), y(t), vx (t), t) =vx (t) = −3dt(x(t)2 + y(t)2 ) 2dh(y(t), vy (t), t) =y(t) = vy (t)dtdy(t)
k(x(t), y(t), vy (t), t) =vy (t) = −3dt(x(t)2 + y(t)2 ) 2
f (x(t), vx (t), t) =とする．

(8)
(9)
(10)
(11)

このスキームが Taylor 展開 (2) の O((∆t)4 ) まで正しいことは式 (7) の第一式を ∆t に関してTaylor 展開しみれば確認できる．このように，ルンゲクッタ法の制度はオイラー法の制度に比べてずっと高く，比較的大きな刻み幅 ∆t を使っても一切り誤差か堆積することなく，長時間の時間発展を計算するのに適している．

7

課題

5 節の課題を，４次のルンゲクッタ法を用いて解き，その精度をオイラー法の精度と比較せよ．特に長時間の振る舞いで両者の違いが明白になるはずである．

8

参考word 等を用いて報告書を書くにあたって，目的の Window 画面を Word に貼り付ける方法について説明する．まず，目的の Window をアクティブにした状態で，「Ctrl」＋「Alt」＋「PrintScreen」を押すことでクリップボードに保存される．この状態で Word の画面に戻り，「Ctrl」＋「v」を押すことで貼り付けることができる．