 # ForwardDiff and custom structure in DifferentialEquations

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 = -0.5*u + u.f1
du = -0.5*u
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])
``````

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

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