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[1]=@. u[2]
    du[2]=@. ((-0.5*u[2]^2)*(3-u[2]/(p[4]))+(1+(1-3*p[7])*u[2]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[9]))*(p[9]/u[1])^(3*p[7])-2*p[1]/(p[2]*u[1])-4*p[3]*u[2]/(p[2]*u[1])-(1+u[2]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[8]*t))/p[2]-p[10]*u[1]*cos(2*pi*p[8]*t)*2*pi*p[8]/(p[2]*p[4]))/((1-u[2]/p[4])*u[1]+4*p[3]/(p[2]*p[4]))
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[4]))+(1+(1-3*p[7])*u[2,:]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[9]))*(p[9]/u[1,:])^(3*p[7])-2*p[1]/(p[2]*u[1,:])-4*p[3]*u[2,:]/(p[2]*u[1,:])-(1+u[2,:]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[8]*t))/p[2]-p[10]*u[1,:]*cos(2*pi*p[8]*t)*2*pi*p[8]/(p[2]*p[4]))/((1-u[2,:]/p[4])*u[1,:]+4*p[3]/(p[2]*p[4]))

where u2 is bascially du[2] 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 dus 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[1]=@. u[2]
    du[2]=@. ((-0.5*u[2]^2)*(3-u[2]/(p[4]))+(1+(1-3*p[7])*u[2]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[9]))*(p[9]/u[1])^(3*p[7])-2*p[1]/(p[2]*u[1])-4*p[3]*u[2]/(p[2]*u[1])-(1+u[2]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[8]*t))/p[2]-p[10]*u[1]*cos(2*pi*p[8]*t)*2*pi*p[8]/(p[2]*p[4]))/((1-u[2]/p[4])*u[1]+4*p[3]/(p[2]*p[4]))
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[1], 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[1]=@. u[2]
    du[2]=@. ((-0.5*u[2]^2)*(3-u[2]/(p[4]))+(1+(1-3*p[7])*u[2]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[9]))*(p[9]/u[1])^(3*p[7])-2*p[1]/(p[2]*u[1])-4*p[3]*u[2]/(p[2]*u[1])-(1+u[2]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[8]*t))/p[2]-p[10]*u[1]*cos(2*pi*p[8]*t)*2*pi*p[8]/(p[2]*p[4]))/((1-u[2]/p[4])*u[1]+4*p[3]/(p[2]*p[4]))
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.

1 Like