Reproducing a hand-drawn plot

Posted on 3 March 2025 by John

The plots in old (i.e. pre-computer) math books are impressive. These plots took a lot of effort to produce, and so they weren’t created lightly. Consequently they tend to be aesthetically and mathematically interesting. A few weeks ago I recreated a plot from A&S using Mathematica, and today I’d like to do the same for the plot below from a different source [1].

Here is my approximate reproduction.

I’ll give the mathematical details below for those who are interested.

The plot shows “normalized associated Legendre functions of the first kind.” There are two families of Legendre polynomials, denoted P and Q; we’re interested in the former. “Associated” means polynomials that are derived the Legendre polynomials by taking derivatives. Normalized means the polynomials are multiplied by constants so that their squares integrate to 1 over [−1, 1].

Mathematica has a function LegendreP[n, m, x] that implements the associated Legendre polynomials P_n^m(x). I didn’t see that Mathematica has a function for the normalized version of these functions, so I rolled by own.

f[n_, m_, x_] := (-1)^n LegendreP[n, m, x] 
        Sqrt[(2 n + 1) Factorial[n - m]/(2 Factorial[n + m])]

I added the alternating sign term up front after discovering that apparently the original plot used a different convention for defining P_n^m than Mathematica uses.

I make my plot by stacking the plots created by

Plot[Table[f[n, n, x], {n, 1, 8}],  {x, 0, 1}]

and

Plot[Table[f[n + 4, n, x], {n, 0, 4}],  {x, 0, 1}]

The original plot shows P₄(x). I used the fact that this equals P₄⁰(x) to simplify the code. I also left out the flat line plotting P₀ because I thought that looked better.

Related post: Duplicating a Hankel function plot from A&S.

[1] Tables Of Functions With Formulae And Curves by Fritz Emde, published in 1945. Available on Archive.org.

Multiple angle asymmetry

Posted on 13 February 2025 by John

The cosine of a multiple of θ can be written as a polynomial in cos θ. For example,

cos 3θ = 4 cos³ θ − 3 cos θ

and

cos 4θ = 8 cos⁴ θ − 8 cos² θ + 1.

But it may or may not be possible to write the sine of a multiple of θ as a polynomial in sin θ. For example,

sin 3θ = −4 sin³ θ + 3 sin θ

but

sin 4θ = − 8 sin³ θ cos θ + 4 sin θ cos θ

It turns out cos nθ can always be written as a polynomial in cos θ, but sin nθ can be written as a polynomial in sin θ if and only if n is odd. We will prove this, say more about sin nθ for even n, then be more specific about the polynomials alluded to.

Proof

We start by writing exp(inθ) two different ways:

cos nθ + i sin nθ = (cos θ + i sin θ)ⁿ

The real part of the left hand side is cos nθ and the real part of the right hand side contains powers of cos θ and even powers of sin θ. We can convert the latter to cosines by replacing sin² θ with 1 − cos² θ.

The imaginary part of the left hand side is sin nθ. If n is odd, the right hand side involves odd powers of sin θ and even powers of cos θ, in which case we can replace the even powers of cos θ with even powers of sin θ. But if n is even, every term in the imaginary part will involve odd powers of sin θ and odd powers of cos θ. Every odd power of cos θ can be turned into terms involving a single cos θ and an odd power of sin θ.

We’ve proven a little more than we set out to prove. When n is even, we cannot write sin nθ as a polynomial in sin θ, but we can write it as cos θ multiplied by an odd degree polynomial in sin θ. Alternatively, we could write sin nθ as sin θ multiplied by an odd degree polynomial in cos θ.

Naming polynomials

The polynomials alluded to above are not arbitrary polynomials. They are well-studied polynomials with many special properties. Yesterday’s post on Chebyshev polynomials defined T_n(x) as the nth degree polynomial for which

T_n(cos θ) = cos nθ.

That post didn’t prove that the right hand side is a polynomial, but this post did. The polynomials T_n(x) are known as Chebyshev polynomials of the first kind, or sometimes simply Chebyshev polynomials since they come up in application more often than the other kinds.

