第4回

常微分方程式の数値解法（ルンゲ・クッタ法）

（2021/05/10 最終修正）

内容

常微分方程式の数値解法（ルンゲ・クッタ法）

目標

ルンゲ・クッタ法を使えるようになる
ルンゲ・クッタ法とオイラー法を比較し，両者の良い点と悪い点を見極める

1階常微分方程式の数値解法（ルンゲ・クッタ法）

前回はオイラー法を用いた数値解法の演習を行いました．オイラー法は単純でそれはそれでよいのですが、実際の解に近い数値解を出すためには、時間刻み幅 h をかなり小さくとらねばならないことが前回の課題3によってわかります．そこで、より精度の高い解法が必要となりますが、今回はそのような解法として、ルンゲ・クッタ法を扱い、それを用いた数値解法の演習を行います．

ルンゲ・クッタ法の導出については、このページの最後を参照してください（詳細は参考書を参照してください）．ルンゲ・クッタ法は次にしめすようなアルゴリズムで、オイラー法に比べれば多少複雑なものとなっています（太字の部分が異なります）．

T を解を求める時間の上限とし、h を時間刻み幅、N を時間刻み数とします．また、変数は y とし、a を初期値とします．r1, r2, r3, r4 はそれぞれ実数型の変数です．

T,N,h を定める

