What exactly is modified about a modified Bessel function?

Posted on 23 August 2026 by John

Special functions often have arcane names that not very helpful without some context. The previous post goes into some reasons for this. This post will expand on a point at the end of the post about “modified” functions.

Things are given their names for reasons. Discovering those reasons may help you understand their motivation and use.

Pure math perspective

For each integer n, the modified Bessel function I_n is essentially the Bessel function J_n evaluated along the imaginary axis. Specifically,

$I_n(x) = i^{-n} J_n(ix)$

From a certain shallow perspective, that’s the end of the story: modified Bessel functions are modified in the sense that the argument is multiplied by i. And there’s a fiddly constant term up front for no apparent reason.

But of course that’s not the end of the story or else this wouldn’t be worth an entire post.

The equation above is analogous to the relationships between circular and hyperbolic functions

$\begin{align*} \sin(ix) &= i \sinh(x) \\ \cos(ix) &= \phantom{i} \cosh(x) \\ \tan(ix) &= \phantom{i} \tanh(x) \end{align*}$

These relationships are interesting because the circular and hyperbolic functions are independently meaningful. If you view these equations merely as definitions you lose their significance. Circular and hyperbolic functions were widely used before Euler discovered the connection between them.

Similarly, there’s a reason the modified Bessel functions were given a name their own. If you were led to Bessel functions and modified Bessel functions separately by different applications, you would regard the equation

$I_n(x) = i^{-n} J_n(ix)$

as a discovery rather than just a definition. The following section explains why someone would be interested in modified Bessel functions.

Before we move on, I’d like to explain the reason for the term i⁻ⁿ term. In general

$I_\nu(x) = \exp(\nu\pi i/2) J_n(ix)$

for all real ν. The reason for the exp(νπi/2) term is that it makes I_ν(x) real for all real x.

Applied math perspective

Bessel functions often arise from solving problems with radial symmetry. Solving the wave equation in cylindrical coordinates using separation of variables leads to Bessel’s differential equation

$x^2 y'' + x y' + (x^2 - \nu^2) y = 0$

and its solutions J_n and Y_n, Bessel functions of the first and second kind.

Solving the heat equation in cylindrical coordinates with separation of variables leads to the modified Bessel equation

$x^2 y^{\prime \prime} + x y^{\prime} - (x^2 + \nu^2) y = 0$

and its solutions I_n and K_n, the modified Bessel functions of the first and second kind.

This is the reason behind the complex analysis perspective above: the change of variables sending x to ix changes the sign of the x² term in Bessel’s equation.

Bessel functions describe radially symmetric oscillations, such as the vibrations of a drum head. Modified Bessel functions describe radially symmetric exponential growth or decay [1], such as in the heat in a cylinder.

Other modified functions

Struve functions are closely related to Bessel functions. The (modified) Struve functions also satisfy Bessel’s (modified) differential equation, but with a non-zero right hand side. The modified Struve functions are proportional to the unmodified Struve functions evaluated along the imaginary axis, with a proportionality constant that makes the modified Struve functions real for real arguments.

There’s a similar relationship between the Mathieu functions and modified Mathieu functions. The general pattern is that “modified” in the context of special functions means “evaluated at ix and multiplied by a constant to make the function real for real arguments.”

[1] The functions I_n grow exponentially and the functions K_n decay exponentially. For this reason, A&S didn’t tabulate I_n and K_n per se. Instead it tabulated e^−xI_n and e^xK_n because these functions varied less over their range.

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

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.

Why polynomial coefficients?

Posted on 1 August 2026 by John

Second order linear differential equations with polynomial coefficients form their own area of study. This seems like a narrow class of equations, but it’s very important in applications.

This class of equations seems like a mathematically natural topic, but why is it so important in applications? I did a PhD in differential equations without ever learning why. The theory of second order linear equations with polynomial coefficients is too complicated for undergraduate courses [0] and too well-established for graduate courses [1].

