Numerical Methods

For Numerical Relativity, we need to

  • evolve the spacetime (hyperbolic PDEs with "smooth" fields);
  • evolve the matter (hyperbolic PDEs with discontinuous fields);
  • solve initial data (elliptic PDEs);
  • extract gravitational waves (interpolation and integration);
  • find and analyse horizons (interpolation, BVPs).

These can be built on some simple foundations.

The general concepts that underpin most numerical methods are

  1. the solution of linear systems $A {\bf x} = {\bf b}$;
  2. the solution of nonlinear root-finding problems ${\bf f} ( {\bf x} ) = {\bf 0}$;
  3. the representation of a function or field $f(x)$ by discrete data $f_i$, by interpolation or other means;
  4. the (discrete) Fast Fourier Transform;
  5. stochastic concepts and methods.

For Numerical Relativity, there has been little need (yet!) for stochastic methods, and the use of FFTs is mostly restricted to analysis. All of these points can be found in standard numerical packages and libraries: the question, however, is

  1. what do we need to understand about these methods before implementing or using them?
  2. when is it faster or better to implement our own version rather than using a library?

Finite differencing

As a first step we'll quickly cover finite differencing: the approximation of derivatives of a function $f$ when the only information about $f$ is its value at a set of points, or nodes, $\{x_i\}$, denoted $\{f_i\}$.

Here we have the "representation of a function" problem. We represent the function $f$ using a piecewise polynomial function $g$. This polynomial must interpolate $f$: that is, $g(x_i) \equiv f(x_i) = f_i$. We then approximate derivatives of $f$ by derivatives of $g$.

As simple examples, let's assume we know three points, $\{f_{i-1}, f_i, f_{i+1}\}$. Then we have the linear polynomial approximations

$$ g_{FD} = \frac{x - x_{i+1}}{x_i - x_{i+1}} f_i + \frac{x - x_{i}}{x_{i+1} - x_{i}} f_{i+1} $$

and

$$ g_{BD} = \frac{x - x_{i}}{x_{i-1} - x_{i}} f_{i-1} + \frac{x - x_{i-1}}{x_i - x_{i-1}} f_i $$

or the quadratic polynomial approximation

$$ g_{CD} = \frac{(x - x_{i})(x - x_{i+1})}{(x_{i-1} - x_{i})(x_{i-1} - x_{i+1})} f_{i-1} + \frac{(x - x_{i-1})(x - x_{i+1})}{(x_{i} - x_{i-1})(x_{i} - x_{i+1})} f_{i} + \frac{(x - x_{i-1})(x - x_{i})}{(x_{i+1} - x_{i-1})(x_{i+1} - x_{i})} f_{i+1}. $$

Note how this Lagrange form is built out of indicator polynomials that take the value $1$ at one node and vanish at all others.

By differentiating these polynomial interpolating functions we get approximations to the derivatives of $f$. Each approximation is different, with different errors.

We'll assume that the nodes are equally spaced, with grid spacing $\Delta x = x_{i+1} - x_i$. The approximations above give the standard forward difference

$$ \left. \frac{\partial g_{FD}}{\partial x} \right|_{x = x_i} \to \left. \frac{\partial f}{\partial x} \right|_{x = x_i} = \frac{1}{\Delta x} \left( f_{i+1} - f_i \right) + {\cal O} \left( \Delta x \right), $$

the standard backward difference

$$ \left. \frac{\partial g_{BD}}{\partial x} \right|_{x = x_i} \to \left. \frac{\partial f}{\partial x} \right|_{x = x_i} = \frac{1}{\Delta x} \left( f_{i} - f_{i-1} \right) + {\cal O} \left( \Delta x \right), $$

and the standard central difference approximations

\begin{align} \left. \frac{\partial g_{CD}}{\partial x} \right|_{x = x_i} & \to \left. \frac{\partial f}{\partial x} \right|_{x = x_i} \\ & = \frac{1}{2 \, \Delta x} \left( f_{i+1} - f_{i-1} \right) + {\cal O} \left( \Delta x^2 \right), \\ \left. \frac{\partial^2 g_{CD}}{\partial x^2} \right|_{x = x_i} & \to \left. \frac{\partial^2 f}{\partial x^2} \right|_{x = x_i} \\ & = \frac{1}{\left( \Delta x \right)^2} \left( f_{i-1} - 2 f_i + f_{i+1} \right) + {\cal O} \left( \Delta x^2 \right). \end{align}

Testing this in code

We'll use finite differencing repeatedly. To test our code we'll be testing the differencing. Let's check the above approximations applied to a simple function,

$$ f(x) = \exp \left[ x \right]. $$

All derivatives match the original function, which evaluated at $x=0$ gives $1$.

First we write the functions, then we test them.


In [ ]:
def backward_differencing(f, x_i, dx):
    """
    Backward differencing of f at x_i with grid spacing dx.
    """
    
    #

In [ ]:
def forward_differencing(f, x_i, dx):
    """
    Forward differencing of f at x_i with grid spacing dx.
    """
    
    #

In [ ]:
def central_differencing(f, x_i, dx):
    """
    Second order central differencing of f at x_i with grid spacing dx.
    """
    
    #

In [ ]:
import numpy

In [ ]:
bd = backward_differencing(numpy.exp, 0.0, dx=1.0)
fd = forward_differencing(numpy.exp, 0.0, dx=1.0)
cd1, cd2 = central_differencing(numpy.exp, 0.0, dx=1.0)

print("Backward difference should be 1, is {}, error {}".format(bd, abs(bd - 1.0)))
print("Forward difference should be 1, is {}, error {}".format(fd, abs(fd - 1.0)))
print("Central difference (1st derivative) should be 1, is {}, error {}".format(cd1, abs(cd1 - 1.0)))
print("Central difference (2nd derivative) should be 1, is {}, error {}".format(cd2, abs(cd2 - 1.0)))

