For Numerical Relativity, we need to
These can be built on some simple foundations.
The general concepts that underpin most numerical methods are
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
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}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$.
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$.
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}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.