Inverse differential equations

Posted on 10 August 2026 by John

In science and engineering classes, you might describe a system using Newton’s laws and end up with a differential equation. You then solve the differential equation, analytically or numerically, to see how the solutions behave.

You might also do the opposite, especially in a mathematics class: look at what differential equation a set of functions satisfy in order to understand those functions.

Bessel functions came out of solving differential equations from astronomy. But then they turned out to be useful, not just in other areas of science, but in pure mathematics as well. So there are Bessel function users who came to the functions first and haven’t seen the differential equation they came from.

You can learn a lot about Bessel functions, and other special functions, by looking at their defining differential equation even if you’re not directly interested in the differential equation or the physical problem that motivated it.

Bessel functions satisfy

x² y″ + x y′ + (x² − n²) y = 0.

You can tell a lot about Bessel functions just by inspecting this equation without solving it. If we divide by x² and write the equation in the form

y″ + (p(x)/x) y′ + (q(x)/x²) y = 0

then p(x) = 1 and q(x) = x² − n². The indicial equation

r(r − 1) + p(0) r + q(0) = 0

reduces to

r² = n²

and so r = ±n. That alone tells us there are two solutions, one analytic at zero and one singular at zero. These are J_n and Y_n respectively. It also tells us the behavior of these functions as x goes to zero and as x goes to infinity. To find out more, look up “method of Frobinius.”

DNA and Bessel functions

Posted on 9 August 2026 by John

I was reading a book on the history of the discovery of the structure of DNA [1] and was surprised by a few passing references to Bessel functions.

According to Claude,

When X-rays are diffracted by a helical structure, the resulting diffraction pattern breaks into a series of horizontal “layer lines.” Cochran, Crick, and Vand showed mathematically that the diffracted amplitude on the n-th layer line is proportional to a Bessel function of the first kind, order n:

J_n(2π r R)

where r is the radius of the helix and R is the distance out from the center (the meridian) in the diffraction pattern.

The citation for this paragraph is a paper from 1952 [2] that amazingly is behind a paywall.

Update: Evan pointed out something fishy in the expression above: “Isn’t there a dimensional problem with J(2π r R)? Maybe 2πr/R instead?”

This was Claude’s response:

r is the radius of the helix, a real-space length (say, in Å).
R is the radial coordinate in reciprocal space — it has units of inverse length (Å⁻¹), because that’s how X-ray diffraction/Fourier space works (spatial frequency, not distance).

So Evan and Claude were both right. There is an implicit division happening, but crystallographers do write r R.

[1] Watson and Crick didn’t “discover DNA” as is commonly said. DNA was discovered in 1878. Watson and Crick discovered the structure of DNA in 1953.

[2] Cochran, W., Crick, F. H. C., & Vand, V. (1952). “The Structure of Synthetic Polypeptides. I. The Transform of Atoms on a Helix.” Acta Crystallographica, 5(5), 581–586.

Holonomic functions

Posted on 2 August 2026 by John

Yesterday I wrote that a lot of the special functions that pop up in mathematical physics are solutions to second order linear differential equations with polynomial coefficients. More generally, holonomic functions are defined to be those functions that are the solutions to linear differential equations, of any order, with polynomial coefficients.

Most special functions are holonomic. To quantify that statement, I went through the special functions covered in Abramowitz and Stegun. The large majority are holonomic, though some common functions like the gamma function are not holonomic.

This report goes through the functions in A&S. For those that are holonomic, it gives the differential equation that the function solves. The large majority of these equations are second order, but not all. And the coefficients are nearly always first or second order polynomials, rarely higher order.

Inverse factorial improved

Posted on 28 July 2026 by John

A couple years ago I wrote about how to compute the inverse of factorial. I used that code in writing the previous post because the post required solving the equation

⌊log₂(n!)⌋ ≥ b

given b. That is, given a number of bits b, find the smallest value of n such that n! ≥ 2^b.

What the code got right

Looking back on the code in that post, there are a few changes I’d like to make. But first of all, I’d like to point out something the post does right: instead of trying to solve

Γ(y) = x

it solves

log Γ(y) = log x.

That’s why the argument to inverse_log_gamma is logarg. That makes the code useful for values of x that would far exceed the maximum floating point value, such as in the calculations for the previous post.

What I’d change

Rounding

The function inverse_factorial from the old post solves finds the closest integer solution. It would be better for it to return the solution without rounding and then let the user round result if they want to. In my calculations in the previous post, I wanted to take the floor, not round.

