Python＋Scipyでルンゲクッタ8次のDOP853(Dormand＆Prince)を使う（その13) KdV方程式の時間発展にDOP853、空間差分はオリジナルのザブスキーとクルスカルのを使ってGIFアニメに。

KdV方程式

∂u/∂t　+u*∂u/∂x+δ^2 ∂^3u/∂x^3 = 0

δ=0.022

ですが、これも空間を差分化すると普通の常微分方程式系と見なせる。そこでルンゲクッタ8次のDOP853を使おう。

空間差分はオリジナルのザブスキーとクルスカルのものと同じとします。

クリックでGIFアニメが始まります。

ソースコードはこんな感じ。

%matplotlib notebook
import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from matplotlib.animation import PillowWriter
import numba

INTERVAL = 0.05 # [s] interval time
FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
FPS = 1000/FRAME_INTERVAL # frames per second

Nt=512
Nx=256

xmax=2.0
xmin=0.0
tmin=0.0
tmax=16/np.pi

t0=0
dt = (tmax-tmin)/Nt
dx = (xmax - xmin) / Nx
d=0.022

t=np.linspace(0,tmax,Nt)
x=np.linspace(0,xmax,Nx)
u0=np.cos(np.pi*x)
sol= np.empty((Nt, Nx))
sol[0] = u0

def kdv(t, u): #odeintのときとt,xの並びが逆
f=np.zeros(Nx)

    for i in range(Nx):
        if i < Nx-1:
            p1 = i + 1
        else:
            p1 = 0

        if i > 0:
            m1 = i - 1
        else:
            m1 = Nx-1

        if i < Nx - 2:
            p2 = i + 2
        elif i == Nx - 2:
            p2 = 0
        else:
            p2 = 1

        if i > 1:
           m2 = i - 2
        elif i == 1:
            m2 = Nx - 1
        else:
            m2 = Nx - 2

        diff3 = (-u[m2] + 2 * u[m1] - 2 * u[p1] + u[p2]) / (2 * (dx ** 3))
        diff1 = (-u[m1] + u[p1]) / (2 * dx)
        f[i] = -(u[m1]+u[i]+u[p1]) * diff1/3 - d*d*diff3

    return f

numba_f = numba.jit(kdv,nopython=True)

solver=ode(kdv)
solver.set_integrator('dop853')
solver.set_initial_value(u0,t0) #なぜか関数と並びが逆

k=1
while solver.successful() and solver.t < tmax:
solver.integrate(t[k])
sol[k] = solver.y
k+= 1

def update(i):
    plt.cla() #現在描写されているグラフを消去
    plt.xlim(0,2)
    plt.ylim(-1,3)
    plt.plot(x,sol[i],c='Red',lw=3)
    plt.xlabel("t")
    plt.ylabel("u")
    plt.title("KdV equation by DOP853")
    plt.grid()

# Plot
fig = plt.figure(figsize=(10,6))
ani = FuncAnimation(fig, update, frames=range(0,Nt,1), interval=FRAME_INTERVAL, blit=True)
plt.show()
ani.save("kdv.gif", writer = 'pillow',fps=FPS)