同様に考えて、Newton 運動方程式のエネルギー保存もうまくやれないだろうか

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

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

として、この系についての保存量である全エネルギー(= 運動エネルギー + 位置エネルギー)の保存について考えて、それを保つように、数値計算式を構成してみよう．

課題: 上の系の保存量を導出しよう

課題: 次の数値計算式を自分で導出しよう

うまい数値計算式

まずは $v = dh/dt$ という従属変数を一つ新たに導入して、問題全体の(見かけの)微分階数を一つ下げよう．

そして、(定数倍を違いを気にせずに) 全エネルギーを例えば

$$ G(h,v) = \frac{1}{2} v^2 + \frac{k}{2m} h^2 + gh $$

としよう(最初の項が運動エネルギー相当、後者二つの項が位置エネルギー相当になっている)．

そして、ある時点 $t = n\Delta t$ での数値解 $h, v$ に対してエネルギーの数値近似値をやはりそっくりな式

$$ G_d(h,v) = \frac{1}{2} v^2 + \frac{k}{2m} h^2 + gh $$

であると定義しよう．

すると例えば、次のような数値計算式で、次のタイムステップ $t = (n+1)\Delta t$ での近似値 $h^+, v^+$ を計算することにすると、$G_d$ は変化しない． つまり、$G_d(h^+, v^+) = G_d(h, v)$ となる．

$$ \left\{\begin{array}{rcl} \displaystyle \frac{h^+ - h}{\Delta t} & = & \displaystyle \frac{v^+ + v}{2} \\ \displaystyle \frac{v^+ - v}{\Delta t} & = & \displaystyle -\frac{k}{m} \left( \frac{h^+ + h}{2} \right) - g\\ \end{array} \right. $$

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

$$ \left\{\begin{array}{rcl} \displaystyle h^+ - \tau v^+ & = & \displaystyle h + \tau v , \\ \displaystyle \frac{k}{m}\tau h^+ + v^+ & = & \displaystyle -\frac{k}{m}\tau h + v -2\tau g\\ \end{array} \right. $$

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

$$ \left(\begin{array}{cc} 1 & -\tau \\ \displaystyle \frac{k}{m}\tau & 1 \end{array}\right) \left(\begin{array}{c} h^+ \\ v^+ \end{array}\right) = \left(\begin{array}{cc} 1 & \tau \\ \displaystyle -\frac{k}{m}\tau & 1 \end{array}\right) \left(\begin{array}{c} h \\ v \end{array}\right) - \left(\begin{array}{c} 0 \\ 2\tau g \end{array}\right) $$

となるので、これを

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

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

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

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

$$ A \left(\begin{array}{c} h^+ \\ v^+ \end{array}\right) = B \left(\begin{array}{c} h \\ v \end{array}\right) - b $$

というシンプルな連立一次方程式で書けることになり、時間ステップを進めるためにはこれを解けばよいことがわかる．

実際にプログラミングして、解いてみよう

m, k, g = 1.0, 1.0, 9.8 # 定数の設定

(1.0, 1.0, 9.8)

Δt = 0.1
τ = Δt/2   # これも設定してしまおう

0.05

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

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

2×2 Array{Float64,2}:
 1.0   -0.05
 0.05   1.0 

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

2×2 Array{Float64,2}:
  1.0   0.05
 -0.05  1.0 

b = [0, 2τ*g]

2-element Array{Float64,1}:
 0.0 
 0.98

すると、数値計算式はいきなり次のように書けてしまう．

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

onestep (generic function with 1 method)

初期値を設定して、

ini_hv = [0, 1.0]

2-element Array{Float64,1}:
 0.0
 1.0

試しに動かしてみよう．ちなみに、この $\Delta t$ はかなり大きくて、値は結構動く．

onestep(ini_hv)

2-element Array{Float64,1}:
 0.0508728
 0.0174564

さて、500ステップほど動かしてみようか．

hv = copy(ini_hv)
sq_hv = copy(hv)

for n in 1:500
    hv = onestep(hv)
    sq_hv = hcat(sq_hv, hv)
end

早速、グラフにしてみよう．

using Plots
gr()

Plots.GRBackend()

plot(sq_hv')

<?xml version="1.0" encoding="utf-8"?>

保存量の式を定義して、

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

invariant (generic function with 1 method)

今の数値解に適用してみる．

sq_invariant = [ invariant(sq_hv[1,n], sq_hv[2,n]) for n in 1:501 ]

501-element Array{Float64,1}:
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 ⋮  
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5
 0.5

ほぼ 0.5 のようなので、それからのずれをプロットしてみよう．

plot(sq_invariant - 0.5)

<?xml version="1.0" encoding="utf-8"?>

ほとんど完璧に保たれていることがよくわかる．

では、相図も描いてみようではないか．

sq_h, sq_v = sq_hv[1,:], sq_hv[2,:]

([0.0, 0.0508728, 0.00348257, -0.141698, -0.383221, -0.718677, -1.14472, -1.6571, -2.2507, -2.91961  …  -7.24117, -6.30503, -5.40375, -4.54633, -3.74131, -2.99673, -2.32, -1.7179, -1.19641, -0.760744], [1.0, 0.0174564, -0.965261, -1.93835, -2.8921, -3.81701, -4.70384, -5.54375, -6.32836, -7.04984  …  9.51275, 9.21006, 8.8155, 8.333, 7.76738, 7.12429, 6.41012, 5.63202, 4.79773, 3.91559])

plot(sq_h, sq_v, aspect_ratio = 1)

<?xml version="1.0" encoding="utf-8"?>

(理論通りの)円になっていることが見てとれる．完璧だ．