Differential equations

Navier-Stokes in the news

Posted on by John

There are rumors that a long-standing math problem, one of the Millennium Prize problems, has been solved.

The problem concerns technical properties of solutions to the Navier-Stokes equations [1], a set of equations that describe the dynamics of fluid flow. Popular accounts of the problem are often oversimplified and misleading.

Some reports will speak of the problem as “solving the Navier-Stokes equations.” The task is not to write down a closed-form solution, which can’t be done, or solve the equations numerically, which has been done for decades. The problem is to prove theoretical properties of solutions which are of little interest in practice.

There has been progress toward settling the Navier-Stokes problem. Terence Tao wrote a post on this yesterday.

What’s also  interesting is the intrigue around the possible solution. A post this morning says

If I am reading this correctly, Tristan Buckmaster is alleging OAI has a resolution of Navier-Stokes … which maybe used info from Buckmaster and Levent Alpöge’s private Codex sessions.

Buckmaster asked OpenAI whether they used his private sessions and they have not responded. Update: Statement from OpenAI.

Personal note

This topic connects parts of my career spanning decades. My graduate work was in PDEs and I had some interest in the Navier-Stokes equations. Here are some notes I wrote back in the day.

Now I work more with privacy than with PDEs. The question of whether OpenAI uses private data, contradicting their stated policy, is more relevant to my current work than whether the Navier-Stokes equations have global regular solutions.

Related posts

[1] I never know whether to say equation or equations. You’ll hear both. You could think of Navier-Stokes as one vector-valued PDE or three scalar-valued equations.

Second solutions

Posted on by John

This post provides a couple examples to go along with two earlier posts.

The pattern we’re illustrating is families of polynomials pn(x) that each satisfy a differential equation and a three-term recurrence. The differential equations have a second solution qn(x) that is the larger solution with respect to x but the smaller solution with respect to n.

In both the examples below pn(x) is a polynomial, and so bounded on the interval [−1, 1], and qn(x) is not a polynomial, with singularities at ±1. This is analogous to the previous examples with Bessel functions Jn(x) and Qn(x) that satisfy the same differential equation but have contrasting behavior with respect to x versus n.

Legendre polynomials

The differential equation

(1-x^2)\,y^{\prime\prime} - 2x\,y^\prime + n(n+1)\,y = 0

has two solutions for each n, Pn(x) and Qn(x).

The solutions Pn(x) are the Legendre polynomials. The solutions Qn(x) are not polynomials but involve a term log((1 + x)/(1 − x)) that blows up at 1 and −1. But for fixed x and increasing n, Pn(x) grows exponentially and Qn(x) decays exponentially, provided |x| > 1.

Chebyshev polynomials

The differential equation

(1-x^2)\,y^{\prime\prime} - x\,y^\prime + n^2\,y = 0

has two solutions for each n, Tn(x) and Vn(x).

The solutions Tn(x) are the Chebyshev polynomials. The solutions Vn(x) are not polynomials but involve a term √(x² — 1) that become vertical up at 1 and −1. But for fixed x with |x| > 1 and increasing n, Tn(x) grows exponentially and Vn(x) decays exponentially.

Junk solutions

Posted on by John

When you’re interested in studying a family of functions, it can be useful to look at a differential equation that the functions solve. This is a theme I’ve written about several times, most recently here and here, but also three years ago here.

Orthogonal polynomials are mathematically elegant as well as very useful in applications [1]. Various families of orthogonal polynomials satisfy various differential equations. These equations have a polynomial and non-polynomial solutions. What use are the latter?

If the differential equation modeled something physical, then the second solution would be necessary to have a complete basis of solutions. But if the differential equation is only instrumental in studying the orthogonal polynomials, what use is a non-polynomial solution?

These non-polynomial solutions turn out to be useful. Just as “junk” DNA turned out not to be junk, these “junk” solutions are important. Junk DNA doesn’t directly code for proteins, but it regulates DNA that does code for proteins and serves other purposes. Similarly, these non-polynomial solutions carry information related to the polynomial solutions.

For example, orthogonal polynomials are used to construct numerical integration methods, such as Gaussian quadrature, and the associated non-polynomial solutions describe the error in these integration methods. Incidentally, Gaussian quadrature is based on Legendre polynomials, mentioned in the previous post. For every family of orthogonal polynomials there is a corresponding integration method. See these notes.

Another tie-in to recent posts is that these non-polynomial solutions are the minimal solution to the polynomial family’s three-term recurrence, the solution that takes extra care to compute numerically.

This post has been very high-level, alluding to ideas without going into details. I’d like to write future posts that go into more depth regarding the ideas introduced here.

 

[1] “Real analysts cannot do without Fourier, complex analysts cannot do without Laurent, and numerical analysts cannot do without Chebyshev [polynomials].” — Lloyd N. Trefethen”

Ultraspherical

Posted on by John

When I hear the term ultraspherical I think of something extremely spherical. For example, a baseball is spherical, but a billiard ball is more spherical. Maybe a highly polished billiard ball is ultraspherical.

Using this line of thought, the term ultraspherical polynomial is inexplicable. This is an example of the arcane terminology I wrote about recently. In this post I’ll explain what it conveys.

A spherical polynomial is a polynomial that naturally falls out of solving Laplace’s equation in spherical coordinates, using separation of variables. Legendre polynomials are spherical polynomials.

Gegenbauer polynomials are so called because a man named Gegenbauer studied them, just as Legendre polynomials take their name from Legendre. Gegenbauer polynomials are also called ultraspherical polynomials. Why is that?

There are two possible reasons. I’m not sure which is the historical reason, but both are plausible and are useful mnemonics.

Ultraspherical polynomials are not extremely spherical, they’re beyond spherical in some sense. More modern terminology uses the hyper- prefix rather than ultra-, which helps a bit.

Ultraspherical polynomials are beyond spherical in two ways. Gegenbauer polynomials are a generalization of Legendre polynomials, so they’re beyond Legendre polynomials in this sense.

More importantly, Gegenbauer polynomials fall out of solving Laplace’s equation on a hypersphere, i.e. a sphere in ℝn for n > 3, just as Legendre polynomials fall out of the case n = 3. It makes sense to call these polynomials hyperspherical because they fall out of solving an equation on a hypersphere. Unfortunately the classical term is ultraspherical rather than hyperspherical.

I think, but I’m not sure, that at one time higher dimensional spheres were called hyperspheres, but the the higher dimensional analog of spherical coordinates was called ultraspherical coordinates. If so, it would be understandable that the adjective modifying coordinates would be applied to the polynomials that result from solving equations in these coordinates.

What exactly is modified about a modified Bessel function?

Posted on 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 In is essentially the Bessel function Jn 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−n 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 Jn and Yn, 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 In and Kn, 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 In grow exponentially and the functions Kn decay exponentially. For this reason, A&S didn’t tabulate In and Kn per se. Instead it tabulated e−xIn and exKn because these functions varied less over their range.

Inverse differential equations

Posted on 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 Jn and Yn 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 Frobenius.”

Holonomic functions

Posted on 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 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 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]

Related posts

[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 / 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 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(hp + 1). Then after solving an ODE over a period of time T with N = T/h steps, the global error is O(hp). 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.