# ODE solver fails with ForwardDiff, but works with central differences for first derivative of an input function

**URL:** <https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033>\
**Category:** Numerics\
**Tags:** ode, autodiff\
**Created:** [June 16, 2021, 4:50pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033 "2021-06-16T16:50:51Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![TS-CUBED](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ts-cubed/32/25292_2.png) [@TS-CUBED](https://discourse.julialang.org/u/TS-CUBED)\
**Post date:** [June 16, 2021, 4:50pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/1 "2021-06-16T16:50:52Z")

</div>

I am modelling the human cardiovascular system and solve equations similar to the windkessel equation in the example:

```julia
        # The 5-element WK formulation, also 4-element parallel for Ls = 0.0
        function wk5(dP, P, params, t)

            #=

            The 5-Element WK with serial L
            set Ls = 0 for 4-Element WK parallel

            Formulation for DifferentialEquations.jl

            P: solution vector (pressures p1 and p2)
            params: parameter array
                        [Rc, Rp, C, Lp, Ls]

            =#

            # Split parameter array:
            Rc, Rp, C, Lp, Ls = params

            dP[1] = (
                -Rc / Lp * P[1]
                + (Rc / Lp - 1 / Rp / C) * P[2]
                + Rc * (1 + Ls / Lp) * ForwardDiff.derivative(I, t)
                + I(t) / C
                )

            dP[2] = -1 / Rp / C * P[2] + I(t) / C

            return dP[1], dP[2]

        end

```

If I try to solve this, I get an error:

```julia
┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq /home/thor/.julia/packages/OrdinaryDiffEq/T2rXL/src/initdt.jl:81
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq /home/thor/.julia/packages/OrdinaryDiffEq/T2rXL/src/solve.jl:510
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ SciMLBase /home/thor/.julia/packages/SciMLBase/grNUR/src/integrator_interface.jl:325

```

If I use a Central Difference (e.g. the derivative from Calculus.jl), the ODE solver works nicely.

In the above example `I` is a function that returns a flow rate, for an argument `t`. `I` is second order differentiable.

When I compare the analytical, the CD, and the AutoDiff results outside the ODE solver, they are absolutely identical up the 10 decimals.

I have to admit, I am not too familiar with AutoDiff (I do understand the principle mathematically, but not the implementation), so I guess (hope) I’m missing something simple here.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 16, 2021, 6:07pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/2 "2021-06-16T18:07:48Z")

</div>

First try ForwardDiff with NaN-safe mode and see if that helps:

