Three Wave Coupling System

M. Kono, 1998, 'Nonlinear Evolution of Instabilities in Plasmas' (in Japanese)

We consider a system described below, $$ \begin{gather} \partial_\tau a_1+V_1\partial_\xi a_1= \beta_1 a_2a_3 \exp(i\epsilon t)+\gamma_1a_1 \\ \partial_\tau a_2+V_2\partial_\xi a_2=-\beta_2 a_1a_3^*\exp(-i\epsilon t)+\gamma_2a_2 \\ \partial_\tau a_3+V_3\partial_\xi a_3=-\beta_3 a_1a_2^*\exp(-i\epsilon t)+\gamma_3a_3-i(\delta+\alpha|a_3|^2)a_3 \end{gather} $$

Here, $V_j$ is a wave propagation speed, $\beta_j$ a coupling parameter, $\epsilon$ frequency mismatch, $\gamma_j$ growth/gamping rate, $\delta$ resonance mismatch, and $\alpha$ frequency shift. Using this system, we consider three cases, 1) soliton
2) chaos
3) soliton+chaos.

Soliton Case

We set $V_1=0.3$, $V_2=-0.254$, $V_3=0.023$, $\beta_{1,2,3}=0.01$, $\epsilon=0$, $\gamma_{1,2,3}=0$, $\delta=0$, and $\alpha=0$.

In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sf
from IPython import display
import math as mt
import matplotlib.animation as animation
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['animation.embed_limit']=60
plt.rcParams['animation.html'] = 'jshtml'
In [2]:
def rk4fft(dx,nx,dt,nt,isav,v1,v2,v3,c1,c2,c3,gam1,gam2,gam3,delt,delt2,E1,E2,E3):
    global kk,kd
    
    jj=np.arange(nx)
    xx=dx*jj; xmax=dx*nx
    tt=dt*np.arange(nt)
    kk=2*mt.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

    #E1=np.zeros(nx,dtype=complex)
    #E2=np.zeros(nx,dtype=complex)
    #E3=np.zeros(nx,dtype=complex)
    #E1=E1+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))
    #E2=E2+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))
    #E3=E3+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))
    E1f=sf.fft(E1)
    E2f=sf.fft(E2)
    E3f=sf.fft(E3)
    
    #v1=0.3;v2=-0.254;v3=0.023
    #v1=0.;v2=-0.;v3=0.0
    #c1=0.01;c2=0.01;c3=0.01
    #gam1=0.001; gam2=-0.015; gam3=-0.015
    #gam1=0.5*gam1; gam2=0.5*gam2; gam3=0.5*gam3
    #delt=0.; delt2=0.002
    
    E1hst=np.zeros((nt//isav,nx),dtype=complex)
    E2hst=np.zeros((nt//isav,nx),dtype=complex)
    E3hst=np.zeros((nt//isav,nx),dtype=complex)    
    
    for it in range(nt):
        t=dt*it
        #---double steps with integrating factor method(4th-order Runge-Kutta)---#
        E1f=np.exp(-v1*1j*kk*dt+gam1*dt)*E1f
        E2f=np.exp(-v2*1j*kk*dt+gam2*dt)*E2f
        E3f=np.exp(-v3*1j*kk*dt+(gam3-1j*delt)*dt)*E3f
        
        gE11,gE21,gE31=f(E1f            ,E2f            ,E3f            ,t       )
        gE12,gE22,gE32=f(E1f+0.5*dt*gE11,E2f+0.5*dt*gE21,E3f+0.5*dt*gE31,t+0.5*dt)
        gE13,gE23,gE33=f(E1f+0.5*dt*gE12,E2f+0.5*dt*gE22,E3f+0.5*dt*gE32,t+0.5*dt)
        gE14,gE24,gE34=f(E1f+    dt*gE13,E2f+    dt*gE23,E3f+    dt*gE33,t+    dt)
        E1f=E1f+dt*(gE11+2*gE12+2*gE13+gE14)/6
        E2f=E2f+dt*(gE21+2*gE22+2*gE23+gE24)/6
        E3f=E3f+dt*(gE31+2*gE32+2*gE33+gE34)/6

        if(it%isav==0):
            E1hst[it//isav,:]=sf.ifft(E1f)
            E2hst[it//isav,:]=sf.ifft(E2f)
            E3hst[it//isav,:]=sf.ifft(E3f)
        
    return locals()
        
def f(E1f,E2f,E3f,t):
    E1=sf.ifft(E1f*kd)
    E2=sf.ifft(E2f*kd)
    E3=sf.ifft(E3f*kd)
    
    E2s=np.conj(E2)
    E3s=np.conj(E3)
    a1=sf.fft(-c1*E2*E3 *np.exp( 1j*delt2*t)) 
    a2=sf.fft(-c2*E1*E3s*np.exp(-1j*delt2*t))
    a3=sf.fft( c3*E1*E2s*np.exp(-1j*delt2*t))
    
    return a1,a2,a3
In [3]:
nt=50000; isav=1000
dt=0.1
dx=1
nx=512

v1=0.3;v2=-0.254;v3=0.023
c1=0.01;c2=0.01;c3=0.01
gam1=0.0; gam2=-0.0; gam3=-0.0
delt=0; delt2=0.00

xx=dx*np.arange(nx)
xmax=xx[-1]
E1=np.zeros(nx,dtype=complex)
E2=np.zeros(nx,dtype=complex)
E3=np.zeros(nx,dtype=complex)
E1=E1+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))
E2=E2+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))
E3=E3+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))

data=rk4fft(dx,nx,dt,nt,isav,v1,v2,v3,c1,c2,c3,gam1,gam2,gam3,delt,delt2,E1,E2,E3)
locals().update(data)
In [4]:
def update_anim(it):
    
    fig.clf()

    ax1 = fig.add_subplot(111)
    
    ax1.clear() 

    xx=dx*np.arange(nx)
    im1=ax1.plot(xx,np.real(E1hst[it,:]),label='$a_1$')
    im2=ax1.plot(xx,np.real(E2hst[it,:]),label='$a_2$')
    im3=ax1.plot(xx,np.real(E3hst[it,:]),label='$a_3$')
    ax1.legend()
    ax1.set_xlabel(r'$\xi$')
    ax1.set_ylabel('$a_j$')
    ax1.set_xlim(0,nx*dx)
    ax1.set_ylim(-4,4)
In [5]:
fig=plt.figure(figsize=(16,5))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)
plt.close()
anim
Out[5]:
In [6]:
xx=dx*np.arange(nx)
plt.figure(figsize=(16,4))
for it in range(nt//isav):
    plt.plot(np.real(E1hst[it,:])+it,xx,'k-')
    
plt.xlabel(r'$\tau$')
plt.ylabel(r'$\xi$')
Out[6]:
Text(0, 0.5, '$\\xi$')

Chaos case

We set $V_j=0$, $\beta_{1,2,3}=0.01$, $\epsilon=0.002$, $\gamma_{1,2}=-0.015$, $\gamma_3=0.001$, $\delta=0$, and $\alpha=0$.

In [9]:
nt=20000; isav=10
dt=1
dx=1
nx=512

v1=0;v2=0;v3=0
c1=0.01;c2=0.01;c3=0.01
gam1=-0.015; gam2=-0.015; gam3=0.001
delt=0.; delt2=0.002

E1=np.ones(nx,dtype=complex)
E2=np.ones(nx,dtype=complex)
E3=np.ones(nx,dtype=complex)

data=rk4fft(dx,nx,dt,nt,isav,v1,v2,v3,c1,c2,c3,gam1,gam2,gam3,delt,delt2,E1,E2,E3)
locals().update(data)
In [10]:
tt=np.linspace(0,nt*dt,nt//isav)
plt.figure(figsize=(16,4))
plt.plot(tt,np.abs(E3hst[:,1]))
plt.xlim(0,nt*dt)
plt.ylim(0,5)
plt.xlabel(r'$\tau$')
plt.ylabel(r'$|a_3|$')
plt.show()

Soliton+Chaos

In [11]:
nt=120000; isav=2400
dt=0.1
dx=1
nx=512

v1=0.3;v2=-0.254;v3=0.023
c1=0.01;c2=0.01;c3=0.01
gam1=-0.015; gam2=-0.015; gam3=0.001
delt=0.; delt2=0.002

E1=np.zeros(nx,dtype=complex)
E2=np.zeros(nx,dtype=complex)
E3=np.zeros(nx,dtype=complex)
E1=E1+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))
E2=E2+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))
E3=E3+1e0*np.exp(1j*(2*np.pi*xx/xmax+2*np.pi*np.random.rand()))

