It’s not obvious when a nonlinear differential equation will have a periodic solution, or asymptotically approach a periodic solution. But there are theorems that give sufficient conditions. This post is about one such theorem.

Consider the equation

x” + f(x) x‘ + g(x) = 0

where f and g are analytic and define F(x) to be the integral of f(t) from 0 to x. The following six hypotheses are sufficient to guarantee that the equation above has a stable limit cycle [1].

1. g is an odd function,
2. g is positive for positive x,
3. f is an even function,
4. f(o) < 0,
5. F(x) → ∞ as x → ∞,
6. F has a unique positive zero.

These hypotheses are satisfied, for example, by Van der Pol’s equation.

Let’s make up coefficient functions f and g that satisfy the conditions above and look for the limit cycle.

We can let g(x) = x³ and f(x) = x² − 1. And here’s what we get:

The graph shows phase portraits for two different initial conditions, given in the legend of the plot. Note that limit cycle appears solid blue because the trajectories of both the dashed line and the dotted line run around the periodic attractor.

Related posts

• Calculating Van der Pol periods
• A stiffening spring
• Spring, rubber band, chaos

[1] You can find this theorem in Topics in Ordinary Differential Equations by William Lakin and David Sanchez.

Suppose we have a differential equation of the form

x” + b(t) x‘ + c(t) x = f(t)

where the functions b and c are periodic with the same period T.

When does the equation above have a solution with period T?

It turns out [1] the equation above will have a T-periodic solution for all T-periodic functions f if and only if the equation

x” + b(t) x‘ + c(t) x = 0

does not have a T-periodic solution. In other words, the non-homogeneous equation has a T-periodic solutions for all f if and only if the homogeneous equation does not [2].

This theorem may not sound right at first. You might think the homogeneous equation having a periodic solution would help if you were looking for a solution to the non-homogeneous equation. An example will show this isn’t the case.

Example

Let’s look at the special case of b = 0, c = a positive constant, and f(t) = cos ωt. As shown here, the natural frequency of the equation is

ω₀ = √c

and the homogeneous equation has general solution

x(t) = R cos(ω₀t – φ)

for some amplitude R and phase φ.

Also, the non-homogeneous equation has a solution

cos(ωt) /(ω₀² – ω²).

So does this contradict the theorem above? Here both the homogeneous and non-homogeneous equations have periodic solutions. However, they don’t have periodic solutions with the same period. The homogeneous solutions have period 2π/ω₀, and the non-homogeneous solution has period 2π/ω.

The non-homogeneous solution above requires the driving frequency ω to be different from the natural frequency ω₀ or else the denominator would be 0.

If the driving frequency does equal the natural frequency, the non-homogeneous equation does have a solution, but the solution is not periodic. The solutions blows up due to resonance. The non-homogeneous equation has a solution

(t/2ω) sin ωt.

and so the amplitudes of the oscillations grow as t grows. With each period, more energy is pumped into the system and the energy increases without bound.

More general periods

The full set of solutions to a non-homogeneous equation is given by the general solution to the homogeneous equation plus any particular solution to the non-homogeneous equation.

If the driving frequency ω does not equal the natural frequency ω₀ then in the example above the homogeneous and non-homogeneous solutions are both periodic, but with different periods. Is their sum periodic? Maybe.

The sum of two periodic functions may or may not be periodic. The image below is from this post where I give an example of two periodic functions, the green and blue lines, whose sum, the red line, is not periodic but rather quasiperiodic.

The sum of two periodic functions is periodic if the ratio of their periods is a rational number. In the example above, the ratio is the golden ratio, an irrational number.

The general solutions to the non-homogeneous equation above are periodic if and only if the ratio of ω to ω₀ is rational.

Contradiction?

But if ω/ω₀ is rational, do we have a contradiction to our theorem above? For example, suppose the natural frequency ω₀ is 2π and the driving frequency ω is 4π. The the homogeneous solutions have period 1 and the particular solution to the non-homogeneous equation has period 1/2. But a function of period 1/2 also has period 1.

