I have an extremely simple case scenario
r(s) = (1 + 0.3*cos(8π*s + 2*sin(2π*s)))*exp(1im*(2π*s - π/2))
r′(s) = ForwardDiff.derivative(r, s)
r′′(s) = ForwardDiff.derivative(r′, s)
h_ref = 1/(5*2^9)
T_ref = collect(0:h_ref:1 - h_ref)
lines(T_ref, real.(r′′.(T_ref)))
It happens that the plot has a weird singularity for a certain value, which I know is wrong (the function is continuous). So not sure what is going on here?
this is a numerical precision issue. specifically, I think the problem is cancellation in exp(1im*(2π*s - π/2)). does the issue go away if you write this as cispi(2s-1/2)?
It looks like the problem arises right at s=0.25:
julia> r′′(0.2499999999999)
65.84713137037353 - 86.15429650286814im
julia> r′′(0.25)
100.39690339578884 - 86.15429650270057im
Indeed this seems to avoid the problem at s=0.25:
julia> r(s) = (1 + 0.3*cos(8π*s + 2*sin(2π*s)))*cispi(2s - 1/2)
r (generic function with 1 method)
julia> r′′(0.2499999999999)
65.84713137037355 - 86.15429650286814im
julia> r′′(0.25)
65.84713137081351 - 86.15429650270057im
But I’m still a little suspicious — it’s not obvious to me why there should be a huge cancellation error in this result right at s=0.25.
In particular, consider:
julia> p(s) = exp(1im*(2π*s - π/2));
julia> p′(s) = ForwardDiff.derivative(p, s);
julia> p″(s) = ForwardDiff.derivative(p′, s);
julia> p″(0.2499999999999)
-39.47841760435743 + 2.4807694081270598e-11im
julia> p″(0.25) # correct answer is -(2π)² ≈ -39.47841760435743
0.0 + 0.0im
This looks like a bug to me, not cancellation error. I filed it as ForwardDiff.jl#653.