Decoupling of time and state types in order to make automatic differentiation possible?

[luchr]

Hello, one often wants to take derivatives of the solution of an ODE w.r.t. the initial-state.

Using ODE.jl and ForwardDiff.jl this is only possible with a very very ugly hack:

using ODE
using ForwardDiff
mytest_rhs = (t,x) -> 5*sin.(x)

# Don't do this: This is ugly
###
import Base.convert;
function convert(::Type{Float64}, x::ForwardDiff.Dual)
    return x.value
end;
###

Sol = x0 -> ode45(mytest_rhs, x0, [0.0, 1.0])[2][end]
Sol(1.0)
ForwardDiff.derivative(Sol,1.0)

Do you see possibilities to get rid of the explicit type casting and inferring done by ODE.jl? Here are some examples:

1. errs = Float64[] in runge_kutta.jl
2. function hinit{T}(F, x0, t0::T, tend, p, reltol, abstol) where T is the type for the time and later with h0 = convert(T,h0) this h0 is converted to T even if h0 is calculated with the help of the state(!).

Why is this hack above bad? Because of this conversion at some places the partial-derivative informations are completely lost which may have negative impact on the accuracy of the values for the computed derivatives.

[ChrisRackauckas]

This is already supported in DifferentialEquations.jl via OrdinaryDiffEq.jl. Your example

using DiffEqBase, OrdinaryDiffEq 
# Or just using DifferentialEquations , but this is the minimal set of deps for this
using ForwardDiff
mytest_rhs = (t,x) -> 5*sin.(x)

function end_from_u0(u0)
  prob = ODEProblem(mytest_rhs,u0,(typeof(u0)(0.0),typeof(u0)(1.0)))
  sol = solve(prob,Tsit5())
  sol[end]
end

end_from_u0(1.0)

ForwardDiff.derivative(end_from_u0,1.0)

so I'm not going to put the work in here, but would review a PR.

The reason why the type for time should be a dual too is nuanced. Two places: initial dt estimation and stepsize control. Intuitively, the each dt actually changes with the initial condition changes, so it needs to be a dual. More concretely, dt needs to match the type of time in order for many things to work (ex: complex independent variables). However, in many cases dt is directly calculated from estimates of the solution (error estimates for adaptivity, slope estimates / error estimates for initial dt). So it needs to be a type (sometimes with units) that is compatible with time, yet be compatible after scaling values from the independent variable. So it needs to be Dual if the indvars are Dual.

[luchr]

OK. Thanks.