[https://juliadiff.org/ForwardDiff.jl/latest/user/advanced.html#Fixing-NaN/Inf-Issues-1](https://juliadiff.org/ForwardDiff.jl/latest/user/advanced.html#Fixing-NaN/Inf-Issues-1)

---

<div class="post-metadata">

**Author:** ![TS-CUBED](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ts-cubed/32/25292_2.png) [@TS-CUBED](https://discourse.julialang.org/u/TS-CUBED)\
**Post date:** [June 16, 2021, 6:56pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/3 "2021-06-16T18:56:14Z")

</div>

That did not fix the issue (will Julia recompile ForwardDiff.jl if I change prelude.jl?), but got me closer to a solution.

The way my function is defined seems to make ForwardDiff.derivative() produce NaN at t=0.0 - need to go back to the whiteboard there and see what causes that.

The function has a `( t / t_norm ) ^ n` term in it which would become complex for negative t, so I have a `rem( t + T, T)` in there for period repetition. Could that cause the derivative to become NaN?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 16, 2021, 7:59pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/4 "2021-06-16T19:59:46Z")

</div>

> [@TS-CUBED](#):
>
> The function has a `( t / t_norm ) ^ n` term in it which would become complex for negative t, so I have a `rem( t + T, T)` in there for period repetition. Could that cause the derivative to become NaN?

Quite possibly.

> [@TS-CUBED](#):
>
> will Julia recompile ForwardDiff.jl if I change prelude.jl?

Since you’re changing a compile time constant, I would restart Julia because there’s no guarantee of recompilation for that (all code may specialize on it)

---

<div class="post-metadata">

**Author:** ![TS-CUBED](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ts-cubed/32/25292_2.png) [@TS-CUBED](https://discourse.julialang.org/u/TS-CUBED)\
**Post date:** [June 16, 2021, 8:10pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/5 "2021-06-16T20:10:35Z")

</div>

I have restarted Julia, but was wondering if it was keeping precompiled libraries or so.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 16, 2021, 8:17pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/6 "2021-06-16T20:17:07Z")

</div>

If you modified it in the library it would trigger recompilation.

---

<div class="post-metadata">

**Author:** ![TS-CUBED](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ts-cubed/32/25292_2.png) [@TS-CUBED](https://discourse.julialang.org/u/TS-CUBED)\
**Post date:** [June 16, 2021, 9:20pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/7 "2021-06-16T21:20:36Z")

</div>

Cheers. I have tested it with another function that is essentially a series of cos and sin, and it works for that. So the problem is with the function.

Since it works for the CD method, I will stick with that for the moment.

Is there a performance benefit with the AutoDiff? Or is it simply the fact that it is essentially an exact solution that is the key difference?

---

<div class="post-metadata">

**Author:** ![TS-CUBED](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ts-cubed/32/25292_2.png) [@TS-CUBED](https://discourse.julialang.org/u/TS-CUBED)\
**Post date:** [June 16, 2021, 9:47pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/8 "2021-06-16T21:47:28Z")

</div>

Yes, it seems to be the modulo periodicity that seems to break ForwardDiff:

If I give it a very simple function that should be differentiable, but use a modulo for time, then this happens:

```julia
#+begin_src julia
function f(t)
    sin(pi * 2 * (t % 1))^2 #* (t % 1 < 0.5)
end
#+end_src

#+RESULTS:
: f (generic function with 1 method)

#+begin_src julia
ForwardDiff.derivative(f, 1)
#+end_src

#+RESULTS:
: NaN

#+begin_src julia
derivative(f, 1)
#+end_src

#+RESULTS:
: -5.736299254432828e-15

```

So the result from the CD is close to zero, as expected. But the ForwardDiff gives NaN. If I take the modulo out, then AutoDiff works.

So while I took care that my functions are continuous, smooth and differentiable at the period boundaries, AutoDiff will not work for them. It does work for cosine series, which are naturally periodic, but not with piecewise continuous functions with modulo periodicity.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 16, 2021, 10:44pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/9 "2021-06-16T22:44:29Z")

</div>

> [@TS-CUBED](#):
>
> Is there a performance benefit with the AutoDiff? Or is it simply the fact that it is essentially an exact solution that is the key difference?

Not really a performance benefit. Forward mode is mostly stability, though you can sometimes get maybe like a 4x if it SIMD’s better.

> [@TS-CUBED](#):
>
> If I give it a very simple function that should be differentiable, but use a modulo for time, then this happens:

Seems like that is worth a ForwardDiff issue.

---

<div class="post-metadata">

**Author:** ![TS-CUBED](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ts-cubed/32/25292_2.png) [@TS-CUBED](https://discourse.julialang.org/u/TS-CUBED)\
**Post date:** [June 17, 2021, 11:09am UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/10 "2021-06-17T11:09:19Z")

</div>

> [@ChrisRackauckas](#):
>
> Seems like that is worth a ForwardDiff issue.

Opened one on github

---

<div class="post-metadata">

**Author:** ![TS-CUBED](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ts-cubed/32/25292_2.png) [@TS-CUBED](https://discourse.julialang.org/u/TS-CUBED)\
**Post date:** [July 6, 2021, 9:17pm UTC](https://discourse.julialang.org/t/ode-solver-fails-with-forwarddiff-but-works-with-central-differences-for-first-derivative-of-an-input-function/63033/11 "2021-07-06T21:17:35Z")

</div>

Based on a reply on github stating that the problem is with how the Dual handles `mod` and `rem`, I tried floor division to get the remainder, and that works:

```julia
#+begin_src julia
function f(t)
    times = t ÷ 1
    t = t - times
    sin(pi * 2 * (t ))^2 * (t < 0.5)
end
#+end_src

#+RESULTS:
: f (generic function with 1 method)

#+begin_src julia
ForwardDiff.derivative(f, 1)
#+end_src

#+RESULTS:
: 0.0

#+begin_src julia
derivative(f, 1)
#+end_src

#+RESULTS:
: 0.00011952987976706194

```