Yesterday’s post also defined Chebyshev polynomials of the second kind by

U_n(cos θ) sin θ = sin (n+1)θ.

So when we say cos nθ can be written as a polynomial in cos θ, we can be more specific: that polynomial is T_n.

And when we say sin nθ can be written as sin θ times a polynomial in cos θ, we can also be more specific:

sin nθ = sin θ U_n−1(cos θ).

Solving trigonometric equations

A couple years ago I wrote about systematically solving trigonometric equations. That post showed that any polynomial involving sines and cosines of multiples of θ could be reduced to a polynomial in sin θ and cos θ. The results in this post let us say more about this polynomial, that we can write it in terms of Chebyshev polynomials. This might allow us to apply some of the numerous identities these polynomials satisfy and find useful structure.

Posthumous Chebyshev Polynomials

Posted on 12 February 2025 by John

Two families of orthogonal polynomials are named after Chebyshev because he explored their properties. These are prosaically named Chebyshev polynomials of the first and second kind.

I recently learned there are Chebyshev polynomials of the third and fourth kind as well. You might call these posthumous Chebyshev polynomials. They were not developed by Mr. Chebyshev, but they bear a family resemblance to the polynomials he did develop.

The four kinds of Chebyshev polynomials may be defined in order as follows.

$\begin{align*} T_n(\cos\theta) &= \cos n\theta \\ U_n(\cos\theta) &= \frac{\sin (n+1)\theta}{\sin \theta} \\ V_n(\cos\theta) &= \frac{\cos \left(n+\frac{1}{2}\right)\theta}{\cos \frac{1}{2}\theta} \\ W_n(\cos\theta) &= \frac{\sin \left(n+\frac{1}{2}\right)\theta}{\sin \frac{1}{2}\theta} \\ \end{align*}$

It’s not obvious that these definitions even make sense, but in each case the right hand side can be expanded into a sum of powers of cos θ, i.e. a polynomial in cos θ. [1]

All four kinds of Chebyshev polynomials satisfy the same recurrence relation

$P_n(x) = 2x\,P_{n-1}(x) - P_{n-2}(x)$

for n ≥ 2 and P₀ = 1 but with different values of P₁, namely x, 2x, 2x − 1, and 2x + 1 respectively [2].

Plots

We can implement Chebyshev polynomials of the third kind using the recurrence relation above.

def V(n, x):
    if n == 0: return 1
    if n == 1: return 2*x - 1
    return 2*x*V(n-1, x) - V(n-2, x)

Here is a plot of V_n(x) for n = 0, 1, 2, 3, 4.

The code for implementing Chebyshev polynomials of the fourth kind is the same, except the middle line becomes

    if n == 1: return 2*x + 1

Here is the corresponding plot.

Square roots

The Chebyshev polynomials of the first and third kind, and polynomials of the second and fourth kind, are related as follows:

$\begin{align*} V_n(x)&=\sqrt\frac{2}{1+x}T_{2n+1}\left(\sqrt\frac{x+1}{2}\right) \\ W_n(x)&=U_{2n}\left(\sqrt\frac{x+1}{2}\right) \end{align*}$

To see that the expressions on the right hand side really are polynomials, note that Chebyshev polynomials of the first and second kinds are odd for odd orders and even for even orders [3]. This means that in the first equation, every term in T_{2n + 1} has a factor of √(1 + x) that is canceled out by the 1/√(1 + x) term up front. In the second equation, there are only even powers of the radical term so all the radicals go away.

You could take the pair of equations above as the definition of Chebyshev polynomials of the third and fourth kind, but the similarity between these polynomials and the original Chebyshev polynomials is more apparent in the definition above using sines and cosines.

The square roots hint at how these polynomials first came up in applications. According to [2], Chebyshev polynomials of the third and fourth kind

have been called “airfoil polynomials”, since they are appropriate for approximating the single square root singularities that occur at the sharp end of an airfoil.

