二階微分はどうする: Newton 運動方程式

時間方向に二階微分がかかっている時間発展常微分方程式

これまでは

\[ \frac{du}{dt} = F(u,t) \]

というタイプの、時間方向の微分が一階しかない時間発展問題の常微分方程式を扱ってきた．

しかしもちろんのことながら、われわれが取り組みたい問題はこういうものばかりではなく、例えば

\[ \frac{d^2 u}{dt^2} = F(u,\frac{du}{dt}, t) \]

といったような、時間方向の微分が二階の問題や、三階、四階といった問題も考えうる．

こういう場合にどうしたらよいか、それを学んでみよう．

いわゆる Newton 運動方程式

時間方向に二階微分が入る時間発展問題の常微分方程式といえば、高校時代に学ぶ Newton の運動方程式がまずは挙げられよう． そこで、ここではそれを扱ってみよう．

具体的には、バネでぶら下げられている重りの挙動を例として考えてみよう． バネ定数を $k$, 重りの重さを $m$, 重力定数を $g$ とすれば、重りの高さ位置 $h(t)$ $(t: $時間$)$ についての Newton の運動方程式は

\[ m \frac{d^2 h}{dt^2} = -kh -mg \]

となる． むろん、解をきちんと考えるためにはさらに初期値等の情報が必要で、例えば

