1

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? I learned that DiffEqSensitivity.jl package is able to provide du/dps to check the sensitivity of the model to p. is there something similar to extract du/dts?

Hossein Haghi
  • 39
  • 3
  • Why not call the `SB` function and extract the wanted components from the returned vector? – Lutz Lehmann Jul 12 '19 at 22:02
  • That is possible as well. However, I did not want to repeat the same operation twice in the process. I am assuming that the ODE solver does calculate du[i]s as it is solving the problem. I wanted to avoid repeating the same process and ask the solver to export the du[i]s in conjunction with u[i]s. – Hossein Haghi Jul 16 '19 at 19:20
  • No, you are not really repeating the evaluation, as the internal steps have nothing to do with the output values, the output is interpolated from the step data. However, it can often be cheaper to compute the derivatives as the derivatives of the interpolating polynomials. – Lutz Lehmann Jul 16 '19 at 20:59

2 Answers2

5

I would use two different components together. First, as you get to really large ODEs, you'll want to only save specific pieces of the solution, or reduced pieces. For this, the SavingCallback is very helpful.

http://diffeq.sciml.ai/latest/features/callback_library#SavingCallback-1

For example, the following solves an ODE and only saves the trace and the norm of the solution at each step:

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

Now you can use that to save what you need. The second piece is to use the integrator to get the derivatives. You can see that get_du! can be used to extract the current (already computed) derivative:

http://diffeq.sciml.ai/latest/basics/integrator#Misc-1

Additionally, you can make use of the interpolation on the integrator. integrator(t,Val{1}) will give the first derivative of the solution at the current t.

Chris Rackauckas
  • 18,645
  • 3
  • 50
  • 81
-2

@ChrisRackauckas I do need every time step that I am defining the solver to solve.

get_du!(out,integrator) gives me an array where all the points have the same value. am I making a mistake somewhere?

prob = ODEProblem(SB,u0,tspan,p) Rdot=zeros(50001,2) u = init(prob,SSPRK22(),dt=1e-9,reltol=1e-8,abstol=1e-8) solve!(u) get_du!(Rdot,u) U=u.sol basically derivative of the second output (du[2]) has to be equal to u2 defined in my previous post.`

Hossein Haghi
  • 39
  • 3
  • For continued discussion, use the [Julia Discourse](https://discourse.julialang.org/) or join the [Slack](https://slackinvite.julialang.org/) (preferred, the DiffEq channel is #diffeq-bridged) – Chris Rackauckas Jul 18 '19 at 14:15