an issue for time-depent lags dde of a common model in laser

[xiaoguo1995]

Dear All,

I'm involved to simulate a result from a dynamic model in laser to match with measured experimental results and I fail to obtain expected result in Matlab (using dde23, ddesd). Therefore, I transfer to use Julia and face a problem as well (cannot even let laser lasing...)

The model describes an optical feedback effect back in laser and the model (rate euqations) are described here, on page 581 for this paper:
https://www.osapublishing.org/aop/abstract.cfm?URI=aop-7-3-570
and this is a screenshot for the rate equation:

The three differential equations above the the photon density, carrier density and phase changing with time.

$$\tau_{ext}$$ is the time-dependent lag in a harmonic motion: $$2L/c + 2lambda/c * sin(2\pi f t)$$

below is the code I write in Julia and refer to official example on documentation:
`

using DifferentialEquations
using Plots
const c = 3e+8;
L_ext_0 = 0.8;
freq_oscillation = 1 / 1e-6; # 1000 ns
lambda = 850e-9; # 850 nm
A = lambda;
 
# time dependent lag
tau(u, p, t) = 2 * L_ext_0 / c + 2 * A / c * sin(2 * pi * freq_oscillation * t)
 
# const
const eta_i = 0.8; I = 1.5 * 30e-3; const q = 1.602e-19; const V = 0.9e-16;
const Gamma = 0.25; const V = 0.9e-16; const tau_n = 1.4e-9;
const n_g = 3.6; const L = 290e-6; const R1 = ((1 - n_g) / (1 + n_g))^2; const R2 = R1;
const alpha_i = 1500; const alpha_m = 1 / (2 * L) * log(1 / (R1 * R2));
const v_g = c / n_g;
const tau_p = 1 / (v_g * (alpha_i + alpha_m));

const beta_sp = 2e-4; const tau_sp = 2.3e-9;

Attenuation = -40; # dB
epsilon_loss = 10^(Attenuation / 10);
R = 0.7;
kappa = epsilon_loss * (1 - R2) * sqrt( R / R2 );
tau_in = (2 * L) / v_g;
kappa_tide = kappa / tau_in;

freq_1 = c / lambda;
omega_th_1 = 2 * pi * freq_1;
freq_2 = freq_1 + c / (2 * n_g * L);
omega_th_2 = 2 * pi * freq_2;

const alpha_Henry = 4.6;
const N_tr = 2e+24;
const a = 0.75e-20;

# define rre function
function rre_multimode(du, u , h, p, t)
    G1 = v_g * a * (u[1] - N_tr)
    du[1] = (eta_i * I) / (q * V) - u[1] / tau_n - (G1 * u[2])
    du[2] = Gamma * G1 * u[2] - u[2] / tau_p + Gamma * beta_sp * u[1] / tau_sp + 2 * kappa_tide * sqrt(u[2] * h(p, t-tau(u, p, t))[2]) * cos(omega_th_1 * tau(u, p, t) + u[3] - h(p, t-tau(u, p, t))[3])
    du[3] = 0.5 * alpha_Henry * Gamma * G1 - 0.5 * alpha_Henry / tau_p - kappa_tide * sqrt(h(p, t-tau(u, p, t))[2] / u[2]) * sin(omega_th_1 * tau(u, p, t) + u[3] - h(p, t-tau(u, p, t))[3])
    #du[2] = Gamma * G1 * u[2] - u[2] / tau_p + Gamma * beta_sp * u[1] / tau_sp + 2 * kappa_tide * sqrt( complex(u[2] * h(p, t-tau(u, p, t)) )[2]) * cos(omega_th_1 * tau(u, p, t) + u[3] - h(p, t-tau(u, p, t))[3])
    #du[3] = 0.5 * alpha_Henry * Gamma * G1 - 0.5 * alpha_Henry / tau_p - kappa_tide * sqrt( complex(h(p, t-tau(u, p, t))[2] / u[2]) ) * sin(omega_th_1 * tau(u, p, t) + u[3] - h(p, t-tau(u, p, t))[3])
end

h(p, t) = [1, 1, pi, 1, pi]
tspan = (0.0, 1000e-9)
u0 = Float64[1, 1, pi, 1, pi]
d_lags = ((u,p,t) -> tau(u,p,t),)
#d_lags = (tau,)

prob = DDEProblem(rre_multimode, u0, h, tspan, dependent_lags = d_lags)
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg, adaptive = false, dt=1e-11, dtmax = 1e-11, reltol = 1e-8, abstol= 1e-9, force_dtmin = true)

plot(sol.t, sol[2,:])

but I cannot get the result even as matlab did from this dynamic model (rate euqation).

Below is the matlab simulation result (simulated by ddesd solver):

Even the above beautiful waveform is not the same as what I simulate from a steady-state equation (derived from rate equation above), which looks like this:

In experiment, we can observe the tilt and sharp fringes like the picture above and this is the expectation.

I appreciate guys who replied to my previous problem and specially for Chris who replied me in email.

I have make my effort during this vacation and try to avoid any trivial mistake and till now I cannot figure out what is the problem behind. I cannot understand why matlab gives different simulation results as well.

Any help is appreciated and for this issue, I think Julia DDE solver for this laser model doesn't work for some reasons.

Thank you very much!
Xiao

[ChrisRackauckas]

Transferred from docs to DelayDiffEq.jl. @devmotion do you think it might be a weird instability issue related to #67 ?

[devmotion]

Hmm.... It seems the provided example errors since S(t)S(t-tau_ext) becomes negative and hence it's not possible to compute the square root of this term anymore. I guess that's also the reason for the lines you commented out (BTW they are not correct since you try to write complex values to a vector of real numbers)?

First of all, it would be interesting to see your MATLAB implementation, if that's possible. Did you consider possibly negative values of S(t)S(t-tau_ext) therein (e.g. by specifying that S should be non-negative)? Or did you just constrain the step size?

Actually, a bit more than one year ago I worked with a DDE that was also not defined for negative values. You can find a Julia implementation of it here, and you can also have a look at simulations with different solvers. As you can see in the Julia implementation, I set the derivatives of variables that should be non-negative to max(0, dxdt), following Shampine's advice. In my case, this implementation together with a suitable algorithm such as Rodas5 was sufficient to be able to compute simulations of my model reliably. However, you might want to have a look at isoutofdomain and the PositiveDomain callback, which are explained in the documentation. Another approach (which might be a bit more elegant) could be to enforce the positivity of S by considering the dynamics of a new variable T = f(S) instead of the dynamics of S, where f is a bijective function and T is unconstrained, contrary to S.

[xiaoguo1995]

Thanks for your reply! @devmotion
Actually, S(t) has a physical meaning - it's photon density. Therefore, it cannot be negative.
the lines you commented out is an naive attempt to let program runs and it fails.
I will try your advice and update how 's it going in Julia.

Thanks,
Peter