Newton’s method

The code in the previous post uses the bisection method. This method is very safe, and fast enough for my purposes, but it could be made faster. Newton’s method is faster, but it can be ill-behaved if you don’t start close enough to the solution.

It’s safe to use Newton’s method to invert log Γ for two reasons. First, you can get a good starting point based on Stirling’s approximation. Second, and more importantly, log Γ is convex. Newton’s method will converge from any starting point when applied to a convex function. A little caution is necessary because log Γ is not convex everywhere, but it is convex on the positive real axis.

Another difficulty with Newton’s method is that you need to supply the derivative of the function whose root you’re trying to find. But this isn’t an issue here because the derivative of log Γ is the digamma function, which is implemented in SciPy.

Tolerance

Finally, the previous code used the default tolerance for deciding when to stop refining the solution. The revised method lets the user specify tolerance. It provides a default value, but that default is visible in the function call, not hidden down in SciPy.

Revised code

Here’s the revised code.

from scipy.special import gammaln, digamma
from scipy.optimize import newton

def inverse_log_gamma(logarg, tol=1e-12):
    assert(logarg > 0)    
    x0 = logarg / log(logarg + 1) + 1 if logarg > 1 else 2.0
    def f(z): return gammaln(z) - logarg
    return newton(f, x0, fprime=digamma, tol=tol)

def inverse_factorial(logarg):
    g = inverse_log_gamma(logarg)
    return g - 1

exp_q

Posted on 26 July 2026 by John

The function exp_q(x) is defined by taking the power series for exp(x) and keeping only the terms whose index is a multiple of q. For example, exp₂(x) keeps only the even-numbered terms in the exponential power series and so equals cosh(x).

$\exp_2(x) = 1 + \frac{x^2}{2!} + \frac{x^4}{4!} + \frac{x^6}{6!} + \cdots = \cosh(x)$

In general,

$\exp_q(x) = \sum_{n=0}^\infty [q \mid n] \frac{x^n}{n!} = \sum_{n=0}^\infty \frac{x^{nq}}{(nq)!}$

The first sum uses Iverson’s bracket notation: a Boolean expression in brackets denotes the function that returns 1 when the expression is true and zero when it is false. Here the bracket equals 1 when q divides n and is zero otherwise.

Closed forms

Let ω = exp(2πi / q). Then

$\exp_q(x) = \frac{1}{q}\sum_{k=0}^{q-1} \exp(\omega^k x)$

This lets us find closed-form expressions for exp_q(x). For example, when q = 4, ω = i and

$\exp_4(x) = \frac{1}{2}\left( \cosh(x) + \cos(x) \right)$

Here’s a proof of the identity above:

$\begin{align*} \frac{1}{q} \sum_{k=0}^{q-1} \exp(\omega^k x) &= \frac{1}{q} \sum_{k=0}^{q-1} \sum_{n=0}^\infty \frac{\omega^{kn}x^n}{n!} \\ &= \sum_{n=0}^\infty \left( \frac{1}{q} \sum_{k=0}^{q-1} \omega^{kn}\right) \frac{x^n}{n!} \\ &= \sum_{n=0}^\infty [q \mid n] \frac{x^n}{n!} \\ &= \exp_q(x) \end{align*}$

In the proof we used the identity

$\frac{1}{q} \sum_{k=0}^{q-1} \omega^{kn} = [q \mid n]$

which is important in deriving the properties of the discrete Fourier transform.

Differential equations

The first time I saw the function exp_q(x) was in differential equations, though I didn’t know at the time the function had a name.

When a course in differential equations gets to power series solutions, a common example or homework problem is to solve

$y^{(k)}(x) = y(x)$

for k = 3 or 4, i.e. to find a function that equals its third or fourth derivative.

If the initial conditions are

and

$y^\prime(0) = y^{\prime\prime}(0) = \cdots = y^{(k-1)}(0) = 0$

the unique solution to

$y^{(k)}(x) = y(x)$

is y(x) = exp_k(x).

Mathematica and Mittag-Leffler

Mathematica does not have a built-in function implementing exp_q(x), but it does have an implementation of the Mittag-Leffler function, and so thanks to a relation between this function and exp_q(x) you can implement the latter as

expq[x_, q_] := MittagLefflerE[q, x^q]

Combinatorics

The first time I saw the notation exp_q(x) was in combinatorics. I had intended to include an application from that book here, but I make that the topic for the next post.

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.

Special functions