Why does scipy.odeint with differing step size in the t array gives qualitatively different results?

Question

I am investigating strange attractors and am currently working on the Duffing attractor (a 2D strange attractor). The solution trajectory makes loops in the positive and negative x direction and I am particularly interested in the sequence in which it switches back and forth between the two (the pattern is chaotic). This is fairly easy to assess since you can plot the x-value of the trajectory vs. time (the independent variable). I use the following code:

a=0.35
b=0.3
w=1

def dX_dt(X, t):
    return np.array([X[1],
                     X[0]-(X[0]**3)-a*X[1]+b*np.cos(w*t)])

tstepmax=5000
tmax=1000
tmin=0
t1 = np.linspace(tmin, tmax+tmin,  tstepmax)
X1=odeint(dX_dt,[0,0],t1)

I was not sure how many steps were required for a given number of time-units, so I experimented with setting tstepmax to 50000 and 500000. To my surprise, after a fairly short time these gave qualitatively different results as seen in this plot image: Three overlaid trajectories I was under the impression that the t1 vector was simply where the trajectory was evaluated for plotting purposes and that the method internally determined the step size to use; indeed, for some methods the step size isn't constant over the trajectory. However, this change gives a fundamentally different trajectory. Why? And, perhaps more importantly, which is the "real" trajectory? How can I reliably obtain the "real" trajectory?

Answer 1

A similar problem was observed and discussed in https://scicomp.stackexchange.com/questions/42613/need-help-to-fully-understand-scipys-odeints-reported-step-sizes-eval-times.

In short, the cause seems to be that the odeint stepper call (in the lsoda Fortran code) does not go a full step but only to the next output value. The stepper state is saved in a workspace array, but apparently not completely. Some control variables seem to be reset at each call, each call virtually constructs a new stepper object mixing old and default information. This leads to slightly different internal step sequences, thus to slightly different error evolutions, which in a chaotic system can lead to visibly diverging solutions.

Using solve_ivp you should not observe this problem, as the time loop there does full internal steps and then interpolates all relevant output values for the current interval. The stepper class continues to exist with all internal information. Thus the internal step sequence should be independent of the t_eval argument.