From Kepler to Bessel

Posted on 6 June 2026 by John

The previous post very briefly said that the integral representation for Bessel functions was motived by solving Kepler’s equation. This post will go into more detail.

Kepler’s equation

There are multiple ways to describe the position of a planet in an elliptical orbit around a star. For historical reasons, these descriptions have arcane names such as mean anomaly, true anomaly, and eccentric anomaly. This post explains how these three are related.

For this post, it is enough to say that often you know mean anomaly M and want to know eccentric anomaly E. These are related via Kepler’s equation

$M = E - e \sin E$
where e is the eccentricity of the orbit. You’d like to solve for E as a function of M and e, but there’s no elementary way to do that.

One way to solve Kepler’s equation is to take a guess at E and plug it into the right hand side of

$E = M + e \sin E$
to get a new E, and keep iterating until the two sides are closer together. I write more about this here.

Another approach to solving Kepler’s equation is to use Newton’s method. I write more about that here.

Still another approach is to expand E in a sine series and find the series coefficients. An advantage to this approach is that once you have the coefficients, you have an expression for E as a function of M, and you can plug in more values of M without having to solve Kepler’s equation for each value of M separately.

Sine series coefficients

Kepler’s equation is easy to solve at E = 0 and at E = π. In both cases, E = M. So the function E − M is zero at both ends of [0, π], which suggests we try to expand E − M in a sine series

$E - M = \sum_{n=1}^\infty a_n \sin nM$

We then calculate the Fourier coefficients a_n as usual.

$\begin{align*} a_n &= \frac{2}{\pi} \int_0^\pi (E-M) \sin(nM) \, dM \\ &= \frac{2}{n \pi} \int_0^\pi (E^\prime - 1) \cos(nM)\, dM \\ &= \frac{2}{n \pi} \int_0^\pi \cos(nM) E^\prime(M) \, dM \\ &= \frac{2}{n \pi} \int_0^\pi \cos\Big(nE - ne\sin(E)\Big) E^\prime(M) \, dM \\ &= \frac{2}{n} \left\{ \frac{1}{\pi} \int_0^\pi \cos\Big(nE - ne \sin(E)\Big)\, dE\right\} \\ &= \frac{2}{n} J_n(ne) \end{align*}$

The second line uses integration by parts. The third line uses Kepler’s equation. The last line uses the definition of the Bessel functions J_n given in the previous post.

Mr. Bessel’s eponymous functions

Posted on 5 June 2026 by John

Yesterday I wrote a post showing that the trapezoid rule evaluates the integral

$\int_{-\pi}^\pi \cos(\sin(x) + x)\, dx$

very efficiently. But how do we know what the exact integral is for comparison? If you ask Mathematica, it will tell you the integral equals −2π J₁(1) where J₁ is a Bessel function. This may seem like rabbit out of a hat, but it’s actually a simple calculation given the integral definition of Bessel functions:

$J_n(z) = \frac{1}{\pi}\int_0^\pi \cos(n\theta - z\sin(\theta))\, d\theta$

Since cosine is even, we can write our integral over [−π, π] as twice the integral over [0, π]. Then a change of variables turns this into the definition of J_n(z) with n = 1 and z = 1.

A deeper question is what have we accomplished by just giving a new name to essentially the same problem we started with. Another question is why in the world are Bessel functions defined as above.

As for what we’ve accomplished, we’ve related out integration problem to a very well-studied function. Bessel functions have been studied for two centuries and it’s easy to find software to evaluate them. Even the usually minimalist command line calculator bc has a function j(x, n) for evaluating J_n(x) for integer values of n. We could calculate our integral to 50 decimal places as follows.

~$ bc -l
>>> scale = 50
>>>  -8*a(1)*j(1,1)
-2.76491937476833705153256665538788207487495025542883

Note that bc doesn’t have a value of π built in, but a(x) evaluates the arctangent function, and π = 4 arctan(1).

There are multiple ways of defining Bessel functions. The three main ways would be in terms of their power series, in terms of the differential equations they solve, and in terms of their integral representation. Friedrich Bessel defined what we now call Bessel functions of the first kind, the J_n functions, in terms of their integral representations.

