今回は、ほぼ理系大学で多分教わるだろう二階微分方程式を中心とした数値解析と解析解の比較記事である。

TL;DR

二階微分方程式はEuler法でも$h^2$オーダーで表現される。※ただし、おそらく$x$の寄与がすくないと$h$のオーダになるはず。
二階線形非斉次微分方程式の解き方知りたければ本記事みてね
FDTD法はnumpy使って書くと比較的簡単にかけるよ

常微分方程式の数値解: Euler法とRunge-Kutta法

微分を差分形式に直して逐次更新することで微分方程式の解となる関数$x(t)$を求めよう。
これには、ある程度小さいステップ刻み$h$を用いて

$$x(t_n=nh)=x_n$$

とすれば、微分は差分として近似できる。

$$\frac{dx(t)}{dt}=\lim_{h\rightarrow 0}\frac{x(t+h)-x(t)}{h}$$

$$\approx \frac{x(nh+h)-x(nh)}{h}$$

$$=\frac{x((n+1)h)-x(nh)}{h}$$

$$=\frac{x_{n+1}-x_n}{h}$$

$$\therefore x_{n+1}\approx x_n+h\frac{dx}{dt}$$

$$\frac{dx}{dt}=v_n$$とすると、

$$x_{n+1}=x_n+hv_n$$

これがオイラー法の更新式である。

一階微分方程式とRunge-Kutta法

ここで、微分方程式が
$$\frac{dx}{dt}+x=6t$$
だったとしよう。これを$v_n$の更新式に使う。

$$v_n=6nh-x_n$$

まとめると、オイラー法の更新式は、

$$v_n=6nh-x_n$$
$$x_{n+1}=x_n+hv_n$$

となる。
ところで、この式は刻み$h$のステップを細かく刻まないと制度がでないことが知られているし、刻みすぎると丸め誤差が影響して制度がでない。

精度を改善するために式を眺めれば、$x_{n+1}$の更新に使う勾配は、実際には以下となるべきであるとわかる。

$$x_{n+1}=x_n+h\int_0^1 v(t+hs)ds$$

ここで、シンプソン則を使えば

$$\int_0^1v(t+hs)ds=\frac{v(t)+4v(t+\frac{h}{2})+v(t+h)}{6}$$

となり、これを利用したものがRunge-Kutta法という。

$$v_n=6nh-x_n$$
$$v_{n+1/2}^{(1)}=6(n+\frac{1}{2})h-(x_n+\frac{h}{2}v_n)$$
$$$$
$$v_{n+1/2}^{(2)}=6(n+\frac{1}{2})h-(x_n+\frac{h}{2}v_{n+1/2}^{(1)})$$
$$v_{n+1}=6(n+1)-(x_n+hv_{n+1/2}^{(2)})$$

$$v_{n+1/2}=\frac{v_{n+1/2}^{(1)}+v_{n+1/2}^{(2)}}{2}$$

つまり、更新式(Runge-Kutra法)は

$$x_{n+1}=x_n+\frac{k_1+2k_2+2k_3+k_4}{6}h$$

$$k_1=f(nh,x_n)$$
$$k_2=f(nh+1/2,x_n+\frac{k_1h}{2})$$
$$k_3=f(nh+1/2,x_n+\frac{k_2h}{2})$$
$$k_4=f(nh+1,x_n+k4h)$$

ここで、$f(t,x_n)$は微分方程式による勾配の評価式である。今回の場合は、微分方程式を勾配について整理した

$$\frac{dx}{dt}=6t-x$$

とかけるので、

$$f(t,x)=6t-x$$

となる。あとは、初期値$x(t=0)$を与えてやれば$t>0$における関数が求まる。

二階微分方程式の数値解

二階微分の際もアプローチは同様である。

微分方程式を以下の形に変形する。

$$\frac{d^2x}{dt^2}=f(t,x,\frac{dx}{dt})$$

初期条件$x_0, v_0$を定める
微分方程式から二階微分係数$a_n$の定義式を求める
$x_{n+1},v_{n+1}$を求める

やってみよう。

$$a_n=f(t,x_n,v_n)$$
$$v_{n+1/2}=v_n+\frac{a_nh}{2}$$
$$v_{n+1/2}=v_n+h\frac{k_1+2k_2+2k_3+k_4}{12}$$
$$v_{n+1}=v_{n+1/2}+h\frac{k_5+2k_6+2k_7+k_8}{12}$$

