An unexpected triangle

Let J(x) be the function plotted below.

This is the Bessel function J₁, 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.

$f(a, b, c) = \int_0^\infty J(ax) \, J(bx) \, J(cx)\, dx$

We can evaluate the integral defining f using Sonine’s formula [1]

$f(a, b, c) = \frac{2\Delta(a, b, c)}{\pi abc}$

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.

[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.

Numerical (mis)application

Related posts

Leave a Reply