In our case, the equation

x” + 4π² x = 0

has solutions of period 1, and the equation

x” + 4π² x = cos 4πt

also has a solution with period 1.

To see why this isn’t a contradiction, we have to go back and read the theorem more carefully. It says the equation

x” + b(t) x‘ + c(t) x = f(t)

has a T-periodic solution for all T-periodic functions f if and only if the equation

x” + b(t) x‘ + c(t) x = 0

does not have a T-periodic solution.

In our case the homogeneous equation has solutions with period 1, and the non-homogeneous equation also has a solution with period 1. But this is just for one particular forcing function f. The theorem makes a claim for all functions f, and so we have not found a contradiction.

As long as there is some function f(t) with period 1 for which the non-homogeneous equation has no solution with period 1 then the theorem is correct. And in fact we’ve already given an example of such a function, namely cos 2πt. In that case the natural frequency and driving frequency are the same, and the solutions are not periodic due to resonance.

Related posts

• Oscillatory solutions to differential equations
• Predator-Prey period
• NASA’s favorite ODE solver

[1] See Topics in Ordinary Differential Equations by William Lakin and David Sanchez.

[2] A differential equation in the form above is called homogeneous if the forcing function on the right side is 0 and it is called non-homogeneous if the forcing function is not zero.

[3] If the driving frequency ω does not equal the natural frequency ω₀ is the solution periodic even if it doesn’t have the same period

The Lotka-Volterra equations are a system of nonlinear differential equations for modeling a predator-prey ecosystem. After a suitable change of units the equations can be written in the form

where ab = 1. Here x(t) is the population of prey at time t and y(t) is the population of predators. For example, maybe x represents rabbits and y represents foxes, or x represents Eloi and y represents Morlocks.

It is well known that the Lotka-Volterra equations have periodic solutions. It is not as well known that you can compute the period of a solution without having to first solve the system of equations.

This post will show how to compute the period of the system. First we’ll find the period by solving the equations for a few different initial conditions, then we’ll show how to directly compute the period from the system parameters.

Phase plot

Here is a plot of (x(t), y(t)) showing that the solutions are periodic.

And here’s the Python code that made the plot above.

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.integrate import solve_ivp
    
    # Lotka-Volterra equations
    def lv(t, z, a, b):
        x, y = z
        return [a*x*(y-1), -b*y*(x-1)]
    
    begin, end = 0, 20
    t = np.linspace(begin, end, 300)
    
    styles = ["-", ":", "--"]
    prey_init = [2, 4, 8]
    a = 1.5
    b = 1/a
    
    for style, p0 in zip(styles, prey_init):
        sol = solve_ivp(lv, [begin, end], [p0, 1], t_eval=t, args=(a, b))
        plt.plot(sol.y[0], sol.y[1], style)
    
    plt.xlabel("prey")
    plt.ylabel("predator")
    plt.legend([f"$p_0$ = {p0}" for p0 in prey_init])    

Note that the derivative of x is zero when y = 1. Since our initial condition sets y(0) = 1, we’re specifying the maximum value of x. (If our initial values of x were less than 1, we’d be specifying the minimum of x.)

Time plot

The components x and y have the same period, and that period depends on a and b, and on the initial conditions. The plot below shows how the period increases with x(0), which as we noted above is the maximum value of x.

And here’s the code that made the plot.

    fig, ax = plt.subplots(3, 1)
    for i in range(3):
        sol = solve_ivp(lv, [begin, end], [prey_init[i], 1], t_eval=t, args=(a, b))
        ax[i].plot(t, sol.y[0], "-")
        ax[i].plot(t, sol.y[1], "--")    
        ax[i].set_xlabel("$t$")
        ax[i].set_ylabel("population")
        ax[i].set_title(f"x(0) = {prey_init[i]}")
    plt.tight_layout()

