`ForwardDiff` fails to compute correct derivative

[nickeltingupta]

It seems I've stumbled onto a case where ForwardDiff.jl gives incorrect result.

I have a function in variable t whose derivative I want to calculate at t=0.0. Using TaylorSeries.jl gives correct result while using ForwardDiff.jl gives incorrect result. This happens only for t=0.0.

I know that ForwardDiff.jl gives incorrect result because I cross-checked with Mathematica (both by Taylor expansion and direct derivative).

The Julia code :

using TaylorSeries,LinearAlgebra,ZChop
using ForwardDiff

ϵ = 1.5+0.0im 
γ = 0.6+0.0im 
ω = 1.0+0.0im

Ω = sqrt(ω^2 - γ^2);

t=Taylor1(BigFloat,5)

RA(w)=((cos(π/4)/(cos(π/8)^2))*2.0)/(sqrt((1.0/(-ϵ^2 + 4.0 * γ^2 * Ω^2)) * ((ϵ^2 + 4.0 * γ^2 * Ω^2) * sec(π/8)^4 + (1/sqrt(γ^2 * Ω^2 - Ω^4 + sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2)))) * cosh((w* sqrt(γ^2 * Ω^2 - Ω^4 - sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))))/Ω) * sec(π/8)^4 
* ((-3.0 * ϵ^2 + 4.0 * γ^2 * Ω^2) * sqrt(γ^2 * Ω^2 - Ω^4 + sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))) * cosh((w* sqrt(γ^2 * Ω^2 - Ω^4 + sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))))/Ω) - im* sqrt(2.0) * ϵ * (ϵ^2 - 4.0 * γ^2 * Ω^2 - 2.0 * sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))) * sinh((w* sqrt(γ^2 * Ω^2 - Ω^4 + sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))))/Ω)) 
+ (sinh((w* sqrt(γ^2 * Ω^2 - Ω^4 - sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))))/Ω) * (32.0 * (ϵ^4 + 8.0 * γ^2 * Ω^4 * (γ^2 + Ω^2) - 2.0 * ϵ^2 * Ω^2 * (γ^2 + 3.0 * Ω^2)) * sin(π/8)^4 * sinh((w* sqrt(γ^2 * Ω^2 - Ω^4 + sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))))/Ω) + 2.0im* ϵ * sqrt(γ^2 * Ω^2 - Ω^4 + sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))) * (ϵ^2 - 4.0 * γ^2 * Ω^2 + 2.0 * sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))) * cosh((w* sqrt(γ^2 * Ω^2 - Ω^4 + sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))))/Ω) * sec(π/8)^2 * (-1.0 + tan(π/8)^2)))
/(sqrt(γ^2 * Ω^2 - Ω^4 - sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2))) * sqrt(γ^2 * Ω^2 - Ω^4 + sqrt(Ω^4 * (ϵ^2 - 4.0 * γ^2 * Ω^2)))))))

println(ForwardDiff.derivative(RA,0.0))
println(TaylorSeries.getcoeff(taylor_expand(RA,0.0,order=1),1))

#Output
#0.0 - 0.012968711710371502im
#0.0 - 1.3258252147247764im

The mathematica code to cross-check :

\[Epsilon] = 1.5;
\[Gamma] = 0.6;
\[Omega] = 1.0;
\[CapitalOmega] = \[Sqrt](\[Omega]^2 - \[Gamma]^2);

RetAmp[t_] := (Cos[\[Pi]/4]/Cos[\[Pi]/8]^2)*2.0/
   Sqrt[(1.0/(-\[Epsilon]^2 + 
         4.0*\[Gamma]^2*\[CapitalOmega]^2))*((\[Epsilon]^2 + 
          4.0*\[Gamma]^2*\[CapitalOmega]^2)*
        Sec[\[Pi]/8]^4 + (1/
          Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \[CapitalOmega]^4 + 
            Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                4.0*\[Gamma]^2*\[CapitalOmega]^2)]])*
        Cosh[(t*Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \[CapitalOmega]^4 \
- Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                  4.0*\[Gamma]^2*\[CapitalOmega]^2)]])/\[CapitalOmega]\
]*Sec[\[Pi]/
           8]^4*((-3.0*\[Epsilon]^2 + 
             4.0*\[Gamma]^2*\[CapitalOmega]^2)*
           Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \[CapitalOmega]^4 + 
             Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                 4.0*\[Gamma]^2*\[CapitalOmega]^2)]]*
           Cosh[(t*Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \
\[CapitalOmega]^4 + 
                 Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                    4.0*\[Gamma]^2*\[CapitalOmega]^2)]])/\
