# Extracting derivatives (du/dt) from an ODE problem

I have multiple ODE problems that I am solving. where I need the solutions (u) and derivative of the solutions (du). For smaller ODEs it is practical for me to do the following

``````using DifferentialEquations

function SB(du,u,p,t)
du=@. u
du=@. ((-0.5*u^2)*(3-u/(p))+(1+(1-3*p)*u/p)*((p-p)/p+2*p/(p*p))*(p/u)^(3*p)-2*p/(p*u)-4*p*u/(p*u)-(1+u/p)*(p-p+p*sin(2*pi*p*t))/p-p*u*cos(2*pi*p*t)*2*pi*p/(p*p))/((1-u/p)*u+4*p/(p*p))
end

R0=2e-6
ps=250e3
f=2e6

u0=([R0 0])
tspan=(0,100/f)

p=[0.0725, 998, 1e-3,1481, 0, 1.01e5,7/5,f, R0, ps]

prob = ODEProblem(SB,u0,tspan,p)
@time u = solve(prob,Tsit5(),alg_hints=[:stiff],saveat=0.01/f,reltol=1e-8,abstol=1e-8)
t=u.t
u2=@. ((-0.5*u[2,:]^2)*(3-u[2,:]/(p))+(1+(1-3*p)*u[2,:]/p)*((p-p)/p+2*p/(p*p))*(p/u[1,:])^(3*p)-2*p/(p*u[1,:])-4*p*u[2,:]/(p*u[1,:])-(1+u[2,:]/p)*(p-p+p*sin(2*pi*p*t))/p-p*u[1,:]*cos(2*pi*p*t)*2*pi*p/(p*p))/((1-u[2,:]/p)*u[1,:]+4*p/(p*p))
``````

where u2 is bascially du in the SB function. This quickly becomes impractical as the size of my ODEs grow (>500 coupled ODEs with >500X500 matrices). Is there way to ask DifferentialEquations.jl package (or any other way) to export du[i]s as it is solving the ODEs? as a side note, it is important to mention that i need every time step that I am asking the solver to save the output at (saveat=0.01/f). I learned that using `integrator` and `get_du!` I should be able to do the job.

``````Rdot=zeros(50001,2)
u = init(prob,SSPRK22(),dt=1e-9,reltol=1e-8,abstol=1e-8)
solve!(u)
get_du!(Rdot,u)
``````

However the obtained results from get_du! do not match u2. Am I missing something? appreciate the help.

Yes, using the `SavingCallback`, as described in the link at StackOverflow:

Here’s the same example, but now saving the `du`s out:

