# Sum of all Spheres

I ran across a video this afternoon that explains that the sum of volumes of all even-dimensional unit spheres equals eπ.

Why is that? Define vol(n) to be the volume of the unit sphere in dimension n. Then

and so the sum of the volumes of all even dimensional spheres is

But what if you wanted to sum the volumes of all odd dimensional unit spheres? Or all dimensions, even and odd?

The answers to all these questions are easy in terms of the Mittag-Leffler function that I blogged about a while back. This function is defined as

and reduces to the exponential function when α = β = 1.

The sum of the volumes of all unit spheres of all dimensions is E1/2, 1(√π). And from the post mentioned above,

where erfc is the complementary error function. So the sum of all volumes of spheres is exp(π) erfc(-√π).

Now erfc(-√π) ≈ 1.9878 and so this says the sum of the volumes of spheres of all dimensions is about twice the sum of the even dimensional spheres alone. And so the sum of the odd dimensional unit sphere volumes is almost the same as the sum of the even dimensional ones.

By the way, you could easily answer other questions about sums of sphere volumes in terms of the Mittag-Leffler function. For example, if you want to add up the volumes in all dimensions that are a multiple of 3, you get E3/2, 13/2).

# Bessel function crossings

The previous looked at the angles that graphs make when they cross. For example, sin(x) and cos(x) always cross with the same angle. The same holds for sin(kx) and cos(kx) since the k simply rescales the x-axis.

The post ended with wondering about functions analogous to sine and cosine, such as Bessel functions. This post will look at that question in more detail. Specifically we’ll look at the functions Jν and Yν.

Because these two Bessel functions satisfy the same second order linear homogeneous differential equation, the Strum separation theorem says that their zeros are interlaced: between each pair of consecutive zeros of Jν is exactly one zero of Yν, and between each pair of consecutive zeros of Yν there is exactly one zero of Jν.

In the following Python code, we find zeros of Jν, then look in between for places where Jν and Yν cross. Next we find the angle the two curves make at each intersection and plot the angles.

    from scipy.special import jn_zeros, jv, yv
from scipy.optimize import bisect
from numpy import empty, linspace, arccos
import matplotlib.pyplot as plt

n = 3 # bessel function order
N = 100 # number of zeros

z = jn_zeros(n, N) # Zeros of J_n
crossings = empty(N-1)

f = lambda x: jv(n,x) - yv(n,x)
for i in range(N-1):
crossings[i] = bisect(f, z[i], z[i+1])

def angle(n, x):
# Derivatives of J_nu and Y_nu
dj = 0.5*(jv(n-1,x) - jv(n+1,x))
dy = 0.5*(yv(n-1,x) - yv(n+1,x))

top = 1 + dj*dy
bottom = ((1 + dj**2)*(1 + dy**2))**0.5
return arccos(top/bottom)

y = angle(n, crossings)
plt.plot(y)
plt.xlabel("Crossing number")
plt.show()


This quadratic behavior is what we should expect from the asymptotics of Jν and Yν: For large arguments they act like shifted and rescaled versions of sin(x)/√x. So if we looked at √xJν and √xYν rather than Jν and Yν we’d expect the angles to reach some positive asymptote, and they do, as shown below.

# 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.)

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.

The A functions are given recursively by

and

for n > 1.

The F functions are given by

where

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.

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,

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.

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

## 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

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.

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

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

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

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

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.

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.

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.

# 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 φ:

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

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

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

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

## Derivatives

The derivatives of the basic Jacobi functions are given below.

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.

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.

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.

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

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

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

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.

# 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
]


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

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]


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


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


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