Why did Bessel care about these integrals? They came out of his calculations in celestial mechanics. One example of this is solving Kepler’s equation with Fourier series; the Fourier coefficients are given by Bessel functions. This is worked out in detail in the next post.

Bessel functions had occurred in applications before Mr. Bessel drew attention to them. The functions are named after him because he was the first to systematically study them.

Mathematics is developed inductively but taught deductively. It’s common for things to be kicked around for years before someone decides they deserve a name and systematic study. See this post on the central limit theorem for another example. The CLT is older than the Gaussian distribution, even older than Gauss.

Couth and uncouth function pairs

Posted on 21 May 2026 by John

“You can’t always get what you want. But sometimes you get what you need.” — The Rolling Stones

Circular functions and hyperbolic functions aren’t invertible, but we invert them anyway. These functions map many points in the domain to each point in the range, and we invert them by mapping a point in the range back to some point in the domain. Often this works as expected, but sometimes it doesn’t.

In the previous post I said

You can relate each trig function “foo” with its hyperbolic counterpart “fooh” by applying one of the functions to iz then multiplying by a constant c that depends on foo: c = i for sin and tan, c = 1 for cos and sec, and c = −i for csc and cot.

In symbols,

c foo(z) = fooh(iz).

Let’s suppose foo and fooh are invertible, ignoring any complications, and solve foo(z) = w for z. We get

i foo⁻¹(w) = fooh⁻¹(cw)

or using “arc” nomenclature for inverse functions

i arcfoo(w) = arcfooh(cw).

When the naive calculation above holds, except possibly at a finite number of points, we say the pair (foo, fooh) is couth. Otherwise we say the pair is uncouth. These term were coined by Robert Corless and his coauthors in their paper [1].

Whether the pair (foo, fooh) is couth depends not only on foo and fooh, but also on the details of how arcfoo and arcfooh are defined.

In Python’s NumPy library, the pairs (sin, sinh) and (tan, tanh) are couth, but the pair (cos, cosh) is uncouth.

Numpy doesn’t define the reciprocal functions sec, sech, csc, csch, cot, and coth. I used to find that annoying, but I’m beginning to think that was wise. These functions cause problems. For example, there may be two reasonable ways to define these functions, one of which forms a couth pair and one of which forms an uncouth pair.

For example, how should you define cot and coth? There would be no disagreement over the definition

cot = lambda x: 1/tan(x)

but there are at least two definitions of inverse coth that you’ll find in practice:

arccot = lambda z: 0.5*pi - arctan(z)
arccot = lambda z: arctan(1/z).

Both definitions have their advantages, but the former is uncouth while the latter is couth. You can verify that both definitions are the same at z = 1 but not at z = −1.

With the following definitions, the pairs (cos, cosh) and (sec, sech) are uncouth and the rest are couth.

from numpy import *

csc     = lambda x: 1/sin(x)
sec     = lambda x: 1/cos(x)
cot     = lambda x: 1/tan(x)
csch    = lambda x: 1/sinh(x)
sech    = lambda x: 1/cosh(x)
coth    = lambda x: 1/tanh(x)

arccot  = lambda z: arctan(1/z)
arcsec  = lambda z: arccos(1/z)
arccsc  = lambda z: arcsin(1/z)
arccoth = lambda z: arctanh(1/z)
arcsech = lambda z: arccosh(1/z)
arccsch = lambda z: arcsinh(1/z)

[1] “According to Abramowitz and Stegun” or arccoth needn’t be uncouth. Robert M. Corless et al. ACM SIGSAM Bulletin, Volume 34, Issue 2, pp 58 – 65 https://doi.org/10.1145/362001.362023

Approximating even functions by powers of cosine

Posted on 30 April 2026 by John

A couple days ago I wrote a post about turning a trick into a technique, finding another use for a clever way to construct simple, accurate approximations. I used as my example approximating the Bessel function J(x) with (1 + cos(x))/2. I learned via a helpful comment on Mathstodon that my approximation was the first-order part of a more general series

