Numerical test for the coupled nonlinear oscillator

Consider the coupled nonlinear oscillators defined by $$ \frac{d\psi_1}{dt} = i\omega_1 \psi_1 + (1- |\psi_1|^2)\psi_1 + D(\psi_2 - \psi_1) $$ $$ \frac{d\psi_2}{dt} = i\omega_2 \psi_2 + (1- |\psi_2|^2)\psi_2 + D(\psi_1 - \psi_2) $$

Our interest is on how these two oscillators' behavior depends on the difference of frequencies $\delta\omega = |\omega_1 - \omega_2 |$ and the interaction strength $D$. An analytical consideration leads to the following phase diagram which means the stationary state is divided into three regions.

Defining functions for numerical integration and plotting results

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def oscillator_deriv(x, t0, omega1=1., omega2= 1., D=1.5):
    """Compute the time-derivative of a coupled oscillator system. array of sets consisting of real part and imaginary part  """
    ampFac1 = 1.0 - (x[0]*x[0] + x[1]*x[1])
    ampFac2 = 1.0 - (x[2]*x[2] + x[3]*x[3])
    return [ampFac1*x[0] - omega1* x[1] + D*(x[2]- x[0]),
            ampFac1*x[1] + omega1* x[0] + D*(x[3]- x[1]),
            ampFac2*x[2] - omega2* x[3] + D*(x[0]- x[2]),
            ampFac2*x[3] + omega2* x[2] + D*(x[1]- x[3])]



def calc_and_plot(deltaOmega, D, max_time=20, run=True) :
    fig =plt.figure(figsize=(14,5))
    
    # show parameter space
    ax1 =fig.add_subplot(1,2,1)   # 1 in (1 row, 2 columns)
    t = np.linspace(0, 2.0, 100)
    x = [0.5*t[i] for i  in range(len(t))]
    y = x    
    t1 = np.linspace(2.0, 4.0, 100)
    x1 = [1.0 for i in range(len(t1))]
    x = np.concatenate((x, x1), axis=0)
    y = np.concatenate((y,[0.5 + 0.125*t1[i]*t1[i] for i in range(len(t1))]), axis=0)
    t = np.concatenate((t,t1),axis=0)
    
    ax1.plot(t,x)
    ax1.plot(t,y)
    ax1.fill_between(t,x,y,color="green")
    ax1.fill_between(t,[0 for i in range(len(t))], x,color="orange")
    ax1.set_xlabel(u"$\delta \omega$")
    ax1.set_ylabel(u"$D$")
    ax1.text(0.5,2.0, u"Syncronized Oscillation")
    ax1.text(3.0, 1.3, u"Arrest of Oscillation")
    ax1.text(2.0, 0.5, u"Drifting Oscillation")
    ax1.plot(deltaOmega, D, 'or') # plot the parameter set (deltaOmega, D)
    
    if not run:
        plt.show()
        return
    
    # solve diff. eq. and plot
    x0 = [0.2, 0.2, 0.2, 0.2 ]  # starting vector
    t = np.linspace(0, max_time, max_time*100) # time step is fixed
    x_t = integrate.odeint(oscillator_deriv, x0, t, args=(1.0, 1.0+deltaOmega, D))
    
    ax2= fig.add_subplot(1,2,2)
    ax2.plot(t,x_t[:,0], label = "$Re(\psi_1)$")
    ax2.plot(t,x_t[:,2], label = "$Re(\psi_2)$")
    ax2.set_xlabel("time")
    ax2.legend(loc=1)
    plt.title('$\omega_1 = 1.0, \delta\omega =' + str(deltaOmega) + '$, $D =' + str(D) + '$')
    
    plt.show()

Inline GUI for setting parameters and controling runs

from ipywidgets import interact, interactive, Checkbox, FloatSlider, Layout
from IPython.display import clear_output, display, HTML

interact(calc_and_plot,
         deltaOmega=FloatSlider(min=0.0, max=4.0, step=0.02, value=2.0, description="$\delta\omega$", layout=Layout(width='80%')),
         D=FloatSlider(min=0, max=2.5, step=0.01,value= 1.2, description="$D$", layout=Layout(width='80%')),
         max_time=FloatSlider(min=10,max=100,step=1, value=20, description="$t_{max}$", layout=Layout(width='80%')),
         run=Checkbox(value=False))