The errors here are significant. What matters is how fast the errors reduce as we change the grid spacing. Try changing from $\Delta x = 1$ to $\Delta x = 0.1$:


In [ ]:
bd = backward_differencing(numpy.exp, 0.0, dx=0.1)
fd = forward_differencing(numpy.exp, 0.0, dx=0.1)
cd1, cd2 = central_differencing(numpy.exp, 0.0, dx=0.1)

print("Backward difference should be 1, is {}, error {}".format(bd, abs(bd - 1.0)))
print("Forward difference should be 1, is {}, error {}".format(fd, abs(fd - 1.0)))
print("Central difference (1st derivative) should be 1, is {}, error {}".format(cd1, abs(cd1 - 1.0)))
print("Central difference (2nd derivative) should be 1, is {}, error {}".format(cd2, abs(cd2 - 1.0)))

We see roughly the expected scaling, with forward and backward differencing errors reducing by roughly $10$, and central differencing errors reducing by roughly $10^2$.

Convergence

The feature that we always want to show is that the error $\cal E$ reduces with the grid spacing $\Delta x$. In particular, for most methods in Numerical Relativity, we expect a power law relationship:

$$ {\cal E} \propto \left( \Delta x \right)^p. $$

If we can measure the error (by knowing the exact solution) then we can measure the convergence rate $p$, by using

$$ \log \left( {\cal E} \right) = p \, \log \left( \Delta x \right) + \text{constant}. $$

This is the slope of the best-fit straight line through the plot of the error against the grid spacing, on a logarithmic scale.

If we do not know the exact solution (the usual case), we can use self convergence to do the same measurement.

We check this for our finite differencing above.


In [ ]:
from matplotlib import pyplot
%matplotlib notebook

dxs = numpy.logspace(-5, 0, 10)
bd_errors = numpy.zeros_like(dxs)
fd_errors = numpy.zeros_like(dxs)
cd1_errors = numpy.zeros_like(dxs)
cd2_errors = numpy.zeros_like(dxs)

for i, dx in enumerate(dxs):
    #

In [ ]:
pyplot.figure()
pyplot.loglog(dxs, bd_errors, 'kx', label='Backwards')
pyplot.loglog(dxs, fd_errors, 'b+', label='Forwards')
pyplot.loglog(dxs, cd1_errors, 'go', label='Central (1st)')
pyplot.loglog(dxs, cd2_errors, 'r^', label='Central (2nd)')
pyplot.loglog(dxs, dxs*(bd_errors[0]/dxs[0]), 'k-', label=r"$p=1$")
pyplot.loglog(dxs, dxs**2*(cd1_errors[0]/dxs[0]**2), 'k--', label=r"$p=2$")
pyplot.xlabel(r"$\Delta x$")
pyplot.ylabel("Error")
pyplot.legend(loc="lower right")
pyplot.show()

Forwards and backwards differencing are converging at first order ($p=1$). Central differencing is converging at second order ($p=2$) until floating point effects start causing problems at small $\Delta x$.

Extension exercises

Higher order

Show, either by Taylor expansion, or by constructing the interpolating polynomial, that the fourth order central differencing approximations are

\begin{align} \left. \frac{\partial f}{\partial x} \right|_{x = x_i} & = \frac{1}{12 \, \Delta x} \left( -f_{i+2} + 8 f_{i+1} - 8 f_{i-1} + f_{i-2} \right) + {\cal O} \left( \Delta x^4 \right), \\ \left. \frac{\partial^2 f}{\partial x^2} \right|_{x = x_i} & = \frac{1}{12 \left( \Delta x \right)^2} \left( -f_{i-2} + 16 f_{i-1} - 30 f_i + 16 f_{i+1} - f_{i+2} \right) + {\cal O} \left( \Delta x^4 \right). \end{align}
Measure the convergence rate

Using numpy.polyfit, directly measure the convergence rate for the algorithms above. Be careful to exclude points where finite differencing effects cause problems. Repeat the test for the fourth order formulas above.

Self convergence

By definition, the error ${\cal E}(\Delta x)$ is a function of the grid spacing, as is our numerical approximation of the thing we're trying to compute $F(\Delta x)$ (above $F$ was the derivative of $f$, evaluated at $0$). This gives

$$ {\cal E}(\Delta x) = F \left( \Delta x \right) - F \left( 0 \right) $$

or

$$ F \left( \Delta x \right) = F \left( 0 \right) + {\cal E}(\Delta x). $$

Of course, $F(0)$ is the exact solution we're trying to compute. However, by subtracting any two approximations we can eliminate the exact solution. Using the power law dependence

$$ {\cal E}(\Delta x) = C \left( \Delta x \right)^p $$

this gives

$$ F \left( 2 \Delta x \right) - F \left( \Delta x \right) = C \left( \Delta x \right)^p \left( 2^p - 1 \right). $$

We still do not know the value of the constant $C$. However, we can use three approximations to eliminate it:

$$ \frac{F \left( 4 \Delta x \right) - F \left( 2 \Delta x \right)}{F \left( 2 \Delta x \right) - F \left( \Delta x \right)} = \frac{\left( 4^p - 2^p \right)}{\left( 2^p - 1 \right)} = 2^p. $$

So the self-convergence rate is

$$ p = \log_2 \left| \frac{F \left( 4 \Delta x \right) - F \left( 2 \Delta x \right)}{F \left( 2 \Delta x \right) - F \left( \Delta x \right)} \right|. $$

Compute this self-convergence rate for all the cases above.