``````using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
prob = ODEProblem((du,u,p,t) -> du .= u, rand(4), (0.0,1.0))
saved_values = SavedValues(Float64, Vector{Float64})
cb = SavingCallback((u,t,integrator)->(get_du(integrator), saved_values)
sol = solve(prob, Tsit5(), callback=cb)
``````

I am trying to run this as an example

``````using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
prob = ODEProblem((du,u,p,t) -> du .= u, rand(4), (0.0,1.0))
saved_values = SavedValues(Float64, Vector{Float64})
cb = SavingCallback((u,t,integrator)->(get_du(integrator), saved_values)
sol = solve(prob, Tsit5(), callback=cb)
``````

and expanding it to my code but I am not entirely sure what the error message means by this

``````UndefRefError: access to undefined reference
in expression starting at untitled:5
``````

line 5 is `sol = solve(prob, Tsit5(), callback=cb)`

``````using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
prob = ODEProblem((du,u,p,t) -> du .= u, rand(4), (0.0,1.0))
saved_values = SavedValues(Float64, Vector{Float64})
cb = SavingCallback((u,t,integrator)->get_du(integrator), saved_values, saveat = 0.1:0.1:1.0)
sol = solve(prob, Tsit5(), callback=cb)
saved_values
``````
1 Like

I tried the above suggested example and also used it in my problem as

``````
using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
function SB(du,u,p,t)
du=@. u
du=@. ((-0.5*u^2)*(3-u/(p))+(1+(1-3*p)*u/p)*((p-p)/p+2*p/(p*p))*(p/u)^(3*p)-2*p/(p*u)-4*p*u/(p*u)-(1+u/p)*(p-p+p*sin(2*pi*p*t))/p-p*u*cos(2*pi*p*t)*2*pi*p/(p*p))/((1-u/p)*u+4*p/(p*p))
end

R0=2e-6
ps=250e3
f=2e6

u0=([R0 0])
tspan=(0,100/f)

p=[0.0725, 998, 1e-3,1481, 0, 1.01e5,7/5,f, R0, ps]

prob = ODEProblem(SB,u0,tspan,p)
saved_values = SavedValues(Float64, Array{Float64})
cb = SavingCallback((u,t,integrator)->get_du(integrator), saved_values, saveat = 1e-9:1e-9:100/f)
sol =solve(prob,Tsit5(),saveat=1e-9,reltol=1e-8,abstol=1e-8,callback=cb)
``````

when I plot `saved_values.saveval`, the values do not change with `t` and are constant in both the provided example and also my code

@ChrisRackauckas Can you please take a look at this?

I am migrating all of my codes from MATLAB to Julia and getting du[i]s is one of my main problems. Appreciate the help.

Is it possible to differentiate the solution interpolation functions? Or…

• non-differentiable interpolation functions
• too inaccurate?

You mean like
(u[i+1]-u[i])/(t[i+1]-t[i]) ?

This usually ends up being too inaccurate (when I used it previously in MATLAB it was very inaccruate), also, the output vector (du[i]s) will have one less data point than the u[i]s.

No, that’s not what I mean. Observe that `u` has two interpretations:

1. An array, where an element is addressed, e.g., as `u`, etc.
2. An interpolation function addressed as, e.g., `u(1.5)`, i.e., evaluated at time 1.5, which does not need to be an element in the time vector.

If the interpolation function has continuous derivative, and is supported by, e.g., ForwardDiff, you can easily compute the derivative.

Alternatively, if the interpolation function supports imaginary time, you can use the complex Taylor series trick, see e.g., Prof Higham’s discussion in a past JuliaCon.

Don’t know if this is possible (support for AD, or admits imaginary time), and don’t know how accurate it will be.

The derivative of the interpolation is already hardcoded with `Val{1}` .

``````using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
function SB(du,u,p,t)
du=@. u
du=@. ((-0.5*u^2)*(3-u/(p))+(1+(1-3*p)*u/p)*((p-p)/p+2*p/(p*p))*(p/u)^(3*p)-2*p/(p*u)-4*p*u/(p*u)-(1+u/p)*(p-p+p*sin(2*pi*p*t))/p-p*u*cos(2*pi*p*t)*2*pi*p/(p*p))/((1-u/p)*u+4*p/(p*p))
end

R0=2e-6
ps=250e3
f=2e6

u0=([R0 0])
tspan=(0,100/f)

p=[0.0725, 998, 1e-3,1481, 0, 1.01e5,7/5,f, R0, ps]

prob = ODEProblem(SB,u0,tspan,p)
saved_values = SavedValues(Float64, Matrix{Float64})
cb = SavingCallback((u,t,integrator)->integrator(t,Val{1}), saved_values, saveat = 1e-9:1e-9:100/f)
sol =solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8,saveat=1e-9,callback=cb)

saved_values.t
saved_values.saveval
``````
2 Likes

Thanks! this is exactly what I needed!
In my applications, Julia has been 40-60 times faster than MATLAB so far! Really happy with DifferentialEquations.jl package.

2 Likes

@ChrisRackauckas,
Is this still that best way to do this? The corresponding documentation pages linked in the StackOverflow page no longer exist.

Thanks,

DS

Thanks, I fixed the links. Yes, still the same.

1 Like

Question regarding ` saved_values.saveval` from callback. I have an example, very similar to the one presented by Hossein_Haghi:

``````function model(du,u,p,t)
du = dx = -(pu + p*u)
du = dy = -p*u - p

end

tspan = (0.0,500.0)

u0 = [20.0;5.0]

p = [0.002,  0.0100, 0.0020,  1.71799639e-03]

prob = ODEForwardSensitivityProblem(model, u0, tspan, p)
saved_values = SavedValues(Float64, Array{Float64})
cb = SavingCallback((u, t, integrator)->integrator(t,Val{1}), saved_values, saveat=0.0:1:500.0)

sol = solve(prob,DP8(),saveat=1,callback=cb)
x,dp = extract_local_sensitivities(sol);
saved_values.t
saved_values.saveval
``````

What values/derivatives are contained in the arrays of, ` saved_values.saveval`? I know that the first two values of each array are the derivatives with respect to time of u1 and u2, but what are the other 8 values? Is it (du/dp)/dt, for each u and p?

Is there a way to just get the derivatives of the u’s w.r.t time and not the other values?

Thanks,

DS

Yup that’s exactly what they are.

like du/dt? That’s just `integrator(t,Val{1})`, or `sol(t,Val{1})` depending on which interface.

2 Likes

Thanks so much.

I am trying to the same thing but for higher orders of d/dt for my du and du in function as follows

``````function SB(du,u,p,t)
du=@. u
du=@. ((-0.5*u^2)*(3-u/(p))+(1+(1-3*p)*u/p)*((p-p)/p+2*p/(p*p))*(p/u)^(3*p)-2*p/(p*u)-4*p*u/(p*u)-(1+u/p)*(p-p+p*sin(2*pi*p*t))/p-p*u*cos(2*pi*p*t)*2*pi*p/(p*p))/((1-u/p)*u+4*p/(p*p))
end

########################
## Control Parameters ##
########################
ps=100e3
R0=2e-6
f=1e6

u0=([R0 0])

tspan=(0,1e-10)
p=[0.0725, 998, 1e-3,1481, 0, 1.01e5,7/5,f, R0, ps]
prob = ODEProblem(SB,u0,tspan,p)

saved_values = SavedValues(Float64, Array{Float64})
saved_values2=SavedValues(Float64, Array{Float64})
saved_values3=SavedValues(Float64, Array{Float64})
saved_values4=SavedValues(Float64, Array{Float64})
saved_values5=SavedValues(Float64, Array{Float64})
saved_values6=SavedValues(Float64, Array{Float64})
saved_values7=SavedValues(Float64, Array{Float64})

cb1 = SavingCallback((u,t,integrator)->integrator(t,Val{1}), saved_values , saveat = 0:span/f:TT) #R2
cb2 = SavingCallback((u,t,integrator)->integrator(t,Val{2}), saved_values2, saveat = 0:span/f:TT) #R3
cb3 = SavingCallback((u,t,integrator)->integrator(t,Val{3}), saved_values3, saveat = 0:span/f:TT) #R4
cb4 = SavingCallback((u,t,integrator)->integrator(t,Val{4}), saved_values4, saveat = 0:span/f:TT) #R5
cb5 = SavingCallback((u,t,integrator)->integrator(t,Val{5}), saved_values5, saveat = 0:span/f:TT) #R6
cb6 = SavingCallback((u,t,integrator)->integrator(t,Val{6}), saved_values6, saveat = 0:span/f:TT) #R7
cb7 = SavingCallback((u,t,integrator)->integrator(t,Val{7}), saved_values7, saveat = 0:span/f:TT) #R8

@time u = solve(prob,Tsit5(),saveat=1e-12,reltol=1e-12,abstol=1e-12,callback=CallbackSet(cb1,cb2,cb3,cb4,cb5))

``````

it can go up to cb4 but crashes when I attempt to go to cb5 and higher returning the following error

``````MethodError: no method matching hermite_interpolant!(::Array{Float64,2}, ::Float64, ::Float64, ::Array{Float64,2}, ::Array{Float64,2}, ::Array{Array{Float64,2},1}, ::Nothing, ::Type{Val{5}})
Closest candidates are:
hermite_interpolant!(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Nothing, !Matched::Type{Val{0}}) at C:\Users\The Hossein\.julia\packages\OrdinaryDiffEq\8Pn99\src\dense\generic_dense.jl:419
hermite_interpolant!(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, !Matched::Type{Val{0}}) at C:\Users\The Hossein\.julia\packages\OrdinaryDiffEq\8Pn99\src\dense\generic_dense.jl:427
hermite_interpolant!(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Nothing, !Matched::Type{Val{1}}) at C:\Users\The Hossein\.julia\packages\OrdinaryDiffEq\8Pn99\src\dense\generic_dense.jl:452
``````

I basically only need these values for higher du/dt at t=0 as initial conditions for another DDE I am trying to solve.

1 Like

You need to use a method with a really high order interpolation. Directly choose `Vern9` for this. `Tsit5` only has a 4th order interpolation, so it’s failing at 5 derivatives of the interpolant. I don’t know if we’ve hardcoded all 9 derivatives though, but if you need it, I could show you what PR to make. It’s not hard, just tedious.

1 Like

I could certainly try if you show me how to do it.

it is still failing at `cb5` with `Vern9`

We’ve only written out the first 4 derivatives of the polynomial, but we put it in a style where it’s clear how to do it:

Just need a soldier to make the next 4

1 Like