Finding the period

In [1] the author develops a way to compute the period of the system as a function of its parameters without solving the differential equations.

First, calculate an invariant h of the system:

Since this is an invariant we can evaluate it anywhere, so we evaluate it at the initial conditions.

Then the period only depends on a and h. (Recall we said we can scale the equations so that ab = 1, so a and b are not independent parameters.)

If h is not too large, we can compute the approximate period using an asymptotic series.

where σ = (a + b)/2 = (a² + 1)/2a.

We use this to find the periods for the example above.

    def find_h(a, x0, y0):
        b = 1/a
        return b*(x0 - np.log(x0) - 1) + a*(y0 - np.log(y0) - 1)

    def P(h, a):
        sigma = 0.5*(a + 1/a)
        s = 1 + sigma*h/6 + sigma**2*h**2/144
        return 2*np.pi*s

    print([P(find_h(1.5, p0, 1), 1.5) for p0 in [2, 4, 8]])

This predicts periods of roughly 6.5, 7.5, and 10.5, which is what we see in the plot above.

When h is larger, the period can be calculated by numerically evaluating an integral given in [1].

Related posts

• Finding the period of Van der Pol oscillations
• The Brusselator
• Lorenz system

[1] Jörg Waldvogel. The Period in the Volterra-Lotka Predator-Prey Model. SIAM Journal on Numerical Analysis, Dec., 1983, Vol. 20, No. 6, pp. 1264-1272

The solution to a differential equation is called oscillatory if its set of zeros is unbounded. This does not necessarily mean that the solution is periodic, but that it crosses the horizontal axis infinitely often.

Fowler [1] studied the following differential equation demonstrates both oscillatory and nonoscillatory behavior.

x” + t^σ |x|^γ sgn x = 0

He proved that

1. If σ + 2 ≥ 0 then all solutions oscillate.
2. If σ + (γ + 3)/2 < 0 then no solutions oscillate.
3. If neither (1) nor (2) holds then some solutions oscillate and some do not.

The edge case is σ = -2 and γ = 1. This satisfies condition (1) above. A slight decrease in σ will push the equation into condition (2). And a simultaneous decrease in σ and increase in γ can put it in condition (3).

It turns out that the edge case can be solved in closed form. The equation is linear because γ = 1 and solutions are given by

c √t sin(φ + √3 log(t) / 2).

The solutions oscillate because the argument to sine is an unbounded function of t, and so it runs across multiples of π infinitely often. The distance between zero crossings increases exponentially with time because the logarithms of the crossings are evenly spaced.

This post took a different turn after I started writing it. My original intention was to solve the differential equation numerically and say “See: when you change the parameters slightly you get what the theorem says.” But that didn’t work out.

First I tried numerically computing the solution above with σ = −2 and γ = 1, but I didn’t see oscillations. I tried making σ a little larger, but still didn’t see oscillations. When I increased σ more I could see the oscillations.

I was able to go back and see the oscillations with σ = −2 and γ = 1 by tweaking the parameters in my ODE solver. It’s best practice to tweak your parameters even if everything looks right, just to make sure your solution is robust to algorithm changes. Maybe I would have done that, but since I already knew the analytic solution, I knew to tweak the parameters until I saw the right behavior.

Without some theory as a guide, it would be difficult to determine from numerical solutions alone whether a solution oscillates. If you see oscillations, then they’re probably real (unless your numerical method is unstable), but if you don’t see oscillations, how do you know you won’t see them if you look further out? How much further out should you look? In a practical application, the context of your problem will tell you how far out to look.

My intention was to recommend the equation above as an illustration for an introductory DE class because it demonstrates the important fact that small changes to an equation can qualitatively change the behavior of the solutions. This is especially true for nonlinear DEs, which the above equation is when γ ≠ 1. It could still be used for that purpose, but it would make a better demonstration than homework exercise. The instructor could tune the DE solver to make sure the solutions are qualitatively correct.

