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

Specific Domains Numerics
,
1

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

        # 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:

┌ 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.

2

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

1 Like
3

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?

4

Quite possibly.

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)

5

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

6

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

7

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?

8

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:

#+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.

1 Like
9

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

Seems like that is worth a ForwardDiff issue.

1 Like
10

Opened one on github

1 Like
11

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:

#+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