New expansions of confluent hypergeometric function

Last week Bujanda et al published a paper [1] that gives new expansions for the confluent hypergeometric function. I’ll back up explain what that means before saying more about the new paper.

Hypergeometric functions

Hypergeometric functions are something of a “grand unified theory” of special functions. Many functions that come up in application are special cases of hypergeometric function, as shown in the diagram below. (Larger version here.)

special function relationships

I give a brief introduction to hypergeometric functions in these notes (4-page pdf).

Confluent hypergeometric functions

The confluent hypergeometric function corresponds to Hypergeometric 1F1 in the chart above [2]. In Mathematica it goes by Hypergeometric1F1 and in Python goes by scipy.special.hyp1f1.

This function is important in probability because you can use it to compute the CDFs for the normal and gamma distributions.

New expansions

Bujanda et al give series expansions of the confluent hypergeometric functions in terms of incomplete gamma functions. That may not sound like progress because the incomplete gamma functions are still special functions, but the confluent hypergeometric function is a function of three variables M(ab; z) whereas the incomplete gamma function γ(sz) is a function of two variables.

They also give expansions in terms of elementary functions which may be of more interest. The authors give the series expansion below.

M_n(a, b; z) = \frac{\Gamma(b)}{\Gamma(a)\,\Gamma(b-a)} \sum_{k=0}^{n-1} A_k(a, b)\, F_k(z)

The A functions are given recursively by

and

A_n(a, b) = \frac{2}{n} \left( (2a - b)A_{n-1}(a, b) + 2(n-b)A_{n-2}(a,b) \right) 

for n > 1.

The F functions are given by

F_n(z) = \frac{n!}{(-z)^{n+1}} \Big( e_n(z/2) - \exp(z) e_n(-z/2) \Big)

where

e_n(z) = \sum_{k=0}^n \frac{z^k}{k!}

The authors give approximation error bounds in [1]. In the plot below we show that n = 3 makes a good approximation for M(3, 4.1, z). The error converges to zero uniformly as n increases.

Related posts

[1] Blanca Bujanda, José L. López, and Pedro J. Pagola. Convergence expansions of the confluent hypergeometric functions in terms of elementary functions. Mathematics of computation. DOI https://doi.org/10.1090/mcom/3389

[2] “Confluent” is a slightly archaic name, predating a more systematic naming for hypergeometric functions. The name mFn means the hypergeometric function has m upper parameters and n lower parameters. The confluent hypergeometric function has one of each, so it corresponds to 1F1. For more details, see these notes.

Gamma gamma gamma!

There are several things in math and statistics named gamma. Three examples are

  • the gamma function
  • the gamma constant
  • the gamma distribution

This post will show how these are related. We’ll also look at the incomplete gamma function which connects with all the above.

The gamma function

The gamma function is the most important function not usually found on a calculator. It’s the first “advanced” function you’re likely to learn about. You might see it in passing in a calculus class, in a homework problem on integration by parts, but usually not there’s not much emphasis on it. But it comes up a lot in application.

\Gamma(x) = \int_0^\infty t^{x-1 } e^{-t}\, dt

You can think of the gamma function as a way to extend factorial to non-integer values. For non-negative integers n, Γ(n + 1) = n!.

(Why is the argument n + 1 to Γ rather than n? There are a number of reasons, historical and practical. Short answer: some formulas turn out simpler if we define Γ that way.)

The gamma constant

The gamma constant, a.k.a. Euler’s constant or the Euler-Mascheroni constant, is defined as the asymptotic difference between harmonic numbers and logarithms. That is,

\gamma = \lim_{n\to\infty} \left( 1 + \frac{1}{2} + \frac{1}{3} + \cdots + \frac{1}{n}- \log n\right)

The constant γ comes up fairly often in applications. But what does it have to do with the gamma function? There’s a reason the constant and the function are both named by the same Greek letter. One is that the gamma constant is part of the product formula for the gamma function.

\frac{1}{\Gamma(x)} = x e^{\gamma x} \prod_{r=1}^\infty \left(1 + \frac{x}{r} \right) e^{-x/r}

If we take the logarithm of this formula and differentiation we find out that

