The Derivative NonLinear Schrödinger (DNLS) Equation

Reference: Nariyuki2006NGP

\begin{gather} \frac{\partial b}{\partial t} + \frac{\partial}{\partial x}(|b|^2b)+i\frac{\partial^2b}{\partial x^2}=0; \hfill \\ b=b_y+ib_z \end{gather}

%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
import scipy.fftpack as sf
import matplotlib.animation as animation
from tqdm import tqdm
import matplotlib as mpl
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5

def DNLS(nx,xmax,nt,dt,b,isav):
    global kk,kkd

    bf=sf.fft(b)

    kk =2*np.pi/xmax*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    kkd=np.zeros(nx)
    kkd[np.abs(kk)<=2/(3+1)*np.max(abs(kk))]=1.0 ### 2/(n+1)-rule for dealiasing see for Caballero2024arXiv: "cuPSS: a package for pseudo-spectral integration of stochastic PDEs"

    bhst=np.zeros((int(nt/isav),nx), dtype=complex)
    bhst[0,:]=sf.ifft(bf)
    for it in tqdm(range(nt)):

        #---4th-order Runge-Kutta---#
        #g1=f(bf          )
        #g2=f(bf+0.5*dt*g1)
        #g3=f(bf+0.5*dt*g2)
        #g4=f(bf+    dt*g3)
        #bf=bf+dt*(g1+2*g2+2*g3+g4)/6
        
        #---double steps with integrating factor method(4th-order Runge-Kutta)---#
        bf=np.exp(1j*kk**2*dt)*bf
        g1=f(bf          )
        g2=f(bf+0.5*dt*g1)
        g3=f(bf+0.5*dt*g2)
        g4=f(bf+    dt*g3)
        bf=bf+dt*(g1+2*g2+2*g3+g4)/6
        
        if it%isav==0: bhst[it//isav,:]=sf.ifft(bf)    
    return locals()

def f(bf):
    b=sf.ifft(bf*kkd)
    b_nl=-1j*kk*sf.fft(abs(b)**2*b)
    return b_nl#+1j*kk**2*bf

nx=2048;xmax=256;dx=xmax/nx;dt=0.1*dx;isav=200
nt=80000
#---initial condition---#
xx=dx*np.arange(nx)
b0=0.4;km=11
b=b0*np.cos(km*2*np.pi*xx/xmax)+1j*b0*np.sin(km*2*np.pi*xx/xmax)
for ik in range(-256,256):
    bw  =1e-5*np.random.randn()
    phiw=0#2*np.pi*np.random.rand()
    b=b+bw*np.cos(ik*2*np.pi*xx/xmax+phiw)+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
data=DNLS(nx,xmax,nt,dt,b,isav)
locals().update(data)

100%|█████████████████████████████| 80000/80000 [00:22<00:00, 3585.26it/s]

Spatiotemporal Evolution

## xx=np.linspace(0,2,nx)
plt.imshow(abs(bhst).T,origin='lower',aspect='auto',cmap='terrain',extent=([0,nt*dt,0,nx*dx]))
plt.title('|b|')
plt.xlabel('t');plt.ylabel(r'$x$')
plt.colorbar()
plt.clim(0,1.4)
plt.show()

Time evolution of the power spectrum

kmin=2*np.pi/(dx*nx)*(-nx/2); kmax=2*np.pi/(dx*nx)*(nx/2)
plt.imshow(np.abs(sf.fftshift(sf.fft(((bhst)).T,axis=0),axes=0)),origin='lower',aspect='auto',cmap='jet', norm=mpl.colors.LogNorm(),extent=([0,nt*dt,kmin,kmax]))
plt.colorbar()
plt.title(r'$\log |b|$')
plt.xlabel('t');plt.ylabel(r'$k_m$')
plt.ylim(kmin/8,kmax/4)
plt.clim(1e-3,1000)
plt.show()

Animation

#---make an animation---#
ims=[]
ntani=int(nt/isav)//5
fig=plt.figure(figsize=(15,5))
flag_legend=True

for it in range(ntani):
    im1=plt.plot(xx,np.real(bhst[it*5,:]),'C0' ,label=r'$b_y$')
    im2=plt.plot(xx,np.imag(bhst[it*5,:]),'C1' ,label=r'$b_z$')
    im3=plt.plot(xx,    abs(bhst[it*5,:]),'k--',label=r'$|b|$')
    im4=plt.plot(xx,   -abs(bhst[it*5,:]),'k--')
    plt.xlabel(r'$x$')
    ims.append(im1+im2+im3+im4)
    if flag_legend:
        plt.legend()
        flag_legend = False
        
plt.rcParams['animation.html'] = 'html5'
ani=animation.ArtistAnimation(fig,ims)
plt.close()

ani

Energy conservation

plt.plot(np.mean(abs(bhst.T)**2,axis=0))
plt.xlabel(r'$t$');plt.ylabel('$|b|^2$')
plt.show()