Soliton KdV_equation¶
We reproduce results of the famous paper of KdV equation (Zabusky&Kruskal1965 PRL) .
The governing equation(the KdV equation) is $${u_t} + u{u_x} + {\delta ^2}{u_{xxx}} = 0.$$
In this simulation, the space is discretized by the spectral method and the time is discretized by a double step method with integrating factor method (+4th-order Jameson method)(Kitagawa2012 in Japanese).
Simulation parameters are $\delta=0.022$, $n_x=512$ where $n_x$ is the grid number of the simulation.
We run the simulation until $\pi t=30.4$.
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
import scipy.fftpack as sf
from tqdm import tqdm
import matplotlib.animation as animation
def KdV(nx,xmax,dx,nt,dt,nu,u,idata):
global kk,kd,ccnv,cdif
#---parameters---#
kk=2*np.pi/xmax*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
kd=np.zeros(nx)
kd[np.abs(kk)<=2/3*np.max(abs(kk))]=1.0 ### 2/3-rule for dealiasing
ccnv=1j*kk*0.5
cdif=nu*(1j*kk)**3
udata=np.zeros((int(nt/idata)+1,nx))
udata[0,:]=u
uf=sf.fft(u)
for it in tqdm(range(nt)):
#---4th-order Runge-Kutta---#
#g1=f(uf)
#g2=f(uf+0.5*dt*g1)
#g3=f(uf+0.5*dt*g2)
#g4=f(uf+ dt*g3)
#uf=uf+dt*(g1+2*g2+2*g3+g4)/6
#---double steps with integrating factor method(4th-order Runge-Kutta)---#
uf=np.exp(-cdif*dt)*uf
g1=f(uf)
g2=f(uf+0.5*dt*g1)
g3=f(uf+0.5*dt*g2)
g4=f(uf+ dt*g3)
uf=uf+dt*(g1+2*g2+2*g3+g4)/6
if it%idata==0:
#print(it)
#display.clear_output(wait=True)
#plt.plot(xx,u,'-')
#plt.show()
udata[it//idata,:]=np.real(sf.ifft(uf))
return locals()
def f(uf):
u=sf.ifft(uf*kd)
u_nl=ccnv*sf.fft(u*u)
return -u_nl#-cdif*uf
1-soliton solution¶
The KdV equation has a solution for one soliton, $$u(x,t) = {u_\infty } + A{{\mathop{\rm sech}\nolimits} ^2}\left( {\frac{{x - {x_0} - Ct}}{\Delta }} \right),$$ where $C=u_\infty+A/3$, $\Delta^2=12\delta^2/A$.
We assume here that $u_\infty=0$, $A=1$, and $x_0=1$.
nx=512;xmax=2;dx=xmax/nx;dt=0.01*dx;nu=0.022**2;idata=500
tb=1/np.pi
tr=30.4*tb
nt=int(tr/dt)
#---initial condition---#
xx=dx*np.arange(nx)
A=1.0
xi=1.0
dd=np.sqrt(12*nu/A)
C=A/3
u=A/np.cosh((xx-xi)/dd)**2
data=KdV(nx,xmax,dx,nt,dt,nu,u,idata)
locals().update(data)
xx=np.linspace(0,2,nx)
plt.imshow(udata,extent=[0,xmax,0,1],origin='lower',aspect=3,cmap='terrain')
plt.xlabel('x');plt.ylabel(r'$t/t_r$')
plt.xlim(0,2)
plt.ylim(0,1)
plt.plot(xx,(xx-0.8)/C/30.4*np.pi,'w--')
plt.show()
2-solition solution¶
The KdV equation can also descrive a interaction between two solitons.
The initial condition is
$$u(x,0) = A_1{{\mathop{\rm sech}\nolimits} ^2}\left( {\frac{{x - {x_1}}}{\Delta_1 }} \right)+A_2{{\mathop{\rm sech}\nolimits} ^2}\left( {\frac{{x - {x_2}}}{\Delta_2 }} \right),$$
where $\Delta_1^2=12\delta^2/A_1$ and $\Delta_2^2=12\delta^2/A_2$.
We use $A_1=0.9$, $A_2=0.4$, $x_1=0.5$, and $x_2=1.25$.
nx=512;xmax=2;dx=xmax/nx;dt=0.01*dx;nu=0.022**2;idata=500
tb=1/np.pi
tr=30.4*tb
nt=int(tr/dt)
#---initial condition---#
xx=dx*np.arange(nx)
A1=0.9
xi1=0.5
C1=A1/3
dd1=np.sqrt(12*nu/A1)
A2=0.4
xi2=1.25
C2=A2/3
dd2=np.sqrt(12*nu/A2)
u=A1/np.cosh((xx-xi1)/dd1)**2+A2/np.cosh((xx-xi2)/dd2)**2
data=KdV(nx,xmax,dx,nt,dt,nu,u,idata)
locals().update(data)
tt=np.linspace(0,1,nt)
xx=np.linspace(0,2,nx)
plt.imshow(udata,extent=[0,xmax,0,1],origin='lower',aspect=3,cmap='terrain')
plt.xlabel('x');plt.ylabel(r'$t/t_r$')
plt.xlim(0,2)
plt.ylim(0,1)
plt.plot(xx,(xx-0.65)/C1/30.4*np.pi,'w--')
plt.plot(xx,(xx-1.15)/C2/30.4*np.pi,'w--')
plt.show()
Zabusky&Kruskal's experiment¶
The initial condition is $$u(x,0)=\cos(\pi x)$$.
nx=512;xmax=2;dx=xmax/nx;dt=0.01*dx;nu=0.022**2;idata=500
tb=1/np.pi
tr=30.4*tb
nt=int(tr/dt)
u=np.cos(np.pi*np.linspace(0,xmax,nx,endpoint=False))
data=KdV(nx,xmax,dx,nt,dt,nu,u,idata)
locals().update(data)
xx=np.linspace(0,2,nx)
it0=0;itb1=int(tb/dt/idata);itb36=int(3.6*tb/dt/idata)
plt.plot(xx,udata[it0, :], 'k:', label='t=0')
plt.plot(xx,udata[itb1,:], 'k--',label=r'$t=1/\pi$')
plt.plot(xx,udata[itb36,:],'k-', label=r'$t=3.6/\pi$')
plt.xlim(0,2);plt.ylim(-1,3)
plt.xlabel('x')
plt.ylabel('u')
plt.legend()
plt.show()
tt=np.linspace(0,1,nt)
xx=np.linspace(0,2,nx)
plt.imshow(udata,extent=[0,xmax,0,1],origin='lower',aspect=3,cmap='terrain')
plt.xlabel('x');plt.ylabel(r'$t/t_r$')
plt.show()
#---make an animation---#
ims=[]
ntani=int(nt/idata)//2
fig=plt.figure(figsize=(15,5))
for it in range(ntani):
#if(it%100==0):
im=plt.plot(xx,udata[it*2,:],'k-')
plt.xlabel(r'$x$')
plt.ylabel(r'$u$')
ims.append(im)
plt.rcParams['animation.html'] = 'html5'
ani=animation.ArtistAnimation(fig,ims)
plt.close()
ani
!jupyter nbconvert --to html --template tmp.tpl KdV.ipynb