\[CapitalOmega]] - 
          I*Sqrt[2.0]*\[Epsilon]*(\[Epsilon]^2 - 
             4.0*\[Gamma]^2*\[CapitalOmega]^2 - 
             2.0*Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                  4.0*\[Gamma]^2*\[CapitalOmega]^2)])*
           Sinh[(t*Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \
\[CapitalOmega]^4 + 
                 Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                    4.0*\[Gamma]^2*\[CapitalOmega]^2)]])/\
\[CapitalOmega]]) + (Sinh[(t*
              Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \[CapitalOmega]^4 - 
                Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                    4.0*\[Gamma]^2*\[CapitalOmega]^2)]])/\
\[CapitalOmega]]*(32.0*(\[Epsilon]^4 + 
               8.0*\[Gamma]^2*\[CapitalOmega]^4*(\[Gamma]^2 + \
\[CapitalOmega]^2) - 
               2.0*\[Epsilon]^2*\[CapitalOmega]^2*(\[Gamma]^2 + 
                  3.0*\[CapitalOmega]^2))*Sin[\[Pi]/8]^4*
             Sinh[(t*Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \
\[CapitalOmega]^4 + 
                   Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                    4.0*\[Gamma]^2*\[CapitalOmega]^2)]])/\
\[CapitalOmega]] + 
            2.0*I*\[Epsilon]*
             Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \[CapitalOmega]^4 + 
               Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                   4.0*\[Gamma]^2*\[CapitalOmega]^2)]]*(\[Epsilon]^2 \
- 4.0*\[Gamma]^2*\[CapitalOmega]^2 + 
               2.0*Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                    4.0*\[Gamma]^2*\[CapitalOmega]^2)])*
             Cosh[(t*Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \
\[CapitalOmega]^4 + 
                   Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                    4.0*\[Gamma]^2*\[CapitalOmega]^2)]])/\
\[CapitalOmega]]*
             Sec[\[Pi]/8]^2*(-1.0 + 
               Tan[\[Pi]/
                  8]^2)))/(Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \
\[CapitalOmega]^4 - 
            Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                4.0*\[Gamma]^2*\[CapitalOmega]^2)]]*
          Sqrt[\[Gamma]^2*\[CapitalOmega]^2 - \[CapitalOmega]^4 + 
            Sqrt[\[CapitalOmega]^4*(\[Epsilon]^2 - 
                4.0*\[Gamma]^2*\[CapitalOmega]^2)]]))]


Series[RetAmp[t], {t, 0.0, 2}]
RetAmp'[t] /. t -> 0.0

[mcabbott]

Not sure why, but appears to be one more bug fixed by #481, as master branch gives the desired answer:

julia> println(ForwardDiff.derivative(RA,0.0))
0.0 - 1.3258252147247764im

(jl_v9rYMr) pkg> st
Status `/private/var/folders/yq/4p2zwd614y59gszh7y9ypyhh0000gn/T/jl_v9rYMr/Project.toml`
  [f6369f11] ForwardDiff v0.11.0-DEV `https://github.com/JuliaDiff/ForwardDiff.jl.git#master`

[nickeltingupta]

Not sure why, but appears to be one more bug fixed by #481, as master branch gives the desired answer:

julia> println(ForwardDiff.derivative(RA,0.0))
0.0 - 1.3258252147247764im

(jl_v9rYMr) pkg> st
Status `/private/var/folders/yq/4p2zwd614y59gszh7y9ypyhh0000gn/T/jl_v9rYMr/Project.toml`
  [f6369f11] ForwardDiff v0.11.0-DEV `https://github.com/JuliaDiff/ForwardDiff.jl.git#master`

I seem to be on v0.10.36 which is the latest version from what I understand. I guess you have a developer version?

[mcabbott]

Yes, ] add ForwardDiff#master will install v0.11.0-DEV.