Extracting derivatives (du/dt) from an ODE problem

Hossein_Haghi · July 18, 2019, 7:23pm

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.

ChrisRackauckas · July 18, 2019, 7:32pm

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)

Hossein_Haghi · July 18, 2019, 8:43pm

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)

ChrisRackauckas · July 19, 2019, 1:45am

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

Hossein_Haghi · July 19, 2019, 6:47pm

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

Hossein_Haghi · July 22, 2019, 5:46pm

@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.

BLI · July 22, 2019, 6:53pm

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

non-differentiable interpolation functions
too inaccurate?

Hossein_Haghi · July 22, 2019, 7:15pm

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.

BLI · July 22, 2019, 8:13pm

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

An array, where an element is addressed, e.g., as u[1], etc.
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.

ChrisRackauckas · July 23, 2019, 12:35am

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

Hossein_Haghi · July 23, 2019, 5:27pm

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.

donkeysaddle · April 2, 2020, 2:05am

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

Thanks,

DS

ChrisRackauckas · April 2, 2020, 2:07am

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

donkeysaddle · April 6, 2020, 2:35am

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[1] = dx = -(p[1]u[1] + p[2]*u[1])
  du[2] = dy = -p[3]*u[1] - p[4]

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

ChrisRackauckas · April 6, 2020, 2:45am

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.

donkeysaddle · April 6, 2020, 3:02am

Thanks so much.

Hossein_Haghi · April 6, 2020, 10:18pm

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

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

########################
## 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

Can you please advice how to go to higher d/dts?
I basically only need these values for higher du/dt at t=0 as initial conditions for another DDE I am trying to solve.

ChrisRackauckas · April 6, 2020, 10:37pm

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.

Hossein_Haghi · April 6, 2020, 10:52pm

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

it is still failing at cb5 with Vern9

ChrisRackauckas · April 6, 2020, 10:52pm

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:

github.com

SciML/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl#L1785-L1842


      
          end
          
          
## Vern7
          @def vern7unpack begin @unpack r011, r012, r013, r014, r015, r016, r017, r042, r043, r044, r045, r046, r047, r052, r053, r054, r055, r056, r057, r062, r063, r064, r065, r066, r067, r072, r073, r074, r075, r076, r077, r082, r083, r084, r085, r086, r087, r092, r093, r094, r095, r096, r097, r112, r113, r114, r115, r116, r117, r122, r123, r124, r125, r126, r127, r132, r133, r134, r135, r136, r137, r142, r143, r144, r145, r146, r147, r152, r153, r154, r155, r156, r157, r162, r163, r164, r165, r166, r167 = cache.tab.interp end
          
          
@def vern7pre0 begin
              @vern7unpack
              b1Θ = @evalpoly(Θ, 0, r011, r012, r013, r014, r015, r016, r017)
              b4Θ = @evalpoly(Θ, 0, 0, r042, r043, r044, r045, r046, r047)
              b5Θ = @evalpoly(Θ, 0, 0, r052, r053, r054, r055, r056, r057)
              b6Θ = @evalpoly(Θ, 0, 0, r062, r063, r064, r065, r066, r067)
              b7Θ = @evalpoly(Θ, 0, 0, r072, r073, r074, r075, r076, r077)
              b8Θ = @evalpoly(Θ, 0, 0, r082, r083, r084, r085, r086, r087)
              b9Θ = @evalpoly(Θ, 0, 0, r092, r093, r094, r095, r096, r097)
              b11Θ = @evalpoly(Θ, 0, 0, r112, r113, r114, r115, r116, r117)
              b12Θ = @evalpoly(Θ, 0, 0, r122, r123, r124, r125, r126, r127)
              b13Θ = @evalpoly(Θ, 0, 0, r132, r133, r134, r135, r136, r137)
              b14Θ = @evalpoly(Θ, 0, 0, r142, r143, r144, r145, r146, r147)
              b15Θ = @evalpoly(Θ, 0, 0, r152, r153, r154, r155, r156, r157)
              b16Θ = @evalpoly(Θ, 0, 0, r162, r163, r164, r165, r166, r167)

This file has been truncated. show original

Just need a soldier to make the next 4