« かつやで全力大人飯(チキンカツ丼)を食す。チキンカツ、焼きそば、から揚げが入った大人様ランチ。 | トップページ | すき家でチーズカルビ牛丼(大盛)を食す。ものすごく甘いタレだが、タバスコと一緒にするとなかなかおいしい。 »

2020年7月14日 (火)

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アニメが始まります。

Kdv

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

%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)

« かつやで全力大人飯(チキンカツ丼)を食す。チキンカツ、焼きそば、から揚げが入った大人様ランチ。 | トップページ | すき家でチーズカルビ牛丼(大盛)を食す。ものすごく甘いタレだが、タバスコと一緒にするとなかなかおいしい。 »

学問・資格」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« かつやで全力大人飯(チキンカツ丼)を食す。チキンカツ、焼きそば、から揚げが入った大人様ランチ。 | トップページ | すき家でチーズカルビ牛丼(大盛)を食す。ものすごく甘いタレだが、タバスコと一緒にするとなかなかおいしい。 »

最近の記事

2021年5月
            1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31          
フォト
無料ブログはココログ