$J_0(x) = 1 + \frac{\cos(x) - 1}{2} - \frac{(\cos(x) - 1)^2}{48} + \frac{7(\cos(x) - 1)^3}{1440} + \cdots$

The first-order approximation has error O(x⁴), as shown in the earlier post. Adding the second-order term makes the error O(x⁶), and adding the third-order term makes it O(x⁸).

I’ve written a few times about cosine approximations to the normal probability density. For example, see this post. We could use the same idea as the series above to approximate the normal density with a series of powers of cosine. This gives us

$\exp(-x^2/2) = 1 + (\cos(x) - 1) + \frac{(\cos(x) - 1)^2}{3} + \frac{2(\cos(x) - 1)^3}{45} + \cdots$

and as before, the first, second, and third order truncated series have error O(x⁴), O(x⁶), and O(x⁸).

The general theory behind what’s going on here is an extension of Bürmann’s theorem. The original version of the theorem relies on a series inversion theorem that in turn relies on the approximating function, in our case cos(x) − 1, not having zero derivative at the center of the series. But there is a more general form of Bürmann’s theorem based on a more general form of series inversion. We will always need a more general version of the theorem when working with even functions because even functions have zero derivative at zero.

Here’s another example, this time using the Bessel function J₁, an odd function, which does use the original version of Bürmann’s theorem to approximate J₁ by powers of sine.

$J_1(x) = \frac{1}{2} \sin(x) + \frac{1}{48} \sin^3(x) + \frac{17}{1920} \sin^5(x) + \cdots$

In this case truncating the series after sin^k(x) gives an error O(x^{k + 2}).

You can find more on Bürmann’s theorem in Whittaker and Watson.

Turning a trick into a technique

Posted on 28 April 2026 by John

Someone said a technique is a trick that works twice.

I wanted to see if I could get anything interesting by turning the trick in the previous post into a technique. The trick created a high-order approximation by subtracting a multiple one even function from another. Even functions only have even-order terms, and by using the right multiple you can cancel out the second-order term as well.

For an example, I’d like to approximate the Bessel function J₀(x) by the better known cosine function. Both are even functions.

J₀(x) = 1 − x²/4 + x⁴/64 + …
cos(x) = 1 − x²/2 + x⁴/24 + …

and so

2 J₀(x) − cos(x) = 1 − x⁴/96 + …

which means

J₀(x) ≈ (1 + cos(x))/2

is an excellent approximation for small x.

Let’s try this for a couple examples.

J₀(0.2) = 0.990025 and (1 + cos(0.2))/2 = 0.990033.

J₀(0.05) = 0.99937510 and (1 + cos(0.05))/2 = 0.99937513.

Closed-form solution to the nonlinear pendulum equation

Posted on 25 April 2026 by John

The previous post looks at the nonlinear pendulum equation and what difference it makes to the solutions if you linearize the equation.

If the initial displacement is small enough, you can simply replace sin θ with θ. If the initial displacement is larger, you can improve the accuracy quite a bit by solving the linearized equation and then adjusting the period.

You can also find an exact solution, but not in terms of elementary functions; you have to use Jacobi elliptic functions. These are functions somewhat analogous to trig functions, though it’s not helpful to try to pin down the analogies. For example, the Jacobi function sn is like the sine function in some ways but very different in others, depending on the range of arguments.

We start with the differential equation

θ″(t) + c² sin( θ(t) ) = 0

where c² = g/L, i.e. the gravitational constant divided by pendulum length, and initial conditions θ(0) = θ₀ and θ′(0) = 0. We assume −π < θ₀ < π.

Then the solution is

θ(t) = 2 arcsin( a cd(ct | m ) )

where a = sin(θ₀/2), m = a², and cd is one of the 12 Jacobi elliptic functions. Note that cd, like all the Jacobi functions, has an argument and a parameter. In the equation above the argument is ct and the parameter is m.

The last plot in the previous post was misleading, showing roughly equal parts genuine difference and error from solving the differential equation numerically. Here’s the code that was used to solve the nonlinear equation.

