Use a time-series as parameter of ODEProblem

New to Julia
1

Hello,

I am quite new to Julia and I would like to use it to solve some ODE problems.
Let’s take the following example:

function test!(du,u,p,t)
    yield = 10
    decay = 5
    du[1] = yield*p - decay*u[1]
end

I can then define some variables and run the solver:

u0 = [0]
tspan = (0,10)
p = 3

prob = ODEProblem(test!,u0,tspan,p)
sol = solve(prob);

That’s fine.
But I then would like to use a time-series for p (p is a time-varying variable which is recorded in some file). For example let’s take:

ts_time = (0:11)
ts_val = [1, 1, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0]

I saw in the examples of DifferentialEquations.jl that one casn define a time dependent function to send to the ODEProblem.
So I tried by linear interpolation:

ts = LinearInterpolation(ts_time, ts_val)

Now I can use ts as a time dependent function. I tried to sent it to the solver:

prob = ODEProblem(test!,u0,tspan,ts)
sol = solve(prob);

But it results to some error:

BoundsError: attempt to access 12-element extrapolate(scale(interpolate(::Vector{Float64}, BSpline(Linear())), (0:11,)), Throw()) with element type Float64 at index [12]

I do not understand what’s the problem here.

I also tried a formulation more similar to the example:

ts2 = t->ts(t)
prob = ODEProblem(test!,u0,tspan,ts2)
sol = solve(prob);

In this case I get the following error:

MethodError: no method matching *(::Int64, ::var"#25#26")
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:560
  *(::T, ::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} at int.jl:88
  *(::Union{Int16, Int32, Int64, Int8}, ::BigInt) at gmp.jl:541

I can I do what I am trying to do (potentially with a different approach…)?
Thanks in advance

1 Like
2

Since it’s an interpolation, you have to use the interpolation.

function test!(du,u,p,t)
    yield = 10
    decay = 5
    du[1] = yield*p(t) - decay*u[1]
end

p is just whatever you passed it. So when you did:

prob = ODEProblem(test!,u0,tspan,ts)
sol = solve(prob);

p is the LinearInterpolation object. So then yield*p tried to do multiplication between a scalar and a LinearInterpolation, which is what then causes the error message you see. But yield*p(t), well, that first interpolates p to get a scalar, then use that scalar, and it etc. does what you want.

From that description you can see how the other error you get is similar, with a similar solution to call the function.

1 Like
3

Thanks for your help…now it looks so obvious :innocent:

I am wondering whether the interpolation technique is a little bit too excessive. Maybe the best (fastest) approach would be to impose the time-step so that t always corresponds to one of the points of my table. So that no interpolation is necessary.
In my case, the phenomena is “slow enough” compared to the input time-series frequency so that no further time-step refinement would be necessary.

But can I impose to the solver the list of time points? Than I could maybe define p as a Dict so that the right p-value would corresponds to each time.
Again, I do not know if this is the best approach…

I tried to develop a similar approach via a for loop, where at each iteration I define and solve an ODEProblem dealing with a single time-step. But I suppose that doing that with a for loop is not the fastest way…

4

Well, for loops in Julia are fast. Nothing wrong with using a for loop…

5

Depends on how nonlinear your equation is.

You can pass tstops as an array of time points for the solver to hit, and then have a DiscreteCallback that updates a parameter value at those points and just use the constant value between them. You could also set adaptive=false and then those would be exactly the time steps, though whether that is okay depends on the required accuracy of the solution.

6

I think I am doing something weird here: I tried to implement the same for loop in python (with odeint) and Julia and I get more or less the same execution time…(0.5 seconds for about 2300 time steps) Julia is even a little bit slower and I am surprised…

However, with the linear interpolation approach the execution time is now reduced to about 0.03 seconds, so a reduction greater than 10x!

But to make things a little bit more confused, I tried to implement the linear interpolation approach also in python with scipy.interpolate.interp1d: in this case, I get by far the worst execution times (even compared to the for loops).

I should probably try to identify the bottleneck of each configuration…

Finally, not knowing how DiscreteCallback works, I tried to set adaptive=false…but I got some kind of errors saying that the time-step does not correspond to any index of my p (which is now a Dict)…I tried to round but with no success…I do not know if this pathway is worth the effort…

In any case thanks to everybody for your support and for making me process in Julia. :slight_smile:

7

Was it in global scope?

8

Not sure to understand what should be in global scope and what that would change. So, let me write down what I tested :slight_smile::

# Define the variable "p" and its interpolation
ts_time = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 100.] #Array(0.:11.)
ts_val = [1., 1., 3., 3., 3., 0., 0., 0., 0., 0., 0., 0.]
ts = LinearInterpolation(ts_time, ts_val);

# The function with "p(t)"
function test_t!(du,u,p,t)
    yield = 10
    decay = 5
    du[1] = yield*p(t) - decay*u[1]
end

#And then a test
for i in 1:5
    start = time()
    prob = ODEProblem(test_t!,u0,tspan,ts)
    global sol_t = solve(prob,saveat=1.);
    println(time()-start)
end

Which results in:

2.2445149421691895
0.0003249645233154297
0.00017905235290527344
0.00015592575073242188
0.0001480579376220703

Let’s test the for-loop approach:

# Function to run the transient as a `for-loop`
function test_transient(time_steps,p,u0)
    
    n_of_steps = length(time_steps)
    u = u0 * ones(n_of_steps)

    for i in range(2,stop=n_of_steps)
        tspan = (time_steps[i-1],time_steps[i])
        pi = p[i]
        u0 = [u[i-1]]
        
        prob = ODEProblem(test!,u0,tspan,pi)
        local sol = solve(prob,save_everystep=false)      
        u[i] = sol.u[end][1]       
    end

    return u
end

# Testing it
for i in 1:5
    start = time()
    global sol = test_transient(ts_time[1:end-1],ts,0)
    println(time()-start)
end

Which leads to:

4.448390007019043
0.0010211467742919922
0.0009670257568359375
0.0009632110595703125
0.0009682178497314453

So 5-10x slower
Surprisingly (at least for me), the results are not really identical. There is some “shift” in the ‘p’ parameter that I have to correct…but I do not think that this affect the speed test…
image

9

Reconstructing a new ODE problem and new solve in a loop is going to be really slow because it’s building and releasing the cache every step. If you want to do that, then using the integrator interface is better:

https://diffeq.sciml.ai/stable/basics/integrator/

And yes, make sure the definition is matching left vs right value choices.

10

Thanks for the suggestion: the integrator interface looks great for what I want to do.

However, I am getting errors since the very start of my code using it…probably due to my inexperience with Julia. I use the init function as solve, as I understood from instructions. However, I am getting a “typical Julia” error message:

MethodError: no method matching __init(::ODEProblem{Float64, Tuple{Float64, Float64}, true, Float64, ODEFunction{true, typeof(test!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem})

Here follows my test:

function test_integrate(time_steps,p,u0)

    n_of_steps = length(time_steps)
    
    u = u0 * ones(n_of_steps)

    println("define prob")
    prob = ODEProblem(test!,u0,(time_steps[1],time_steps[end]),p[1])
    println("init integrator")
    integrator = init(prob)
 
    for i in range(2,stop=n_of_steps)
        println("step",str(i))
        integrator.p = power[i]
        step!(integrator,time_steps[i-1]-time_steps[i],stop_at_tdt=true)        
    end
    
    return integrator
end

sol = test_integrate(ts_time[1:end-1],ts,0)

Thanks in advance for helping me with my newbie-issues…

11

Similar to solve, you should pass an algorithm if you’re directly using a solver package instead of DifferentialEquations.jl.

Open an issue and I can have it catch this one. The typical one is solve which looks like this:

julia> using OrdinaryDiffEq

julia> prob = ODEProblem((u,p,t)->1.0,1.0,(0.0,1.0))
ODEProblem with uType Float64 and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 1.0

julia> solve(prob)
ERROR: Default algorithm choices require DifferentialEquations.jl.
Please specify an algorithm (e.g., `solve(prob, Tsit5())` for an ODE)
or import DifferentialEquations directly.

You can find the list of available solvers at https://diffeq.sciml.ai/stable/solvers/ode_solve/
and its associated pages.

Stacktrace:
  [1] __solve(prob::ODEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Nothing; default_set::Bool, second_time::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
12

Nevermind, no need for the issue, the better error message was quick to add for the init case: Handle no defaults on init case by ChrisRackauckas · Pull Request #780 · SciML/DiffEqBase.jl · GitHub thanks for the report.

1 Like
13

Happy of having being somehow useful. :slight_smile:

So, I specified an algorithm (‘Tsit5()’) and I could go on.

(FYI, I also tried to directly import DifferentialEquations, then to add DifferentialEquations. prefix to any command (solve, ìnit…) otherwise they were no more ‘defined’ : but I still got the error with DifferentialEquations.init).

Now the problem is that the step function does not accept any argument exception made for the integrator.

step!(integrator,dt=time_steps[i]-time_steps[i-1],stop_at_tdt=true)

…results in:

MethodError: no method matching step!(::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, Vector{Float64}, Nothing, Float64, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Float64, ODEFunction{true, typeof(xenon!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(xenon!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(xenon!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}; dt=1.0, stop_at_tdt=true)

So, it is not accepting my “dt”, which is necessary in order to change the p value at the proper time…

14

You wrote those as keyword arguments, not positional arguments. In Julia, keyword arguments are documented (and implemented) by being after a ;, so step!(integrator;dt[,stop_at_tdt=false]) would mean keyword. But

https://diffeq.sciml.ai/stable/basics/integrator/#SciMLBase.step!

These are positional (because the second is required for the third), so the syntax is:

step!(integrator,time_steps[i]-time_steps[i-1],true)