How to simulate the famous Lang-Kobayashi semiconductor laser model with delayed optical feedback

[YouknowIdoknow]

Hi, sir. When I use the Jitcdde to solve the the famous Lang-Kobayashi semiconductor laser model with delayed optical feedback.

It doesn't work. Can you help me solve it? The script is following:

from jitcdde import t, y, jitcdde
import numpy as np
import pylab as plt
from symengine import cos, sin, pi, Symbol
#--------------------------the system parameters----------------------------
Gn = 8.4e-13
N0 = 1.4e24
taup = 1.927e-12
k = 6e9
τ = 1.5e-9
taus = 2.04e-9
w0 = 0
a = 3.0
J = 1.098e33

f = [
    0.5 * ( Gn * ( y(2) - N0 ) - 1/taup) * y(0) + k * y(0, t - τ) * cos(w0*τ + y(1) - y(1, t - τ)),
    a * 0.5 * ( Gn * ( y(2) - N0) - 1/taup) - k * y(0, t - τ) / y(0) * sin(w0*τ + y(1) - y(1, t - τ)),
    J - y(2)/taus - Gn * ( y(2) - N0)*(y(0)**2),
]

DDE = jitcdde(f)
DDE.constant_past( [1.3e10, 0.0, 1.9e24], time=0.0 )
DDE.step_on_discontinuities()

dt = 5e-12
times = np.arange(dt, 200e-9, dt)
Y = np.vstack([ DDE.integrate(T) for T in times ])

plt.figure()
plt.plot(Y[10000:50000,0],Y[10000:50000,1], 'r-', lw=0.5)
plt.xlabel('x', fontsize=15)
plt.ylabel('y', fontsize=15)
plt.tick_params(labelsize=15)
# plt.xlim(0, 7)
# plt.ylim(-1.5, 1.5)
plt.title('Phase portrait')
plt.show()

[Wrzlprmft]

The problem with your implementation is that your time and dynamical variables have extreme orders of magnitude and – on top – these orders also differ between the dynamical variables. Adding/changing the following lines will make the integration run:

DDE.set_integration_parameters(first_step=1e-14,min_step=1e-15,max_step=1e-11,atol=1e-5)
DDE.step_on_discontinuities(min_distance=1e-14)

However, I cannot tell whether the result is any good and you certainly do not have proper error handling (as different error tolerances for different dynamical variables are not implemented). I strongly recommend that you rescale parameters such that the dynamical variables and timescale of the dynamics are roughly in the order of magnitude of 1.

Also see the third point in this FAQ and this Q&A on Stack Overflow.