numerical – Calculus VII

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

Unreasonable effectiveness of the left endpoint rule

The left endpoint rule, and its twin right endpoint rule, are ugly ducklings of integration methods. The left endpoint rule is just the average of the values of the integrand over left endpoints of equal subintervals:

$\displaystyle \int_a^b f(x)\,dx \approx \frac{b-a}{n} \sum_{i=0}^{n-1} f(a+i/n)$

Here is its graphical representation with ${n=10}$ on the interval ${[-1,1]}$ : the sample points are marked with vertical lines, the length of each line representing the weight given to that point. Every point has the same weight, actually.

Primitive, ineffective, with error ${O(1/n)}$ for ${n}$ points used.

Simpson’s rule is more sophisticated, with error ${O(1/n^4)}$ . It uses weights of three sizes:

Gaussian quadrature uses specially designed (and difficult to compute) sample points and weights: more points toward the edges, larger weights in the middle.

Let’s compare these quadrature rules on the integral ${\int_{-1}^1 e^x \cos \pi x\,dx}$ , using ${10}$ points as above. Here is the function:

The exact value of the integral is ${\dfrac{e^{-1}-e}{1+\pi^2}}$ , about -0.216.
Simpson’s rule gets within 0.0007 of the exact value. Well done!
Gaussian quadrature gets within 0.000000000000003 of the exact value. Amazing!
And the lame left endpoint rule outputs… a positive number, getting even the sign wrong! This is ridiculous. The error is more than 0.22.

Let’s try another integral: ${\int_{-1}^1 \sqrt{4+\cos \pi x +\sin \pi x} \,dx}$ , again using ${10}$ points. The function looks like this:

The integral can be evaluated exactly… sort of. In terms of elliptic integrals. And preferably not by hand:

Maple output, "simplified" — Maple output, “simplified”

Simpson’s rule is within 0.00001 of the exact value, even better than the first time.
Gaussian quadrature is within 0.00000003, not as spectacular as in the first example.
And the stupid left endpoint rule is … accurate within 0.00000000000000005. What?

The integral of a smooth periodic function over its period amounts to integration over a circle. When translated to the circle, the elaborate placement of Gaussian sample points is… simply illogical. There is no reason to bunch them up at any particular point: there is nothing special about (-1,0) or any other point of the circle.

The only natural approach here is the simplest one: equally spaced points, equal weights. Left endpoint rule uses it and wins.

Unreasonable effectiveness of Runge-Kutta 2

I used to give this example of an initial-value problem whenever numerical solutions of ODE came up in class:

$\displaystyle y' = -y^2, \quad y(0)=1$

The exact solution is ${y(t)=1/(1+t)}$ . Using Euler’s method with step ${h=1}$ (admittedly large) we get, starting with ${y_0=1}$ :

$\displaystyle y_1 = y_0 + h (-y_0^2) = 0,\quad y_2 = y_1 + h (-y_1^2) = 0, \dots$

The numerical solution drops down to zero at once and stays there. Not good at all.

The trapezoidal version of RK2 (also known as improved Euler, or Heun’s method) treats the above value of ${y_1}$ merely as prediction ${\tilde y_1}$ , and computes the actual ${y_1}$ using the average of the slopes at ${y_0}$ and ${\tilde y_1}$ :

$\displaystyle y_1 = y_0 + \frac{h}{2} (-y_0^2 - \tilde y_1^2) = 1 + \frac12 (-1 - 0 ) = \frac12$

… which is right on target. Let’s do it again: the prediction ${\tilde y_2 = y_1 + h (-y_1^2) = 1/4}$ and correction

$\displaystyle y_2 = y_1 + \frac{h}{2} (-y_1^2 - \tilde y_2^2) = 1 + \frac12 \left(-\frac{1}{4}-\frac{1}{16} \right) = \frac{11}{32}$

… which is ${\approx 0.34}$ , compared to the exact solution ${1/3}$ . Nice!

I like the example, because all numbers are manageable with mental arithmetics for two steps. And it’s clear that ${11/32}$ is pretty close to ${1/3 = 11/33}$ .

But the fact that the relative error stays below ${3.7\%}$ throughout the computation despite the large step size ${h=1}$ seems… too good to be true. Consider that for the exact solution, the expression ${y^2}$ can differ by the factor of ${4}$ between two ${t}$ -values at distance ${1}$ .

The midpoint version of RK2, which generally performs about as well as the trapezoidal version, is nowhere as close in this example. Indeed, using the same data as above, we get

$\displaystyle y_1 = y_0 - h \left(\frac{y_0+\tilde y_1}{2}\right)^2 = \frac{3}{4}$

and so forth: the relative error reaches ${59\%}$ at the second step. That’s 16 times worse than the trapezoidal version.

What happens here has less to do with numerical analysis than with algebra of rational functions. Using ${h=1}$ in trapezoidal RK2, we are in effect iterating the function

$\displaystyle \varphi(y) = y + \frac12 (-y^2-(y-y^2)^2) = y-y^2+y^3-\frac12 y^4$

The exact solution would be obtained by iterating

$\displaystyle \Phi (y) = \frac{y}{1+y} = y-y^2+ y^3 - y^4 + \dots$

Two functions just happen to coincide at ${y=1}$ , which is our starting point here.

From there we get to ${0.5}$ , and on ${[0,0.5]}$ they are really close anyway, due to another coincidence: the truncation error ${|\varphi(y)-\Phi(y)|}$ is ${O(y^4)}$ instead of ${O(y^3)}$ as it is normally for second-order methods.

The midpoint method with ${h=1}$ amounts to iterating

$\displaystyle \psi(y) = y - \left(\frac{y+y-y^2}{2} \right)^2 = y-y^2+y^3-\frac14 y^4$

which is not substantially further away from ${\Phi}$ , but does not enjoy the lucky coincidence at ${y=1}$ .

The tables are turned with the initial value ${y_0=3}$ . The exact solution is ${y(t) = 3/(1+3t)}$ , which drops sharply from ${t=0}$ to ${t=1}$ ; its slope decreases by the factor of ${16}$ during one step. Still, midpoint-RK2 does a decent job with ${h=1}$ :

while trapezoidal RK2 outputs ${y_1=-19.5}$ , ${y_2=-80110}$ , ${y_3 = -2\cdot 10^{19}}$ and promptly overflows.

With a reasonable step size, like ${h=0.1}$ , normality is restored: both methods perform about equally well, with ${0.13\%}$ error for trapezoidal and ${0.21\%}$ for midpoint.

Calculating 2*2 with high precision

The definition of derivative,

$\displaystyle f'(x) = \lim_{h\rightarrow 0}\frac{f(x+h) - f(x)}{h} \ \ \ \ \ \ \ \ \ \ (1)$

is not such a great way to actually find the derivative numerically. Its symmetric version,

$\displaystyle f'(x) = \lim_{h\rightarrow 0}\frac{f(x+h) - f(x-h)}{2h} \ \ \ \ \ \ \ \ \ \ (2)$

performs much better in computations. For example, consider the derivative ${f(x)=e^x }$ at the point ${x=1}$ . We know that ${f'(1)=2.718281828\dots}$ . Numerically, with ${h=0.001}$ , we get

$\displaystyle \frac{f(1+h) - f(1)}{h} \approx \mathbf{2.71}9 \dots$

(error ${>10^{-3}}$ ) versus

$\displaystyle \frac{f(1+h) - f(1-h)}{2h} \approx \mathbf{2.71828}2 \dots$

(error ${<10^{-6}}$ ).

Considering this, why don’t we ditch (1) altogether and adopt (2) as the definition of derivative? Just say that by definition,

$\displaystyle f'(x) = \lim_{h\rightarrow 0}\frac{f(x+h)-f(x-h)}{2h}$

whenever the limit exists.

This expands the class of differentiable functions: for example, ${f(x)=|x|}$ becomes differentiable with ${f'(0)=0}$ . Which looks more like a feature than a bug: after all, ${f}$ has a minimum at ${0}$ , and the horizontal line through the minimum is the closest thing to the tangent line that it has.

Another example: the function

$\displaystyle f(x) = \begin{cases} x ,\quad & x\le 0 \\ 3x,\quad & x>0 \end{cases}$

has ${f'(0)=2}$ under this definition, because

$\displaystyle \lim_{h\rightarrow 0^+}\frac{f(x+h)-f(x-h)}{2h} = \lim_{h\rightarrow 0^+}\frac{3h-(-h)}{2h} = 2$

and

$\displaystyle \lim_{h\rightarrow 0^-}\frac{f(x+h)-f(x-h)}{2h} = \lim_{h\rightarrow 0^-}\frac{h-3(-h)}{2h} = 2$

This example also makes sense: since ${f(x)=|x|+2x}$ , getting ${f'(0)=0+2}$ is expected. In fact, with the new definition we still have basic derivative rules: if ${f,g}$ are differentiable, then ${f+g}$ , ${f-g}$ , ${fg}$ , ${f/g}$ are also differentiable (with the usual caveat about ${g\ne 0}$ ) and the familiar formulas hold.