How to extract the right-hand side function from a differential equation solution?

Specific Domains Numerics
1

Hey,

Basically, I want to extract the right hand side function f from the solution of a differential equation \dot{x} = f(x). For a minimum working example consider the following


julia> using DifferentialEquations

julia> myfunc(dx, x, u, t) = (dx .= -x)

myfunc (generic function with 1 method)

julia> prob = ODEProblem(myfunc, rand(1), (0.0, 1.0))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true

timespan: (0.0, 1.0)

u0: 1-element Vector{Float64}:

0.848128719065019

julia> sol = solve(prob)

retcode: Success

Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation

t: 5-element Vector{Float64}:

0.0

0.10002357021250141

0.34209643917644617

0.6554220498434524

1.0

u: 5-element Vector{Vector{Float64}}:

[0.848128719065019]

[0.7674005123680626]

[0.6024086282879929]

[0.4403679037389535]

[0.3120092483176733]

I remember that it was once possible to extract the function f from sol as


f = sol.prob.f.f

(maybe wrong, not hundred percent sure!). However, for DifferentialEquations@v7.3.0, the type of sol.prob.f.f is


julia> typeof(sol.prob.f.f)

FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}

not typeof(myufunc) as I expected. So, I defined the following function


extract_f(sol) = first(sol.prob.f.f.fw).obj[].f.f

to extract f from sol.


julia> extract_f(sol) = first(sol.prob.f.f.fw).obj[].f.f

extract_f (generic function with 1 method)

julia> extract_f(sol)

myfunc (generic function with 1 method)

What I want to ask is this: Is there a simpler and maybe tidier way of getting f from sol?

2

You could use the Integrator interface

julia> integrator = init(sol.prob, sol.alg)

integrator.f is the rhs.

1 Like
3

It is definitely a little bit tidier approach, it initializes a new integrator though. :frowning_face:

4

SciMLBase.unwrapped_f(sol.prob.f), though using undocumented internals is always bound to change.

1 Like
5

Sure. That is why I am having problems now :smile:

6

I think I should take some more steps to access f from here

f = SciMLBase.unwrapped_f(sol.prob.f)  # ODEFunction
f.f.f # myfunc
7

That will generate an unwrapped ODEFunction. SciMLBase.unwrapped_f(sol.prob.f.f) is probably as far as you’d ever want to go.

8

Thank you @ChrisRackauckas for your help. But, for the completeness of the thread, I just want to point out this.

Actually, SciMLBase.unwrapped_f(sol.prob.f.f) is not quite the final point that I wanted to reach. Being different than the MWE given above, in my use case the right-hand function f was actually a functor, e.g.

mutable struct MyFunc{T}
    a::T 
end

and I wanted to reach an instance myfunc of MyFunc from the solution sol for the mutation of its parameter myfunc.a . So, what I wanted was not the ODEFunction, but myfunc instead.

1 Like
9

Ahh okay, and for completeness of the thread we found an issue with the setup before, fixed it (Fix and better test wrapping levels for AutoSpecialize by ChrisRackauckas · Pull Request #820 · SciML/DiffEqBase.jl · GitHub), and now have tests and documentation on such unwrapping functionality.

2 Likes