Back to 498CQM Home
Back to 498CQM Calendar
Previous Lecture
Next Lecture

Phys 498CQM Lecture 1

Author: Erik Koch, Modified by R. Martin
HW 1 (due 1/31)

Outline

Numerical Differentiation

Numerical differentiation approximates derivatives by finite differences. All formulae assume that the function being differentiated can be approximated by an analytic function, usually a polynomial.

We will use shorthand notation in all of the following formulae:

fi = f(xi)
fi+1 = f(xi+1), etc.
Often we will take the points to be equally spaced on a grid with spacing h, so that xn = x0 + n h.

Taylor Expansion

The most common method for numerical differentiation is to derive formulae based on the Taylor expansion of f about xi.
fi = fi (1)
fi+1 = fi + h f'i + h2/2 f''i + h3/6 f'''i + . . . (2)
fi-1 = fi - h f'i + h2/2 f''i - h3/6 f'''i + . . . (3)
fi+2 = fi + 2h f'i + 2h2 f''i + 4h3/3 f'''i + . . . (4)
fi-2 = fi - 2h f'i + 2h2 f''i - 4h3/3 f'''i + . . . (5)
. . .
Combining Eqs. (1) and (2) gives a derivative formula to O(h).
Combining Eqs. (2) and (3) gives a derivative formula to O(h2).
Combining Eqs. (1), (2), (3), and (4) gives a derivative formula to O(h3).
Combining Eqs. (2), (3), (4), and (5) gives a derivative formula to O(h4).

These numerical differentiation formulae are:
f' = ( fi+1 - fi ) / h + O( h )
f' = ( fi+1 - fi-1 ) / 2h + O( h2 )
f' = ( -fi+2 + 6 fi+1 - 3 fi - 2 fi-1 ) / 6h + O( h3 )
f' = ( -fi+2 + 8 fi+1 - 8 fi-1 + fi-2 ) / 12h + O( h4 )
A few comments about the formulae:

Error Analysis

Since the error in the result is a power of the discretization length h, the error should always be plotted as on log-log plots. The slope should agree with the analytic expected slope!

Examples in class.

The Method of Undetermined Coefficients

When deriving differentation formulae, it is useful to have a systematic approach. The method described here is quite general and can be used on non-uniform grids. We will also use it in the next section to derive an integration formula.

This method is very simple. First make an ansatz for the form of the formula you want, and leave the coefficients as variables. For example,

f'i = ( a-1 fi-1 + a0 fi + a1 fi+1 + a2 fi+2 )

This example has four undetermined coefficents. Now choose four functions you want this formula to treat exactly. For our example we will choose f(x) = 1, f(x) = x - xi, f(x) = (x - xi)2, and f(x) = (x - xi)3.

Plugging these functions into our ansatz gives a linear system of equations,

0 = a-1 + a0 + a1 + a2
1 = - a-1 + a1 + 2 a2
0 = + a-1 + a1 + 4 a2
0 = - a-1 + a1 + 8 a2
The solution to these equations is a-1 = -1/6, a-1 = 1, a-1 = -1/2, and a-1 = -1/3.

Numerical Quadrature

The term "quadrature" refers to evaluating the value of an integral.

We can derive a formula for numerical quadrature by using the method of undetermined coefficients.

Suppose we want to numerically integrate f(x) over a finite interval, and use the function values at the end and the midpoint. If we let h be half the length of the interval, we can expect our integral to be approximated by :

ai-1 h fi-1 + ai h fi + ai+1 h fi+1

To determine the coefficients, we require that our formula be exact for the functions f(x) = 1, f(x) = x - xi, and f(x) = (x - xi)2. These fuctions integrate to 2h, 0, and 2/3 h, respectively.

The linear system of equations that determine the coefficients is

2 = ai-1 + ai + ai+1
0 = - ai-1 + ai+1
2/3 = + ai-1 + ai+1
The solution to these equations is ai-1 = 1/3, ai = 4/3, and ai+1 = 1/3. Thus, the integration formula (know as Simpon's formula) is

Integral f(x) = ( fi-1 + 4 fi + fi+1 ) h / 3 + O(h3)

The usual way to apply this formula is to break the integral into a sum of small integrals, and evaluate each small integral using Simpson's rule. Thus, for grid spacing h = L / n, where L is the range of integration and n is an even integer, Simpon's rule applied to the n / 2 small integrals yields the formula:

Integral f = h / 3 × ( f0 + 4 f1 + 2 f2 + 4 f3 + 2 f4 + . . . 4 fn-1 + fn )

Using the method of undetermined coefficients, we can derive elementary integration formulas.
Integration Range Integration Formula Elementary Error Composite Error
0 to h h / 2 × ( f0 + f1 ) O( h3 ) O( h2 )
0 to 2 h h / 3 × ( f0 + 4 f1 + f2 ) O( h5 ) O( h4 )
0 to 3 h 3 h / 8 × ( f0 + 3 f1 + 3 f2 + f3 ) O( h5 ) O( h4 )
0 to 3 h 2 h / 45 × ( 7 f0 + 32 f1 + 12 f2 + 12 f3 + 32 f4 ) O( h7 ) O( h6 )
Observations:

Low order methods ( i. e. Simpson's ) vs. higher order methods.

Error Analysis

Examples in class of log-log plots.

Richardson Extrapolation

Extrapolate the results for finite h to approximate the h = 0 result.

Romberg Integration (see Chapter 4.3 of Numerical Recipies) applies Richardson extrapolation to numerical quadrature. Let Fexact be the exact value of the integral, and let F(h) be the numerical approximation with integration step h. For composite integration with Simpson's rule, which is O( h2 ), the relation between F(h) and Fexact for small h should be:

F( h ) = Fexact + c h2,
where c is a constant for the particular integral. After evaluating the integral for several different values of h, the constant can be estimated, and the equation can be solved to give a much better estimate of Fexact.
Last Modified January 1
Email question/comments/corrections to rmartin@uiuc.edu