The explanation that I was missing can be found in the first chapter of [2]. The PDEs that are common in physics are separable in various coordinate systems, meaning that in these coordinate systems the PDEs reduce to ODEs. These ODEs either have polynomial coefficients, or there is a change of variables which makes the ODEs have polynomial coefficients.

See this writeup that looks at the Helmholtz and Laplace equations in 11 coordinate systems.

[0] You may see the simplest parts of the theory in a section on solving ODEs with power series. But textbooks don’t go very far for good reasons.

[1] Unfortunately, a lot of really useful topics are left out of the graduate curriculum because they’re too well understood to provide thesis topics. Or the problems that are still open have been open for so long that they’re likely too hard to be cracked by a graduate student.

[2] Gerhard Kristensson. Second Order Differential Equations: Special Functions and their Classification. Springer, 2010.

Counting rooted trees

Posted on 1 August 2026 by John

Combinatorial problems can be interesting for their own sake, but they are more interesting when there is a connection to a problem outside combinatorics, and the more unexpected the connection the better.

Counting the number of unlabeled rooted trees [1] with n nodes is a pure mathematics problem. Designing numerical methods for solving differential equations is an applied mathematics problem. And yet the two are closely linked.

Let t(n) be the number of distinct unlabeled rooted trees with n nodes. The diagram below shows that the first few terms of this sequence are 1, 1, 2, and 4.

Recursive calculation

The values of t(n) can be computed recursively using

$\begin{align*} g_k &= \sum_{d\mid k} d t_d \\ t_1 &= 1 \\ t_n &= \frac{1}{n-1} \sum_{k=1}^{n-1} g_k t_{n-k} \text{ for } n > 1<br /> \end{align*}<br />$

You can implement this in Python as follows.

from sympy import divisors

def t(n):
    if n <= 1:
        return 1 if n == 1 else 0
    return sum(g(k) * t(n - k) for k in range(1, n)) // (n - 1)

def g(k):
    return sum(d * t(d) for d in divisors(k))

This code is correct, but it will run more efficiently if you cache function values to avoid calculating the same values over and over. You can do this by adding

from functools import lru_cache

and writing @lru_cache(maxsize=None) above both function definitions.

Connection to Runge-Kutta

In an earlier post I showed that designing a 4-stage explicit Runge-Kutta method required solving a system of 8 equations in 10 unknowns, leaving two degrees of freedom in the solutions.

The number of constraints c(s) needed to design an s-stage explicit RK method is equal to the number of rooted trees with up to s nodes:

c(s) = t(1) + t(2) + t(3) + … + t(s)

This is because there is a one-to-one correspondence between constraints on the nth derivative of an RK formula and rooted trees, and an s stage method has to satisfy the constraints of all stages up to s. In the example of the 4th order RK method, we have

c(4) = t(1) + t(2) + t(3) + t(4) = 1 + 1 + 2 + 4 = 8.

The first few values [2] of t(n) are

1, 1, 2, 4, 9, 20, 48, 115, 286, 719, 1842, 4766, 12486, 32973, …

and so you can see that t(n) grows quickly. In fact, it grows exponentially [3].

However, the number of parameters in an s stage RK method is s(s + 1)/2. The number of equations grows exponentially and the number of variables grows only quadratically, so at some point you have more equations than variables. That’s already the case for s = 5 because you have 17 constraints on 15 variables. The system has a solution because symmetry considerations render some of the equations redundant.

A 10th order RK method requires 17 stages. (See the previous post for why the number of stages exceeds the order when the order is greater than 4.) Designing such a method would require solving over a million equations in 153 variables, and yet it can be done. [4]

[1] This is a slightly contradictory term. Unlabeled means the we don’t distinguish the nodes. But we do distinguish one node, namely the root.

[2] See OEIS A000081.

[2] Richard Otter proved in 1948 that the number of unlabeled rooted trees with n nodes is asymptotically C αⁿ / n^−3/2 where C = 0.4399… and α = 2.9557…. The cumulative sum is at least this large since Otter’s estimate gives the size of the last term in the sum.

[3] E. Hairer. A Runge-Kutta Method of Order 10. J. Inst. Maths Applics (1978) 21, 47-59

