`ForwardDiff` fails to compute correct derivative
nickeltingupta opened this issue · 3 comments
nickeltingupta commented
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.3258252147247764imThe 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.0mcabbott commented
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 commented
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 commented
Yes, ] add ForwardDiff#master will install v0.11.0-DEV.