Sage – Calculus VII

Rarefaction colliding with two shocks

The Burgers equation ${v_t+vv_x =0}$ is a great illustration of shock formation and propagation, with ${v(x,t)}$ representing velocity at time ${t}$ at position ${x}$ . Let’s consider it with piecewise constant initial condition

$\displaystyle u(x,0) = \begin{cases} 1,\quad & x<0 \\ 0,\quad & 0<x<1 \\ 2,\quad & 1<x<2 \\ 0,\quad & x>2 \end{cases}$

The equation, rewritten as ${(v_t,v_x)\cdot (1,v)=0}$ , states that the function ${v}$ stays constant on each line of the form ${x= x_0 + v(x_0,0)t}$ . Since we know ${v(x_0,0)}$ , we can go ahead and plot these lines (characteristics). The vertical axis is ${t}$ , horizontal ${x}$ .

burger0

Clearly, this picture isn’t complete. The gap near ${1}$ is due to the jump of velocity from ${0}$ to ${2}$ : the particles in front move faster than those in back. This creates a rarefaction wave: the gap gets filled by characteristics emanating from ${(1,0)}$ . They have different slopes, and so the velocity ${v}$ also varies within the rarefaction: it is ${v(x,t)=(x-1)/t}$ , namely the velocity with which one has to travel from ${(1,0)}$ to ${(x,t)}$ .

burger1

The intersecting families of characteristics indicate a shock wave. Its propagation speed is the average of two values of ${v}$ to the left and to the right. Let’s draw the shock waves in red.

burger2

Characteristics terminate when they run into shock. Truncating the constant-speed characteristics clears up the picture:

burger3

This is already accurate up to the time ${t=1}$ . But after that we encounter a complication: the shock wave to the right separates velocity fields of varying intensity, due to rarefaction to the left of the shock. Its propagation is now described by the ODE

$\displaystyle \frac{dx}{dt} = \frac12 \left( \frac{x-1}{t} + 0 \right) = \frac{x-1}{2t}$

which can be solved easily: ${x(t) = 2\sqrt{t} +1}$ .

Similarly, the shock on the left catches up with rarefaction at ${t=2}$ , and its position after that is found from the ODE

$\displaystyle \frac{dx}{dt} = \frac12 \left( 1 + \frac{x-1}{t} \right) = \frac{x-1+t}{2t}$

Namely, ${x(t) = t-\sqrt{2t}+1}$ for ${t>2}$ . Let’s plot:

burger4

The space-time trajectories of both shocks became curved. Although initially, the shock on the right moved twice as fast as the shock on the left, the rarefaction wave pulls them toward each other. What is this leading to? Yet another collision.

The final collision occurs at the time ${t=6+4\sqrt{2} \approx 11.5}$ when the two shock waves meet. At this moment, rarefaction ceases to exist. The single shock wave forms, separating two constant velocity fields (velocities ${1}$ and ${0}$ ). It propagates to the right at constant speed ${1/2}$ . Here is the complete picture of the process:

burger5

I don’t know where the holes in the spacetime came from; they aren’t predicted by the Burgers equation.

Oscillatory explosion

