Finite Difference Approximations

Columbia University
APAM 1601

Introduction

Finite difference approximations to derivatives are the basic tools to the numerical solution to the dynamical equations of physics. As you know, by the very nature of these approximations, numerical "solutions" to the continuous motions we observe in the world can never be exact solutions. By choosing the correct method for making these approximations, we try to minimize the errors between the real solution and the numerical solution. Whenever possible, we try to show that (with a arbitrary precise computer) our method converges to the "real" solution in the limit of vanishing step size.

In this notebook, we'll derive low-order finite difference approximations to derivatives. 
We use Mathematica to perform simple algebra (and this is not the usual activity of computational physicists!) In order to understand how these derivations are performed, you need to understand how to use the following Mathematica commands: Series[…] and Solve[…]. (Notice also the use of "post-fix" commands. You will find this form useful since Mathematica is an interpretive language, and you can more easily apply additional commands to an input statement at the end than at the beginning.)

?Series

?Solve

First Derivative

In this section, we examine approximations to the first derivative of the function f[x]. We assume that we know the function only at fixed intervals, δ, or x ∈ δ, 2δ, 3δ,…

This section introduces "approximation" errors. These are errors due to the mathematical algorithms used by the scientist to evaluate functions.

Forward Difference

The first-order forward difference approximates the derivative using f[x0] and f[x0 + δ]. (A backward difference approximation uses the values, f[x0] and f[x0 - δ].)

Series[f[x], {x, x0, 4}]

t1 = f[x0+δ] == (Series[f[x], {x, x0, 4}] // Normal)/. (x - x0) -> δ

(f'[x0] /. Solve[t1, f'[x0]]) // First // Collect[#,δ, Simplify]&

The second-order approximation also uses f[x0 + 2δ].

t2 = f[x0+2δ] == (Series[f[x], {x, x0, 4}] // Normal)/. (x - x0) -> 2δ

(f'[x0] /. Solve[{t1, t2}, f'[x0], f''[x0]]) // First // Collect[#,δ, Simplify]&

Central Difference

The central difference approximates the derivative using f[x0 + δ] and f[x0 - δ]. The truncation error is of order (δ^3). Therefore, a central difference is very accurate while requiring little numerical computation.

t3 = f[x0-δ] == (Series[f[x], {x, x0, 4}] // Normal)/. (x - x0) -> -δ

(f'[x0] /. Solve[{t1, t3}, f'[x0], f''[x0]]  ) // First//Collect[#,δ, Simplify]&

Higher-Orders

In some cases, higher-order approximations to derivatives are needed. Mathematica can derive these approximations easily. For example, below is the next-order central difference approximation to the first derivative.

t4 = f[x0-2δ] == (Series[f[x], {x, x0, 4}] // Normal)/. (x - x0) -> -2δ

(f'[x0] /. Solve[{t1, t2, t3, t4}, f'[x0], {f''[x0],f'''[x0], f''''[x0]}]  ) // First//Collect[#,δ, Simplify]&

which is accurate to order (δ^4).

Second Derivative

Here, we find finite difference approximations to the second derivative.

Forward Difference

(f''[x0] /. Solve[{t1, t2}, f''[x0], f'[x0]]) // First // Collect[#,δ, Simplify]&

Central Difference

(f''[x0] /. Solve[{t1, t3}, f''[x0], f'[x0]]  ) // First// Collect[#,δ, Simplify]&

Subtractive Cancellation Example (from 1-Numbers.nb)

Another practical and useful example of a major source of numerical error in computational physics is subtractive cancellation (a type of "round-off" error.) When two numbers are subtracted from each other, the result is less precise as the original two numbers! This problem becomes severe when the difference between the two numbers become small.

We will illustrate this by computing a first-order and second-order approximation to a derivative and compare this approximation with the "exact" value.

A first-order approximation to the derivative of the function f[x] is
	(∂_x) f ~ ((f[x + h] - f[x])/h)
(It is called "first-order" because it's accuracy improves linearly with decreasing h.)

A second-order approximation to the derivative of the function f[x] is
	(∂_x) f ~ ((f[x + h] - f[x-h])/(2h))
(It is called "second-order" because it's accuracy improves linearly with decreasing (h^2). This is called a "central difference" because the difference of f[x] on both sides of x are evaluated.)

f[x_] := x^3;

First Order Derivative

error1[h_,x_:1] := Abs[3 x^2 - (f[x+h] - f[x])/h]

err1List = Table[{10.0^h,error1[10.0^h]},{h,-14,-1,0.5}];

err1List // Shallow

ListLogLogPlot[err1List,PlotStyle->PointSize[0.02`]]

Often, I set Graphics options to make the plot easier to read.

myPlot[items__]:=ListLogLogPlot[items,Frame->True,GridLines->Automatic,PlotStyle->{PointSize[0.02`],RGBColor[1,0,0]},BaseStyle->{FontWeight->"Bold",FontSize->12}]

p1 = myPlot[err1List,PlotLabel->"1st Order Difference Approximation"];

Second Order Derivative

error2[h_,x_:1] := Abs[3 x^2 - (f[x+h] - f[x-h])/(2h)]

err2List = Table[{10.0^h,error2[10.0^h]},{h,-14,-1,0.5}];

p2 = myPlot[err2List,PlotLabel->"2nd Order Difference Approximation"];

Show[p2,p1,PlotLabel->"Difference Error"]

Arbitrary-Precision

From Section 3.1.4:
  When you enter an approximate real number,Mathematica to decide whether to treat it as machine precision or arbitrary precision.Unless you specify otherwise,then if you give less than$MachinePrecisiondigits,Mathematica will treat the number as machine precision,and if you give more digits,it will treat the number as arbitrary precision.

10.0^-16

10.000000000000000000000000000000^-16

10.0`32^-16

10.0`200^-16

err3List = Table[{10.0^h,error2[10.0`32^h]},{h,-14,-1, 0.5`32}];

err3List // Shallow

p3 = myPlot[err3List,PlotLabel->"2nd Order Arbitrary-Precision"];

Show[p3,p2,p1,PlotLabel->"Finite  Difference  Error"]

Fixed Grid Size Example

In this example, we compare the central difference approximation to the first derivative of Sin[x]. We will use the built-in function called ListCorrelate[…] to easily perform the finite-difference over the entire list (think "vector" or "array") of data. As in many numerical problems, we are required to perform special calculations at the ends or boundaries of the data. For this example, we select our data to be periodic.

gridSize=100;
dx=((2 π)/gridSize);sinList=Table[N[Sin[x]],{x,0,2 π-dx,dx}];
dsinList=Table[N[Cos[x]],{x,0,2 π-dx,dx}];
n=Length[sinList];
ListPlot[sinList]

?ListCorrelate

ListCorrelate[{-1,0,1},{f[1],f[2],f[3],f[4],f[5]},2]/(2 d)

dsinListN=((ListCorrelate[{-1,0,1},sinList,2])/(2 dx));
ListPlot[dsinListN]

ListPlot[dsinListN-dsinList]

(dx^2 /6)// N (* dx^2 is the approximate size of trunction error *)

Summary

We have derived the formula for the finite difference approximations for first and second derivatives. This begins with a Taylor series expansion at x = x0 and used to relate the derivatives to nearby, known values of f[x].

We find that central difference approximations are significantly more accurate than equivalent forward (and backward) differences. The error associated with the approximation is called "truncation error" since we are truncating the Taylor series expansion. Reduced truncation error gives a more accurate solution to the differential equation and allows for longer integration times.