Runge-Kutta order versus stages

Posted on 1 August 2026 by John

The textbook version of the Runge-Kutta method for solving differential equations has 4 stages and has 4th order error. For lower order versions of RK the number of stages s also matches the order of the error p. But in order to achieve error on the order of p ≥ 5, you need more than p stages. This is known as the Butcher barrier.

Before going any further, let’s back up and say what we mean by stages and by order.

Stages

The number of stages in an RK method to solve the equation

y' = f(t, y)

is the number of evaluations of the function f on the right-hand side. For example, the textbook RK4 method estimates the solution at each step by

$y_{n+1} = y_n + \frac{h}{6}\left( k_{n1} + 2k_{n2} + 2k_{n3} + k_{n4}\right)$

where

$k_{n1} &=& f(t_n, y_n) \\ k_{n2} &=& f(t_n + 0.5h, y_n + 0.5hk_{n1}) \\ k_{n3} &=& f(t_n + 0.5h, y_n + 0.5hk_{n2}) \\ k_{n4} &=& f(t_n + h, y_n + hk_{n3}) \\$

which requires four stages, i.e. four evaluations of f.

Order

A differential equation solver is said to have order p if the local error, the error after one step of size h, is O(h^{p + 1}). Then after solving an ODE over a period of time T with N = T/h steps, the global error is O(h^p). So, for example, if p = 4, you would expect that cutting your step size h in half would cut your error at T by a factor of 16.

More stages than the order

John C. Butcher proved that an explicit RK method of order p requires s stages where s > p if p > 4.

An important example is the Dormand-Prince method. It is a version of RK that has order 5 and 7 stages. The clever thing about this method is that you can make a 4th order solver out of a subset of its function evaluations.

That means that after you’ve evaluated one step of the 5th order method, you can also evaluate a 4th order method essentially for free. And by comparing them, you can get a sense of the error. If the solutions given by the two methods are substantially different, you have probably taken too big a step and need to back up. If the two solutions essentially agree, you’re probably good to take the next step.

For an explict RK method to have order 5, 6, or 7 you need at least 6, 7, or 9 stages respectively.

Solving the RK4 design equations

Posted on 31 July 2026 by John

I was digging into the Runge-Kutta method for solving differential equations and a line from [1] piqued my curiosity.

These calculations, which are not reproduced in Kutta’s paper (they are however in Huen (1900)), are very tedious.

The calculations are a set of eight constraints that the parameters of a fourth-order Runge-Kutta method must satisfy. I wondered how well Mathematica might have done at assisting Mr. Huen in his “very tedious” calculations if it had been available in 1900.

I go into Runge-Kutta methods in this post. Here I’d like to concentrate on a step in the design of the methods, namely solving the set of equations alluded in the quote above.

$\begin{align*} b_1 + b_2 + b_3 + b_4 &= 1 \\ b_2 c_2 + b_3 c_3 + b_4 c_4 &= \frac{1}{2} \\ b_2 c_2^2 + b_3 c_3^2 + b_4 c_4^2 &= \frac{1}{3} \\ b_3 a_{32} c_2 + b_4(a_{42} c_2 + a_{43} c_3) &= \frac{1}{6} \\ b_2 c_2^3 + b_3 c_3^3 + b_4 c_4^3 &= \frac{1}{4} \\ b_3 c_3 a_{32} c_2 + b_4 c_4(a_{42} c_2 + a_{43} c_3) &= \frac{1}{8} \\ b_3 a_{32} c_2^2 + b_4(a_{42} c_2^2 + a_{43} c_3^2) &= \frac{1}{12} \\ b_4 a_{43} a_{32} c_2 &= \frac{1}{24} \end{align*}$

The first thing to note is that there are 10 variables and only 8 equations, and so the solution is not fully determined. What we think of as the fourth order Runge-Kutta method is in fact a fourth order Runge-Kutta method.

One could argue that we should have b₂ = b₃ and c₂ = c₃. With these additional equations, the system of equations has a unique solution, and Mathematic finds it easily.

