常微分方程式を用いたシミュレーション

放物運動や古典力学などの物理シミュレーションは、常微分方程式 (ordinary differential equation, O.D.E.) で記述できるものが非常に多いです。 常微分方程式とは変数が1つだけの関数の微分を求める方程式で、これを数値的に解く方法を数値微分や数値解法と呼びます。 常微分方程式の数値解法には、差分方程式をもとにしたオイラー法やホイン法、ルンゲ・クッタ法などがあります。

本講義では、常微分方程式の数値解法を述べたあとに、いくつかの物理法則を実際にシミュレーションしていきます。

常微分方程式とは

微分とは、ある独立変数を持つ関数の変化率を求める数学的な作業です。 つまり、点xから点(x + h)までの平均増加率を求める作業ですので、数式で表すと以下の様になり、これを差分方程式と呼びます。 $h$は刻み幅で、$\frac{y(x + h) - y(x)}{h}$の部分は引き算（差分）したものを割ったもの（商）ですので、差分商と呼びます。 左辺の$f^{\prime}(x)$は、$f(x)$を微分して得られる導関数と呼ばれるものです。

$$\displaystyle f^{\prime}(x) = \lim_{h \to \infty} \frac{y(x + h) - y(x)}{h}$$

導関数に、点$x = a$などの実際の値を入れたものを微分係数と呼びます。

$$\displaystyle f^{\prime}(a) = \lim_{h \to \infty} \frac{y(a + h) - y(a)}{h}$$

例として、簡単な関数の微分を解析的に行ってみます。 $y = 12 x - 8 e^x + 9$という関数の導関数とx=0.7における微分係数は、それぞれ以下の式と値になります。 また、原始関数と微分係数は以下の図の様になり、微分係数はある点での接線の傾きと言えます。

• 原始関数：$y = 12 x - 8 e^x + 9$
• 導関数：$y^{\prime} = y - 12 - 8 e^x$
• 微分係数 （x=0.7）：$y^{\prime}(0.7) = -4.11002166$

図：$y = 12 x - 8 e^x + 9$のグラフと$x=0.7$における微分係数

微分方程式は、微分を含んだ方程式のことです。 例えば、$y$が$x$の関数のとき、$x$と$y$及び$y^{\prime}$の導関数を含むものなどです。 $y^{\prime}$が変数$x,y$のみで与えられるので、以下の様な式として表せます。(1)の方は1階常微分方程式、(2)の方は2階常微分方程式です。 これらの式の様に、変数が$x$しかないものを常微分方程式といい、本講義資料ではこの形式の常微分方程式のみを取り扱います。

$$\begin{align} y^{\prime}(x) &= \frac{dy}{dx} \\ y^{\prime\prime}(x) &= \frac{d^2y}{dx^2} \end{align}$$

常微分方程式を解く

常微分方程式の一般的な話がわかったところで、ここからは、運動方程式を例に、常微分方程式を解くという意味とその用語について見て行きます。 具体的には、質量$m$の剛体の時刻$t$におけるx軸上の座標$x(t)$を解析的に求めてみます。 ここまでは、まだ数値解法の話ではなく、常微分方程式を解くことで何が出来るかを述べています。

まず、$F$を力、時刻$t$の加速度を$a(t)$として、この剛体の運動方程式を立ててみます。

$$F = m a(t)$$

加速度は速度$v(t)$の1階常微分や位置$x(t)$の2階常微分ですので、以下の様に書けます。 また、時刻$t$とx軸上の座標$x$だけが変数ですので、この方程式は常微分方程式です。 x,y,z軸の三次元空間での運動方程式を解きたい場合は、変数が3つになるため、偏微分方程式となります。

$$\displaystyle a(t) = \frac{dv(t)}{dt} = \frac{d^2 x}{d t^2}$$

次に、位置を求めたいので$x(t)$を運動方程式に代入し、$x(t)$について解くべく両辺を積分します。

