JuliaMath / Calculus.jl

Calculus functions in Julia
Other
278 stars 76 forks source link

Hessian low precision when specific values are small #150

Open tautomer opened 4 years ago

tautomer commented 4 years ago

I have a 1D system x coupled to a harmonic oscillator y with a frequency of omega, so f_yy should be omega^2.

However, when the frequency of the harmonic oscillator is fairly small, like 5e-5, H[2, 2] returned from Calculus.hessian is 0, though Calculus.second_derivative still give the right result.

Simple code to reproduce the problem

using Calculus

# parameters
omega = 5e-5
chi = 0.02
a = 1.9657
b = 25.2139
c = 9.04337

# potential
v(x) = 4.254987695360661e-5x^4 - 0.00278189309952x^2
# coupling
inter(x, q) = sqrt(2omega) * chi * (a*x-b*tanh(x/c)) * q
# harmonic oscillator 
vp(q) = 0.5 * omega^2 * q^2
# total potential
u(x, q) = inter(x, q) + v(x) + vp(q)
# hessian
h = Calculus.hessian(x -> u(x[1], x[2]))
# print
println(h([0.0, 0.0]), " hessian")
println(Calculus.second_derivative(vp, 0.0), " 2nd derivative")

The output is

[-0.005563827244039653 -0.0001644434349728066; -0.0001644434349728066 0.0] hessian
2.5e-9 2nd derivative

The output of the first println is [-0.005563827244039653 -0.0001644434349728066; -0.0001644434349728066 0.0] (H[2, 2] is 0), while the second one gives right value 2.5e-9.

For the time being, I just replace the element h[2, 2] = omega^2 as a work-around.

I understand the value vp gets really really small, so it's kind of difficult to obtain nuermical derivatives, but why would second_derivative work, yet hessian doesn't?