Dirichlet kernel

There’s an interesting connection between Chebyshev polynomials of the fourth kind and Fourier series.

The right hand side of the definition of W_n is known in Fourier analysis as D_n, the Dirichlet kernel of order n.

$D_n(\theta) = \frac{\sin \left(n+\frac{1}{2}\right)\theta}{\sin \frac{1}{2}\theta}$

The nth order Fourier series approximation of f, i.e. the sum of terms −n through n in the Fourier series for f is the convolution of f with D_n, times 2π.

$(D_n * f)(\theta) = 2\pi \sum_{k=-n}^n \hat{f}(k) \exp(ik\theta)$

Note that D_n(θ) is a function of θ, not of x. The equation W_n(cos θ) = D_n(θ) defines W_n(x) where x = cos θ. To put it another way, D_n(θ) is not a polynomial, but it can be expanded into a polynomial in cos θ.

[1] Each function on the right hand side is an even function, which implies it’s at least plausible that each can be written as powers of cos θ. In fact you can apply multiple angle trig identities to work out the polynomials in cos θ.

[2] J.C. Mason and G.H. Elliott. Near-minimax complex approximation by four kinds of Chebyshev polynomial expansion. Journal of Computational and Applied Mathematics 46 (1993) 291–300

[3] This is not true of Chebyshev polynomials of the third and fourth kind. To see this note that V₁(x) = 2x − 1, and W₁(x) = 2x + 1, neither of which is an odd function.

Duplicating Hankel plot from A&S

Posted on 22 January 2025 by John

Abramowitz and Stegun has quite a few intriguing plots. The post will focus on the follow plot, Figure 9.4, available here.

A&S figure 9.4

We will explain what the plot is and approximately reproduce it.

The plot comes from the chapter on Bessel functions, but the caption says it is a plot of the Hankel function H₀⁽¹⁾. Why a plot of a Hankel function and not a Bessel function? The Hankel functions are linear combinations of the Bessel functions of the first and second kind:

H_0⁽¹⁾ = J₀ + i Y₀

More on that Hankel functions and their relations to Bessel functions here.

The plot is the overlay of two kinds of contour plots: one for lines of constant magnitude and one for lines of constant phase. That is, if the function values are written in the form re^iθ then one plot shows lines of constant r and one plot shows lines of constant θ.

We can roughly reproduce the plot of magnitude contours with the following Mathematica command:

ContourPlot[Abs[HankelH1[0, x + I y]], {x, -4, 2 }, {y, -1.5 , 1.5 }, 
 Contours -> 20, ContourShading -> None, AspectRatio -> 1/2]

This produces the following plot.

Absolute value contour

Similarly, we can replace Abs with Arg in the Mathematica command and increase Contours to 30 to obtain the following phase contour plot.

Phase contour

Finally, we can stack the two plots on top of each other using Mathematica’s Show command.

Magnitude and phase contours

By the way, you can clearly see the branch cut in the middle. The Hankel function is continuous (even analytic) as you move clockwise from the second quadrant around to the third, but it is discontinuous across the negative real axis because of the branch cut.

An unexpected triangle

Posted on 18 January 2025 by John

Let J(x) be the function plotted below.

This is the Bessel function J₁, but we drop the subscript because it’s the only Bessel function we’re interested in for this post. You can think of J as a sort of damped sine.

We can create versions of J with different frequencies by multiplying the argument x by different constants. Suppose we create versions with three different frequencies — J(ax), J(bx), and J(cx) — and integrate their product. This defines a function f of the frequencies.

$f(a, b, c) = \int_0^\infty J(ax) \, J(bx) \, J(cx)\, dx$

We can evaluate the integral defining f using Sonine’s formula [1]

$f(a, b, c) = \frac{2\Delta(a, b, c)}{\pi abc}$

where Δ(a, b, c) is the area of a triangle with sides a, b, c, if such a triangle exists, and zero otherwise.

It’s amazing that this formula takes three parameters with no intrinsic geometric meaning and out pops the area of a triangle with such sides.