$$\displaystyle \begin{aligned} F = m a(t) &= m \frac{d^2 x(t)}{dt^2} \\ \frac{d^2 x(t)}{dt^2} &= \frac{F}{m} \\ \frac{d x(t)}{dt} &= \int \left[ \frac{d^2 x}{dt^2} dt \right] = \int \frac{F}{m} dt = \frac{F}{m} + C_1 \\ x(t) &= \int \left[ \frac{d x(t)}{dt} \right] = \int \left( \frac{F}{m} + C_1 \right) dt = \frac{F}{2m} t^2 + C_1 t + C_2 \\ \end{aligned}$$

積分定数$C_1, C_2$を含む微分方程式の解を、一般解と呼びます。 積分定数は任意の値ですので、この微分方程式には無数に解が存在することになります。 今回の運動方程式の場合は、それぞれ以下の物理的な意味を持ちます。

• $C_2$：初期位置。これは$t=0$のとき、$x(t=0) = C_2$となるため
• $C_1$：初期速度。これは$t=0$のとき、$\displaystyle \frac{dx}{dt}(t=0) = C_1$となるため

積分定数$C_1,C_2$の具体的な値、つまり、関数の初期値を指定して微分方程式を解く問題のことを初期値問題と呼びます。 初期値が与えられるとこの微分方程式には一意に決まるため、これを特殊解と呼びます。

例えば、以下のグラフに示すように、初期位置$C_2$、初期速度$C_1$を実際に指定（初期値問題）すると具体的な位置（特殊解）を求めることができます。なお、微分方程式としては、どれも正しい答え（一般階）です。

図：2つの積分定数$(C_1, C_2)$を変化させた時の運動方程式の結果

数値解法のイメージ

ここでは、実際に常微分方程式の数値解法の（差分法の）イメージについて述べていきます。

その大まかなイメージは以下の3ステップになります。

1. 常微分方程式を差分法的式で近似して
2. 現在の点での値を計算し
3. 差分方程式で得られた解から本当の常微分方程式の解を推定する、というものです

これを物理的な意味として言い換えると、現在の状態からなんらかの式に従って次の状態を求める作業に相当します。 なので、常微分方程式を数値的に解く考え方は、解析的に解くのとは少し違うものになります。

---

ここからは、上記の3ステップについて詳しく見ていきます。

微分の定義は、差分商の刻み幅を$h$を$\infty$にするものですが、計算機では実際には$h \to \infty$にすることはできません。 なので、十分に小さな$h$を利用して、これを微分の近似値とします。

$$\displaystyle y^{\prime}(x) \approx \frac{y(x + h) - y(x)}{h}$$

これを、常微分方程式$y^{\prime} = f(x, y)$に代入して$y(x+h)$について解くと、以下の式が得られます。 つまり、常微分方程式は有限の差分方程式で近似できる、ということになります。 また、得られた式は、$h$だけ先の$y$の値$y(x+h)$を現在の$x,y$から近似的に求められることを示しています。 これが、前述した常微分方程式の数値解法が、現在の状態からなんらかの式に従って次の状態を求めるもの、という意味になります。

$$\displaystyle \begin{aligned} \frac{y(x + h) - y(x)}{h} &\approx f(x, y) \\ y(x + h) &\approx y(x) + h f(x, y) = y(x) + h y^{\prime} \end{aligned}$$

この式を図示したものが以下になります。 点$x_{j-1}$での傾き$y^{\prime}$から次の点$x_{j}$を予想するのが、常備分方程式の数値的解法の基本的なイメージです。 しかし、単純に点$x_{j-1}$での傾きを使って刻み幅$h$だけ離れた次の点を予想するだけでは、2次以上の曲線ではうまくいかなそうなことが分かるかと思います。 何らかの方法で$x_{j-1}$から$x_{j}$の間の平均的な傾きを推定して、実際の$y(x_{j})$の値を精度良く求める必要があります。 つまり、常微分方程式の数値解法のポイントは、傾きの取扱いと言えます。

図：常微分方程式を数値的に解くイメージ。何らかの方法で平均的な傾きを推定し、実際の$y(x_{j})$の値を求める