eqs = {
    b1 + b2 + b3 + b4 == 1,
    b2*c2 + b3*c3 + b4*c4 == 1/2,
    b2*c2^2 + b3*c3^2 + b4*c4^2 == 1/3,
    b3*a32*c2 + b4*(a42*c2 + a43*c3) == 1/6,
    b2*c2^3 + b3*c3^3 + b4*c4^3 == 1/4,
    b3*c3*a32*c2 + b4*c4*(a42*c2 + a43*c3) == 1/8,
    b3*a32*c2^2 + b4*(a42*c2^2 + a43*c3^2) == 1/12,
    b4*a43*a32*c2 == 1/24,
    b2 == b3,
    c2 == c3
};

vars = {b1, b2, b3, b4, c2, c3, c4, a32, a42, a43};

solution = Solve[eqs, vars]

This returns the parameters used for the version of Runge-Kutta presented in every textbook.

If you keep the requirement b₂ = b₃ but substitute the requirement 2c₂ = c₃ for c‘s Mathematica will return the coefficients for the so-called Runge-Kutta 3/8 rule. This method has some slight advantages by some criteria.

In 1951 Gill [2] discovered a fourth order Runge-Kutta rule optimized for running in extremely constrained computer hardware. It’s a strange method, with irrational parameters, but one that was a very clever response to the limitations of its time.

Update: See this post for a discussion of the parameters and constriants for higher-ordered RK methods.

[1] Hairer, Nørsett, and Wanner. Solving Ordinary Differential Equations I: Nonstiff Problems. Springer-Verlag 1987.

[2] A. Gill. A process for the step-by-step integration of differential equations in an automatic digital computing machine. Proc. Cambridge Philos. Soc., vol 27, pp 95–108.

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.

Derivative equals inverse

Posted on 29 June 2026 by John

Here’s kind of a strange problem with an interesting solution: find a function f such that the derivative of f equals the inverse of f for all positive x.

f ′(x) = f⁻¹(x)

This is a differential equation, but a very unusual one, one that cannot be solved using any of the techniques taught in a class on differential equations.

The unique solution is

f(x) = φ(x / φ)^φ

where φ is the golden ratio. What an unexpected appearance of the golden ratio!

The problem was proposed by H. L. Nelson and solved by A. C. Hindmarsh. See The American Mathematical Monthly, Vol. 76, No. 6 p. 696.

Real and imaginary parts

Posted on 23 May 2026 by John

The previous post announced some notes I wrote up based on an article by Henry Baker implementing functions of a complex variable in terms of functions of a real variable. That is, it finds functions u(x, y) and v(x, y) such that

f(x + iy) = u(x, y) + i v(x, y)

where x, y, u, and v are all real-valued. Not only that, but if f is an elementary function, so are u and v. Here “elementary” has a technical meaning, but essentially it means functions that you could evaluate on a scientific calculator. A couple of the functions might be unfamiliar, such as sgn and atan2, but there are no functions like the gamma function that are defined in terms of integrals.

One application of Baker’s equations would be to bootstrap a math library that doesn’t support complex numbers into one that does. But the equations could be useful in pure math when you’d like to have a convenient expression for the real or imaginary part of a function.

The real and imaginary parts of a complex analytic function are harmonic functions. So the functions on the right hand side of Baker’s equations satisfy Laplace’s equation:

u_xx + u_yy = 0

and

v_xx + v_yy = 0.

Furthermore, the functions u and v form harmonic conjugate pairs, meaning each is the Hilbert transform of the other.

Differential equations

What exactly is modified about a modified Bessel function?

Pure math perspective

Applied math perspective

Other modified functions

Inverse differential equations

Holonomic functions

Why polynomial coefficients?

Counting rooted trees

Recursive calculation

Connection to Runge-Kutta

Related posts

Runge-Kutta order versus stages

Stages

Order

More stages than the order

Solving the RK4 design equations

Related posts

exp_q

Closed forms

Differential equations

Mathematica and Mittag-Leffler

Combinatorics

Derivative equals inverse

Real and imaginary parts