JuliaMath / JuliaMath/Calculus.jl
Hessian low precision when specific values are small
- 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.