常微分方程式の数値解法には、平均的な傾きの取り扱いの違いで、以下に示す3つの代表的な方法があります。 3つの方法それぞれ名前が付いているのですが、実はルンゲ・クッタ法として一般化することができます。 このため、オイラー法は1次のルンゲ・クッタ法、ホイン法は2次のルンゲ・クッタ法とも呼ばれます。

数値解法	使用する式	意味	次数
オイラー法	$f(x,y)$	点$x_{j-1}$での傾き	1
ホイン法	$\frac{1}{2}(k_1 + k_2), ~~~~k_1 = h f(x, y), ~~~~k_2 = h f(x + h, y + k_1)$	点$x_{j-1}$と$x_{j}$での傾きの平均	2
4次のルンゲ・クッタ法	$\frac{1}{6}(k_1 + 2 k_2 + 2 k_3 + k_4), ~~~~k_1 = h f(x, u), ~~~~k_2 = h f(x + \frac{h}{2}, y + \frac{k_1}{2}), ~~~~k_3 = h f(x + \frac{h}{2}, y + \frac{k_2}{h}), ~~~~k_4 = h f(x + h, y + k_3)$	点$x_{j-1}$と$x_{j}$、その間での2点の合計4点での傾きの重み付け平均	4

これらの方法は、$x_{j}$の値の計算に$x_{j-1}$の値のみを用いる1段階法に分類されるものです。 $x_{j-1}$だけでなく$x_{j-2}$などの複数の値を用いるものは、多段法に分類されます。

運動方程式

常微分方程式を使った最もわかりやすく簡単な例が自由落下のシミュレーションです。 まずはこの自由落下のシミュレーションプログラムを示します。 なお、常微分方程式は変数が1つですので、x座標のみしか取り扱えません。 このプログラムで利用しているオイラー法と連立常微分方程式の解法については、後述します。

以下に、教科書にある自由落下のシミュレーションをJuliaに移植し、アニメーションを出力するものを示します。

自由落下のシミュレーションを行うJuliaプログラム

using Plots

# 定数
const G = 9.80665  # 重力加速度

# メイン実行部
t = 0.0   # 時刻t
h = 0.01  # 時刻の刻み幅

# 係数の入力
println("初速度v0を入力してください:")
v = parse(Float32, readline())
println("初期高度x0を入力してください:")
x = parse(Float32, readline())

println("時刻t, 位置x, 速度v")
println(t, ", ", x, ", ", v)
# グラフデータに現在位置などを追加
tlist = [t]
xlist = [x]
vlist = [v]

# グラフの表示
plot(tlist, xlist, title="Freefall", aspect_ratio=1.0)
# aspect_ratio=1.0, xlim=(-25.0, 25.0), ylim=(-25.0, 25.0)

# 自由落下の計算（オイラー法）
println("Computing freefall")
while x >= 0.0   # 地面に達するまで計算
    global t += h      # 時刻の更新
    global v += G * h  # 速度の計算
    global x -= v * h  # 位置の更新
    # println(t, ", ", x, ", ", v)  # 現在時刻と現在の位置
    # グラフデータに現在位置を追加
    push!(tlist, t)
    push!(xlist, x)
    push!(vlist, v)
end

simulation_period = length(xlist)
println("Simulation period: ", simulation_period)

# アニメーション
println("Making animation")
anim = Animation()
# yの座標は刻み幅hと一緒
ylist = repeat([h], simulation_period)
anim = @animate for t in 1:simulation_period
        plot(aspect_ratio=1.0, ylim=(0, findmax(xlist)[1]), xlim=(-5, 5))
        my_title = "Timestep: $(t), Velocity: $(round(vlist[t]; digits=1))"
        plot!(ylist[1:t], xlist[1:t], label="", title=my_title)
        # plot!(ylist[t:t], xlist[t:t], label="", seriestype=:scatter)
        plot!(ylist[t:t], xlist[t:t], label="", markershape=:circle, markersize=5)
end
println("Generating animation")
gif(anim, "freefall.gif", fps=25)

savefig("freefall.png")

実行結果