from scipy.special import ellipj, ellipk
from numpy import sin, cos, pi, linspace, arcsin
from scipy.integrate import solve_ivp

def exact_period(θ):
    return 2*ellipk(sin(θ/2)**2)/pi

def nonlinear_ode(t, z):
    x, y = z
    return [y, -sin(x)]    

theta0 = pi/3
b = 2*pi*exact_period(theta0)
t = linspace(0, 2*b, 2000)

sol = solve_ivp(nonlinear_ode, [0, 2*b], [theta0, 0], t_eval=t)

The solution is contained in sol.y[0].

Let’s compare the numerical solution to the exact solution.

def f(t, c, theta0):
    a = sin(theta0/2)
    m = a**2
    sn, cn, dn, ph = ellipj(c*t, m)
    return 2*arcsin(a*cn/dn)

There are a couple things to note about the code. First,SciPy doesn’t implement the cd function, but it can be computed as cn/dn. Second, the function ellipj returns four functions at once because it takes about as much time to calculate all four as it does to compute one of them.

Here is a plot of the error in solving the differential equation.

And here is the difference between the exact solution to the nonlinear pendulum equation and the stretched solution to the linear equation.

A lesser-known characterization of the gamma function

Posted on 18 March 2026 by John

The gamma function Γ(z) extends the factorial function from integers to complex numbers. (Technically, Γ(z + 1) extends factorial.) There are other ways to extend the factorial function, so what makes the gamma function the right choice?

The most common answer is the Bohr-Mollerup theorem. This theorem says that if f: (0, ∞) → (0, ∞) satisfies

f(x + 1) = x f(x)
f(1) = 1
log f is convex

then f(x) = Γ(x). The theorem applies on the positive real axis, and there is a unique holomorphic continuation of this function to the complex plane.

But the Bohr-Mollerup theorem is not the only theorem characterizing the gamma function. Another theorem was by Helmut Wielandt. His theorem says that if f is holomorphic in the right half-plane and

f(z + 1) = z f(z)
f(1) = 1
f(z) is bounded for {z: 1 ≤ Re z ≤ 2}

then f(x) = Γ(x). In short, Wielandt replaces the log-convexity for positive reals with the requirement that f is bounded in a strip of the complex plane.

You might wonder what is the bound alluded to in Wielandt’s theorem. You can show from the integral definition of Γ(z) that

|Γ(z)| ≤ |Γ(Re z)|

for z in the right half-plane. So the bound on the complex strip {z: 1 ≤ Re z ≤ 2} equals the bound on the real interval [1, 2], which is 1.

Γ(1/n)

Posted on 4 February 2026 by John

If n is a positive integer, then rounding Γ(1/n) up to the nearest integer gives n. In symbols,

$\left\lceil \Gamma\left( \tfrac{1}{n}\right) \right\rceil = n$

We an illustrate this with the following Python code.

>>> from scipy.special import gamma
>>> from math import ceil
>>> for n in range(1, 101):
    ... assert(ceil(gamma(1/n)) == n)

You can find a full proof in [1]. I’ll give a partial proof that may be more informative than the full proof.

The asymptotic expansion of the gamma function near zero is

$\Gamma(z) = \frac{1}{z} - \gamma + {\cal O}(z^2)$

where γ is the Euler-Mascheroni constant.

So when we set z = 1/n we find Γ(1/n) ≈ n − γ + O(1/n²). Since 0 < γ < 1, the theorem above is true for sufficiently large n. And it turns out “sufficiently large” can be replaced with n ≥ 1.

[1] Gamma at reciprocals of integers: 12225. American Mathematical Monthly. October 2022. pp 789–790.

Trigamma

Posted on 1 September 2025 by John

Guy wearing a ΓΓΓ shirt

The most important mathematical function after the basics is the gamma function. If I could add one function to a calculator that has trig functions, log, and exponential, it would be the gamma function. Or maybe the log of the gamma function; it’s often more useful than the gamma function itself because it doesn’t overflow as easily.

The derivative of the log gamma function is the digamma function, denoted ψ. It comes up often in application. I just did a quick search and found I’ve written six posts containing the word “digamma.”