\[ \left\{ \begin{array}{ccl} h(0) & = & h_{ini}, \cr \displaystyle \frac{dh}{dt}(0) & = & v_{ini} \end{array} \right. \]

などという情報を付与しないといけない．

さて、これを例に、しばらく考察を進めてみよう．

(その前に) 保存量

計算について考える前に、この系の保存量を考えてみよう． これは高校の物理で学んだ通りで、全エネルギー = 運動エネルギー ＋ 位置エネルギー が保存されるはずだ．

そして、運動エネルギーは

\[ \frac{1}{2} m \left( \frac{dh}{dt} \right)^2 \]

と表され、位置エネルギー = 力の積分 on 逆らって移動した距離 は

\[ \int_0^h (kx + mg) dx = \frac{1}{2} kh^2 + mgh \]

となることは覚えているだろうか？ まあ思い出してもらうとして、全エネルギーは

\[ E(h(t)) = \frac{1}{2} m \left( \frac{dh}{dt} \right)^2 + \frac{1}{2} kh^2 + mgh \]

となるわけなので、これがこの系の保存量である…はずだな．

高階微分が含まれる時間発展問題へのアプローチ 1a

高階微分が含まれる時間発展問題へのアプローチはいくつか考えられるが、いちばんストレートなのは「方程式はそのままに、そのまま離散近似する」という手法だろう．

例えば、先の問題に対して 近似解を $h_n \cong h(n\Delta t)$, $n = 0,1,2, \cdots$ として

\[ m \left( \frac{h_{n+1} - 2h_n + h_{n-1} }{\Delta t^2} \right) = -k h_n + mg \]

という近似式を提案する、などだ． もちろん初期値等の付与情報の離散化も必要で、この場合の安直な離散化は例えば次のようになるだろう．

\[ \left\{ \begin{array}{ccl} h_0 & = & h_{ini}, \cr \displaystyle \frac{h_1 - h_0}{\Delta t} & = & v_{ini} \end{array} \right. \]

tips: これはあくまで「例」であって、良い離散化かどうかは全く考慮していないことに注意しよう．実際は、たぶん、相当「悪い」離散化になっているだろう.

高階微分が含まれる時間発展問題へのアプローチ 1b

アプローチ 1a の改善版として考えられるアプローチは、「方程式はそのままに、離散化は工夫して」という手法だろう． これは悪くないアプローチだが、説明が少し必要だ．そこで詳しくは次回以降にするとしよう．

高階微分が含まれる時間発展問題へのアプローチ 2a

次に考えられるアプローチは、「従属変数を新たに必要なだけ導入して方程式を一階常微分方程式に変形し、その常微分方程式に対する汎用解法を適用する」という手法で王道といえるものだ．

話は簡単で、従属変数を新たに一つ導入すると全体の微分階数が 1階下げられるので(方程式は連立になる), これを必要なだけ繰り返せば良い． どういうことか、上の例で示してみよう．

まず、従属変数として、

\[ v(t) := \frac{dh}{dt} \]

を新たに導入する．すると元の方程式そのものは

\[ m\frac{d v}{dt} = -kh -mg \]

と書き換えることができて微分階数を一つ下げられる． よって全体としては(連立ではあるが)一階常微分方程式である次の形

\[ \left\{ \begin{array}{rcl} \displaystyle \frac{d h}{dt} & = & v, \cr\cr \displaystyle m \frac{d v}{dt} & = & -kh -mg \end{array} \right. \]

に帰着させられる、というわけだ．

あとはこの問題に対して、Euler 法でも時間対称離散化でも Runge-Kutta 法でもなんでも好きな汎用解法を適用すれば良い、というわけだ．

tips: この場合の全エネルギーは

\[ E(h,v) = \frac{1}{2} m v^2 + \frac{1}{2} kh^2 + mgh \]

と表現される．

実際にやってみよう: Euler 法

Euler 法だと、近似値を $h_n \cong h(n\Delta t)$, $v_n \cong v(n\Delta t)$ として, まあ普通は近似式は次のようになるかな．

\[ \left\{ \begin{array}{rcl} \displaystyle \frac{h_{n+1} - h_n}{\Delta t} & = & v_n, \cr\cr \displaystyle m \left( \frac{v_{n+1} - v_n}{\Delta t} \right) & = & -k h_n -mg . \end{array} \right. \]

書き換えると次のような感じだ．

\[ \left\{ \begin{array}{rcl} h_{n+1} & = & h_n + \Delta t \, v_n, \cr\cr v_{n+1} & = & \displaystyle v_n + \Delta t \, (- \frac{k}{m} h_n -g), \end{array} \right. \]

あとはプログラミングだ．

定数等をまず定義して…

1
2
3

m, k, g = 1.0, 1.0, 9.8 # 定数の設定
Δt = 0.1
h_ini, v_ini = 0.0, 1.0

Euler 法の 1step を関数としていつものように定義しよう．

1
2
3
4
5
6

function euler(h,v)
    h_new = h + Δt * v  
    v_new = v + Δt * (-(k/m)*h - g)

    return h_new, v_new
end

エネルギーの計算式も定義しておいて、

1

energy(h,v) = (m*v^2)/2 + (k*h^2)/2 + m*g*h

早速動かしてみようではないか！

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

h,v = h_ini, v_ini
h_sq = [h]
v_sq = [v]
E_sq = [ energy(h,v) ]

for i in 1:500
    h,v = euler(h,v)
    push!(h_sq, h)
    push!(v_sq, v)
    push!(E_sq, energy(h,v))
end

近似解をプロットしてみよう． まず, いつもとまったく同じ準備をする．

1
2

using Plots
gr()

そして、時間軸方向もデータを丁寧に用意して、

1

t_sq =  Δt * [ 0:500 ]

1-element Array{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},1}: 0.0:0.1:50.0

早速プロットしよう．

1

plot(t_sq, h_sq, xaxis = "time t", yaxis = "height", legend = false)

う～ん、バネにぶら下がった重りがびよんびよんしてるだけのはずなのに、振幅が大きくなっているってのはおかしいな…

エネルギーが保存されているか、プロットしてみよう．

1

plot(t_sq, E_sq, xaxis = "time t", yaxis = "Energy", legend = false)

う～ん、エネルギーがすごい勢いで「増えている」．こりゃだめだ．

念のために相図も描いておこう．

1
2
3

plot(h_sq, v_sq)

plot!((h_ini, v_ini), marker = :circle, label = "initial pt")

本来は保存系であることが明白なのに、数値解は一周してもとに戻ってこない．やはりまずいねえ…

Runge-Kutta 法なら？

というわけで、アプローチ 2a の「汎用解法」の部分を Runge-Kutta 法に切り替えてみよう． つまり、もとの微分方程式を、左辺 = 一階微分 の形式に書き直した

\[ \left\{ \begin{array}{rcl} \displaystyle \frac{d h}{dt} & = & v, \cr\cr \displaystyle \frac{d v}{dt} & = & \displaystyle -\frac{k}{m} h -g \end{array} \right. \]

に対して Runge-Kutta 法を適用するわけだ．

あとはいつもどおりなので、さくさくやってみよう．

まず、微分方程式の右辺を定義する．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

function r(hv) 

    h,v = hv

    eq1 = v
    eq2 = - (k/m)*h - g

    return [ eq1, eq2 ]

end

そしてそれを元に Runge-Kutta 法を適用したものを作ってしまう．

1
2
3
4
5
6
7

function RK(hv)
    r1 = r(hv)
    r2 = r(hv + Δt/2 * r1)
    r3 = r(hv + Δt/2 * r2)
    r4 = r(hv + Δt * r3)
    return hv + (Δt/6) * (r1 + 2*r2 + 2*r3 + r4)
end

エネルギー計算式も呼び出しやすい形式にしておこう．

1

energy(hv) = energy(hv[1], hv[2])

あとは実際に計算するだけだ．

1
2
3
4
5
6
7
8
9

hv = [ h_ini, v_ini ]
hv_sq = copy(hv)
E_sq = [ energy(hv) ]

for i in 1:500
    hv = RK(hv)
    hv_sq = hcat(hv_sq, hv)
    push!(E_sq, energy(hv))
end

プロットしやすいように結果を切り出して…

1
2

h_sq = hv_sq[1, :]
v_sq = hv_sq[2, :]

プロットする．

1

plot(t_sq, h_sq, xaxis = "time t", yaxis = "height", legend = false)

おお、振動が一定になったぞ！

エネルギーもプロットしよう．

1

plot(t_sq, E_sq, xaxis = "time t", yaxis = "Energy", legend = false)

グラフとしては下がっているが、縦軸の数字を読む限り、だいぶ一定っぽい．

初期値を引いてプロットして、ズレの形で視覚化しよう．

1

plot(t_sq, E_sq - 0.5, xaxis = "time t", yaxis = "Energy", legend = false)

計算の全過程を通じて、エネルギーがおおよそ 3桁程度の精度で保たれていることがわかる． ん～、まあ悪くない．まあ、もう少し良くてもいいかなあとは思うが．

相図も描いてみよう．

1
2
3

plot(h_sq, v_sq)

plot!((h_ini, v_ini), marker = :circle, label = "initial pt")

ふむ、図の上では一周して戻ってきた、と言ってよさそうだ．

高階微分が含まれる時間発展問題へのアプローチ 2b

次に考えられるアプローチは、「従属変数を新たに必要なだけ導入して方程式を一階常微分方程式に変形し、その常微分方程式に対する 特別な 解法を設計、適用する」という手法で、これもまあ王道といえよう．

この場合の 特別 というのはもちろん文脈依存なので、以下に「エネルギー保存性を再現するという意味で特別な」例を示そう．

まずは、先に示したように、 $v = dh/dt$ という従属変数を一つ新たに導入したときの全エネルギーは

\[ E(h,v) = \frac{1}{2} m v^2 + \frac{1}{2} k h^2 + mgh \]

となるはずだ．

そして、離散化版のエネルギーについては、ある時点 $t = n\Delta t$ での数値解 $h_n, v_n$ に対してやはりそっくりな式

\[ E_d(h_n, v_n) = \frac{1}{2} m v_n^2 + \frac{1}{2} k h_n^2 + mg h_n \]

であると定義しよう．

そして、天下りであるが、例えば、次のような数値計算式で近似値 $h_{n+1}, v_{n+1}$ を計算すると、$E_d$ は変化しない( Why? は課題へ)． つまり、$E_d(h_{n+1}, v_{n+1}) = E_d(h_n, v_n)$ となる．

\[ \left\{ \begin{array}{rcl} \displaystyle \frac{h_{n+1} - h_n}{\Delta t} & = & \displaystyle \frac{v_{n+1} + v_n}{2}, \cr\cr \displaystyle \frac{v_{n+1} - v_n}{\Delta t} & = & \displaystyle -\frac{k}{m} \left( \frac{h_{n+1} + h_n}{2} \right) - g \end{array} \right. \]

こいつは幸いにも線形な式なので、キレイに変形できるはずだ． やってみると、まず、簡単のために $\tau = \Delta/2$ とおいて

\[ \left\{ \begin{array}{rcl} \displaystyle h_{n+1} - \tau v_{n+1} & = & \displaystyle h_n + \tau v , \cr\cr \displaystyle \frac{k}{m}\tau h_{n+1} + v_{n+1} & = & \displaystyle -\frac{k}{m}\tau h + v_n -2\tau g \end{array} \right. \]

となり、これは線形代数の記号では

\[ \left( \begin{array}{cc} 1 & -\tau \cr \displaystyle \frac{k}{m}\tau & 1 \end{array} \right) % \left( \begin{array}{c} h_{n+1} \cr v_{n+1} \end{array} \right) = \left( \begin{array}{cc} 1 & \tau \cr \displaystyle -\frac{k}{m}\tau & 1 \end{array} \right) \left( \begin{array}{c} h_n \cr v_n \end{array} \right) - \left( \begin{array}{c} 0 \cr 2\tau g \end{array} \right) \]

となる. よって、定数部分をわかりやすく

\[ A = \left(\begin{array}{cc} 1 & -\tau \cr \displaystyle \frac{k}{m}\tau & 1 \end{array}\right), \]

\[ B = \left(\begin{array}{cc} 1 & \tau \cr \displaystyle -\frac{k}{m}\tau & 1 \end{array}\right), \]

\[ b = \left(\begin{array}{c} 0 \cr 2\tau g \end{array}\right) \]

とおくと、この数値計算式は

\[ \left( \begin{array}{c} h_{n+1} \cr v_{n+1} \end{array} \right) = A^{-1} \left\{ B \left(\begin{array}{c} h_n \cr v_n \end{array}\right) - b \right\} \]

と書けることになり、シンプルに解けることがわかる．

さてあとはプログラムだ．

まず適当に定義をして、

1

τ = Δt/2

続けて、$A, B, b$ を定義してしまおう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

A = [
  1 -τ
  (k*τ)/m 1
    ]

B = [
  1 τ
  -(k*τ)/m 1
    ]

b = [0, 2τ*g]

すると、この近似計算式の 1step 分はいきなり次のように書けてしまう．

1

onestep(hv) = A \ (B*hv - b)

onestep (generic function with 1 method)

あとはいつものように計算するだけだ．

1
2
3
4
5
6
7
8
9

hv = [ h_ini, v_ini ]
hv_sq = copy(hv)
E_sq = [ energy(hv) ]

for i in 1:500
    hv = onestep(hv)
    hv_sq = hcat(hv_sq, hv)
    push!(E_sq, energy(hv))
end

プロット用にデータを切り出して、

1
2

h_sq = hv_sq[1, :]
v_sq = hv_sq[2, :]

まずは近似解そのものをプロット．

1

plot(t_sq, h_sq, xaxis = "time t", yaxis = "height", legend = false)

うむ．うまくいっているようだな．

ではエネルギーも．

1

plot(t_sq, E_sq, xaxis = "time t", yaxis = "Energy", legend = false)

GKS: Rectangle definition is invalid in routine SET_WINDOW GKS: Rectangle definition is invalid in routine SET_WINDOW GKS: Rectangle definition is invalid in routine SET_WINDOW GKS: Rectangle definition is invalid in routine SET_WINDOW

値の変化が乏しすぎてプロットに失敗するようだ． というわけで初期値からのズレ、をプロットさせよう．

1

plot(t_sq, E_sq -0.5, xaxis = "time t", yaxis = "Energy", legend = false)

エネルギー値は 11桁ほど保たれていることがわかる．ふむ、これは良い．

一応、相図もプロットしておこう．目で見てずれがわかるとも思えないが．

1
2
3

plot(h_sq, v_sq)

plot!((h_ini, v_ini), marker = :circle, label = "initial pt")

というわけで、最後の「謎」の近似式だとだいぶ良いことがわかる． さて、これはいったいどういう仕組なのだろう？

Report No.07 Newton 運動方程式に対して…

1. 自力でエネルギー $E(h)$ がもとの運動方程式に従うと保存されることを確認しよう．
   
   $h(t)$ が元の二階常微分方程式である運動方程式の厳密解であるとき、 $h$ のみで書かれた全エネルギー $E(h)$ の時間微分$\displaystyle \frac{dE}{dt}$ を自分で手で計算し、確かにこれがゼロになることを確認しよう
   
   
2. 最後の解法が、なぜこんなによくエネルギー $E(h,v)$ を保つのか、自分なりに説明してみよ．