a の入力を受ける (input 関数）

y = a

t = 0 での y の値の記録

i = 0,1,2,....,N-1の順に

t = i*h

r1 = f(t, y)

r2 = f(t + h/2, y + (h/2)*r1)

r3 = f(t + h/2, y + (h/2)*r2)

r4 = f(t + h, y + h*r3)

y_new = y + (h/6)*(r1 + 2*r2 + 2*r3 + r4)

y = y_new

(一定間隔毎に)t+h と y_new を記録

を繰り返す

課題1（演習）

上記のルンゲ･クッタ法を用いて、次の問題を実際に解き、ルンゲ･クッタ法がオイラー法と比べてどの程度性能が良いか比べよ．次の初期値問題を 0 ≦ t ≦ 5 の範囲で h = 0.5 とし、実際にオイラー法とルンゲ･クッタ法を用いて解け．（動画では、上記アルゴリズムをそのままコーディングする。それだけのことなのだが、これまでの経験で、なかなかうまくいかないようで、摩訶不思議なコードができあがることが多いようなので。）

結果は大体次のようになります．h = 0.5とかなり粗く時間刻み幅をとっているにも関わらず、ルンゲ･クッタ法は良く微分解を近似できていることがわかります．それに比べてオイラー法は微分解から大きく離れてしまっています．これでは正しい計算とは言えないでしょう．

Code4-1.py : ルンゲ・クッタ法（1階常微分方程式）

# Code4-1.py

import math

import matplotlib.pyplot as plt

Y0 = 1.0

T = 5.0

Ns = 100

hs = T/float(Ns)

N = 10

h = T/float(N)

INTV = 1

t, y = [0.0]*(Ns + 1), [0.0]*(Ns + 1)

def f(t, y):

return y

#Y0 = float(input("Input Y0: "))

for i in range(0, Ns + 1):

t[i] = hs*i

y[i] = math.exp(t[i])

plt.plot(t, y, color="red", linewidth=2, label="$y(t)=\mathrm{e}^{t}$")

tg, Yg = [], []

Y = Y0

tg.append(0.0)

Yg.append(Y)

for i in range(0, N):

tt = i*h

Y = Y + h*f(tt, Y)

if i%INTV == 0:

tg.append(h*(i + 1))

Yg.append(Y)

plt.plot(tg, Yg, color="green", linewidth=2, label="Euler")

tg, Yg = [], []

Y = Y0

tg.append(0.0)

Yg.append(Y)

for i in range(0, N):

tt = i*h

r1 = f(tt, Y)

r2 = f(tt + h/2.0, Y + (h/2.0)*r1)

r3 = f(tt + h/2.0, Y + (h/2.0)*r2)

r4 = f(tt + h, Y + h*r3)

Y_new = Y + h*(r1 + 2*r2 + 2*r3 + r4)/6.0

Y = Y_new

if i%INTV == 0:

tg.append(h*(i + 1))

Yg.append(Y)

plt.plot(tg, Yg, color="blue", linewidth=2, label="Runge-Kutta")

plt.title("4-1")

plt.xlabel("t", fontsize=20)

plt.ylabel("y(t)", fontsize=20)

plt.legend()

plt.show()

#plt.savefig("fig4-1.png")

課題2

第3回の課題7の問題3（誤差の検討）について、ルンゲ・クッタ法を加えて検討せよ。

結果を言えば、オイラー法は誤差は h に比例して減少し、ルンゲ・クッタ法では誤差は h⁴ に比例して減少する。つまり、オイラー法では h を 1/10 にすれば誤差も 1/10 になるが、ルンゲ・クッタ法では、h を 1/10 にすれば誤差は 1/10000 になる。もしも、本課題でそのような関係が出ない場合には、プログラムが間違っていると思われる。但し、極端に小さな h については、有限桁での計算をしている関係から、おかしな結果が出る事があるので注意が必要である。

おまけ演習

課題2を自分で確かめて欲しいが、ここでは先の Code4-1.py を変更し、課題1の問題について、オイラー法とルンゲ・クッタ法の誤差の検討を行う。おまけとはいえ、結構大事な内容だ。

Code4-1-a.py : 誤差の検討

# Code4-1-a.py

import math

import matplotlib.pyplot as plt

Y0 = 1.0

T = 5.0

NN = [10, 100, 1000]

N = 10

err_e = []

err_r = []

def f(t, y):

return y

for N in NN:

h = T/float(N)

Y = Y0

for i in range(0, N):

tt = i*h

Y = Y + h*f(tt, Y)

err_e.append(abs(math.exp(5.0) - Y))

Y = Y0

for i in range(0, N):

tt = i*h

r1 = f(tt, Y)

r2 = f(tt + h/2.0, Y + (h/2.0)*r1)

r3 = f(tt + h/2.0, Y + (h/2.0)*r2)

r4 = f(tt + h, Y + h*r3)

Y_new = Y + h*(r1 + 2*r2 + 2*r3 + r4)/6.0

Y = Y_new

err_r.append(abs(math.exp(5.0) - Y))

print(err_e, err_r)

plt.plot(NN, err_e, color="green", linewidth=2, label="Euler")

plt.plot(NN, err_r, color="blue", linewidth=2, label="Runge-Kutta")

plt.xscale("log")

plt.yscale("log")

plt.title("4-1-a")

plt.xlabel("N", fontsize=20)

plt.ylabel("err", fontsize=20)

plt.legend()

plt.show()

#plt.savefig("fig4-1-a.png")

2階常微分方程式の数値解法（ルンゲ・クッタ法）

第3回の課題で扱った

をy1(t)=y(t), y2(t)=y'(t) とおいて次のように連立の１階の常微分方程式に帰着したもの

を再び考えます．このような連立の常微分方程式をルンゲ･クッタ法を用いて解く場合のアルゴリズムを次に示します．ただし、 T を解を求める時間の上限とし、h を時間刻み幅、N を時間刻み数とします．また、変数は y1, y2 とし、それぞれの初期値を a1, a2 とします．r11, r21, r31, r41, r12, r22, r32, r42 および y1_new, y2_new はそれぞれ実数型の変数です．また、f1(t, y1, y2) および f2(t, y1, y2) はそれぞれ実数型の関数で（２）式の場合には f1 が第一式の右辺、 f2 が第二式の右辺となります．

T,N,h を定める

a1, a2 の入力を受ける (input 関数）

y1 = a1

y2 = a2

t = 0 での y1, y2 を記録

i = 0,1,2,....,N-1の順に

t = i*h

r11 = f1(t, y1, y2)

r12 = f2(t, y1, y2)

r21 = f1(t + h/2, y1 + (h/2)*r11, y2 + (h/2)*r12)

r22 = f2(t + h/2, y1 + (h/2)*r11, y2 + (h/2)*r12)

r31 = f1(t + h/2, y1 + (h/2)*r21, y2 + (h/2)*r22)

r32 = f2(t + h/2, y1 + (h/2)*r21, y2 + (h/2)*r22)

r41 = f1(t + h, y1 + h*r31, y2 + h*r32)

r42 = f2(t + h, y1 + h*r31, y2 + h*r32)

y1_new = y1 + (h/6)*(r11 + 2*r21 + 2*r31 + r41)

y2_new = y2 + (h/6)*(r12 + 2*r22 + 2*r32 + r42)

y1 = y1_new

y2 = y2_new

(一定間隔毎に) t+h と y1_new, y2_new を画面を記録

を繰り返す

課題3（演習）

ルンゲ・クッタ法を用いて(2)を解け．前回同様、0 ≦ t ≦ 20の範囲において、h=0.0001, 0.001, 0.01, 0.1それぞれの場合について微分解 y(t)=cos(t) とどの程度あっているか各自調べよ．(第3回の課題3のルンゲ・クッタ法版)

Code4-3.py : 単振動・2階常微分方程式（ルンゲ・クッタ法）

# Code4-3.py

import math

import matplotlib.pyplot as plt

y10 = 1.0

y20 = 0.0

T = 20.0

NN = [200, 2000, 20000, 200000]

def f1(t, y1, y2):

return y2

def f2(t, y1, y2):

return -y1

for N in NN:

h = T/float(N)

INTV = N/200

y1 = y10

y2 = y20

tg, y1g, y2g = [], [], []

tg.append(0.0)

y1g.append(y1)

y2g.append(y2)

for i in range(0, N):

t = i*h

r11 = f1(t, y1, y2)

r12 = f2(t, y1, y2)

r21 = f1(t + h/2.0, y1 + (h/2.0)*r11, y2 + (h/2.0)*r12)

r22 = f2(t + h/2.0, y1 + (h/2.0)*r11, y2 + (h/2.0)*r12)

r31 = f1(t + h/2.0, y1 + (h/2.0)*r21, y2 + (h/2.0)*r22)

r32 = f2(t + h/2.0, y1 + (h/2.0)*r21, y2 + (h/2.0)*r22)

r41 = f1(t + h, y1 + h*r31, y2 + h*r32)

r42 = f2(t + h, y1 + h*r31, y2 + h*r32)

y1_new = y1 + h*(r11 + 2*r21 + 2*r31 + r41)/6.0

y2_new = y2 + h*(r12 + 2*r22 + 2*r32 + r42)/6.0

y1 = y1_new

y2 = y2_new

if (i + 1)%INTV == 0:

tg.append(h*(i + 1))

y1g.append(y1)

y2g.append(y2)

l = "h=" + str(h)

plt.plot(tg, y1g, linewidth=2, label=l)

plt.title("4-3")

plt.xlabel("t", fontsize=20)

plt.ylabel("y(t)", fontsize=20)

plt.legend()

plt.show()

#plt.savefig("fig4-3.png")

課題4（演習）

これまでは、横軸を t 、縦軸を y としたグラフを描いてもらいました．例えば横軸を y、縦軸を y' としたグラフを描いてみて下さい．どのようにすればよいでしょうか？またどのようなグラフになりますか？やってみて下さい．

Code4-4.py : 相平面表示

# Code4-4.py

import math

import matplotlib.pyplot as plt

y10 = 1.0

y20 = 0.0

T = 20.0

NN = [200]

def f1(t, y1, y2):

return y2

def f2(t, y1, y2):

return -y1

for N in NN:

h = T/float(N)

INTV = N/200

y1 = y10

y2 = y20

tg, y1g, y2g = [], [], []

tg.append(0.0)

y1g.append(y1)

y2g.append(y2)

for i in range(0, N):

t = i*h

r11 = f1(t, y1, y2)

r12 = f2(t, y1, y2)

r21 = f1(t + h/2.0, y1 + (h/2.0)*r11, y2 + (h/2.0)*r12)

r22 = f2(t + h/2.0, y1 + (h/2.0)*r11, y2 + (h/2.0)*r12)

r31 = f1(t + h/2.0, y1 + (h/2.0)*r21, y2 + (h/2.0)*r22)

r32 = f2(t + h/2.0, y1 + (h/2.0)*r21, y2 + (h/2.0)*r22)

r41 = f1(t + h, y1 + h*r31, y2 + h*r32)

r42 = f2(t + h, y1 + h*r31, y2 + h*r32)

y1_new = y1 + h*(r11 + 2*r21 + 2*r31 + r41)/6.0

y2_new = y2 + h*(r12 + 2*r22 + 2*r32 + r42)/6.0

#y1_new = y1 + h*r11

#y2_new = y2 + h*r12

y1 = y1_new

y2 = y2_new

if (i + 1)%INTV == 0:

tg.append(h*(i + 1))

y1g.append(y1)

y2g.append(y2)

l = "h=" + str(h)

plt.plot(y1g, y2g, linewidth=2, label=l)

plt.title("4-4")

plt.xlabel("y(t)", fontsize=20)

plt.ylabel("y'(t)", fontsize=20)

plt.legend()

plt.show()

#plt.savefig("fig4-4.png")

減衰振動

課題5

下の図のように、壁に一端を固定されたバネの他端におもりがつけられている状況を考えます。

このおもりを少しだけ右に引っ張ってから手をはなす。おもりの中心の時刻 t におけるつり合いの位置からのずれを y(t) とします。また、時刻 0 でのおもりの位置を y(0)=1、おもりの速度を y'(0)=0 とするとつぎの単振動の問題となります。（フックの法則を用いた）

これではおもりに働く力はバネの復元力のみで、ある意味非現実的です。今回は復元力以外におもりの速度 y' に比例する抵抗が働く場合を考えましょう（ダンパを加えたと考えよう）。改良されたモデル方程式は次のようになります。

ここで k,ω はそれぞれ抵抗力、復元力の強さを表す正の定数です。上式の問題を連立１階常微分方程式に変換し、ルンゲ・クッタ法を用いて解き、matplotlib を用いて解を可視化（グラフ）してください。

以下の３通りについて，それぞれグラフを描いてください．（例えば，ωをある値に固定して，k のみを変化させ，下記３通りを満たすような k を選ぶとよい．３通りのグラフを色を変えて重ね描きするとよい．）

k < ω の場合：減衰振動
k > ω の場合：過減衰
k = ω の場合：臨界減衰

ところで，ドアは一種のバネ振り子＋ダンパ系と見る事ができる．設計する際，ダンパの強さを臨界減衰になるよう設計するのが一番良いとされる．なぜであろうか？考えてみてください．（良いドアである為には次の二つが重要である．素早くしまる事，しまるときの音が小さい事．）

自動車やバイクのサスペンションも一種のバネ振り子＋ダンパ系であるが，おそらくはやはり臨界減衰に近いような設定が良いのであろうと想像はするが，実際にはより複雑な状況もあり，単純ではないだろう．

強制振動

ここからは強制振動系を考えよう．先の減衰振動系に外から周期的な外力を加える事を考える．

扱う方程式は，以下のようになる．

見ての通り，右辺が０であれば先の減衰振動系となります．つまりは，先の減衰振動系（放っておけば減衰して静止する）に対して，外から周期的な外力を加えているのが上式であります．

課題6（演習）

上式をルンゲ・クッタ法で解くプログラムを作成せよ．例えばパラメータ，初期値としては，

とし，時間刻み幅 h は h=0.1 とする．解 y(t) を赤色，外力 Acosω_et を緑色で表示することにする．

念のため、y1(t) = y(t), y2(t) = y'(t) とおいて上式を連立の１階の常微分方程式に帰着すると、以下になる。これを、ルンゲ・クッタ法のアルゴリズムで解く。

y1(0) = 1, y2(0) = 0

A = 0 の場合（つまり外力が無い場合）には，

となり，当然ながら減衰振動の様子がみられる．

A = 0.1 の場合には，

となり，過渡応答，定常応答が見られる．

A = 0.1, ω_e = 1.0 の場合には，次のようになる．

こちら、所謂共振であるが、減衰振動の初期値が大きいために面白みが少なく見えるが、例えば静止状態に外力を加えたとしても、同等の結果になる。僅かな外力から大きな振動がでるということろが面白い。

下にある、 Code4-6.py について、 y10 = 0.0, OmegaE = 1.0 とした結果が次の図になる。

Code4-6.py : 強制振動・2階常微分方程式（ルンゲ・クッタ法）

# Code4-6.py

import math

import matplotlib.pyplot as plt

y10 = 1.0

y20 = 0.0

T = 100.0

NN = [1000]

A = 0.1

k = 0.1

Omega = 1.0

OmegaE = 1.3

def f1(t, y1, y2):

return y2

def f2(t, y1, y2):

return -Omega*Omega*y1 - 2*k*y2 + A*math.cos(OmegaE*t)

for N in NN:

h = T/float(N)

INTV = N/200

y1 = y10

y2 = y20

tg, y1g, y2g, eg = [], [], [], []

tg.append(0.0)

y1g.append(y1)

y2g.append(y2)

eg.append(A*math.cos(OmegaE*0.0))

for i in range(0, N):

t = i*h

r11 = f1(t, y1, y2)

r12 = f2(t, y1, y2)

r21 = f1(t + h/2.0, y1 + (h/2.0)*r11, y2 + (h/2.0)*r12)

r22 = f2(t + h/2.0, y1 + (h/2.0)*r11, y2 + (h/2.0)*r12)

r31 = f1(t + h/2.0, y1 + (h/2.0)*r21, y2 + (h/2.0)*r22)

r32 = f2(t + h/2.0, y1 + (h/2.0)*r21, y2 + (h/2.0)*r22)

r41 = f1(t + h, y1 + h*r31, y2 + h*r32)

r42 = f2(t + h, y1 + h*r31, y2 + h*r32)

y1_new = y1 + h*(r11 + 2*r21 + 2*r31 + r41)/6.0

y2_new = y2 + h*(r12 + 2*r22 + 2*r32 + r42)/6.0

y1 = y1_new

y2 = y2_new

if (i + 1)%INTV == 0:

tg.append(h*(i + 1))

y1g.append(y1)

y2g.append(y2)

eg.append(A*math.cos(OmegaE*h*(i + 1)))

plt.plot(tg, y1g, color="red", linewidth=2, label="y(t)")

plt.plot(tg, eg, color="green", linewidth=2, label="Cos")

plt.title("4-6")

plt.xlabel("t", fontsize=20)

plt.ylabel("y(t)", fontsize=20)

plt.legend()

plt.show()

#plt.savefig("fig4-6.png")

課題7（演習）

1960年代初頭に気象学者エドワード・ローレンツ (Edward N. Lorenz) は流体力学に基づく次のような単純な３変数系の常微分方程式を示した．

このモデルはいわゆるカオスを生じるモデルとして有名である．パラメータは P, R, b の３種類があり、P=16, b = 4, R=35として、u-w 平面に解の軌道をプロットすると次のような図があらわれる．（初期値 (u0, v0, w0) = (5, 5, 5), t=0 から t=100 までの計算結果を線でつないだもの）

これはローレンツアトラクター (Lorenz attractor) と呼ばれるもので、ストレンジアトラクター（奇妙なアトラクター）とも呼ばれる。大変有名なので目にしたことがある人も多いと思う．このようなカオス的な解の挙動が天気予報があたらない理由となっているともいわれている．例えば P=16, b=4 として R を色々と変えてみると、解の挙動があれこれ変わる．ということで、次の動画でコーディングして色々と試してみる．但し、動画ではオイラー法で解くことにする。ルンゲ・クッタ法への変更はお任せする。

Code4-7.py : Lorenz Attractor (ローレンツ・アトラクター)

# Code4-7.py Lorenz Attractor

import math

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

u0 = 5.0

v0 = 5.0

w0 = 5.0

T = 100.0

N = 50000

h = T/float(N)

INTV = 1

P = 16.0

b = 4.0

R = 35.0

def fu(u, v, w):

return P*(v - u)

def fv(u, v, w):

return -u*w + R*u - v

def fw(u, v, w):

return u*v - b*w

ug, vg, wg = [], [], []

u = u0

v = v0

w = w0

ug.append(u)

vg.append(v)

wg.append(w)

for i in range(0, N):

t = i*h

u_new = u + h*fu(u, v, w)

v_new = v + h*fv(u, v, w)

w_new = w + h*fw(u, v, w)

u = u_new

v = v_new

w = w_new

if (i + 1)%INTV == 0:

ug.append(u)

vg.append(v)

wg.append(w)

fig = plt.figure()

ax = Axes3D(fig)

ax.plot(ug, vg, wg, linewidth=0.5)

ax.set_xlabel("u")

ax.set_ylabel("v")

ax.set_zlabel("w")

#plt.plot(ug, wg, linewidth=0.5)

#plt.title("Lorenz Attractor")

#plt.xlabel("u(t)")

#plt.ylabel("w(t)")

plt.show()

これより下の課題は、各自で取り組んでみるとよい。

課題8

重力の作用のみを受けて運動する物体を考えます（真空中でボールを投げるような状況）．ニュートンの運動の第二法則から物体の時刻 t での鉛直方向の位置を y(t) とするとその物体の鉛直方向での運動は次のような常微分方程式で表されます．（下向きを正とする）

m は物体の質量、g は重力加速度です．色々な初期値 y(0), y'(0) を与えてシミュレーションしてみて下さい．

課題9

さらに、上の問題に空気の抵抗のようなものを加味してみましょう．物体にかかる抵抗がそのときの物体の速度に比例するとすると、運動方程式は、

となります．ただし、r は抵抗の強さを表します．色々な初期値および m, r で計算してみてください．空気抵抗を加味する前と様子は変わりましたか？

なお、両問題とも解析的に解ける問題ですので、実際に解いてみて、今回のシミュレーション結果と見比べて理解を深めて下さい．

課題10

上の両問題では物体の鉛直方向のみの運動を考えましたが、水平方向の運動も加味するとどのような運動方程式になり、その解はどのようになりますか？運動方程式を立て、シミュレーションしてみて下さい．（抵抗が速度に比例する場合と速度の2乗に比例する場合があるようだ。自動車などの話をするときには、速度の2乗に比例と言ったことを良く言うように思う。それぞれ、粘性抵抗と慣性抵抗と呼ぶそうだ。速度が遅い場合には粘性抵抗が強く働き、速度が速いと慣性抵抗が強く働くようではあるが、どの程度わかっているのかよくわからない。）

課題10について、上向きを正とした場合（課題8，9では下向きを正としていました）には、次の方程式となります。

これを x = x1, x' = x2, y = y1, y' = y2 として連立の一階微分方程式に変換すると、次式となります。

こちらを解けば良い事になります。

次のコードは、ここまでの内容がわかっていればわかると思うので、コーディングのデモはしない。コードの解説と、このコードを使って遊ぶデモ動画を近々作成する。こちらもオイラー法であるが、ルンゲ・クッタ法に変更したい人はするとよい。

Code4-10.py : 放物運動のシミュレーション（打ち出し角度と飛距離の関係）

# Code4-10.py

# -*- coding: utf-8 -*-

import math

import matplotlib.pyplot as plt

#初期速さ

V0 = 10.0

T = 3.0

N = 300

h = T/float(N)

INTV = 1

#パラメータ

g = 9.8

m = 1.0

#空気抵抗

r = 1.0

def fx1(x1, x2):

return x2

def fx2(x1, x2):

return -r*x2/m

def fy1(y1, y2):

return y2

def fy2(y1, y2):

return -g - r*y2/m

#角度と飛距離の関係を記録

distg, thetag = [], []

#打ち出し角度を10度〜80度まで変化させる

for k in range(10, 80):

theta = k

#初期値

x1 = 0.0

x2 = V0*math.cos(math.radians(theta))

y1 = 0.0

y2 = V0*math.sin(math.radians(theta))

#グラフ用リスト

xg, yg = [], []

xg.append(x1)

yg.append(y1)

for i in range(0, N):

t = i*h

#オイラー法

x1_new = x1 + h*fx1(x1, x2)

x2_new = x2 + h*fx2(x1, x2)

y1_new = y1 + h*fy1(y1, y2)

y2_new = y2 + h*fy2(y1, y2)

# y=0 となる地点を線形補間で求める（飛距離）

# 接地したところで break でループを抜ける

if y1 >= 0 and y1_new < 0:

dy = y1_new - y1

dx = x1_new - x1

dist = -(dx/dy)*y1 + x1

xg.append(dist)

yg.append(0)

distg.append(dist)

thetag.append(theta)

break

x1 = x1_new

x2 = x2_new

y1 = y1_new

y2 = y2_new

if i%INTV == 0:

xg.append(x1)

yg.append(y1)

# x-y グラフを描く

#plt.plot(xg, yg, linewidth=0.5)

# theta-dist グラフを描く

plt.plot(thetag, distg)

plt.title("Code4-10")

plt.xlabel("theta", fontsize=20)

plt.ylabel("max dist.", fontsize=20)

#plt.legend()

plt.show()

#plt.savefig("fig4-10.png")

課題11

惑星の運動を考えます．万有引力の法則から２物体（質量はそれぞれ M,m）に働く力は

となります．rは２物体間の距離、Gは万有引力定数で

です．原点に太陽（質量 M）があり、惑星（質量 m）の位置を

とします．惑星に働く力は

となり、

より

とおくと、２階の常微分方程式の初期値問題

が得られます．これを１階の常微分方程式系に直し、数値計算にて解をもとめて下さい．

ルンゲ・クッタ法の導出について

簡単にルンゲ・クッタ法（Runge-Kutta method）について説明します。演習で取り扱ったルンゲ・クッタ法は4段公式とよばれるものですが、ここでは簡単の為、2段公式を取り扱います。(Runge と Kutta はカール・ルンゲおよびマルティン・クッタという数学者の名前由来)

次の常微分方程式の初期値問題の数値解を求めたいとしましょう。

(R1)

2段公式とは、t = t_n で（R1）の近似解が求められたとして、時間ステップ幅 h だけ進んだ t_n+1 = t_n + h における値 x_n+1 を次の形で計算する公式です。

(R2)

つまり、 x_n から x_n+1 へ進む１ステップにおいて方程式（R1）の第一式の右辺の f(t, x) を2回計算するので、2段公式と呼びます。この2段公式には、パラメータ α, β, γ₁, γ₂ が含まれますが、これらは方程式（R1）の真の解 x(t) のテイラー展開

(R3)

と（R2）式の第一式の同様の展開とが h のなるべく高い次数の項まで一致するように定めます。つまり、x_n = x(t_n) と仮定して（R2）式の k₂ の式の右辺を f(t, x) 周りでテイラー展開すると、

(R4)

となります。これと（R3）が h² の項まで一致するように

(R5)

とおく。公式（R2）で定められる近似解のテイラー展開が h の2次の項まで真の解のテイラー展開に一致するという意味で、この公式を2次の公式と呼びます。この公式は2段公式でもあるので、2段2次の公式とも呼ばれます。パラメータを定める方程式（R5）には自由度が2残っているのでその定め方によって色々な公式が導かれます。（改良オイラー法、ホイン法、修正オイラー法など）Wikipediaの次ページを見ると良い。ルンゲ＝クッタ法のリスト

4段4次のルンゲ・クッタ法の公式も同様にして求めることが出来ます。

目次に戻る