How do I setup the timestep when using DifferentialEquations.jl in Julia for an irregular time series?

Question

Playing with the harmonic oscillator, the differential equation is driven by a regular time series w_i in the millisecond range.

ζ = 1/4pi                                           # damped ratio

function oscillator!(du,u,p,t)

      
           du[1] = u[2]                             # y'(t) = z(t)
           du[2] = -2*ζ*p(t)*u[2] - p(t)^2*u[1]     # z'(t) = -2ζw(t)z(t) -w(t)^2y(t)
       end
y0 = 0.0                                            # initial position
z0 = 0.0002                                         # initial speed
u0 = [y0, z0]                                       # initial state vector
tspan = (0.0,10)                                    # time interval

dt = 0.001                                          # timestep

w = t -> freq[Int(floor(t/dt))+1]                   # time series

prob = ODEProblem(oscillator!,u0,tspan,w)           # define ODEProblem

sol = solve(prob,DP5(),adaptive=false,dt=0.001)

How do I setup the timestep when the parameter w_i is an irregular time series in the millisecond range.

date                    │  w  
────────────────────────┼───────
2022-09-26T00:00:00.023 │  4.3354
2022-09-26T00:00:00.125 │  2.34225
2022-09-26T00:00:00.383 │ -2.0312
2022-09-26T00:00:00.587 │ -0.280142 
2022-09-26T00:00:00.590 │  6.28319
2022-09-26T00:00:00.802 │  9.82271
2022-09-26T00:00:00.906 │ -5.21289
....................... |  ........

Answer 1

While it's possible to disable adaptivity, and even if it was possible to force arbitrary step sizes, this isn't in general what you want to do, as it limits the accuracy of the solution greatly.

Instead, interpolate the parameter to let it take any value of t. Fortunately, it's really simple to do!

using Interpolations

...

ts = [0, 0.1, 0.4, 1.0]
ws = [1.0, 2.0, 3.0, 4.0]
w = linear_interpolation(ts, ws)
tspan = first(ts), last(ts)

prob = ODEProblem(oscillator!, u0, tspan, w)
sol = solve(prob, DP5(), dt=0.001)

Of course, it doesn't need to be a linear interpolation.

If you still need the solution saved at particular time points, have a look at saveat for solve. E.g. saving the solution using ts used in the interpolation:

sol = solve(prob, DP5(), dt=0.001, saveat=ts)

Edit: Follow up on comment:

Mathematically, you are always making some assumption about the w(t) over the entire domain tspan. There is no such as "driven by a time series".

For example, the standard Runge-Kutta method you have chosen here will require that the ODE function is at h/2. For the better DP5() it is evaluated at several more sub-steps. This is of course unavoidable, regardless of adaptivity is used or not. Try adding println(t) into your ODE function and you will see this.

In case someone comes from matlab's ode45, not that it simply still uses adaptivity, and just treats explicit time steps the same as saveat does. And, of course, it will evaluate the function at various t outside of the explicit steps as well.

So even in your first example, you are interpolating your w. You are making a strange type of constant_interpolation (but with floor, which combined with floats will cause other issues, since floor(n*dt/dt) might evaluate to n or n-1.).

And even if you were to pick a method that only will try to evaluate at exactly the predetermined time steps, say e.g. ExplicitEuler(), you are still implicitly making the same assumption that w(t) is constant up until the next time step. Only now, you are also getting a much worse solution from just the ODE integration.

If a constant-previous type interpolation really is how w(t) is defined over the entire domain (which is what you did with floor(t/dt)) here, then what we have is:

w = extrapolate(interpolate((ts,), ws, Gridded(Constant{Previous}())), Flat())

There is simply mathematically no way we get to ignore what happens across the time-step, and there is no reason to limit the time-stepping to the sample points of our "load" function. It's not natural is correct in any mathematical sense. u'(t) has to be defined on the entire domain we integrate over.