I recommend the equation as an illustration, but for a course covering numerical solutions to DEs. It illustrates a different point than I first intended, namely the potentially finicky behavior of software for solving differential equations.

There are a couple morals to this story. One is that numerical methods have not eliminated the need for theory. The ideal is to combine analytic and numerical methods. Analytic methods point out qualitative behavior that might be hard to discover numerically. And analytic methods are not likely to produce closed-form solutions for nonlinear DEs.

Another moral is that it’s best to twiddle the parameters to your numerical method to make sure the solution doesn’t qualitatively change. If you’ve adequately computed a solution, computing it again with smaller integration steps shouldn’t make much difference. But if you do see a substantial difference, your first solution probably wasn’t as good as you thought.

Related posts

• A stiffening spring
• Driving vibrations with a sawtooth wave
• Symplectic Euler method

[1] R. H. Fowler. Further studies of Emden’s and similar differential equations. Quarterly Journal of Mathematics. 2 (1931), pp. 259–288.

This post will tie together many things I’ve blogged about before. The previous post justified separation of variables. This post will illustrate separation of variables.

Also, this post will show why you might care about Bessel functions and their zeros. I’ve written about Bessel functions before, and said that Bessel functions are to polar coordinates what sines and cosines are to rectangular coordinates. This post will make this analogy more concrete.

Separation of variables

The separation of variables technique is typically presented in three contexts in introductory courses on differential equations:

1. One-dimensional heat equation
2. One-dimensional wave equation
3. Two-dimensional (rectangular) Laplace equation

My initial motivation for writing this post was to illustrate separation of variables outside the most common examples. Separation of variables requires PDEs to have a special form, but not as special as the examples above might imply. A secondary motivation was to show Bessel functions in action.

Radially symmetric wave equation

Suppose you have a thin membrane, like a drum head, and you want to model its vibrations. By “thin” I mean that the membrane is sufficiently thin that we can adequately model it as a two-dimensional surface bobbing up and down in three dimensions. It’s not so thick that we need to model the material in more detail.

Let u(x, y, t) be the height of the membrane at location (x, y) and time t. The wave equation is a PDE modeling the motion of the membrane by

where Δ is the Laplacian operator. In rectangular coordinates the Laplacian is given by

We’re interested in a circular membrane, and so things will be much easier if we work in polar coordinates. In polar coordinates the Laplacian is given by

We will assume that our boundary conditions and initial conditions are radially symmetric, and so our solution will be radially symmetric, i.e. derivatives with respect to θ are zero. And in our case the wave equation simplifies to

Boundary and initial conditions

Let a be the radius of our membrane. We will assume our membrane is clamped down on its boundary, like a drum head, which means we have the boundary condition

for all t ≥ 0. We assume the initial displacement and initial velocity of the membrane are given by

for all r between 0 and a.

Separating variables

Now we get down to separation of variables. Because we’re assuming radial symmetry, we’re down to a function of two variables: r for the distance from the center and t for time. We assume

and stick it into the wave equation. A little calculation shows that

The left side is a function of t alone, and the right side is a function of r alone. The only way this can be true for all t and all r is for both sides to be constant. Call this constant -λ². Why? Because looking ahead a little we find that this will make things easier shortly.

Separation of variables allowed us to reduce our PDE to the following pair of ODEs.

The solutions to the equation for R are linear combinations of the Bessel functions J₀(λr) and Y₀(λr) [1].

And the solutions to the equation for T are linear combinations of the trig functions cos(cλt) and sin(cλt).

The boundary condition u(a, t) = 0 implies the boundary condition R(a) = 0. This implies that λa must be a zero of our Bessel function, and that all the Y₀ terms drop out. This means that our solution is

where

Here α_n are the zeros of the Bessel function J₀.

The coefficients A_n and B_n are determined by the initial conditions. Specifically, you can show that

The function J₁ in the expression for the coefficients is another Bessel function.

