Finite difference methods

Questions

  • How do I use a finite difference method to calculate derivatives?
  • What is the Laplacian operator?

Objectives

  • Use the finite difference method to calculate the derivative of an unknown function
  • Express the Laplacian as a differential operator

Finite difference methods are a family of techniques used to calculate derivatives

Finite-difference methods are a class of numerical techniques for solving differential equations by approximating derivatives with finite differences. They are widely used for solving ordinary and partial differential equations, as they can convert equations that are unsolvable analytically into a set of linear equations that can be solved on a computer.

They rely on the idea of discretization: breaking a domain (for example, the space domain) into a finite number of discrete elements.

Finite Differences

A numerical derivative can be calculated using the forward, backward or central difference methods

Finite difference method

The standard definition of a derivative is

To calculate the derivative numerically we make \(h\) very small and calculate

This is the forward difference because it is measured in the forward direction from \(x\).

The backward difference is measured in the backward direction from \(x\):

and the central difference uses both the forwards and backwards directions around \(x\):

Let’s start with a simple example - let’s use the forward difference method to calculate the derivative of \(x^2\) at \(x=5\) with \(h=0.01\).

def x_squared(x):
    return x**2

def forward_difference(f_x, x, h):
    return (f_x(x+h) - f_x(x)) / h
    
forward_difference(x_squared,5, 0.01)
10.009999999999764

Note that there is a more concise method of defining short functions which is to use a Lambda statement. For example, instead of defining the x_squared function we could use a lambda statement:

x_2 = lambda x: x**2
forward_difference(x_2, 5, 0.01)
10.009999999999764

Or we can even pass it directly into the function:

forward_difference(lambda x: x**2, 5, 0.01)
10.009999999999764

We need to converge with respect to the step size \(h\)

Our expressions above are approximations as they are only exactly equal to the derivative when the step size \(h\) is zero. Whether using forwards, backwards or central differences it is important to converge with respect to a decreasing step size \(h\).

Note that in the next tutorial we will see that it is also possible to make \(h\) too small!

As we can see from the example in the image at the top of the page, the central difference is (in general) more accurate than the forward or backward differences. In fact, the error is order \(h\) for the forwards and backwards methods, and order \(h^2\) for the central difference.

Let’s test this idea using our simple \(x^2\) example that we started above:

exact_answer = 20

def calculate_x2_error(h):
    error = exact_answer - forward_difference(x_squared,5, h)
    print ("error for h={} is {}".format(h,round(error,10)))
calculate_x2_error(0.01)
calculate_x2_error(0.005)
calculate_x2_error(0.0025)
error for h=0.01 is 9.99
error for h=0.005 is 9.995
error for h=0.0025 is 9.9975

We can see that as the step size \(h\) is halved, the error halves.

Second-order derivatives can be calculated using finite differences

The second derivative is a derivative of a derivative, and so we can calculate it be applying the first derivative formulas twice. The resulting expression (after application of central differences) is:

Let’s test this out using the \(x^2\) example that we started above:

def second_order_forward_difference(f_x, x, h):
    return (f_x(x+h) - 2*f_x(x) + f_x(x-h)) / (h**2)
second_order_forward_difference(x_squared, 5, 0.01)
1.9999999999953388

The second derivative of \(x^2\) is 2, so the implementation appears correct.

Partial derivatives can be calculated using finite differences

The extension to partial derivatives is also relatively straight-forward. At this point we also consider a second dependent variable, \(y\).

Let’s consider another example, where we calculate the \(x\)-component of a force \(F\) in a potential energy \(U = x^2+y^2\), at \(x=5,y=10\). We know that the force and potential energy are calculated as follows:

\[F_x = \frac{\partial U}{\partial x}\]

def potential_energy(x,y):
    return x**2 + y**2

def partial_dfdx(f_xy, x, y, h):
    return (f_xy(x+(h/2),y) - f_xy(x-(h/2),y)) / h
partial_dfdx(potential_energy, 5, 10, 0.01)
10.000000000000853

Which is close to the analytic answer of 10.

The Laplacian operator corresponds to an average rate of change

The Laplacian operator \(\nabla^2\) is a very important differential operator in physics. We will see it later in the course, when studying partial differential equations. It is used to mathematically describe a variety of physical processes, including diffusion, gravitational potentials, wave propogation and fluid flow.

When applied to \(f\) and written in full for a two dimensional cartesian coordinate system with dependent variables \(x\) and \(y\) it takes the following form:

With equivalent expressions for a single dimension or extension to higher dimensions.

We can think of the laplacian as encoding an average rate of change. To develop an intuition for how the laplacian can be interpreted physically, we need to understand two related operators - div and curl. We will not explore these operators further in this lesson, but a related video is listed under external resources.

The Laplacian can be calculated using finite differences

By adding the two equations for second order partial derivatives (Equations 8 and 9), we find that the Laplacian in two dimensions is:

The expression above is known as a five-point stencil as it uses five points to calculate the Laplacian.

Keypoints

  • Finite difference methods are a family of techniques used to calculate derivatives
  • A numerical derivative can be calculated using the forward, backward or central finite difference
  • We need to converge with respect to the step size \(h\)
  • Second-order derivatives, partial derivatives and the Laplacian can also be calculated using finite differences
  • The Laplacian operator corresponds to an average rate of change

Test your understanding

Implementing the Laplacian

  1. Use a five-point stencil to calculate

\[\nabla^2\phi = \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} + \frac{\partial^2\phi}{\partial z^2}\]

for \(\phi = 6\cos(x)+7\sin(y)\) at \(x=\pi\) and \(y=\pi\), using a step-size of your choice.

  1. Compare to the exact answer.

Show answer

import math

def integrand(x,y):
    return 6*math.cos(x) + 7*math.sin(y)

def laplacian(f_xy, x, y, h):
    return (f_xy(x+h,y) + f_xy(x-h,y) + f_xy(x,y+h) + f_xy(x,y-h) - 4*(f_xy(x,y))) / (h**2)    
  
laplacian(integrand, math.pi, math.pi, 1E-2)
5.999950000159515
  1. This is within 1e-5 to the exact answer of 6.