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.
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 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.
The post An unexpected triangle first appeared on John D. Cook.