JuliaMath / JuliaMath/Calculus.jl

Hessian low precision when specific values are small

Open
#150 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
285
Forks
78
PR merge metrics
No merged PRs in 30d

Description

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

```julia
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

```text
[-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?

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.