ForwardDiff and custom structure in DifferentialEquations

New to Julia
1

Hi, I am using DifferentialEquations with callbacks to solve a differential equation with parameters that change at some point in time. The solution works fine, it is taken directly from the documentation for DifferentialEquations. Now I want to calculate the derivative of the obtained solution with respect to the initial conditions. Following the example from https://diffeq.sciml.ai/dev/analysis/sensitivity/#concrete_solve-Examples-1 I made a function solutionSim that takes only the initial condition. The issue is with the SimType structure:

LoadError: MethodError: no method matching SimType(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solutionSim),Float64},Float64,2},2}, ::Float64)
Closest candidates are:
SimType(!Matched::Array{T,1}, ::T) where T

How can I adjust SimType so it is able to accept dual numbers?

MWE:

using DiffEqSensitivity, DifferentialEquations, ForwardDiff

mutable struct SimType{T} <: DEDataVector{T}
    x::Array{T,1}
    f1::T
end

function f(du,u,p,t)
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
end

const tstop1 = [5.]
const tstop2 = [8.]


function condition(u,t,integrator)
  t in tstop1
end

function condition2(u,t,integrator)
  t in tstop2
end

function affect!(integrator)
  for c in full_cache(integrator)
    c.f1 = 1.5
  end
end

function affect2!(integrator)
  for c in full_cache(integrator)
    c.f1 = -1.5
  end
end

save_positions = (true,true)

cb = DiscreteCallback(condition, affect!, save_positions=save_positions)

save_positions = (false,true)

cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

cbs = CallbackSet(cb,cb2)

u0 = SimType([10.5;1.0], 0.0)
prob = ODEProblem(f,u0,(0.0,10.0))

const tstop = [5.;8.]
sol = DifferentialEquations.solve(prob,Tsit5(),callback = cbs, tstops=tstop)

#####Does not work from here
function solutionSim(x)
  u0=SimType(x,0.0)
  _prob = remake(prob,[u0])
  DifferentialEquations.solve(_prob,Tsit5(),callback = cbs, tstops=tstop)[end]
end
# #
out=zeros(2,2)
dxJac = ForwardDiff.jacobian!(out,solutionSim,[10.0 10.0])
2

OK, not sure whether I should reply to my own post, but here is something that works.

Redefinition of the structure:

mutable struct SimType11111{Real} <: DEDataVector{Real}
    x::Array{Real,1}
    f1::Float64
end

Redefinition of solutionSim:

function solutionSim(x)
  u0=SimType11111(x,0.0)
  _prob = ODEProblem(f,u0,(0.0,10.0))
  DifferentialEquations.solve(_prob,Tsit5(),callback = cbs, tstops=tstop)[end]
end
3

Yeah that works. It’s not the most efficient and you should instead define a conversion so that if x is made dual f1 auto-converts to a trivial dual, but this will also get the job done.

1 Like