julia> include("freefall.jl")
初速度v0を入力してください:
0
初期高度x0を入力してください:
10
時刻t, 位置x, 速度v
0.0, 10.0, 0.0
Computing freefall
Simulation period: 144
Making animation
Generating animation
[ Info: Saved animation to /root/freefall.gif
"/root/freefall.png"

左図：初速度v0=0, 初期高度x0=10のときの結果。右図：初速度v0=-10, 初期高度x0=10のときの結果。

オイラー法

オイラー法 (Euler method) は、以下に示す様な式を用いるもので、数値解法のイメージで解説したものと同一です。 1次のルンゲ・クッタ法とも呼ばれます。 オイラー法のイメージを以下に示します。

$$\begin{aligned} x_{i+1} &= x_i + h \\ y_{i+1} &= y_i + h f(x_i,y_i) \\ \end{aligned}$$

これを図示したものが以下の図になります。

図：オイラー法のイメージ。点$x_{j-1}$での傾きのみを使って、刻み幅$h$だけ離れた次の点を求める。

積分の精度を高めるには、刻み幅を小さくするしかありません。 h=0.1よりも、h=0.01で計算する方が、高い精度で積分が行えます。 しかし、大きくしすぎると足し算の回数が増えすぎてしまい、丸め誤差が累積して数値的な精度が落ちてしまいます。 これは数値積分の区分求積法と同様です。

また、誤差は刻み幅$h$に比例します。つまり、刻み幅を$h$から$h/2$へと半分にすると、誤差も半分になります。 なお、常微分方程式の数値解法において、1ステップで現れる誤差を局所離散化誤差、最終的なステップ数に至った時の誤差を大域離散化誤差と呼びます。

プログラム例

djs-office-hours/07-solving_odes.jl at master · da-james/djs-office-hoursにオイラー法のJuliaコードが実装されていますので、こちら少し改変して紹介します。

オイラー法のjuliaコード

using Plots

"""
eulers_method(f::Function, α::Real, a::Real, b::Real, N::Int64)
Approximate the solution of the initial-value problem (IVP)
    y' = f(t,y), a <= t <= b, y(a) = α
at (N + 1) equally spaced numbers in the interval [a,b].
# Arguments
- `f::Function` : the ODE
- `α::Real` : the initial condition such that y(a) = α
- `a::Real` : the left-sided endpoint
- `b::Real` : the right-sided endpoint
- `N:Int64` : the spacing on the interval of [a,b]
"""
function eulers_method(f::Function, α::Real, a::Real, b::Real, N::Int64)

    n1 = N + 1
    u = zeros(n1,3)

    h = (b - a) / N
    u[1,1] = a
    u[1,2] = α

    for i in 2:n1
        u[i,3] = h * f(u[i-1,1], u[i-1,2])
        u[i,2] = u[i-1,2] + u[i,3]
        u[i,1] = a + (i - 1) * h
    end

    return u
end

df(t,y) = y - t^2 + 1
f(t) = -0.5*exp(t) + t^2 + 2t + 1
u0 = 0.5
a = 0
b = 2
n = 10
tlin = 0:0.2:2
sol1 = eulers_method(df, u0, a, b, n)
plot(tlin, f.(tlin), label="exact", legend=:bottomright, dpi=150)
plot!(sol1[:,1], sol1[:,2], markershape=:o, label="eulers")

実行結果

図：オイラー法のJuliaコードの実行結果

原始関数が二次関数なので、あまりうまく行っていません。 局所離散化誤差が累積し、最終的な10ステップ目では大域離散化誤差が大きくなってしまうことが分かります。

ルンゲ・クッタ法

ルンゲ・クッタ法 (Runge-Kutta method) は、微分方程式の数値解法で最も有名なものの1つです。 この方法では仮のステップを4つに増やして、より多数の点で傾きを計算し、平均的な傾きの推定精度を高めます。

まずは、公式を以下に示します。 4つの点の傾きを加重平均するため、正確には4次のルンゲ・クッタ法と呼ばれます。 ホイン法と同様に、$k_1 \to k_2 \to k_3 \to k_4$の順に計算をしていく必要があります。

$$\begin{aligned} x_{i+1} &= x_{i} + h, \\ k_1 &= h f(x_i, y_i),\\ k_2 &= h f(x_i + \frac{h}{2}, y_i + \frac{k_1}{2}), \\ k_3 &= h f(x_i + \frac{h}{2}, y_i + \frac{k_2}{2}), \\ k_4 &= h f(x_i + h, y_i + k_3), \\ y_{i+1} &= y_i + \frac{1}{6}(k_1 + 2 k_2 + 2 k_3 + k_4)\\ \end{aligned}$$

$y = 12 x - 8 e^x + 9$をホイン法で解いた場合のイメージを図示したものが以下になります。 注意するポイントとしては、$k_3$は、$k_2$での傾きを始点から$h/2$だけ進んだものである、という事です。 点$(x_i, y_i)$から$k_1$、$k_2$と$h/2$ずつ進んだ後は、一旦、点$(x_i, y_i)$に戻ってから$k_3$、$k_4$を計算する感じです。 図中では、それぞれ実線と点で表されています。

図：ルンゲ・クッタ法のイメージ。$k_1$と$k_4$に加え、それぞれで得られた傾きで$h/2$だけ進んだ点$k_2, k_3$での傾きの加重平均を傾きとします。

オイラー法の公式と図の対応関係は以下の様なものになります。

1. オイラー法で増分$k_1$を求め、その傾きで仮のステップを半ステップ$h/2$だけ進み、点$(x_i + h/2, y_i + k_1/2)$に到達します。
2. その点$(x_i + h/2, y_i + k_1/2)$での傾き$k_2$を計算し、その傾きで仮のステップを半ステップだけ進んで次の点$(x_i + h/2, y_i + k_2/2)$に到達します。
3. 点$(x_i + h/2, y_i + k_2/2)$での傾き$k_3$を計算し、その傾きで仮のステップをフルに進んで、点$(x_i + h, y_i + k_3)$に到達します。
4. 点$(x_i + h, y_i + k_3)$での傾きを$k_4$を計算します。
5. 4つの増分$k_1, k_2, k_3, k_4$を、$1:2:2:1$で荷重平均したものを最終的な増分の推定値として真のステップを進みます。

プログラム例

ルンゲクッタ法を行うjuliaプログラム

using Plots

function rk4_method(f::Function, α::Real, a::Real, b::Real, N::Int64)
    n1 = N + 1
    u = zeros(n1,2)
    h = (b - a) / N
    u[1,1] = a
    u[1,2] = α
    for i in 2:n1
        t = u[i-1,1]
        w = u[i-1,2]
        k1 = h * f(t, w)
        k2 = h * f(t + h / 2, w + k1 / 2)
        k3 = h * f(t + h / 2, w + k2 / 2)
        k4 = h * f(t + h, w + k3)
        u[i,2] = w + (k1 + 2 * k2 + 2 * k3 + k4) / 6
        u[i,1] = a + (i - 1) * h
        println(i, u[i,:])
    end
    return u
end

df(t,y) = y - t^2 + 1
f(t) = -0.5*exp(t) + t^2 + 2t + 1
u0 = 0.5
a = 0
b = 2
#n = 10
println("Please input a total step")
input_n = readline()
n = parse(Int64, input_n)
tlin = 0:0.2:2
sol1 = rk4_method(df, u0, a, b, n)

println("computing exact result...")
exact = plot(tlin, f.(tlin), label="exact", legend=:bottomright, dpi=150)
println("Done!\n")

println("Generating result animation...\n")
n1 = n + 1
anim = Animation()
anim = @animate for i in 2:n1
    plot!(sol1[1:i,1], sol1[1:i,2], markershape=:o, markercolor=:blue, linecolor=:blue, label="")
end

file_name = "rk4_method_" * input_n * ".gif"
@time gif(anim, file_name, fps = 2)
println("Done!\n")

# plot!(sol1[:,1], sol1[:,2], markershape=:o, label="rk4")
# savefig("./rk4_method.png")

実行結果

図：ルンゲ・クッタ法のJuliaコードの実行結果

かなりうまく行っていることがわかるかと思います。

連立常微分方程式

ここまでは、従属変数と方程式の数が1つだけの$y^{\prime}=f(x,y)$の様な常微分方程式を扱ってきました。 ここからは、従属変数と方程式が2つ以上の連立常微分方程式の解法について見ていきます。 複雑に思えるかもしれませんが、基本的な考え方だけでなく利用できる解法も、方程式が1つだけのものと一緒です。

例えば、以下の様な$y^{\prime}_0, y^{\prime}_1, y^{\prime}_2, ..., y^{\prime}_n$についての複数の連立微分方程式を解くことに なります。 解く範囲は[a, b]として、初期の値 $y^{\prime}_0 (a), y^{\prime}_1 (a), ..., y^{\prime}_n (a)$が既知である様なものです。

$$\begin{aligned} y^{\prime}_0 (x) &= f_0 (x, y_0, y_1, y_2, ..., y_n) \\ y^{\prime}_1 (x) &= f_1 (x, y_0, y_1, y_2, ..., y_n) \\ ... \\ y^{\prime}_n (x) &= f_n (x, y_0, y_1, y_2, ..., y_n) \end{aligned}$$

$y_0, y_1, ..., y_0$と$f_0, f_1, ..., f_n$にベクトル表記${\bf y }(x), {\bf f }(x, {\bf y})$を用いると、これまでに学んだ（連立していない）常微分方程式と見た目が一緒になります。 また、初期値は${\bf y}(x_0) = {\bf y}_0$の様になります。

$${\bf y }(x) = \left\{ \begin{array}{ll} y_0 (x) \\ y_1 (x) \\ ... \\ y_n (x) \\ \end{array} \right. , ~~~ {\bf f }(x, {\bf y}) = \left\{ \begin{array}{ll} f_0 (x, {\bf y }) \\ f_1 (x, {\bf y }) \\ ... \\ f_n (x, {\bf y }) \\ \end{array} \right.$$

$${\bf y }^{\prime}(x) = {\bf f }(x, {\bf y})$$

例えば、オイラー法は以下の様に表せます。 ここに、初期値を代入して次々と計算していくことで実際に連立常微分方程式を解くことができます。

$${\bf y_{i+1} } = {\bf y_{i} } + h {\bf f }(x_i, {\bf y_i})$$

運動方程式の数値解法

前述の様に、自由落下のシミュレーションに用いる運動方程式は、以下のような2階の常微分方程式となっています。 これを連立常微分方程式の解法と、オイラー法を使ってシミュレーションしてみましょう。

まずは、運動方程式を連立常微分方程式として扱います。

$$F = m a(t) = m \frac{d^2 x(t)}{dt^2}$$

この式を、速度$v_f$と位置$x_f$に関する連立1次常微分方程式とみなし、以下の2つの式を導きます。

$$\begin{aligned} \frac{d v_f}{dt} &= g \\ \frac{d x_f}{dt} &= v_f \end{aligned}$$

これを、時刻$t$を刻み幅 $h$で変化させながら、順に求めていきます。 ここでは、常微分方程式を解くのにオイラー法を用いています。

$$\begin{aligned} v_{f(i+1)} &= v_{f(i)} + g \times h \\ x_{f(i+1)} &= x_{f(i)} + v_{f1} \times h \end{aligned}$$

この様にして導き出した連立1次常微分方程式をプログラムとして実装したものが、前述の自由落下のシミュレーションになります。 具体的には、ハイライトされた部分が運動方程式をオイラー法で実装した部分になります。

ポテンシャルに基づく2次元運動

ここからは2次元平面を移動する質点の運動をシュミレーションします¹。 具体的には、惑星の重力によって与えられたポテンシャルの中を運動する質点を扱います。

• 軌道設計への挑戦 - 宇宙科学研究所
• 高校数学・物理で学ぶスイングバイ軌道の作り方①　〜軌道伝播をしてみよう〜
• Watch how the James Webb Space Telescope will orbit the sun _ Science News - YouTube

質点であれば、2次元運動もニュートンの運動方程式で記述できます。 先程の式を2次元運動に拡張し、${\bf F, \alpha, v, x}$をベクトル量とします。 2次元運動シュミレーションするためには、この式を数値的に解くことになります。 それぞれ、X軸およびY軸の各成分について数値計算します。

$${\bf F} = m {\bf \alpha} = m \frac{d{\bf v(t)}}{dt} = m \frac{d^2 {\bf x(t)}}{dt^2}$$

シュミレーション対象とする2次元平面は、以下の図のように設定します。 この平面内には、質点に対して力を及ぼすものとしてQ1とQ2が置いてあります。 Q1とQ2は質点に対し引力や斥力を及ぼすもので、惑星のようなものだと考えてください。 ここでは、Q1とQ2はこの2次元平面の中で固定されており、質点は動き回れるものとします。 また、質点は単位となる質量を持っており、Q1とQ2から力を受けるものとします。

図：2次元平面に、2つの惑星と1つの質点が置かれた空間

この平面の中を質点が運動すると、質点はQ1とQ2から引力や斥力を受けます。 受ける力は、Q1とQ2からの距離 $r$の2乗に反比例するとします。 Q1とQ2の質点に対する影響力の強さを$q_1, q_2$とすれば、質点に働く力の大きさ$|F_{Q1}|, |F_{Q2}|$は、以下のようになります。ただし、kは適当な定数です。また、質点の持っている力は-1であるとします。

$$\begin{aligned} |F_{Q1}| &= \frac{k q_1}{r^2} \\ |F_{Q2}| &= \frac{k q_2}{r^2} \end{aligned}$$

この設定は、固定された複数の惑星が存在する平面で、ある単位質量を持った質点が平面内で運動を行う際のシミュレーションを行うのと同等になります。 これを簡単な質点の軌道シミュレーションとします。

2次元運動シミュレーション

前述の2次元に拡張したニュートンの運動方程式を、数値的に解くことを考えます。 まず1つの惑星から質点に働く力を、x軸方向の力とy軸方向の力に分けて、次のように考えます。

$$\begin{aligned} F_{q_x} &= \frac{x_x - q_x}{r} \times {\bf |F_{q}|} = \frac{x_x - q_x}{r} \times \frac{kq}{r^2} \\ F_{q_y} &= \frac{x_y - q_y}{r} \times {\bf |F_{q}|} = \frac{x_y - q_y}{r} \times \frac{kq}{r^2} \\ r^2 &= (x_x - q_x)^2 + (x_y - q_y)^2 \end{aligned}$$

ここで、係数k=1、質点の質量m=1とすると、加速度${\bf \alpha = (\alpha_x, \alpha_y)}$となります。 惑星が2個以上存在する場合は、以下の式を惑星の数だけ重ね合わせて、最終的な加速度を求めます。

$$\begin{aligned} \alpha_x &= \frac{x_x - q_x}{r^3} \times q \\ \alpha_y &= \frac{x_y - q_y}{r^3} \times q \end{aligned}$$

ここまでで、2次元運動シミュレーションを数値的に解く準備ができました。

ここからは、シミュレーションを行うステップについて考えます。

1. 以下の各変数について、適当な初期値を決める。

   • 時刻の刻み幅$h$
   • 質点の初速度$(v_{x_0}, v_{y_0})$
   • 質点の初期位置$(x_{x_0}, x_{y_0})$
   • 惑星の個数$no_{fq}$
   • すべての惑星$q_i$の位置$(q_{x_i}, q_{y_i})$とその大きさ$q_{i_q}$

2. オイラー法により、次のステップの速度$v_{next}$を求める。 ただし$(\alpha_x, \alpha_y)$は、上記で求めた値を全ての惑星について加算する。

$$\begin{aligned} v_{next_x} &= v_x + \alpha_x \times h \\ v_{next_y} &= v_y + \alpha_y \times h \end{aligned}$$

1. オイラー法により、次のステップの位置$x_{next}$を求める。

$$\begin{aligned} x_{next_x} &= x_x + v_{next_x} \times h \\ x_{next_y} &= x_y + v_{next_y} \times h \end{aligned}$$

1. 2と3を何かしらの終了条件を満たすまで繰り返す。

サンプルコード

2次元運動シミュレーションを行うJuliaプログラム

using Plots

# 定数
const Q = [((0.0, 0.0), 10.0), ((5.0, -5.0), 5.0)]  # 惑星の位置と値
const TIMELIMIT = 30.0                              # シミュレーション打ち切り時刻
const RLIMIT = 0.1                                  # 距離rの最低値
const H = 0.01                                      # 時刻の刻み幅

# メイン実行部
t = 0.0  # 時刻t

# 係数の入力
println("初速度v0xを入力してください:")
vx = parse(Float32, readline())
println("初速度v0yを入力してください:")
vy = parse(Float32, readline())
println("初期位置xを入力してください:")
x  = parse(Float32, readline())
println("初期位置yを入力してください:")
y  = parse(Float32, readline())

# 現在時刻と現在の位置
# println(t, ", ", x, ", ", y, ", ", vx, ", ", vy)

# 現在位置などを追加
xlist = [x]
ylist = [y]

# 2次元運動の計算
while t < TIMELIMIT   # 打ち切り時間まで計算
    global t += H            # 時刻の更新
    rmin = Inf        # 距離の最小値を初期化
    for qi = 1:length(Q)
        rx = Q[qi][1][1] - x  # 距離rxの計算
        ry = Q[qi][1][2] - y  # 距離ryの計算
        r = sqrt(rx * rx + ry * ry)    # 距離rの計算
        if r < rmin
            rmin = r  # 距離の最小値を更新
                end
        global vx += (rx / r / r / r * Q[qi][2]) * H  # 速度vxの計算
        global vy += (ry / r / r / r * Q[qi][2]) * H  # 速度vyの計算
        end
        global x += vx * H  # 位置xの計算
        global y += vy * H  # 位置yの計算
        # 現在時刻と現在の位置
        # println(t, ", ", x, ", ", y, ", ", vx, ", ", vy)
        push!(xlist, x)
        push!(ylist, y)
        # 現在時刻と現在の位置
        if rmin < RLIMIT
        break  # 近づきすぎたら終了
        end
end

simulation_period = length(xlist)
println("Simulation period: ", simulation_period)

# 計算結果の可視化
# グラフ
margin=5
x_graph_range = (findmin(xlist)[1]-margin, findmax(xlist)[1]+margin)
y_graph_range = (findmin(ylist)[1]-margin, findmax(ylist)[1]+margin)
plot(xlist[:], ylist[:], aspect_ratio=1.0, ylim=y_graph_range, xlim=x_graph_range)
for qi = 1:length(Q)
    plot!(Q[:][qi][1], seriestype=:scatter, label="", markersize=Q[:][qi][2])
end
savefig("orbit_calculation.png")

# アニメーション
println("Making animation")
anim = Animation()
animation_step = trunc(Int, simulation_period / 200)
anim = @animate for t in 1:animation_step:simulation_period
        my_title = "Initial position ($(xlist[1]), $(ylist[1]))"
        plot(aspect_ratio=1.0, ylim=y_graph_range, xlim=x_graph_range)
        for qi = 1:length(Q)
                plot!(Q[:][qi][1], seriestype=:scatter, label="", markersize=Q[:][qi][2])
        end
        plot!(xlist[1:t], ylist[1:t], label="Time: $(t)", title=my_title)
        plot!(xlist[t:t], ylist[t:t], label="", markershape=:circle, markersize=5)
end
println("Generating animation")
gif(anim, "orbit_calculation.gif", fps=25)

実行例

julia> include("electoric_field.jl")
初速度v0xを入力してください:
-2
初速度v0yを入力してください:
1
初期位置xを入力してください:
2
初期位置yを入力してください:
2

Simulation period: 3001
Making animation
Generating animation
[ Info: Saved animation to /root/electoric_field.gif
Plots.AnimatedGif("/root/electoric_field.gif")

図：初速度v0x=-2, 初速度v0y=1, 初期位置x=2, 初期位置y=2としたときの計算結果

参考文献

---

1. [教科書P.35 2.2節 「ポテンシャルに基づく2次元運動」] ↩