data=rk4fft(dx,nx,dt,nt,isav,v1,v2,v3,c1,c2,c3,gam1,gam2,gam3,delt,delt2,E1,E2,E3)
locals().update(data)
In [12]:
def update_anim(it):
    
    fig.clf()

    ax1 = fig.add_subplot(111)
    
    ax1.clear() 

    xx=dx*np.arange(nx)
    im1=ax1.plot(xx,np.real(E1hst[it,:]),label='$a_1$')
    im2=ax1.plot(xx,np.real(E2hst[it,:]),label='$a_2$')
    im3=ax1.plot(xx,np.real(E3hst[it,:]),label='$a_3$')

    ax1.legend()
    ax1.set_xlabel(r'$\xi$')
    ax1.set_ylabel('$a_j$')
    ax1.set_xlim(0,nx*dx)
    ax1.set_ylim(-5,5)
In [13]:
fig=plt.figure(figsize=(16,4))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav)
plt.close()
anim
Out[13]:
In [14]:
xx=dx*np.arange(nx)
plt.figure(figsize=(16,4))
for it in range(nt//isav):
    plt.plot(np.real(E1hst[it,:])+it,xx,'k-')
    
plt.xlabel(r'$\tau$')
plt.ylabel(r'$\xi$')
plt.show()
In [3]:
!jupyter nbconvert --to html --template tmp.tpl three_wave_resonance.ipynb

[NbConvertApp] Converting notebook three_wave_resonance.ipynb to html
[NbConvertApp] Writing 38279223 bytes to three_wave_resonance.html