\frac{\Gamma'(1)}{\Gamma(1)} = -\gamma

The gamma distribution

If you take the integrand defining the gamma function and turn it into a probability distribution by normalizing it to integrate to 1, you get the gamma distribution. That is, a gamma random variable with shape k has probability density function (PDF) given by

f(x) = \frac{1}{\Gamma(k)} x^{k-1} e^{-x}

More generally you could add a scaling parameter to the gamma distribution in the usual way. You could imaging the scaling parameter present here but set to 1 to make things simpler.

The incomplete gamma function

The incomplete gamma function relates to everything above. It’s like the (complete) gamma function, except the range of integration is finite. So it’s now a function of two variables, the extra variable being the limit of integration.

\gamma(s, x)= \int_0^x t^{s-1} e^{-t} \, dt

(Note that now x appears in the limit of integration, not the exponent of t. This notation is inconsistent with the definition of the (complete) gamma function but it’s conventional.)

It uses a lower case gamma for its notation, like the gamma constant, and is a generalization of the gamma function. It’s also essentially the cumulative distribution function of the gamma distribution. That is, the CDF of a gamma random variable with shape s is γ(sx) / Γ(s).

The function γ(sx) / Γ(s) is called the regularized incomplete gamma function. Sometimes the distinction between the regularized and unregularized versions is not explicit. For example, in Python, the function gammainc does not compute the incomplete gamma function per se but the regularized incomplete gamma function. This makes sense because the latter is often more convenient to work with numerically.

Related posts

Fourier series for Jacobi functions

I’ve mentioned a couple times that Jacobi elliptic functions are generalizations of sines and cosines. In an earlier post I showed how you can make sn and cn from sine and cosine by a nonlinear rescaling of the input. In this post I’ll look at a linear scaling of the input and a sum of sines or cosines.

This post described how the Jacobi functions depend on a parameter m. As this parameter goes to zero, the functions sn and cn converge to sine and cosine respectively. Here we’ll look more quantitatively at what happens as m varies between 0 and 1.

The Jacobi functions are elliptic functions, so they’re periodic in two directions in the complex plane. The quarter period along the real axis is called K and is given as a function of m by the following integral.

K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}

You can see that as m goes to zero, K goes to π/2, and so we have a function with period 2π.

Along the imaginary axis, the quarter period is given by K‘ = K(1-m). So as m goes to zero, K‘ goes to infinity. Said another way, the sn and cn functions lose periodicity in the imaginary direction as they converge to sine and cosine. Their fundamental rectangle gets stretched into an infinite vertical strip.

We said we’d get more quantitative, and to do that we introduce a function with the quaint name nome. The nome of a Jacobi elliptic function is defined as

q = \exp(-\pi K'/K)

Since K and K‘ are both functions of m, q is a function of m, which we plot below.

nome plot

The nome is important because it quantifies the size of Fourier coefficients.

We mentioned a linear scaling at the top of the post. That scaling is v = πu/(2K). With that change of variables, we can write down Fourier series for sn, cn, and dn in terms of q.

\begin{align*} \mathrm{sn}(u) &=\frac{2\pi}{K\sqrt{m}} \sum_{n=0}^\infty \frac{q^{n+1/2}}{1-q^{2n+1}} \sin ((2n+1)v) \\ \mathrm{cn}(u)&=\frac{2\pi}{K\sqrt{m}} \sum_{n=0}^\infty \frac{q^{n+1/2}}{1+q^{2n+1}} \cos ((2n+1)v)\\ \mathrm{dn}(u)&=\frac{\pi}{2K} + \frac{2\pi}{K} \sum_{n=1}^\infty \frac{q^{n}}{1+q^{2n}} \cos (2nv). \end{align*}

These equations are complicated, but we can begin to understand them by simply looking at the size of the coefficients.

If q is small, the denominators of the coefficients are approximately 1, so we can concentrate on the numerators, which are a geometric sequence. The coefficients go to zero rapidly, and so a small number of terms in the Fourier series are enough to approximate the Jacobi functions well.

As m goes to 1, so does q(m), but the function q(m) waits until the last minute to catch up with m because it has such a steep slope near 1. This means the function q(m) can be small even when m is large.

In the plot below, m = 0.99, and yet sn is approximated fairly well with only two sine terms. This is because q(0.99) = 0.262. For smaller m, say 0.8, it’s hard to distinguish sn and its two-term sine series approximation visually.

Two-term Fourier approximation to sn

Note that Fourier series for sn and cn only contain odd multiples of v. It’s not obvious why this is so. Because sn is an odd function, its Fourier series only contains sine terms. But this doesn’t explain why it should only contain odd frequencies. This gives our approximation extra accuracy. Our two-term approximation is effectively a three-term approximation, with the middle term being zero.

Related posts

Comparing trig functions and Jacobi functions

My previous post looked at Jacobi functions from a reference perspective: given a Jacobi function defined one way, how do I relate that to the same function defined another way, and how would you compute it?

This post explores the analogy between trigonometric functions and Jacobi elliptic functions.

Related basic Jacobi functions to trig functions

In the previous post we mentioned a connection between the argument u of a Jacobi function and the amplitude φ:

u = \int_0^{\varphi} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}

We can use this to define the functions sn and cn. Leaving the dependence on m implicit, we have

\begin{align*} \mathrm{sn}(u) &= \sin(\varphi) \\ \mathrm{cn}(u) &= \cos(\varphi) \end{align*}