Numerical (mis)application

It would be ill-advised, but possible, to invert Sonine’s formula and use it to find the area of a triangle; Heron’s formula would be a far better choice. But just for fun, we’ll take the ill-advised route.

from numpy import linspace, pi, sqrt, inf
from scipy.special import j1
from scipy.integrate import quad

def heron(a, b, c):
    s = (a + b + c)/2
    return sqrt(s*(s-a)*(s-b)*(s-c))

def g(a, b, c):
    integrand = lambda x: j1(a*x) * j1(b*x) * j1(c*x)
    i, _ = quad(integrand, 0, inf, limit=500)
    return i

def area(a, b, c):
    return g(a, b, c)*pi*a*b*c/2

print(area(3,4,5), heron(3,4,5))

SciPy’s quad function has difficulty with the integration, and rightfully issues a warning. The code increases the limit parameter from the default value of 50 to 500, improving the accuracy but not eliminating the warning. The area function computes the error of a 3-4-5 triangle to be 5.9984 and the heron function computes the exact value 6.

Update: I tried the numerical integration above in Mathematica and it correctly computed the integral to 10 decimal places with no help. I suspect it is detecting the oscillatory nature of the integral and using Levin’s integration method; when I explicit set the integration to be Levin’s method, I get the same result.

Impossible triangles

You could calculate the area of a triangle from Heron’s theorem or from Sonine’s theorem. The results are identical when the parameters a, b, and c form a valid triangle. But the two approaches diverge when a, b, and c do not form a triangle. If you pass in parameters like (3, 1, 1) then Heron’s theorem gives a complex number and Sonine’s theorem yields zero.

[1] Discovered by Nikolay Yakovlevich Sonin (1849–1915). The formula is usually referred to as Sonine’s formula rather than Sonin’s formula, as far as I’ve seen. This variation is minor compared to what transliteration does to other Russian names.

Sonine’s formula is more general, extending to Bessel functions J_ν with Re(ν) > ½. I chose ν = 1 for this post because the Sonin’s formula is simplest in this case.

Dimensional analysis for gamma function values

Posted on 16 January 2025 by John

Sometimes it’s useful to apply dimensional analysis where it doesn’t belong, to imagine things having physical dimension when they don’t. This post will look at artificially injecting dimensions into equations involving factorials and related functions.

Factorials

The factorial of n is defined as the product of n terms. If each of these terms had units of length, the factorial would have units of n-dimensional volume. It occurred to me this morning that a lot of identities involving factorials make sense dimensionally if you pretend the terms in the identities have units. The expressions on both sides of an equation have the same dimension, and only terms of the same dimension are added together. This isn’t always the case, so caveat emptor.

We could also think of an nth rising power or nth falling power as an n-dimensional volume. If we do, then the dimensions in a Newton series cancel out, for example.

Dimensional analysis for factorials and related functions could make it easier to remember identities, and easier to spot errors. And it could suggest the correct form of a result before the details of the result are filled in. In other words, artificial dimensional analysis can provide the same benefits of physically meaningful dimensional analysis.

Gamma function

For integers n, Γ(n) = (n − 1)!, and so we could assign dimension n − 1 to Γ(n), and more generally assign dimension z − 1 to Γ(z) for any complex z.

Now for some examples to show this isn’t as crazy as it sounds. For starters, take the identity Γ(z + 1) = z Γ(z). We imagine the left hand side is a z-dimensional volume, and the right hand side is the product of a term with unit length and a term that represents a (z − 1) dimensional volume, so the units check out.

Next let’s take something more complicated, the Legendre duplication formula:

$\Gamma (z)\Gamma \left(z+\frac{1}{2}\right)=2^{1-2z}\sqrt{\pi} \; \Gamma (2z)$
The left hand side has dimension (z − 1) + (z − ½) = 2z − ½. The right hand side has dimension 2z − 1, apparently contradicting our dimensional analysis scheme. But √π = Γ(½), and if we rewrite √π as Γ(½) in the equation above, both sides have dimension 2z − ½. The dimensions in other identities, like the reflection formula, also balance when you replace π with Γ(½)².

