I am solving the dynamics of a system when it interacts with a pulse, which basically is solving a time-dependent differential equation. In general it works fine, but whenever I take the bandwidth of the pulse small, i.e. around unity, the solver depends on where the pulse starts t0. Let me give you the code and some pictures
%matplotlib inline import matplotlib.pyplot as plt import numpy as np import matplotlib as mpl from scipy.integrate import solve_ivp import scipy as scipy ## Pulses def Box(t): return (abs(t-t0) < tau) #return np.exp(-(t - t0)**2.0/(tau ** 2.0)) ## Differential eq. def Leq(t,u,pulse): v=u[:16].reshape(4,4) a0=u da0=-a0+pulse(t)*E_0 M=np.array([[-1,0,a0,0], [0,-1,0,-a0], [a0,0,-1,0], [0,-a0,0,-1]])*kappa Dr=np.array([[1,1j,0,0], [1j,1,0,0], [0,0,1,1j], [0,0,1j,1]])*kappa/2.0 dv=M.dot(v)+v.dot(M)+Dr return np.concatenate([dv.flatten(), [da0]]) ## Covariance matrix cov0=np.array([[1,1j,0,0], [1j,1,0,0], [0,0,1,1j], [0,0,1j,1]])/4 ##initial vector cov0 = cov0.reshape(-1); ## vectorize initial vector a0_in=np.zeros((1,1),dtype=np.complex64) ##initial vector a0_in = a0_in.reshape(-1); ## vectorize initial vector u_0=np.concatenate([cov0, a0_in]) ### Parameters kappa=1.0 ##adimenstional kappa: kappa0/kappa tau=1.0 ##bandwidth pump #Eth=kappa0*kappa/g ##threshold intensity E_0=0.8 # pump intensity normalized to threshold t0=2.0 ##off set Tmax=10 ##max value for time Nmax=100000 ##number of steps dt=Tmax/Nmax ##increment of time t=np.linspace(0.0,Tmax,Nmax+1) Box_sol=solve_ivp(Leq, [min(t),max(t)] , u_0, t_eval= t, args=(Box,)) print(Box_sol.y[0:16,int(Nmax/2)].reshape(4,4))
and then just some expectation values and the plotting.
As you can see from the pictures, the parameters are the same, the only difference is the time shift t0. It seems as whenever the pulse’s bandwidth is small, I have to start really close to t=0, which I don’t fully understand why. It should,’t be like that.
It is a problem of the solver? If so, what can I do to fix it? I’m now worried that this problem can manifest in other simulations I am running and I don’t realise.
This is a well-known behavior, there have been several questions on this topic.
In short, it is the step size controller. The assumption behind it is that the ODE function is smooth to a high order, and that local behavior informs the global behavior in a medium range. Thus if you start flat, with vanishing higher derivatives, the step size is quickly increased. If the situation is unfortunate, this will jump over the non-smooth bump in all stages of the integrator step.
There are two strategies to avoid that
2*tau so that the bump is surely encountered, the step size selection will make sure that the sampling is dense around the jumps.
use separate integration for the pieces of the jump, so that the control input inside each segment is a constant.
The first variant is in a small way an abuse of the step size controller, as it works outside specifications in some heuristic emergency mode around the jumps. The second variant requires a little more coding effort, but will have less internal steps, as the ODE function for each segment is again completely smooth.