Subsections

Numerical differentiation

To understand the basics of numerical differentiation, we merely need to recall some high school calculus. To keep stuff nice and simple here, we shall only look at the numerical approximations to first and second derivatives.

First derivative

Let's think about calculating the first derivative to start with. Figure 3.1 may be of some help in understanding the processes I am about to describe.

Figure 3.1: Diagram illustrating the ideas of numerical differentiation. To approximate the derivative of the function $ f$ at the point $ x$, one must either use a point to the left or right of $ x$ and calculate the slope of the line joining the two points, or use both points to the left and right of $ x$ and calculate the slope of the line joining the two points.
\begin{figure}\centerline{\includegraphics[width=100mm]{diffDiag}}\end{figure}

Imagine that we want to find the first derivative (i.e. slope) of some function $ f(x)$ at the point $ x$. If we know the analytical form of the function then we could calculate the derivative and use that. In practice in computational physics, we don't know the exact analytical form of the function we are attempting to find the derivative of, hence we have to do what any good physicist would do: have a guess (some physicists prefer to call this process approximation, but we all know the truth). How do we get this guess of the slope? Well, we can be guided by remembering how to calculate the derivative of a curve from first principles. Remember this from high school calculus?

$\displaystyle f'(x) = \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}.$ (3.1)

This equation says that if we choose a point a distance $ h$ away on the $ x$-axis from our point $ x$, we can get a good idea of the slope of the curve at $ f(x)$ by taking the slope of the line going through $ f(x+h)$ and $ f(x)$. The limit tells us that we can get a better and better idea just by letting $ h$ get really small. On a computer we are limited by just how small we can let $ h$ become, and so Equation (3.1) is only ever able to be an approximation (guess) of the actual derivative of $ f$ at $ x$. Just to let you know, there are more funky ways of improving the accuracy and/or precision of this guess, but this is the basic idea used in many cases. It is possible to show3.1 using a Taylor expansion about the point $ f(x+h)$ that Equation (3.1) has an error (in technical jargon) ``of the order of $ h$''. This statement we will write as $ O(h)$ and therefore, to be really explicit about the fact that we are only approximating the derivative, we will write the first derivative as:

$\displaystyle f'(x) = \frac{f(x+h) - f(x)}{h} + O(h).$ (3.2)

This technique of finding the derivative is quite naïve in the sense that we really aren't using a huge amount of brain power to work out the derivative. So, how can we do a better job? We can use a centred formula.

A centred formula uses the idea that the line joining the points $ f(x-h)$ and $ f(x+h)$ is also an approximation to the slope at the point $ f(x)$. The first derivative from first principles then becomes

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

Note how this equation is ``centred'' in $ x$. Even though this equation looks very similar to Equation (3.2), there is a big difference when one considers finite $ h$. Again, using a Taylor series expansion3.2 it can be shown that this centred form has error of order $ h^2$, hence we write Equation (3.3) as

$\displaystyle f'(x) = \frac{f(x+h) - f(x-h)}{2h} + O(h^2),$ (3.4)

which is really good as we will choose $ h$ to be nice and small, and so $ h^2$ will be really small. Just as a jargon-type note, if we have an equation that has truncation error / has error / is accurate to order $ h$, then we say that it is linear in $ h$. If the equation is accurate to order $ h^2$, we say that it is quadratic in $ h$.

Second derivative

The second derivative $ f''(x)$ of the function $ f(x)$ is just an extension of the concepts in obtaining the first derivative. By using a Taylor series expansion3.3 we can show that the centred second derivative is

$\displaystyle f''(x) = \frac{f(x+h) + f(x-h) - 2 f(x)}{h^2} + O(h^2),$ (3.5)

which again is quadratic in $ h$, and therefore nice and accurate.

Just knowing Equations (3.4) and (3.5) is sufficient understanding to be able to cope with the remaining sections, in which we discuss various forms of algorithm for solving PDEs.



Footnotes

... show3.1
see Numerical Methods in Physics by A. L. Garcia for details.
... expansion3.2
and again, this is discussed in more detail in Garcia.
... expansion3.3
guess whose book you should look at for details.
Paul Cochrane 2002-04-18