Hypergeometric functions

The hypergeometric function F(a, b; c; z) is defined by its power series representation. If we assign dimensions to the coefficients in the series as we’ve done here, then the numerator and denominator of each term have the same dimensions, and so the hypergeometric function should be dimensionless in our sense. This implies we should expect that identities for hypergeometric functions to be dimensionless as well. And indeed they are. I’ll give two examples.

First, let’s look at Gauss’ summation identity

$F (a,b;c;1)= \frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}$

provided the real part of c is greater than the real part of a + b. Notice that the numerator and denominator both have dimension 2c − a − b − 2.

Next, let’s look at Kummer’s identity

$F (a,b;1+a-b;-1)= \frac{\Gamma(1+a-b)\Gamma(1+\tfrac12a)}{\Gamma(1+a)\Gamma(1+\tfrac12a-b)}$

Both the numerator and denominator have dimension 3a/2 − b.

Finally, let’s look at a more complicated formula.

$\begin{align*} F(a, b; c; z) &= \frac{\Gamma(c) \Gamma(b-a)}{\Gamma(b) \Gamma(c-a)} \,(-z)^{-a\phantom{b}} F\left(a, 1-c+a; 1-b+a; \frac{1}{z}\right) \\ &+ \frac{\Gamma(c) \Gamma(a-b)}{\Gamma(a) \Gamma(c-b)} \,(-z)^{-b\phantom{a}} F\left(\,b, 1-c+b; 1-a+b; \,\frac{1}{z}\right) \\ \end{align*}$

In both the terms involving gamma functions, the dimensions of the numerator and denominator cancel out.

Series for the reciprocal of the gamma function

Posted on 24 December 2024 by John

Stirling’s asymptotic series for the gamma function is

$\Gamma(z) \sim (2\pi)^{1/2} z^{z - 1/2} e^{-z} \sum_{n=0}^\infty (-1)^n \frac{\gamma_n}{z^n}$

Now suppose you’d like to find an asymptotic series for the function 1/Γ(z).

Since the series for Γ has the form f(z) times an infinite sum, it would make sense to look for a series for 1/Γ of the form 1/f(z) times an infinite sum. The hard part would be finding the new infinite sum. In general the series for function and the series for its reciprocal look nothing alike.

Here’s where we have a pleasant surprise: the coefficients in the series for 1/Γ are exactly the same as the coefficients in the series for Γ, except the signs don’t alternate.

$\frac{1}{\Gamma(z)} \sim (2\pi)^{-1/2} z^{-z +1/2} e^{z} \sum_{n=0}^\infty \frac{\gamma_n}{z^n}$

Illustration

The following is not a proof, but it shows that the result is at least plausible.

Define Γ* to be Γ divided by the term in front of the infinite series:

$\Gamma^*(z) = (2\pi)^{-1/2} z^{-z +1/2} e^{z} \Gamma(z)$

Then the discussion above claims that Γ* and 1/Γ* have the same asymptotic series, except with alternating signs on the coefficients. So if we multiply the first few terms of the series for Γ* and 1/Γ* we expect to get something approximately equal to 1.

Now

$\Gamma^*(z) = 1 + \frac{1}{12z} + \frac{1}{288z^2} - \frac{139}{51840z^3} - \cdots$

and we claim

$\frac{1}{\Gamma^*(z)} = 1 - \frac{1}{12z} + \frac{1}{288z^2} + \frac{139}{51840z^3} - \cdots$

So if we multiply the terms up to third order we expect to get 1 and some terms involving powers of z in the denominator with exponent greater than 3. In fact the product equals

$1 + \frac{571}{1244160 z^4} -\frac{19321}{2687385600 z^6}$

which aligns with our expectations.

Simple error function approximation

Posted on 18 December 2024 by John

I recently ran across the fact that