If m = 0, then u = φ and so sn and cn are exactly sine and cosine.

There’s a third Jacobi function we didn’t discuss much last time, and that’s dn. It would be understandable to expect dn might be analogous to tangent, but it’s not. The function dn is the derivative of φ with respect to u, or equivalently

\mathrm{dn}(u) = \sqrt{1 - m \sin^2\varphi}

The rest of the Jacobi functions

Just as there are more trig functions than just sine and cosine, there are more Jacobi functions than sn, cn, and dn. There are two ways to define the rest of the Jacobi functions: in terms of ratios, and in terms of zeros and poles.

Ratios

I wrote a blog post one time asking how many trig functions there are. The answer depends on your perspective, and I gave arguments for 1, 3, 6, and 12. For this post, lets say there are six: sin, cos, tan, sec, csc, and ctn.

One way to look at this is to say there are as many trig functions as there are ways to select distinct numerators and denominators from the set {sin, cos, 1}. So we have tan = sin/cos, csc = 1/sin, etc.

There are 12 Jacobi elliptic functions, one for each choice of distinct numerator and denominator from {sn, cn, dn, 1}. The name of a Jacobi function is two letters, denoting the numerator and denominator, where we have {s, c, d, n} abbreviating {sn, cn, dn, 1}.

For example, cs(u) = sn(u) / cn(u) and ns(u) = 1 / sn(u).

Note that to take the reciprocal of a Jacobi function, you just reverse the two letters in its name.

Zeros and poles

The twelve Jacobi functions can be classified [1] in terms of their zeros and poles over a rectangle whose sides have length equal to quarter periods. Let’s look at an analogous situation for trig functions before exploring Jacobi functions further.

Trig functions are periodic in one direction, while elliptic functions are periodic in two directions in the complex plane. A quarter period for the basic trig functions is π/2. The six trig functions take one value out of {0, 1, ∞} at 0 and different value at π/2. So we have one trig function for each of the six ways to chose an permutation of length 2 from a set of 3 elements.

In the previous post we defined the two quarter periods K and K‘. Imagine a rectangle whose corners are labeled

s = (0, 0)
c = (K, 0)
d = (KK‘)
n = (0, K‘)

Each Jacobi function has a zero at one corner and a pole at another. The 12 Jacobi function correspond to the 12 ways to chose a permutation of two items from a set of four.

The name of a Jacobi function is two letters, the first letter corresponding to where the zero is, and the second letter corresponding to the pole. So, for example, sn has a zero at s and a pole at n. These give same names as the ratio convention above.

Identities

The Jacobi functions satisfy many identities analogous to trigonometric identities. For example, sn and cn satisfy a Pythagorean identity just like sine and cosine.

\mathrm{sn}^2 u + \mathrm{cn}^2 u = 1

Also, the Jacobi functions have addition theorems, though they’re more complicated than their trigonometric counterparts.