The functions J₀ and J₁ are so important in applications that even the otherwise minimalist Unix calculator bc includes these functions. (As much as I appreciate Bessel functions, this still seems strange to me.) And you can find functions for zeros of Bessel functions in many libraries, such as scipy.special.jn_zeros in Python.

Related posts

• The acoustics of kettle drums
• bc surprises
• Fourier-Bessel series

[1] This is why introductory courses are unlikely to include an example in polar coordinates. Separation of variables itself is no harder in polar coordinates than in rectangular coordinates, and it shows the versatility of the method to apply it in a different setting.

But the resulting ODEs have Bessel functions for solutions, and it’s understandable that an introductory course might not want to go down this rabbit trail, especially since PDEs are usually introduced at the end of a differential equation class when professors are rushed and students are tired.

[2] Snare drum image above “Snare drum” by YannickWhee is licensed under CC BY 2.0

The separation of variables technique for solving partial differential equations looks like a magic trick the first time you see it. The lecturer, or author if you’re more self-taught, makes an audacious assumption, like pulling a rabbit out of a hat, and it works.

For example, you might first see the heat equation

u_t = c² u_xx.

The professor asks you to assume the solution has the form

u(x, t) = X(x) T(t).

i.e. the solution can be separated into the product of a function of x alone and a function of t alone.

Following that you might see Laplace’s equation on a rectangle

u_xx + u_yy = 0

with the analogous assumption that

u(x, y) = X(x) Y(y),

i.e. the product of a function of x alone and a function of y alone.

There are several possible responses to this assumption.

1. Whatever you say, doc.
2. How can you assume that?
3. How do you know you’re not missing any possibilities?
4. What made someone think to try this?

As with many things, separation of variables causes the most consternation for the moderately sophisticated students. The least sophisticated students are untroubled, and the most sophisticated student can supply their own justification (at least after the fact).

One response to question (2) is “Bear with me. I’ll show that this works.”

Another response would be “OK, how about assuming the solution is a sum of such functions. That’s a much larger space to look in. And besides, we are going to take sums of such solutions in a few minutes.” One could argue from functional analysis or approximation theory that the sums of separable functions are dense in reasonable space of functions [1].

This is a solid explanation, but it’s kind of anachronistic: most students see separation of variables long before they see functional analysis or approximation theory. But it would be a satisfying response for someone who is seeing all this for the second time. Maybe they were exposed to separation of variables as an undergraduate and now they’re taking a graduate course in PDEs. In an undergraduate class a professor could do a little foreshadowing, giving the students a taste of approximation theory.

Existence of solutions is easier to prove than uniqueness in this case because you can concretely construct a solution. This goes back to the “it works” justification. This argument deserves more respect than a sophomoric student might give it. Mathematics research is not nearly as deductive and mathematics education. You often have to make inspired guesses and then show that they work.

Addressing question (3) requires saying something about uniqueness. A professor could simply assert that there are uniqueness theorems that allow you to go from “I’ve found something that works” to “and so it must be the only thing that works.” Or one could sketch a uniqueness theorem. For example, you might apply a maximum principle to show that the difference between any two solutions is zero.

Question (4) is in some sense the most interesting question. It’s not a mathematical question per se but a question about how people do mathematics. I don’t know what was going through the mind of the first person to try separation of variables, or even who this person was. But a plausible line of thinking is that ordinary differential equations are easier than partial differential equations. How might you reduce a PDE to an ODE? Well, if the solution could be factored into functions of one variable, …

The next post will illustrate using separation of variables by solving the wave equation on a disk.

Related posts

• Hilbert space methods for PDEs
• Structure in jazz and math
• Boundary conditions are the hard part

[1] Also, there’s the mind-blowing Kolmogorov-Arnol’d theorem. This theorem says any continuous function of several variables can be written as a sum of continuous separable functions. It doesn’t say you can make the functions in your sum smooth, but it suggests that sums of separable functions are more expressive than you might have imagined.