$\int_0^x \exp(-t^2)\, dt \approx \sin(\sin(x))$

is a remarkably good approximation for −1 ≤ x ≤ 1.

Since the integral above defines the error function erf(x), modulo a constant, this says we have a good approximation for the error function

$\text{erf}(x) \approx \frac{2}{\sqrt{\pi}} \sin( \sin(x) )$

again provided −1 ≤ x ≤ 1.

The error function is closely related to the Gaussian integral, i.e. the normal probability distribution CDF Φ. The relation between erf and Φ is simple but error-prone. I wrote up some a page notes for myself a few years ago so I wouldn’t make a mistake again moving between these functions and their inverses.

Update: This post makes the connection to probability explicit.

You can derive the approximation by writing out the power series for exp(t), substituting −t² for t, and integrating term-by-term from 0 to x. You’ll see that the result is the same as the power series for sine until you get to the x⁵ term, so the error is on the order of x⁵. Here’s a plot of the error.

The error is extremely small near 0, which is what you’d expect since the error is on the order of x⁵.

Chebyshev polynomials as distorted cosines

Posted on 15 August 2024 by John

Forman Acton’s book Numerical Methods that Work describes Chebyschev polynomials as

cosine curves with a somewhat disturbed horizontal scale, but the vertical scale has not been touched.

The relation between Chebyshev polynomials and cosines is

T_n(cos θ) = cos(nθ).

Some sources take this as the definition of Chebyshev polynomials. Other sources define the polynomials differently and prove this equation as a theorem.

It follows that if we let x = cos θ then

T_n(x) = cos(n arccos x).

Now sin x = cos(π/2 − x) and for small x, sin x ≈ x. This means

arccos(x) ≈ π/2 − x

for x near 0, and so we should expect the approximation

T_n(x) ≈ cos(n(π/2 − x)).

to be accurate near the middle of the interval [−1, 1] though not at the ends. A couple plots show that this is the case.

Mote Chebyshev posts

Interpolating the gamma function

Posted on 8 August 2024 by John

Suppose you wanted to approximate Γ(10.3). You know it’s somewhere between Γ(10) = 9! and Γ(11) = 10!, and linear interpolation would give you

Γ(10.3) ≈ 0.7 × 9! + 0.3 × 10! = 1342656.

But the exact value is closer to 716430.69, and so our estimate is 53% too high. Not a very good approximation.

Now let’s try again, applying linear interpolation to the log of the gamma function. Our approximation is

log Γ(10.3) ≈ 0.7 × log 9! + 0.3 × log 10! = 13.4926

while the actual value is 13.4820, an error of about 0.08%. If we take exponentials to get an approximation of Γ(10.3), not log Γ(10.3), the error is larger, about 1%, but still much better than 53% error.

The gamma function grows very quickly, and so the log gamma function is usually easier to work with.

As a bonus, the Bohr–Mollerup theorem says that log gamma is a convex function. This tells us that not only does linear interpolation give an approximation, it gives us an upper bound.

The Bohr–Mollerup theorem essentially says that the gamma function is the only function that extends factorial from a function on the integers to a log-convex function on the real numbers. This isn’t quite true since it’s actually &Gamma(x + 1) that extends factorial. Showing the gamma function is unique is the hard part. In the preceding paragraph we used the easy direction of the theorem, saying that gamma is log-convex.

Special functions

Reproducing a hand-drawn plot

Multiple angle asymmetry

Proof

Naming polynomials

Solving trigonometric equations

Related posts

Posthumous Chebyshev Polynomials

Plots

Square roots

Dirichlet kernel

Related posts

Duplicating Hankel plot from A&S

Related posts

An unexpected triangle

Numerical (mis)application

Impossible triangles

Related posts

Dimensional analysis for gamma function values

Factorials

Gamma function

Hypergeometric functions

Series for the reciprocal of the gamma function

Illustration

Simple error function approximation

Chebyshev polynomials as distorted cosines

Mote Chebyshev posts

Interpolating the gamma function

Related posts