$\psi(z) = \frac{d}{dz} \log \Gamma(z)$

The derivative of the digamma function ψ′ is the trigamma function.

$\psi^\prime(z) = \frac{d}{dz} \psi(z) = \frac{d^2}{dz^2} \log \Gamma(z)$

The trigamma function, and higher derivatives of the digamma function, appear in applications. I remember, for example, a researcher asking me to add the trigamma function to the mathematical library I wrote for the biostatistics department at MD Anderson.

I was thinking about the trigamma function because I ran across a the following series for the function [1].

$\psi^\prime (z+1) = \frac{1}{(z + 1)} + \frac{1!}{2} \frac{1}{(z+1)^{\bar{2}}} + \frac{2!}{3}\frac{1}{(z+1)^{\bar{3}}} + \cdots$

Note the bars on top of the exponents: the denominators are rising powers of z + 1, not ordinary powers.

The series converges uniformly for Re(z) > −1 + δ for δ > 0 [2]. It series converges quickly for large z.

When I saw the title of the paper I thought it sounded like a Greek fraternity. There is a Tri-Delta fraternity, but as far as I know there is no Tri-Gamma fraternity.

[1] Harold Ruben. A Note on the Trigamma Function. The American Mathematical Monthly. Vol 83, No. 8. p. 622.

[2] It may seem unnecessary to say Re(z) > −1 + δ for δ > 0. Couldn’t you just say for Re(z) > −1? Pointwise, yes, but uniform convergence requires the real part of z to be bounded away from −1 by a fixed amount, regardless of the imaginary part of z.

An integral theorem of Gauss

Posted on 31 August 2025 by John

Gauss proved in 1818 that the value of integral

$\int_0^{\pi/2} \left( x^2 \sin^2\theta + y^2 \cos^2 \theta \right)^{-1/2} \, d\theta$

is unchanged if x and y are replaced by (x + y)/2 and √(xy), i.e. if you replaced x and y with their arithmetic mean and geometric mean [1].

So, for example, if you wanted to compute

$\int_0^{\pi/2} \left( 9 \sin^2\theta + 49 \cos^2 \theta \right)^{-1/2} \, d\theta$

you could instead compute

$\int_0^{\pi/2} \left( 25 \sin^2\theta + 21 \cos^2 \theta \right)^{-1/2} \, d\theta$

Notice that the coefficients of sin² θ and cos² θ are more similar in the second integral. It would be nice if the two coefficients were equal because then the integrand would be a constant, independent of θ, and we could evaluate the integral. Maybe if we apply Gauss’ theorem again and again, the coefficients will become more nearly equal.

We started with x = 3 and y = 7. Then we had x = 5 and y = √21 = 4.5826. If we compute the arithmetic and geometric means again, we get x = 4.7913 and y = 4.7874. If we do this one more time we get x = y = 4.789013. The values of x and y still differ, but only after the sixth decimal place.

It would seem that if we keep replacing x and y by their arithmetic and geometric means, we will converge to a constant. And indeed we do. This constant is called the arithmetic-geometric mean of x and y, denoted AGM(x, y). This means that the value of the integral above is π/ (2 AGM(x, y)). Because the iteration leading to the AGM converges quickly, this provides an efficient numerical algorithm for computing the integral at the top of the post.

This is something I’ve written about several times, though in less concrete terms. See, for example, here. Using more advanced terminology, the AGM gives an efficient way to evaluate elliptic integrals.

[1] B. C. Carlson, Invariance of an Integral Average of a Logarithm. The American Mathematical Monthly, Vol. 82, No. 4 (April 1975), pp. 379–382.

Special functions

From Kepler to Bessel

Kepler’s equation

Sine series coefficients

Mr. Bessel’s eponymous functions

Related posts

Couth and uncouth function pairs

Approximating even functions by powers of cosine

Turning a trick into a technique

Closed-form solution to the nonlinear pendulum equation

A lesser-known characterization of the gamma function

Γ(1/n)

Trigamma

Related posts

An integral theorem of Gauss

Related posts