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$.

In [1]:
%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
In [2]:
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.r_[np.ones(nx//3),np.zeros(nx//3+nx%3),np.ones(nx//3)]   #de-aliased
    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$.

In [3]:
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)

100%|██████████| 247721/247721 [00:48<00:00, 5127.34it/s]
In [4]:
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$.

In [5]:
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)

100%|██████████| 247721/247721 [00:46<00:00, 5323.08it/s]
In [6]:
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)$$.

In [7]:
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)

100%|██████████| 247721/247721 [00:56<00:00, 4388.42it/s]
In [8]:
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()
In [9]:
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()
In [12]:
#---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()
In [13]:
ani
Out[13]: