Friday, October 11, 2013

Numerical differentiation

1Numerical differentation¶

Recall that earlier we looked at the limit definition of derivatives:

$$\lim_{h \to 0} \frac{\partial f}{\partial x}$$

There is a range of $h$ for which you get a good approximation. However, as $h$ gets too small, then you get wildly unstable (oscillating) results due to truncation¹ errors, and as $h$ gets too large, the approximation slowly loses accuracy.

If we take a small but finite value of $h$, $\displaystyle \frac{\partial f}{\partial x}$ is called a finite difference. Then we can interpolate $f(x)$ at various points and calculate $\displaystyle \frac{d}{dx}$ for the interpolating function. Suppose that

$$f(x) = \sum_{k=0}^n f(x_k)L_k(x) + \frac{(x-x_0)(x-x_1)\cdots(x-x_n)}{(n+1)!} f^{(n+1)}(\xi)$$

Then, differentiating the above, we get:

$$\frac{df}{dx} = \sum_{k=0}^n f(x_n)L'_k(x) + \frac{d}{dx} \left [ \frac{(x-x_0)\cdots(x-x_n)}{(n+1)!} f^{(n+1)}(\xi)\right ]$$

Then, if we want $\displaystyle \frac{df}{dx} (x_j)$ (i.e., the derivative at one of the points):

[NEEDS VERIFICATION]

$$\frac{df}{dx}(x_j) = \sum_{k=0}^n f(x_j) L'_k(x_j) + \frac{f^{(n+1)}(\xi)}{(n+1)!} \tag{called the $(n+1)$ point formula}$$

Usually, we use either 3 or 5 points.

1.13-point formula¶

Consider the case where $n=2$, and so we have 3 equally spaced points, $x_0$, $x_1 = x_0 + h$, and $x_2 = x_0 + 2h = x_1 + h$.

Then, the interpolating function (Lagrange) looks like:

$$f(x) = \frac{(x-x_1)(x-x_2)}{(x_0-x_2)(x_0-x_2)} f(x_0) + \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1) + \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} f(x_2) + \,\text{error}$$

Now let's calculate the derivative:

$$f'(x) = \frac{(x-x_2) + (x-x_1)}{(x_0-x_1)(x_0-x_2)} f(x_0) + \frac{(x-x_2) + (x-x_0)}{(x_1-x_0)(x_1-x_2)} f(x_1) + \frac{(x-x_1) + (x-x_0)}{(x_2-x_0)(x_2-x_1)}f(x_2) + \,\text{error}$$

Let's say we want the derivative at $x_0$. Recall that the difference between any two adjacent points is some number $h$, due to equidistant points. Then:

$$\begin{align} f'(x_0) & = \frac{(-2h) + (-h)}{(-h)(-2h)} f(x_0) + \frac{(-2h) + 0}{(h)(-h)} f(x_1) + \frac{(-h) + 0}{(2h)(h)}f(x_2) \\ & = \frac{1}{h} \left [ -\frac{3}{2}f(x_0) + 2f(x_1) - \frac{1}{2} f(x_0) \right ] \end{align}$$

The error function, at $n=2$, is given by

$$\frac{f'''(\xi)}{6} \frac{d}{dx} (x-x_0)(x-x_1)(x-x_2) = \frac{f'''(\xi)}{6} [(x-x_1)(x-x_2) + (x-x_0)(x-x_2) + (x-x_2)(x-x_2)]$$

So at $x = x_0$, the error function is

$$\frac{f'''(\xi)}{6}2h^2$$

(He wrote down some other stuff at this point that I omitted .)

Think of $x_0$ and $x_2$ as endpoints for an interval, with $x_1$ as the midpoint.

Now let's do a change of variables in the definition of $f'(x_1)$, such that

• $x_1$, which is equivalent to $x_0 + h$, is now called $x_0$
• $x_0$ is now called $x_0 - h$ (referring to the new $x_0$)
• $x_2$ is now called $x_0 + h$ (referring to the new $x_0$)

Then:

$$f'(x_0) = \frac{1}{2h} [-f(x_0 - h) + f(x_0 + h)] = \frac{h^2}{6} f'''(\xi)$$

Next, we perform a change of variables to get $f'(x_2)$:

• $x_2$, which is equivalent to $x_0 + 2h$, is now called $x_0$
• $x_1$, which is equivalent to $x_0 + h$, is now called called $x_0 - h$
• $x_0$ is now called $x_0 - 2h$

Then:

$$f'(x_2) = \frac{1}{2h} [f(x_0 -2h) - 4f(x_0-h) + 3f(x_0) ] + \frac{h^2}{3} f'''(\xi)$$

So $f'(x_0)$ and $f'(x_2)$ are essentially equivalent (technically, symmetric: $h$ versus $-h$). This makes sense, since they're both endpoints. Thus there are only two independent approximations here: one for the endpoints (called the 3-point endpoint rule), and one for the midpoint (called the 3-point endpoint rule). For both rules, the error terms are proportional to $h^2$; notably, however, the error for the midpoint is that half that for the endpoints.

The analogous derivations for the 5-point rule are left as an exercise for the reader. Note that we use an odd number of points to ensure that there is a midpoint.

1.2Taylor series approach¶

Here's an alternate approach using Taylor series. First, we note the following equations:

$$f(x_0 + h) = f(x_0) + f'(x_0)h + f''(x_0) \frac{h^2}{2} + f'''(x_0)\frac{h^3}{6} + \cdots \tag{1}$$

$$f(x_0 - h) = f(x_0) - f'(x_0)h + f''(x_0)\frac{h^2}{2} - f'''(x_0)\frac{h^3}{6} + \cdots \tag{2}$$

Let's solve for the derivative by rewriting $(1)$:

$$f'(x_0) = \underbrace{\frac{f(x_0 + h) - f(x_0)}{h}}_{\text{finite-difference approximation}} - \underbrace{\frac{f''(x_0)h^2}{2} + \cdots}_{\text{error}}$$

This scheme is called forward difference, or forward Euler. Similarly, for (2), we can write:

$$f'(x_0) = \frac{f(x_0) -f(x_0 - h)}{h} + \frac{f''(x_0)h^2}{2} + \cdots$$

This scheme is called backward difference of backward Euler.

Now consider $(1) - (2)$. All the odd derivatives are doubled, and all the even derivatives vanish. So we get:

$$(1) - (2) = f(x_0 + h) - f(x_0 - h) = 2f'(x_0)h + \frac{f'''(x_0)h^3}{3} + \cdots$$

Thus we can solve for $f'(x_0)$:

$$f'(x_0) = \frac{f(x_0 + h) - f(x_0 - h)}{2h} - \frac{f'''(x_0)h^2}{6} + \cdots$$

This results in a much better scheme, since here the error is proportional to $h^2$. Thus, convergence is faster. This method is called centered difference, or, "leapfrog" (to get the derivative at one point, you use the points on either side). Incidentally, this is exactly the same equation as the midpoint rule that we derived using Lagrange interpolation above. However, this method is easier to derive.

1.2.1Notes on roundoff errors¶

Note that at small values of $h$, we get roundoff errors. Say we have an analytic function $f(x)$, our approximation for it is $\tilde{f}(x)$, and our rounding error is given by $e(x) = f(x) - \tilde{f}(x)$. Then, $f(x_0 \pm h) = \tilde{f}(x_0 \pm h) + e(x_0 \pm h)$. Using centered difference, we get:

$$f'(x_0) = \frac{1}{2h}(f(x_0 \pm h) - f(x_0 - h)) - \frac{h^2f'''(\xi)}{6} \tag{no idea what's going on here}$$

This is derived by replacing $f$ with its definition in terms of $\tilde{f}$ and $e$?????

So the error in the first derivative is given by (something).

Suppose that $|e(x_0+h) - e(x_0 - h)|$ is bounded above by $\epsilon$, and $f'''(\xi)$ is bounded above by $M$. Then:

$$|\text{error}| \leq \underbrace{\left | \frac{\epsilon}{2h} \right |}_{\text{round-off error}} + \underbrace{\left | \frac{h^2M}{6} \right |}_{\text{truncation error}}$$

Round-off error dominates when $h$ is small, and truncation error dominates when $h$ is large.