# How to debug failing DifferentialEquations.solve when called with Dual numbers

**URL:** <https://discourse.julialang.org/t/how-to-debug-failing-differentialequations-solve-when-called-with-dual-numbers/79895>\
**Category:** Numerics\
**Tags:** diffeq, turing, forwarddiff, modelingtoolkit\
**Created:** [April 23, 2022, 2:43pm UTC](https://discourse.julialang.org/t/how-to-debug-failing-differentialequations-solve-when-called-with-dual-numbers/79895 "2022-04-23T14:43:07Z")\
**Posts on this page:** 1\
**Showing post:** 6

<div class="post-metadata">

**Author:** ![bgctw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bgctw/32/22050_2.png) [@bgctw](https://discourse.julialang.org/u/bgctw)\
**Post date:** [April 27, 2022, 7:47am UTC](https://discourse.julialang.org/t/how-to-debug-failing-differentialequations-solve-when-called-with-dual-numbers/79895/6 "2022-04-27T07:47:20Z")

</div>

In addition to advise given in the [PSA](https://discourse.julialang.org/t/psa-how-to-help-yourself-debug-differential-equation-solving-issues/62489) the following strategies helped me debugging.

**Inspecting how the derivate function is called by the solver by wrapping it.**  
Although I was not able to peek into the derivative function generated by MTK, I was able to replace it by a wrapper in ODEProblem.

```julia
df = probos.f
df2 = (du,u,p,t) -> begin
    df(du,u,p,t)
    @show t,u, du
end
probos_d = remake(probos, f = df2)
sol = solve(probos_d, solver);

```

This allowed me to get to the actual states that are tried by the solver - which I cannot see in the solution that ends earlier when problems occur.

**Directly calling the derivative function with Dual numbers at specified state and parameters**

```julia
dud = PreallocationTools.dualcache(probos.u0)
fo = (p) -> begin
    _du = PreallocationTools.get_tmp(dud, first(p))
    probos.f(_du,probos.u0,p,0)
    @show _du
    sum(abs2, _du)
end 
dfop = ForwardDiff.gradient(fo, probos.p);

```

This allowed me to explore the problem in a more agile manner.

**Overriding equations of an assembled ODESystem**  
Rather than changing the source code of the different components and reloading (Revise) and reassembling the system, I was able [modify the ODESystem](https://github.com/SciML/ModelingToolkit.jl/issues/1518) and experiment by replacing specific equations.

```julia
function sesam_constI2(;name, s = sesam3(;name=:s))
    @variables t
    D = Differential(t)
    @unpack α_RP, return_RP, revenue_RP, syn_Enz_w = s
    eqs = [
        revenue_RP ~ IfElse.ifelse(α_RP > 1e-16, return_RP / (α_RP * syn_Enz_w), 0.0),
    ]
    basesys = s
    sys = ODESystem(eqs, t, [], []; name)
    # extend_overwrite so far only in extend_overwrite.jl in project dir
    sys_ext = extend_overide(sys, basesys) 
end
s2 = sesam3(name=:s) 
sci2 = sesam_constI2(;s=s2,name=:s); equations(sci2)
@named spci2 = plant_sesam_system(sci2,pl) # couples to drivers and calls simplify
probo2 = ODEProblem(spci2, u0site, tspan, psite)
probos = remake(probo2, u0 = ForwardDiff.value.(u)) # u from steps above

```

---

_[View the full topic](https://discourse.julialang.org/t/how-to-debug-failing-differentialequations-solve-when-called-with-dual-numbers/79895)._