$$k_1=f(nh,x_n,v_n)$$
$$k_2=f(nh+h/2,x_n+v_nh/2+k_1h^2/8,v_n+k_1h/2)$$
$$k_3=f(nh+h/2,x_n+v_nh/2+k_2h^2/8,x_n+k_2h/2)$$
$$k_4=f(nh+h,x_n+v_nh+k_3h^2/2,v_n+k_3h)$$

$$a_{n+1/2}=\frac{k_2+k_3}{2}$$

あとは、これにシンプソン則を使って$v_{n+1},x_{n+1}$を求める流れになる。
ただし、$x_{n+1}$を求めるために$v_{n+1/2$を見積もる必要がある。これもシンプソン則を使って

$$\int_n^{n+1/2} P(x)dx=\frac{(b-a)(5k_1+4k_2+4k_3-k_4)}{24}$$

とできるので、

$$v_{n+1/2}=v_n+h\frac{5k_1+4k_2+4k_3-k_4}{24}$$

$$v_{n+1}=v_n+h\frac{k_1+2k_2+2k_3+k_4}{6}$$

であり、

$$x_{n+1}=x_n+h\frac{v_n+4v_{n+1/2}+v_{n+1}}{6}$$

$$=x_n+v_nh+h^2\frac{k_1+k_2+k_3}{6}$$

まとめると、二階微分方程式の数値解をRunge-Kutta法利用すると以下の更新式になる。

$$x_{n+1}=x_n+v_nh+\frac{k_1+k_2+k_3}{6}h^2$$

$$v_{n+1}=v_n+h\frac{k_1+2k_2+2k_3+k_4}{6}$$

$$k_1=f(t_n,x_n,v_n)$$
$$k_2=f(t_{n+1/2},x_n+\frac{v_nh}{2}+\frac{k_1h^2}{8},v_n+\frac{k_1h}{2})$$
$$k_3=f(t_{n+1/2},x_n+\frac{v_nh}{2}+\frac{k_2h^2}{8},v_n+\frac{k_2h}{2})$$
$$k_4=f(t_{n+1},x_n+v_nh+\frac{k_3h^2}{2},v_n+k_3h)$$

数値解と解析解の相互検証(検算)

ためしに、簡単な二階微分方程式の解を解析解と比較して精度を確認してみる。

$$\frac{d^2x}{dt^2}=-x(1+\lambda^2)-\lambda\frac{dx}{dt}$$

この解は、

$$x=\exp(-\frac{\lambda}{2}t)(A\cos t+B\sin t)$$

である。
初期条件$x_0,\dot x_0$($\dot x=\frac{dt}{dt}$)が与えられたとき

$$A=x_0$$
$$B=\dot x_0+\frac{x}{2}$$

あとは、これを数値解と比較することで、ここまでの計算の検算としつつ、数値解の精度も見ていこう。

pythonではnumpyを使えば比較的簡単に組めるが、逐次法は並列演算的ではないためnumpyの高速な計算を利用できないことは残念だ。

import numpy as np

# Euler.
def solver_Euler(t, x, dx, f, h):
	ddx=f(t, x, dx)
	dx += ddx*h
	x+=dx*h
	return x,dx

# Runge-Kutta.
def stepf(f, t, x, dx, ddx, dt):
	return f(t+dt, x+dx*dt+ddx*dt*dt/2, dx+ddx*dt)

def solver_rungeKutta(t, x, dx, f, h):
	k1=f(t, x, dx)
	k2=stepf(f, t, x, dx, k1, h/2)
	k3=stepf(f, t, x, dx, k2, h/2)
	k4=stepf(f, t, x, dx, k3, h)
	x=x+dx*h+h*h*(k1+k2+k3)/6
	dx=dx+h*(k1+2*k2+2*k3+k4)/6
	return x, dx

# define lambda const.
LAM=0*-1/20

# define differential equation.
def test_fn(t, x, dx):
	return -x*(1+LAM**2/4)-dx*LAM

# define analytical solution.
def test_gt(t, x0, dx0):
	A=x0
	B=dx0+A*LAM/2
	return np.exp(-t*LAM/2)*(A*np.cos(t)+B*np.sin(t))

# calc solutions.
def get_curve(fn, gt_fn, x0, dx0, h=0.01, maxt=12):
	x=x0
	dx=dx0
	t=np.arange(0, maxt, 0.1)
	gt_x=gt_fn(t, x0, dx0)
	et=np.arange(0, maxt, h)
	# numerical solutions.
	eux=[(x,dx)]
	rkx=[(x,dx)]
	for ti in et[1:]:
		x,dx=solver_Euler(ti, x,dx, fn, h)
		eux.append([x,dx])
		x,dx=solver_rungeKutta(ti, *rkx[-1], fn, h)
		rkx.append([x,dx])
	eux=np.array(eux)
	rkx=np.array(rkx)
	x_eu=eux[:,0]
	x_rk=rkx[:,0]
	return [(t, gt_x), (et, x_eu), (et, x_rk)]

初期条件と数値解のステップを変えて比較してみると以下になる。
なお、青線が解析解(真値)、橙線がEuler法による解、緑線がRunge-Kutta法(二階に拡張)である。なお、以降$\dot x=\frac{dx}{dt}$とする。

$$\lambda=0.0$$

$$\lambda=0.1$$

$$\lambda=-0.05$$

ステップが0.01程度ならeuler法でもほぼ再現性よく、0.8くらいになるとeuler法ではつらい一方Runge-Kutta法は頑張っている。euler法が思った以上に健闘しているが、今回導出した式自体はうまくできてそうだ。
一方、下段のようにステップが1を越えてくると誤差が目立つ。両者とも$h$の高次項が0に収束する前提なので、前提が一致していないことが原因と思える。

Euler法とRunge-Kutta法の性能差が出ているか明確でないので、$\lambda=0.1$のサンプル
で真値からの誤差を$\log_{10}$のスケールで
プロットしてみよう。

だいたいEuler法は$h^2$、Rung-Kutta法は$$h^4$のオーダー付近が誤差になっており、Runge-Kutta法の側は理論値に近いといえる。
一方Euler法は1階微分方程式なら$h$のオーダーの誤差のはずだが、2階微分方程式だと$x$の変化のうち$h$の寄与は二乗であることが関係しているのかもしれない。

二階線形非斉次微分方程式の一般解

もうちょっと複雑な微分方程式の解析解も見てみよう。

$$\frac{dx^2}{dt^2}+k_1\frac{dx}{dt}+k_2x+f(t)=0$$

という形状の方程式は、

$$r^2+k_1r+k2=0$$

とくに、上の二次方程式の解を$r,\bar r$とすると
$$r+\bar r=-k_1, r-\bar r=\sqrt{k_1^2-4k_2}, r\bar r=k2$$

なる$r$を用いて

$$x=C(t)\exp(rt)$$

とかける。C(t)の満たすべき式は、微分方程式の$x$を$C(t)$の式に変形すれば

$$\frac{d^2C}{dt^2}+(2r+k_1)\frac{dC}{dt}+f(t)\exp(-rt)=0$$

$\frac{dC}{dt}$の一階線形非一斉次式になったので、同様に繰り返すと、

$$\frac{dC}{dt}=C_2(t)\exp(-(r-\bar r)t)$$

$$C_2(t)=-\int f(t)\exp(-\bar rt)dt+A$$
$$C(t)=-\int \exp(-(r-\bar r)t_1)\int f(t_2)\exp(-\bar rt_2)dt_2dt_1+A\exp(-(r-\bar r)t)+B$$
$$x=-\exp(rt)\int \exp(-(r-\bar r)t_1)\int f(t_2)\exp(-\bar rt_2)dt_2dt_1+A\exp(\bar rt)+B\exp(rt)$$

第一項の二重積分は以下の部分積分を用いて変形できる。

$$\int \exp(-wt)\int g\exp(wt)dtdt=\frac{1}{w}\left(\int gdt-\int g\exp(wt)dt\exp(-wt)\right)$$

$w=r-\bar r$、$g(t)=f(t)\exp(-rt)$を代入すると二重積分の展開に適用できるので、

$$x=-\frac{1}{r-\bar r}\left(\exp(rt)\int f(t)\exp(-rt)dt-\exp(\bar rt)\int f(t)\exp(-\bar rt)dt\right)+A\exp(\bar rt)+B\exp(rt)$$

検算

AとBの項は斉次方程式の解なので明らかに微分方程式を満たす。確かめたいのは前2項の積分であるので、AとBを一旦0として$x=S-T$として計算していこう。ただし、
$$\frac{dS}{dt}-\frac{dT}{dt}=rS-\bar rT$$
を利用する。

$$\frac{d^2x}{dt^2}+k_1\frac{dx}{dt}=(r^2S-\bar r^2T)-(r+\bar r)(rS-\bar rT)-f(t)$$

ただし、最後の$f(t)$は二階積分項から得られる。

$$\therefore =-r\bar r(S-T)-f(t)=-k_2x-f(t)$$

つまり、最初の微分方程式に戻る。

減衰つき強制振動の運動方程式

二階線形非斉次微分方程式の解析解が分かったので、これを使って振動する外力があり、かつ抵抗もある状態の運動方程式を立てよう。

時間微分は$\dot x$と表現する。

$$\ddot x=-2\sigma\dot x-k^2x+\alpha\exp(wt)$$

$$r,\bar r=-\sigma\pm \sqrt{\sigma^2-k^2}$$

さきほどの解析解に代入していくと、

$$x(t)=\frac{\alpha\exp(wt)}{\bar r-r}(\frac{1}{w-r}-\frac{1}{w-\bar r})+A\exp(rt)+B\exp(\bar rt)$$
$$x(t)=\frac{\alpha\exp(wt)}{w^2+2w\sigma+k^2}+\exp(-\sigma t)\left(A\exp(\sqrt{\sigma^2-k^2}t)+B\exp(-\sqrt{\sigma^2-k^2}t)\right)$$

とくに、減衰なし$\sigma=0$のとき、

$$x(t)=\frac{\alpha\exp(wt)}{w^2+k^2}+A\exp(ikt)+B\exp(-ikt)$$

となり、とくに$w=\pm ik$のときに第一項の分母が無限大になる(共振).

数値解との比較(係数選び)

いま、数値解が計算しやすいように微分方程式内に虚数が出てこないようにしたい。

準備として、先ほどの解を変形しておこう。

$$x(t)=\frac{(w^2+k^2)-2w\sigma}{(w^2+k^2)^2-4w^2\sigma^2}\alpha\exp(wt)+A\exp(rt)+B\exp(\bar rt)$$

さて、外力を以下で考えよう。

$$f(t)=\alpha\cos(wt)=\frac{\alpha}{2}(\exp(iwt)+\exp(-iwt))$$

このとき先ほど準備した式を使えば、

$$\ddot x=-2\sigma\dot x-k^2x+\alpha(4k^2-3\sigma^2)\cos(wt)$$

の一般解は

$$x(t)=\frac{\alpha(4k^2-3\sigma^2)}{(k^2-w^2)^2+4w^2\sigma^2}\left\lbrace(k^2-w^2)\cos(wt)+2w\sigma\sin(wt)\right\rbrace+\exp(-\sigma t)\left\lbrace A\cos(\sqrt{k^2-\sigma^2}t)+B\sin(\sqrt{k^2-\sigma^2}t)\right\rbrace$$

だとわかる。
なお、$w=\sqrt{k^2-\sigma^2}$とするなら

$$(k^2-w^2)^2+4w^2\sigma^2=\sigma^2(4k^2-3\sigma^2)$$
$$k^2-w^2=\sigma^2$$

となる。結局、比較的一般解が分かりやすい外力つき減衰振動の式を逆算的に求めると以下

$$\ddot x=-2\sigma\dot x-(w^2+\sigma^2)x+\alpha(4w^2+\sigma^2)\cos wt$$

とわかったので、これを数値解比較のターゲットにしよう。一般解は

$$x(t)=(\alpha+A\exp(-\sigma t))\cos w t+\left(\frac{2\alpha w}{\sigma}+B\exp(-\sigma t)\right)\sin wt$$

となる。

比較開始

$t=0$を初期条件とすれば、

$$\dot x_{t=0}=w\left(\frac{2\alpha w}{\sigma}+B\right)-\sigma A$$
$$x_{t=0}=\alpha+A$$

つまり、

$$A=x_0-\alpha$$
$$B=\frac{\dot x_0+\sigma A}{w}-\frac{2\alpha w}{\sigma}$$

さて、解析解と数値解と比較しよう。コードはほぼ流用できて、方程式と解析解部分だけ差し替えればいい。

SGM=0.2
OMG=1
ALP=0.12
def test_fn(t, x, dx):
	OMG2=OMG**2
	SGM2=SGM**2
	return -2*SGM*dx-(OMG2+SGM2)*x+ALP*(4*OMG2+SGM2)*np.cos(t)

def test_gt(t, x0, dx0):
	A=x0-ALP
	C=2*ALP*OMG/SGM
	B=(dx0-SGM*A)/OMG-C
	return (ALP+A*np.exp(-SGM*t))*np.cos(OMG*t)+(C+B*np.exp(-SGM*t))*np.sin(OMG*t)

結果。初期条件が外力と位相がずれていれば位相が修正され、位相がそろっているが大きい振動からだと外力の振動幅まで減衰する。大体理屈通りである。

数値解は$h=0.1$くらいまでほぼ解析解の特徴を捉えており、細かい誤差を気にしなければEuler法でも十分使えそうだといえる。

FDTD法で場の運動のシミュレーション

最後に、電磁場など「場の運動」のシミュレーションを行うFDTD法についても軽く触れる。
電磁場の基礎方程式であるMaxwell方程式からスタートする。

$$\mathrm{div} D=\rho$$
$$\mathrm{rot} E=-\frac{\partial B}{\partial t}$$
$$\mathrm{div}B=0$$
$$\mathrm{rot}H=j+\frac{\partial D}{\partial t}$$

いま電流がない条件を考えると、rotに関する式は

$$\mathrm{rot} H=\epsilon\frac{\partial E}{\partial t}$$
$$\mathrm{rot} E=-\mu\frac{\partial H}{\partial t}$$

となる。左辺は空間に関する微分、右辺は時間微分で解の形状は$$E=\mathbf{E}(x,y,z,t)$$
という時間と3次元の座標に依存してベクトルを返す関数である。こういう表現を場(field)と呼称し、返却する値がスカラー、ベクトルである場合はスカラー場、ベクトル場と呼ぶ。

あとは時間発展をEuler法と同様に記述すればいい。
空間微分を計算するためにシミュレーションする空間を細かいセルに分ければ、微分は以下で近似
できる。

$$\phi^\prime(x+\frac{h}{2})=\frac{\phi(x+h)-\phi(x)}{h}$$

これを用いれば、

$$E(x+h_x,t+h_t)=E(x,t)+\frac{h_t}{\epsilon}rotH(h_x;x+\frac{h_x}{2},t+\frac{h_t}{2})$$
$$H(x+\frac{h_x}{2},t+\frac{h_t}{2})=H(x+\frac{h_x}{2},t-\frac{h_t}{2})-\frac{h_t}{2\mu}rotE(h_x;x,t)$$

となり、EもHも$h_x,h_t$のステップごとに飛び飛びの値がわかれば時間発展が計算できることがわかる。EとHで時間、空間とも半ステップずれた位置の値が必要なので、計算上は半分ずらしたセルを定義してやればいい。

実装してみよう。三次元ではなく、z方向には並進対称であるという条件でx,yの2dとして時間発展を求めよう。

import numpy as np
import matplotlib.pyplot as plt

H, E=np.zeros((2,96,96,3))

def mkgrad(M, isX, pLT):
	Mr=M*0
	if pLT:
		s=1
		e=None
	else:
		s=0
		e=-1
	if isX:
		Mr[:,s:e]=M[:,1:]-M[:,:-1]
		if pLT:
			Mr[:,0]=Mr[:,1]*2-Mr[:,2]
		else:
			Mr[:,-1]=Mr[:,-2]*2-Mr[:,-3]
	else:
		Mr[s:e]=M[1:]-M[:-1]
		if pLT:
			Mr[0]=Mr[1]*2-Mr[2]
		else:
			Mr[-1]=Mr[-2]*2-Mr[-3]
	return Mr


def mkrot(M, isLT):
	Mdx=mkgrad(M, True, not isLT)
	Mdy=mkgrad(M, False, not isLT)
	return np.stack((Mdy[...,2], -Mdx[...,2], Mdx[...,1]-Mdy[...,0])).transpose(1,2,0)

mu=1
eps=1

for t in range(8001):
	E[48:49,48:49,2]=np.sin(t/400)*7
	rotE=mkrot(E, True)
	rotH=mkrot(H, False)
	H+=eps*rotE/100
	E-=mu*rotH/100
	if t % 1000==0:
		plt.imshow(E[...,2], cmap="jet", vmin=-3, vmax=3)
		plt.pause(0.2)
		plt.clf()

境界で反射する前の結果。電場のz方向について描画する。

吸収境界条件が雑なので反射がキャンセルできてないが、ひとまずこのままにする(反射を抑えるにはMurの吸収境界条件があるらしい)。