Let J(x) be the function plotted below.
This is the Bessel function J1, but we drop the subscript because it’s the only Bessel function we’re interested in for this post. You can think of J as a sort of damped sine.
We can create versions of J with different frequencies by multiplying the argument x by different constants. Suppose we create versions with three different frequencies — J(ax), J(bx), and J(cx) — and integrate their product. This defines a function f of the frequencies.
We can evaluate the integral defining f using Sonine’s formula [1]
where Δ(a, b, c) is the area of a triangle with sides a, b, c, if such a triangle exists, and zero otherwise.
It’s amazing that this formula takes three parameters with no intrinsic geometric meaning and out pops the area of a triangle with such sides.
Numerical (mis)application
It would be ill-advised, but possible, to invert Sonine’s formula and use it to find the area of a triangle; Heron’s formula would be a far better choice. But just for fun, we’ll take the ill-advised route.
from numpy import linspace, pi, sqrt, inf from scipy.special import j1 from scipy.integrate import quad def heron(a, b, c): s = (a + b + c)/2 return sqrt(s*(s-a)*(s-b)*(s-c)) def g(a, b, c): integrand = lambda x: j1(a*x) * j1(b*x) * j1(c*x) i, _ = quad(integrand, 0, inf, limit=500) return i def area(a, b, c): return g(a, b, c)*pi*a*b*c/2 print(area(3,4,5), heron(3,4,5))
SciPy’s quad
function has difficulty with the integration, and rightfully issues a warning. The code increases the limit
parameter from the default value of 50 to 500, improving the accuracy but not eliminating the warning. The area
function computes the error of a 3-4-5 triangle to be 5.9984 and the heron
function computes the exact value 6.
Update: I tried the numerical integration above in Mathematica and it correctly computed the integral to 10 decimal places with no help. I suspect it is detecting the oscillatory nature of the integral and using Levin’s integration method; when I explicit set the integration to be Levin’s method, I get the same result.
Impossible triangles
You could calculate the area of a triangle from Heron’s theorem or from Sonine’s theorem. The results are identical when the parameters a, b, and c form a valid triangle. But the two approaches diverge when a, b, and c do not form a triangle. If you pass in parameters like (3, 1, 1) then Heron’s theorem gives a complex number and Sonine’s theorem yields zero.
Related posts
[1] Discovered by Nikolay Yakovlevich Sonin (1849–1915). The formula is usually referred to as Sonine’s formula rather than Sonin’s formula, as far as I’ve seen. This variation is minor compared to what transliteration does to other Russian names.
Sonine’s formula is more general, extending to Bessel functions Jν with Re(ν) > ½. I chose ν = 1 for this post because the Sonin’s formula is simplest in this case.
Great post! The numerical error with SciPy’s quad function is a good reminder of how sensitive numerical methods can be. Have you explored using higher precision libraries like mpmath for such calculations?
I’d like to look into this further. The problem isn’t with precision, however, but rather with oscillations. The integrand oscillates forever, and the oscillations decay slowly, like x−3/2.