\begin{align*} \mathrm{sn}(u + v) &= \frac{\mathrm{sn}\,u\, \mathrm{cn}\, v\, \mathrm{dn}\,v\, + \mathrm{sn}\,v\, \mathrm{cn}\, u\, \mathrm{dn}\,u\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \\ \\ \mathrm{\mathrm{cn}\,}(u + v) &= \frac{\mathrm{cn}\, u\, \mathrm{cn}\, v\, - \mathrm{sn}\,u\, \mathrm{dn}\,u\, \mathrm{sn}\,v\, \mathrm{dn}\,v\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \\ \\ \mathrm{dn}(u + v) &= \frac{\mathrm{dn}\,u\, \mathrm{dn}\,v\, - m\, \mathrm{sn}\,u\, \mathrm{cn}\, u\, \mathrm{sn}\,v\, \mathrm{cn}\, v\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \end{align*}

Derivatives

The derivatives of the basic Jacobi functions are given below.

\begin{align*} \mathrm{sn}'(u) &= \mathrm{cn}(u)\, \mathrm{dn}(u) \\ \\ \mathrm{cn}'(u) &= -\mathrm{sn}(u) \,\mathrm{dn}(u) \\ \\ \mathrm{dn}'(u) &= -m\,\mathrm{sn}(u)\, \mathrm{cn}(u) \\ \end{align*}

Note that the implicit parameter m makes an appearance in the derivative of dn. We will also need the complementary parameter m‘ = 1 – m.

The derivatives of all Jacobi functions are summarized in the table below.

\begin{table} \centering \begin{tabular}{l|rrrr} \multicolumn{1}{l}{} & s & n & d & c \\ \cline{2-5} s & & dn cn & nd cd & nc dc \\ n & $-$ ds cs & & $m$ sd cd & sc dc \\ d & $-$ ns cs & $-m$ sn cn & & $m'$ sc nc \\ c & $-$ ns ds & $-$ sn dn & $-m'$ sd nd & \end{tabular} \end{table}

The derivatives of the basic Jacobi functions resemble those of trig functions. They may look more complicated at first, but they’re actually more regular. You could remember them all by observing the patterns pointed out below.

Let wx, yz be any permutation of {s, n, d, c}. Then the derivative of wx is proportional to yx zx. That is, the derivative of every Jacobi function f is a constant times two other Jacobi functions. The names of these two functions both end in the same letter as the name of f, and the initial letters are the two letters not in the name of f.

The proportionality constants also follow a pattern. The sign is positive if and only if the letters in the name of f appear in order in the sequence s, n, d, c. Here’s a table of just the constants.

\begin{table} \centering \begin{tabular}{l|rrrr} \multicolumn{1}{l}{} & s & n & d & c \\ \cline{2-5} s & & 1 & 1 & 1 \\ n & $-1$ & & $m$ & 1 \\ d & $-1$ & $-m$ & & $m'$ \\ c & $-1$ & $-1$ & $-m'$ & \end{tabular} \end{table}

Note that the table is skew symmetric, i.e. its transpose is its negative.

[1] An elliptic function is determined, up to a constant multiple, by its periods, zeros, and poles. So not only do the Jacobi functions have the pattern of zeros and poles described here, these patterns uniquely determine the Jacob functions, modulo a constant. For (singly) periodic functions, the period, zeros, and poles do not uniquely determine the function. So the discussion of zeros and poles of trig functions is included for comparison, but it does not define the trig functions.

Clearing up the confusion around Jacobi functions

The Jacobi elliptic functions sn and cn are analogous to the trigonometric functions sine and cosine. The come up in applications such as nonlinear oscillations and conformal mapping. Unfortunately there are multiple conventions for defining these functions. The purpose of this post is to clear up the confusion around these different conventions.

Plot of Jacobi sn

The image above is a plot of the function sn [1].

Modulus, parameter, and modular angle

Jacobi functions take two inputs. We typically think of a Jacobi function as being a function of its first input, the second input being fixed. This second input is a “dial” you can turn that changes their behavior.

There are several ways to specify this dial. I started with the word “dial” rather than “parameter” because in this context parameter takes on a technical meaning, one way of describing the dial. In addition to the “parameter,” you could describe a Jacobi function in terms of its modulus or modular angle. This post will be a Rosetta stone of sorts, showing how each of these ways of describing a Jacobi elliptic function are related.

The parameter m is a real number in [0, 1]. The complementary parameter is m‘ = 1 – m. Abramowitz and Stegun, for example, write the Jacobi functions sn and cn as sn(um) and cn(um). They also use m1 = rather than m‘ to denote the complementary parameter.

The modulus k is the square root of m. It would be easier to remember if m stood for modulus, but that’s not conventional. Instead, m is for parameter and k is for modulus. The complementary modulus k‘ is the square root of the complementary parameter.

The modular angle α is defined by m = sin² α.

Note that as the parameter m goes to zero, so does the modulus k and the modular angle α. As any one of these three goes to zero, the Jacobi functions sn and cn converge to their counterparts sine and cosine. So whether your dial is the parameter, modulus, or modular angle, sn converges to sine and cn converges to cosine as you turn the dial toward zero.

As the parameter m goes to 1, so does the modulus k, but the modular angle α goes to π/2. So if your dial is the parameter or the modulus, it goes to 1. But if you think of your dial as modular angle, it goes to π/2. In either case, as you turn the dial to the right as far as it will go, sn converges to hyperbolic secant, and cn converges to the constant function 1.

Quarter periods

In addition to parameter, modulus, and modular angle, you’ll also see Jacobi function described in terms of K and K‘. These are called the quarter periods for good reason. The functions sn and cn have period 4K as you move along the real axis, or move horizontally anywhere in the complex plane. They also have period 4iK‘. That is, the functions repeat when you move a distance 4K‘ vertically [2].

The quarter periods are a function of the modulus. The quarter period K along the real axis is

K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}

and the quarter period K‘ along the imaginary axis is given by K‘(m) = K(m‘) = K(1 – m).

The function K(m) is known as “the complete elliptic integral of the first kind.”

Amplitude

So far we’ve focused on the second input to the Jacobi functions, and three conventions for specifying it.

There are two conventions for specifying the first argument, written either as φ or as u. These are related by

u = \int_0^{\varphi} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}

The angle φ is called the amplitude. (Yes, it’s an angle, but it’s called an amplitude.)

When we said above that the Jacobi functions had period 4K, this was in terms of the variable u. Note that when φ = π/2, uK.

Jacobi elliptic functions in Mathematica

Mathematica uses the u convention for the first argument and the parameter convention for the second argument.

The Mathematica function JacobiSN[u, m] computes the function sn with argument u and parameter m. In the notation of A&S, sn(um).

Similarly, JacobiCN[u, m] computes the function cn with argument u and parameter m. In the notation of A&S, cn(um).

We haven’t talked about the Jacobi function dn up to this point, but it is implemented in Mathematica as JacobiDN[u, m].

The function that solves for the amplitude φ as a function of u is JacobiAmplitude[um m].

The function that computes the quarter period K from the parameter m is EllipticK[m].

Jacobi elliptic functions in Python

The SciPy library has one Python function that computes four mathematical functions at once. The function scipy.special.ellipj takes two arguments, u and m, just like Mathematica, and returns sn(um), cn(um), dn(um), and the amplitude φ(um).

The function K(m) is implemented in Python as scipy.special.ellipk.

Related posts

[1] The plot was made using JacobiSN[0.5, z] and the function ComplexPlot described here.

[2] Strictly speaking, 4iK‘ is a period. It’s the smallest vertical period for cn, but 2iK‘ is the smallest vertical period for sn.

Unit speed curve parameterization

A friend asked me a question the other day that came out of a graphics application. He needed to trace out an ellipse in such a way that the length of curve traced out each second was constant. For a circle, the problem is simple: (cos(t), sin(t)) will trace out a circle covering a constant amount of arc length per unit time. The analogous parameterization for an ellipse, (a cos(t), b sin(t)) will move faster near the longer semi-axis and slower near the shorter one.

There’s a general solution to the problem for any curve. Given a parameterization p(t), where p is vector-valued, the length covered from time 0 to time t is

s(t) = \int_0^t || p'(\tau) || \,d\tau

If you change the time parameterization by inverting this function, solving for t as a function of s, then the total length of curve traversed by p(t(s)) up to time s is s. This is called either the unit speed parameterization or parameterization by arc length.

The hard part is inverting s(t). If you had to find a unit speed parameterization in a calculus class, the problem would be carefully designed so the function s(t) is easy to invert. Real applications don’t usually work out so easily.

Digression on elliptic integrals and elliptic functions

My response to my friend’s question was that there probably is a closed-form unit speed parameterization, but it would probably involve elliptic functions. He didn’t need much resolution, and decided to do something ad hoc.

Starting with the parameterization  p(t) = (a cos(t), b sin(t)), the arc length s(t) is given by a special function known as an “incomplete elliptic integral of the second kind.” I remembered that the Jacobi elliptic functions were related to the inverses of elliptic integrals, so my hunch was that you could make a unit speed parameterization using Jacobi elliptic functions. Maybe you can, but it’s not as easy as I thought because the Jacobi functions are related to the inverses of elliptic integrals of the first kind.

Elliptic integrals are so named because they were first identified by computing arc length for a (portion of) an ellipse. Elliptic functions were discovered by inverting elliptic integrals, but not the same class of elliptic integrals that give the arc length of an ellipse. (There may well be a transformation that makes these more directly related, but I’m not aware of it.)

Incidentally, elliptic curves are related to elliptic functions, but they are not ellipses! There is a connection from ellipses to elliptic curves, but it’s historical and indirect.

What if we had a more general curve than an ellipse, say something parameterized by cubic splines? Cubic splines are piecewise cubic polynomials, patched together in such a way that the first and second derivatives are continuous across the patches. We can find length of a spline by finding the length of each polynomial patch.

If p(t) is the parameterization of a curve in 2 or 3 dimensions (or really any number of dimensions) and each component of p is a cubic polynomial, then each component of the derivative of p is a quadratic polynomial, and so the sum of the squares of the components is a fourth degree polynomial. So finding the arc length involves integrating the square root of a fourth degree polynomial. This makes it an elliptic integral!

Unfortunately, knowing that the arc length of a cubic spline corresponds to an elliptic integral is not so useful because it could be any type of elliptic integral, depending on its parameters. You’d have to do some work first to put it into a form where you could call on elliptic integrals to finish your problem.

Numerically computing arc length and unit speed parameterization

The elliptic integral path is something of a dead end. It could still be useful if you needed high accuracy, or if you had some restriction on the class of curves you’re interested in. But in general, you’d need to use numerical integration to find the arc length.

You could also find unit-speed parameterizations numerically, using root-finding to invert s(t) at particular points. Since s is an increasing function of t, you could use a bisection method, which is not the most efficient approach but very robust.

It takes a fair amount of computation to carry root finding where each function evaluation requires computing a numerical integral. But this would work, and depending on your context it could be efficient enough.

If you needed more efficiency, say for a real-time embedded system, you could take a different approach. Your spline is probably an approximation to something else, and so your solution only needs to be as accurate as the spline approximation. This gives you the wiggle room to do something more efficient. You might change your parameterization slightly to make the arc length calculations easier, or come up with a way to approximately invert the arc length function, something that takes advantage of your particular problem and error tolerance.

Related posts

Surprising curves with sine and sn

In the previous post I said that the Jacobi functions are like trig functions. That’s true if you look along the real axis. If you look at the rest of the complex plane you’ll see how they can be very different.

The sine function is periodic along the real axis, but it grows exponentially along the imaginary axis. I mentioned parenthetically last time that the Jacobi functions are not just periodic but doubly periodic. That means that not only are they periodic as you move along the real axis, they’re periodic when you move along any line in the complex plane [1].We’ll illustrate this with some plots.

It will make our code more compact if we first define a function to split a complex number into its real and imaginary parts.

    pair[z_] := {Re[z], Im[z]}

Here’s what it looks like when you plot the real and imaginary parts of the sine function along the line (10 + 0.5i)t.

    
    ParametricPlot[
        pair[ Sin[(0.5 I + 10)t] ], 
        {t, 0, 10}, 
        PlotRange -> All
    ]

plot of sine with complex argument

By contrast, here’s a plot of the sn function along the line 1.25 + it.

    ParametricPlot[
        pair[ JacobiSN[1.25 + I t, 0.5] ], 
        {t, 0, 10}
    ]

It’s a closed loop because sn is periodic in every direction. (By the way, the curve looks like an egg. More posts about egg-shaped curves starting here.)

When you plot sn along various lines you’ll always get closed curves, but you won’t always get such round curves. Here’s an example that doesn’t look like a closed loop because the curve turns around sharply at each end.

 
    ParametricPlot[
        pair[ JacobiCN[ (I + 0.5) t, 0.5] ], 
        {t, 0, 10}
    ]

To show that it really is periodic, we’ll add a vertical time dimension to visualize how the curve is traced out over time.

    triple[z_, t_] = {Re[z], Im[z], t}
    ParametricPlot3D[
        triple[JacobiSN[(I + 0.5) t, 0.5], 0.1 t], 
        {t, 0, 20}
    ]

Here’s another curve plotting sn along a different line.

    ParametricPlot[
        pair[ JacobiSN[(0.9 + I) t, 0.5] ], 
        {t, 0, 55},
        PlotRange -> All, 
        PlotPoints -> 100
    ]

As before, we can add a time dimension to imagine how the curve is traced out.

    ParametricPlot3D[
        triple[JacobiSN[(0.9 + I) t, 0.5], 0.1 t], 
        {t, 0, 150},
        PlotRange -> All,
        PlotPoints -> 200
]

Here’s a variation on the plot above that stretches the vertical axis so you can better see what’s going on.

   
    ParametricPlot3D[
        triple[JacobiSN[(0.9 + I) t, 0.5], 0.2 t], 
        {t, 0, 150},
        PlotRange -> All,
        PlotPoints -> 200
    ]

[1] To be more precise, elliptic functions are periodic in two linearly independent directions, and thus in any direction that’s an integer linear combination of those two directions. So they’re exactly periodic in a countable number of directions almost periodic in every other direction.

System of Jacobi elliptic functions

Jacobi’s elliptic functions are sorta like trig functions. His functions sn and cn have names that reminiscent of sine and cosine for good reason. These functions come up in applications such as the nonlinear pendulum (i.e. when θ is too
large to assume θ is a good enough approximation to sin θ) and in conformal mapping.

I ran across an article [1] yesterday that shows how Jacobi’s three elliptic functions—sn, cn, and dn—could be defined by one dynamical system

\begin{eqnarray*} x' &=& yz \\ y' &=& -zx \\ z' &=& -k^2 xy \end{eqnarray*}

with initial conditions x(0) = 0, y(0) = 1, and z(0) = 1.

The parameter k is the modulus. (In Mathematica’s notation below, k² is the parameter. See this post on parameterization conventions.) As k decreases to 0, sn converges to sine, cn to cosine, and dn to 1. As k increases to 1, sn converges tanh, and cn and dn converge to sech. So you could think of k as a knob you turn to go from being more like circular functions (ordinary trig functions) to more like hyperbolic functions.

Since we have a dynamical system, let’s plot the solution, varying the modulus each time.The Jacobi functions are periodic (in fact they’re doubly periodic) and so the plots will be closed loops.

f[t_, m_] = {JacobiSN[t, m], JacobiCN[t, m], JacobiDN[t, m]}
ParametricPlot3D[f[t, 0.707], {t, 0, 10}, AspectRatio -> 1]

phase portrait k = 1/2

ParametricPlot3D[f[t, 0.99], {t, 0, 20}, AspectRatio -> 1]

phase portrait k = 1/2

ParametricPlot3D[f[t, 0.01], {t, 0, 20}, AspectRatio -> 1]

phase portrait k = 1/2

Note that this last plot is nearly flat because the modulus is small and so z is nearly constant. The small modulus also makes the phase portrait nearly circular because x is approximately sine and y is approximately cosine.

[1] Kenneth R. Meyer. Jacobi Elliptic Functions from a Dynamical Systems Point of View. The American Mathematical Monthly, Vol. 108, No. 8 (Oct., 2001), pp. 729-737

Approximating gamma ratios

Ratios of gamma functions come up often in applications. If the two gamma function arguments differ by an integer, then it’s easy to calculate their ratio exactly by using (repeatedly if necessary) the fact at Γ(x + 1) = x Γ(x).

If the arguments differ by 1/2, there is no closed formula, but the there are useful approximations. I’ve needed something like this a few times lately.

The simplest approximation is

\frac{\Gamma\left(x + 1 \right) }{\Gamma\left(x + \frac{1}{2} \right) } \sim x^{1/2}

You could motivate or interpret this as saying Γ(x + 1/2) is approximately the geometric mean between Γ(x + 1) and Γ(x). As we’ll see in the plot below, this approximation is good to a couple significant figures for moderate values of x.

There is another approximation that is a little more complicated but much more accurate.

\frac{\Gamma\left(x + 1 \right) }{\Gamma\left(x + \frac{1}{2} \right) } \sim \left(x^2 + \frac{x}{2} + \frac{1}{8}\right)^{1/4}

The following plot shows the relative error in both approximations.

gamma ratio approximation errors

By the way, the first approximation above is a special case of the more general approximation

\frac{\Gamma(x+a)}{\Gamma(x)} \sim x^a

Source:  J. S. Frame. An Approximation to the Quotient of Gamma Function. The American Mathematical Monthly, Vol. 56, No. 8 (Oct., 1949), pp. 529-535

 

Hypergeometric functions are key

From Orthogonal Polynomials and Special Functions by Richard Askey:

At first the results we needed were in the literature but after a while we ran out of known results and had to learn something about special functions. This was a very unsettling experience for there were very few places to go to really learn about special functions. At least that is what we thought. Actually there were many, but the typical American graduate education which we had did not include anything about hypergeometric functions. And hypergeometric functions are the key to this subject, as I have found out after many years of fighting them.

Emphasis added.

Askey’s book was written in 1975, and he was describing his experience from ten years before that. Special functions, and in particular hypergeometric functions, went from being common knowledge among mathematicians at the beginning of the 20th century to being arcane by mid century.

I learned little about special functions and nothing about hypergeometric functions as a graduate student. I first ran into hypergeometric functions reading in Concrete Mathematics how they are used in combinatorics and in calculating sums in closed form. Then when I started working in statistics I found that they are everywhere.

Hypergeometric functions are very useful, but not often taught anymore. Like a lot of useful mathematics, they fall between two stools. They’re considered too advanced or arcane for the undergraduate curriculum, and not a hot enough area of research to be part of the graduate curriculum.

Related posts:

Orthogonal polynomials and the beta distribution

This post shows a connection between three families of orthogonal polynomials—Legendre, Chebyshev, and Jacobi—and the beta distribution.

Legendre, Chebyshev, and Jacobi polynomials

A family of polynomials Pk is orthogonal over the interval [-1, 1] with respect to a weight w(x) if

\int_{-1}^1 P_m(x) P_n(x) w(x) \, dx = 0

whenever mn.

If w(x) = 1, we get the Legendre polynomials.

If w(x) = (1 – x²)-1/2 we get the Chebyshev polynomials.

These are both special cases of the Jacobi polynomials which have weight w(x) = (1- x)α (1 + x)β. Legendre polynomials correspond to α = β = 0, and Chebyshev polynomials correspond to α = β = -1/2.

Connection to beta distribution

The weight function for Jacobi polynomials is a rescaling of the density function of a beta distribution. The change of variables x = 1 – 2u shows

\int_{-1}^1 f(x) (1-x)^\alpha (1+x)^\beta \, dx = 2^{\alpha + \beta + 1}\int_0^1 f(1-2u) u^\alpha (1-u)^\beta \,du

The right side is proportional to the expected value of f(1 – 2X) where X is a random variable with a beta(α + 1, β+1) distribution. So for fixed α and β, if mn and X has a beta(α + 1, β+1) distribution, then the expected value of Pm(1 – 2X) Pn(1 – 2X) is zero.

While we’re at it, we’ll briefly mention two other connections between orthogonal polynomials and probability: Laguerre polynomials and Hermite polynomials.

Laguerre polynomials

The Laguerre polynomials are orthogonal over the interval [0, ∞) with weight w(x) = xα exp(-x), which is proportional to the density of a gamma random variable with shape α+1 and scale 1.

Hermite polynomials

There are two minor variations on the Hermite polynomials, depending on whether you take the weight to be exp(-x²) or exp(-x²/2). These are sometimes known as the physicist’s Hermite polynomials and the probabilist’s Hermite polynomials. Naturally we’re interested in the latter. The probabilist’s Hermite polynomials are orthogonal over (-∞, ∞) with the standard normal (Gaussian) density as the weight.

Related posts

Length of a rose

The polar graph of r = cos(kθ) is called a rose. If k is even, the curve will trace out 2k petals as θ runs between 0 and 2π. If k is odd, it will trace out k petals, tracing each one twice. For example, here’s a rose with k = 5.

(I rotated the graph 36° so it would be symmetric about the vertical axis rather than the horizontal axis.)

The arc length of a curve in polar coordinates is given by

\int_a^b \sqrt{r^2 + \left(\frac{dr}{d\theta}\right)^2}\, d\theta

and so we can use this find the length. The integral doesn’t have a closed form in terms of elementary functions. Instead, the result turns out to use a special function E(x), the “complete elliptic integral of the second kind,” defined by

E(m) = \int_0^{\pi/2} \sqrt{1 - m \sin^2 x} \, dx

Here’s the calculation for the length of a rose:

\int_0^{2\pi} \sqrt{r^2 + (r')^2}\, d\theta &=& \int_0^{2\pi} \sqrt{\cos^2 k\theta + k^2 \sin^2 k\theta} \, d\theta \\ &=& \int_0^{2\pi} \sqrt{\cos^2 u + k^2 \sin^2 u} \, du \\ &=& 4 \int_0^{\pi/2} \sqrt{\cos^2 u + k^2 \sin^2 u} \, du \\ &=& 4 \int_0^{\pi/2} \sqrt{1 + (k^2-1) \sin^2 u} \, du \\ &=& 4 E(-k^2 + 1)

So the arc length of the rose r = cos(kθ) with θ running from 0 to 2π is 4 E(-k² + 1). You can calculate E in SciPy with scipy.special.ellipe.

If we compute the length of the rose at the top of the post, we get 4 E(-24) = 21.01. Does that pass the sniff test? Each petal goes from r = 0 out to r = 1 and back. If the petal were a straight line, this would have length 2. Since the petals are curved, the length of each is a little more than 2. There are five petals, so the result should be a little more than 10. But we got a little more than 20. How can that be? Since 5 is odd, the rose with k = 5 traces each petal twice, so we should expect a value of a little more than 20, which is what we got.

As k gets larger, the petals come closer to being straight lines. So we should expect that 4E(-k² + 1) approaches 4k as k gets large. The following plot of E(-k² + 1) – k provides empirical support for this conjecture by showing that the difference approaches 0, and gives an idea of the rate of convergence. It should be possible to prove that, say, that E(-k²) asymptotically approaches k, but I haven’t done this.

Related posts:

Denver airport, Weierstrass, and A&S

Last night I was driving toward the Denver airport and the airport reminded me of the cover of Abramowitz and Stegun’s Handbook of Mathematical Functions.

Here’s the airport:

Denver airport

And here’s the book cover:

I’ve written about the image on book cover before. Someone asked me what function it graphed and I decided it was probably the Weierstrass ℘ function.

For more on Weierstrass’ elliptic function and why I think that’s what’s on the cover of A&S, see this post.

Photo of Denver airport via Wikipedia.

Function on cover of Abramowitz & Stegun

Someone mailed me this afternoon asking if I knew what function was graphed on the cover of Abramowitz and Stegun’s famous Handbook of Mathematical Functions.

Here’s a close-up of the graph from a photo of my copy of A&S.

This was my reply.

It looks like a complex function of a complex variable. I assume the height is the magnitude and the markings on the graph are the phase. That would make it an elliptic function because it’s periodic in two directions.

It has one pole and one zero in each period. I think elliptic functions are determined, up to a constant, by their periods, zeros, and poles, so it should be possible to identify the function.

In fact, I expect it’s the Weierstrass p function. More properly, the Weierstrass ℘ function, sometimes called Weierstass’ elliptic function. (Some readers will have a font installed that will properly render ℘ and some not. More on the symbol ℘ here.)

Related posts:

Bessel series for a constant

Fourier series express functions as a sum of sines and cosines of different frequencies. Bessel series are analogous, expressing functions as a sum of Bessel functions of different orders.

Fourier series arise naturally when working in rectangular coordinates. Bessel series arise naturally when working in polar coordinates.

The Fourier series for a constant is trivial. You can think of a constant as a cosine with frequency zero.

The Bessel series for a constant is not as simple, but more interesting. Here we have

1 = J_0(x) + 2J_2(x) + 2J_4(x) + 2J_6(x) + \cdots

Since

J_{-n}(x) = (-1)^n J_n(x)

we can write the series above more symmetrically as

1 = \sum_{n=-\infty}^\infty J_{2n}(x)

Related posts: