Track all variables within Julia DifferentialEquations.jl

Hi All,

I want to track all variables within Julia’s DifferentialEquation’s efficiently (easy and fast).
My current solution, on a test function, is:

function test2(du,u,p,t)
    
    global i
    global Out_values
    global sampletimes
    
    a,b,c = p
    
    d=a^0.1*(t+1)
    e=u[1]/a
    f=u[2]/d
    
    if t in sampletimes
        
        Out_values[1,i] = d
        Out_values[2,i] = e
        Out_values[3,i] = f
        i=i+1
    
    end
    
    du[1] = a*u[1]
    du[2] = d*u[2]
    du[3] = b*u[2] - c*u[3]
end

sampletimes = tspan[1]:0.3:tspan[2]
Out_values = Array{Float64}(undef, 3, 2*11)
i=1

prob = ODEProblem(test2,u0,tspan,p)
sol = solve(prob,saveat=0.3,tstops=sampletimes)

However, there are two problems:

  1. Somehow I get duplicated results (not sure why) with slightly different results, and
  2. I can’t track variables if I decide to use all results, i.e. sol = solve(prob).

Any help is appreciated.

  1. You could use this: Callback Library · DifferentialEquations.jl
    It kind of does the same as you do, but it works with adaptive step-sizes and avoids using globals.

  2. A second approach could also be to compute the additional values afterwards. For example you could define d_values = @. p.a^(0.1)*(sol.t), e_values = @. sol[1,:] ./ p.a and f_values = @. sol[2,:] ./ p.d… If you don’t have too many values to track, this might be easier.

To answer the questions:

  • You get duplicated results because solvers (especially implicit solvers) evaluate the function many times, to approximate the Jacobian, or for adaptive stepsizes. These additional evaluations are tracked in your approach. The callback library takes care of these things!
  • I don’t understand the second problem fully, but with the saving callbacks you could call the solve function with any arguments (i.e. also without the saveat). With the second approach it directly uses the sol array, so it shouldbe compatible in all cases.
1 Like

I don’t have experience with it, but this might fit the bill: https://github.com/jonniedie/SimulationLogs.jl

3 Likes

@SteffenPL, thanks for the answer. Is there a tutorial or an example of a similar problem where callbacks are applied?

@DrPapa, thank you very much for answering my question. I took a look at the package, and it seems it serves the purpose. I will investigate more later and check how it performs against the callback library that @SteffenPL proposed.

@SteffenPL, thanks for the answer. Is there a tutorial or an example of a similar problem where callbacks are applied?

It’s for example in the documentation when you scroll down.

For you case it looks like this:


function saving_callback(u, t, integrator)
    a, b, c = integrator.p
    d=a^0.1*(t+1)
    e=u[1]/a
    f=u[2]/d
    return (d, e, f)
end

saved_values = SavedValues(Float64, Float64, Float64)
cb = SavingCallback(saving_callback, saved_values)

sol = solve(prob, callback=cb)

now saved_values.saveval contains the values corresponding to the times saved_values.t.

2 Likes

Thank you very much, @SteffenPL. This helps me a lot to understand the callback. Best regards!