Nonlinear differential equations don’t necessarily have globally defined solutions, even if the equation is autonomous (no time is involved explicitly). The simplest example is ${y'=y^2}$ with initial condition ${y(0)=1}$ : the solution is ${y(t) = 1/(1-t )}$ , which ceases to exist at ${t=1}$ , escaping to infinity.

This kind of behavior can often be demonstrated without solving the ODE. Consider ${y''=y^3}$ with ${y(0)=1}$ , ${y'(0)=0}$ . I’ve no idea what the explicit solution is, but it’s clear that the quantity ${E(t) = \frac12 (y')^2 - \frac14 y^4}$ remains constant: indeed, ${E'(t) = y''y' - y^3 y' = 0}$ . (Here, ${E}$ is analogous to the sum of kinetic and potential energy in physics.)

The solution ${y}$ is increasing and convex. Let ${t_n}$ be such that ${y(t_n)=n}$ , e.g., ${t_1=0}$ . The invariance of ${E}$ yields ${y'(t_n) = 2^{-1/2} (n^2 - 1)}$ . By the mean value theorem, ${t_{n+1}-t_n \le \sqrt{2}(n^2-1)^{-1}}$ for ${n=2,3,\dots}$ . Since the series on the right converges, the limit ${T = \lim_{n\rightarrow\infty }t_n}$ is finite; this is the blow-up time.

But no blow-up occurs for the equation ${y''=-y^3}$ , where the nonlinear term pushes back toward equilibrium. Indeed, the invariant energy is now ${E= \frac12 (y')^2 + \frac14 y^4}$ , which implies that ${y}$ and ${y'}$ stay bounded. In the phase space ${(y,y')}$ the solution with initial values ${y(0)=1}$ , ${y'(0)=0}$ traces out this curve:

Closed trajectory in the phase space (y,y') — Closed trajectory in the phase space (y,y’)

Accordingly, the solution is a periodic function of ${t}$ (although this is not a trigonometric function):

Everything so far has been pretty straightforward. But here is a stranger animal: ${y'''=-y^3}$ . As in the previous example, nonlinear term pushes toward equilibrium. Using initial conditions ${y(0)=1}$ , ${y'(0)=y''(0)=0}$ , I get this numerical solution up to time ${T=5}$ :

As in the previous example, ${y}$ oscillates. But the speed and amplitude of oscillation are increasing.

Rapidly increasing:

In the plane ${(y,y')}$ the solution spirals out:

The plots make it clear that the solution ceases to exist in finite time, but I don’t have a proof. The issue is that the function ${y}$ does not escape to infinity in just one direction. Indeed, if ${y(t_0)>0}$ , then regardless of the values of ${y'}$ and ${y''}$ at ${t_0}$ , the negative third derivative will eventually make the function ${y}$ decrease, and so ${y}$ will attain the value ${0}$ at some ${t_1>t_0}$ . After that, the third derivative is positive, guaranteeing the existence of time ${t_2>t_1}$ when ${y}$ returns to ${0}$ again. I haven’t succeeded in proving that the limit ${\lim t_n}$ is finite, giving the time when oscillatory explosion occurs.

The plots were made in SageMathCloud:

var('t y y1 y2')
P = desolve_system_rk4([y1,y2,-y^3],[y,y1,y2], ics=[0,1,0,0],ivar=t,end_points=5,step=0.01)
Q=[ [i,j] for i,j,k,l in P]
list_plot(Q, plotjoined=True)

Implicit Möbius strip

Reading the Wolfram MathWorld™ article on the Möbius strip I came across an equation that was news to me.

After the familiar parametric form

$\displaystyle x=(1+s\cos(t/2))\cos t,\quad y=(1+s\cos(t/2))\sin t, \quad z = s\sin(t/2) \qquad (1)$

it is said that “In this parametrization, the Möbius strip is therefore a cubic surface with equation…”

$\displaystyle y(x^2+y^2+z^2-1)-2z(x^2+y^2+x) =0 \qquad\qquad\qquad(2)$

That’s a neat cubic equation for sure. But… wouldn’t the sign of the left hand side of (2) (call it $F$ ) distinguish the “positive” and “negative” sides of the surface, contrary to its non-orientability?

I went to SageMathCloud™ to check and got this plot of (2):

That’s… not exactly the Möbius strip as I know it. But it’s true that (1) implies (2): the cubic surface contains the Möbius strip. Here is the plot of (1) superimposed on the plot of (2).

The self-intersection arises where the “positive side” $F > 0$ suddenly transitions to negative $F< 0$ .

Trigonometric approximation and the Clenshaw-Curtis quadrature

Trying to approximate a generic continuous function on ${[-1,1]}$ with the Fourier trigonometric series of the form ${\sum_n (A_n\cos \pi nx+B_n\sin \pi nx)}$ is in general not very fruitful. Here’s such an approximation to ${f(x)=\exp(x)}$ , with the sum over ${n\le 4}$ :

It’s better to make a linear change of variable: consider ${f(2x-1)}$ on the interval ${[0,1]}$ , and use the formula for the cosine series. This results in ${\exp(2x-1)}$ , which is approximated by the partial sum ${\sum_{n=0}^4 A_n\cos \pi nx}$ of its cosine Fourier series as follows.

Better approximation after a change of variable

But one can do much better with a different, nonlinear change of variable. Consider ${f(\cos x)}$ on the interval ${[0,\pi]}$ , and again use the formula for the cosine series. This results in ${\exp(\cos x)}$ , which is approximated by the partial sum ${\sum_{n=0}^4 A_n\cos nx}$ of its cosine Fourier series as follows.

Excellent approximation after nonlinear change of variable

Yes, I can’t see any difference either: the error is less than ${10^{-3}}$ .

The composition with cosine improves approximation because ${f(\cos t)}$ is naturally a periodic function, with no jumps or corners in its graph. Fourier series, which are periodic by nature, follow such functions more easily.

A practical implication of this approximation of ${f(\cos t)}$ is the Clenshaw-Curtis integration method. It can be expressed in one line:

$\displaystyle \int_{-1}^1 f(x)\,dx = \int_0^\pi f(\cos t)\sin t\,dt \approx \int_0^\pi \sum_{n=0}^N a_n \cos nt \sin t\,dt = \sum_{n=0}^N \frac{(1+(-1)^n) a_n}{1-n^2}$

The integral ${\int_{-1}^1 f(x)\,dx}$ is approximated by summing ${2a_{2k}/(1-4k^2)}$ , where ${a_{2k}}$ are even-numbered cosine coefficients of ${f(\cos t)}$ . In the example with ${f(x)=\exp(x)}$ using just three coefficients yields

$\displaystyle \frac{2a_0}{1}+\frac{2a_2}{-3}+\frac{2a_4}{-15} \approx 2.350405$

while the exact integral is ${\approx 2.350402}$ .

At first this doesn’t look practical at all: after all, the Fourier coefficients are themselves found by integration. But since ${f(\cos t)}$ is so close to a trigonometric polynomial, one can sample it at equally spaced points and apply the Fast Fourier transform to the result, quickly obtaining many coefficients at once. This is what the Clenshaw-Curtis quadrature does (at least in principle, the practical implementation may streamline these steps.)

Completely monotone imitation of 1/x

I wanted an example of a function ${f}$ that behaves mostly like ${1/x}$ (the product ${xf(x)}$ is bounded between two positive constants), but such that ${xf(x)}$ does not have a limit as ${x\rightarrow 0}$ .

The first thing that comes to mind is ${(2+\sin(1/x))/x}$ , but this function does not look very much like ${1/x}$ .

Then I tried ${f(x)=(2+\sin\log x)/x}$ , recalling an example from Linear Approximation and Differentiability. It worked well:

I can't believe it's not a hyperbola! — I can’t believe it’s not a hyperbola!

In fact, it worked much better than I expected. Not only if ${f'}$ of constant sign, but so are ${f''}$ and ${f'''}$ . Indeed,

$\displaystyle f'(x) = \frac{\cos \ln x - \sin \log x - 2}{x^2}$

is always negative,

$\displaystyle f''(x) = \frac{4 -3\cos \log x + \sin \log x}{x^3}$

is always positive,

$\displaystyle f'''(x) = \frac{10\cos \log x -12}{x^4}$

is always negative. The sign becomes less obvious with the fourth derivative,

$\displaystyle f^{(4)}(x) = \frac{48-40\cos\log x - 10\sin \cos \ln x}{x^5}$

because the triangle inequality isn’t conclusive now. But the amplitude of ${A\cos t+B\sin t}$ is ${\sqrt{A^2+B^2}}$ , and ${\sqrt{40^2+10^2}<48}$ .

So, it seems that ${f}$ is completely monotone, meaning that ${(-1)^n f^{(n)}(x)\ge 0}$ for all ${x>0}$ and for all ${n=0,1,2,\dots}$ . But we already saw that this sign pattern can break after many steps. So let’s check carefully.

Direct calculation yields the neat identity

$\displaystyle \left(\frac{1+a\cos \log x+b\sin\log x}{x^n}\right)' = -n\,\frac{1+(a-b/n)\cos\log x+(b+a/n) \sin\log x}{x^{n+1}}$

With its help, the process of differentiating the function ${f(x) = (1+a\cos \log x+b\sin\log x)/x}$ can be encoded as follows: ${a_1=a}$ , ${b_1=b}$ , then ${a_{n+1}=a_n-b_n/n}$ and ${b_{n+1} = b_n+a_n/n}$ . The presence of ${1/n}$ is disconcerting because the harmonic series diverges. But orthogonality helps: the added vector ${(-b_n/n, a_n/n)}$ is orthogonal to ${(a_n,b_n)}$ .

The above example, rewritten as ${f(x)=(1+\frac12\sin\log x)/x}$ , corresponds to starting with ${(a,b) = (0,1/2)}$ . I calculated and plotted ${10000}$ iterations: the points ${(a_n,b_n)}$ are joined by piecewise linear curve.

The total length of this curve is infinite, since the harmonic series diverges. The question is, does it stay within the unit disk? Let’s find out. By the above recursion,

$\displaystyle a_{n+1}^2 + b_{n+1}^2 = \left(1+\frac{1}{n^2}\right) (a_n^2+b_n^2)$

Hence, the squared magnitude of ${(a_n,b_n)}$ will always be less than

$\displaystyle \frac14 \prod_{n=1}^\infty \left(1+\frac{1}{n^2}\right)$

with ${1/4}$ being ${a^2+b^2}$ . The infinite product evaluates to ${\frac{\sinh \pi}{\pi}a\approx 3.7}$ (explained here), and thus the polygonal spiral stays within the disk of radius ${\frac12 \sqrt{\frac{\sinh \pi}{\pi}}\approx 0.96}$ . In conclusion,

$\displaystyle (-1)^{n} \left(\frac{1+(1/2)\sin\log x}{x}\right)^{(n)} = n!\,\frac{1+a_{n+1}\cos\log x+b_{n+1} \sin\log x}{x^{n+1}}$

where the trigonometric function ${a_{n+1}\cos\log x+b_{n+1} \sin\log x}$ has amplitude strictly less than ${1}$ . Since the expression on the right is positive, ${f}$ is completely monotone.

The plot was generated in Sage using the code below.

a,b,c,d = var('a b c d')
a = 0
b = 1/2
l = [(a,b)]
for k in range(1,10000):
    c = a-b/k
    d = b+a/k
    l.append((c,d))
    a = c
    b = d
show(line(l),aspect_ratio=1)

Linear approximation and differentiability

If a function ${f\colon \mathbb R\rightarrow \mathbb R}$ is differentiable at ${a\in \mathbb R}$ , then it admits good linear approximation at small scales. Precisely: for every ${\epsilon>0}$ there is ${\delta>0}$ and a linear function ${\ell(x)}$ such that ${|f(x)-\ell(x)|<\epsilon \,\delta}$ for all ${|x|<\delta}$ . Having ${\delta}$ multiplied by ${\epsilon}$ means that the deviation from linearity is small compared to the (already small) scale ${\delta}$ on which the function is considered.

For example, this is a linear approximation to ${f(x)=e^x}$ near ${0}$ at scale ${\delta=0.1}$ .

Linear approximation to exponential function

As is done on this graph, we can always take ${\ell}$ to be the secant line to the graph of ${f}$ based on the endpoints of the interval of consideration. This is because if ${L}$ is another line for which ${|f(x)-L(x)|<\epsilon \,\delta}$ holds, then ${|\ell-L|\le \epsilon \,\delta}$ at the endpoints, and therefore on all of the interval (the function ${x\mapsto |\ell(x)-L(x)|}$ is convex).

Here is a non-differentiable function that obviously fails the linear approximation property at ${0}$ .

(By the way, this post is mostly about me trying out SageMathCloud.) A nice thing about ${f(x)=x\sin \log |x|}$ is self-similarity: ${f(rx)=rf(x)}$ with the similarity factor ${r=e^{2\pi}}$ . This implies that no matter how far we zoom in on the graph at ${x=0}$ , the graph will not get any closer to linear.

I like ${x\sin \log |x|}$ more than its famous, but not self-similar, cousin ${x\sin(1/x)}$ , pictured below.

Standard example from intro to real analysis

Interestingly, linear approximation property does not imply differentiability. The function ${f(x)=x\sin \sqrt{-\log|x|}}$ has this property at ${0}$ , but it lacks derivative there since ${f(x)/x}$ does not have a limit as ${x\rightarrow 0}$ . Here is how it looks.

Let’s look at the scale ${\delta=0.1}$

and compare to the scale ${\delta=0.001}$

Well, that was disappointing. Let’s use math instead. Fix ${\epsilon>0}$ and consider the function ${\phi(\delta)=\sqrt{-\log \delta}-\sqrt{-\log (\epsilon \delta)}}$ . Rewriting it as

$\displaystyle \frac{\log \epsilon}{\sqrt{-\log \delta}+\sqrt{-\log (\epsilon \delta)}}$

shows ${\phi(\delta)\rightarrow 0}$ as ${\delta\rightarrow 0}$ . Choose ${\delta}$ so that ${|\phi(\delta)|<\epsilon}$ and define ${\ell(x)=x\sqrt{-\log \delta}}$ . Then for ${\epsilon \,\delta\le |x|< \delta}$ we have ${|f(x)-\ell(x)|\le \epsilon |x|<\epsilon\,\delta}$ , and for ${|x|<\epsilon \delta}$ the trivial bound ${|f(x)-\ell(x)|\le |f(x)|+|\ell(x)|}$ suffices.

Thus, ${f}$ can be well approximated by linear functions near ${0}$ ; it’s just that the linear function has to depend on the scale on which approximation is made: its slope ${\sqrt{-\log \delta}}$ does not have a limit as ${\delta\to0}$ .

The linear approximation property does not become apparent until extremely small scales. Here is ${\delta = 10^{-30}}$ .

Using a paraboloid to cover points with a disk

Find the equation of tangent line to parabola ${y=x^2}$ … borrring calculus drill.

Okay. Draw two tangent lines to the parabola, then. Where do they intersect?

If the points of tangency are ${a}$ and ${b}$ , then the tangent lines are
${y=2a(x-a)+a^2}$ and ${y=2b(x-b)+b^2}$ . Equate and solve:

$\displaystyle 2a(x-a)+a^2 = 2b(x-b)+b^2 \implies x = \frac{a+b}{2}$

Neat! The ${x}$ -coordinate of the intersection point is midway between ${a}$ and ${b}$ .

What does the ${y}$ -coordinate of the intersection tell us? It simplifies to

$\displaystyle 2a(b-a)/2+a^2 = ab$

the geometric meaning of which is not immediately clear. But maybe we should look at the vertical distance from intersection to the parabola itself. That would be

$\displaystyle x^2 - y = \left(\frac{a+b}{2}\right)^2 -ab = \left(\frac{a-b}{2}\right)^2$

This is the square of the distance from the midpoint to ${a}$ and ${b}$ . In other words, the squared radius of the smallest “disk” covering the set ${\{a,b\}}$ .

Same happens in higher dimensions, where parabola is replaced with the paraboloid ${z=|\mathbf x|^2}$ , ${\mathbf x = (x_1,\dots x_n)}$ .

Indeed, the tangent planes at ${\mathbf a}$ and ${\mathbf b}$ are
${z=2\mathbf a\cdot (\mathbf x-\mathbf a)+|\mathbf a|^2}$ and ${z=2\mathbf b\cdot (\mathbf x-\mathbf b)+|\mathbf b|^2}$ . Equate and solve:

$\displaystyle 2(\mathbf a-\mathbf b)\cdot \mathbf x = |\mathbf a|^2-|\mathbf b|^2 \implies \left(\mathbf x-\frac{\mathbf a+\mathbf b}{2}\right)\cdot (\mathbf a-\mathbf b) =0$

So, ${\mathbf x}$ lies on the equidistant plane from ${\mathbf a}$ and ${\mathbf b}$ . And, as above,

$\displaystyle |\mathbf x|^2 -z = \left|\frac{\mathbf a-\mathbf b}{2}\right|^2$

is the square of the radius of smallest disk covering both ${\mathbf a}$ and ${\mathbf b}$ .

The above observations are useful for finding the smallest disk (or ball) covering given points. For simplicity, I stick to two dimensions: covering points on a plane with the smallest disk possible. The algorithm is:

Given points ${(x_i,y_i)}$ , ${i=1,\dots,n}$ , write down the equations of tangent planes to paraboloid ${z=x^2+y^2}$ . These are ${z=2(x_i x+y_i y)-(x_i^2+y_i^2)}$ .
Find the point ${(x,y,z)}$ that minimizes the vertical distance to paraboloid, that is ${x^2+y^2-z}$ , and lies (non-strictly) below all of these tangent planes.
The ${x,y}$ coordinates of this point is the center of the smallest disk covering the points. (Known as the Chebyshev center of the set). Also, ${\sqrt{x^2+y^2-z}}$ is the radius of this disk; known as the Chebyshev radius.

The advantage conferred by the paraboloid model is that at step 2 we are minimizing a quadratic function subject to linear constraints. Implementation in Sage:


points = [[1,3], [1.5,2], [3,2], [2,-1], [-1,0.5], [-1,1]] 
constraints = [lambda x, p=q: 2*x[0]*p[0]+2*x[1]*p[1]-p[0]^2-p[1]^2-x[2] for q in points]
target = lambda x: x[0]^2+x[1]^2-x[2]
m = minimize_constrained(target,constraints,[0,0,0]) 
circle((m[0],m[1]),sqrt(m[0]^2+m[1]^2-m[2]),color='red') + point(points)

Credit: this post is an expanded version of a comment by David Speyer on last year’s post Covering points with caps, where I considered the same problem on a sphere.