<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step3:  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
the solution of linear systems $A {\bf x} = {\bf b}$;
the solution of nonlinear root-finding problems ${\bf f} ( {\bf x} ) = {\bf 0}$;
the representation of a function or field $f(x)$ by discrete data $f_i$, by interpolation or other means;
the (discrete) Fast Fourier Transform;
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
Step4:  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$
Step5:  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
<ASSISTANT_TASK:>
Python Code:
def backward_differencing(f, x_i, dx):
    
    Backward differencing of f at x_i with grid spacing dx.
    
    f_i = f(x_i)
    f_i_minus_1 = f(x_i - dx)
    
    return (f_i - f_i_minus_1) / dx
def forward_differencing(f, x_i, dx):
    
    Forward differencing of f at x_i with grid spacing dx.
    
    f_i = f(x_i)
    f_i_plus_1 = f(x_i + dx)
    
    return (f_i_plus_1 - f_i) / dx
def central_differencing(f, x_i, dx):
    
    Second order central differencing of f at x_i with grid spacing dx.
    
    f_i = f(x_i)
    f_i_minus_1 = f(x_i - dx)
    f_i_plus_1 = f(x_i + dx)
    
    first_derivative = (f_i_plus_1 - f_i_minus_1) / (2.0 * dx)
    second_derivative = (f_i_minus_1 - 2.0 * f_i + f_i_plus_1) / (dx**2)
    
    return first_derivative, second_derivative
import numpy
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)))
Explanation: 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
the solution of linear systems $A {\bf x} = {\bf b}$;
the solution of nonlinear root-finding problems ${\bf f} ( {\bf x} ) = {\bf 0}$;
the representation of a function or field $f(x)$ by discrete data $f_i$, by interpolation or other means;
the (discrete) Fast Fourier Transform;
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
what do we need to understand about these methods before implementing or using them?
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.
End of explanation
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)))
Explanation: 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$:
End of explanation
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):
    bd_errors[i] = abs(backward_differencing(numpy.exp, 0.0, dx) - 1.0)
    fd_errors[i] = abs(forward_differencing(numpy.exp, 0.0, dx) - 1.0)
    cd1, cd2 = central_differencing(numpy.exp, 0.0, dx)
    cd1_errors[i] = abs(cd1 - 1.0)
    cd2_errors[i] = abs(cd2 - 1.0)
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()
Explanation: 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.
End of explanation
<END_TASK>