# Integrals and derivatives

One of the most basic but also most important application of computers in physics is the evaluation of integrals and derivatives. Numerical evaluation of integrals is a particularly crucial topic, because integrals occur widely in physics and, while some integrals can be done analytically in closed form, most cannot. They can, however, almost always be solved on a computer. Conversely, derivatives can almost always be performed analytically, if a close analytic expression of the expression is available. However, we may not always be interested in an analytic expression of the derivative, since the expression could be too complicated or because its evaluation is more involved than a numerical derivative. In this chapter we examine a number of different techniques for evaluating integrals and derivatives, as well as taking a brief look at the related operation of interpolation.

## Integration

In integration our objective is to find the value of a definite integral with fixed bounds from $x=a$ to $x=b$ $$ I(a,b)=\int_a^b f(x)dx.$$ Note that this is a much simpler problem than finding the antiderivative (or the indefinite integral) $F(y)=\int_0^y f(x)dx$. The simplest and most naive way to solve the definite integral is to discretize the sum and to assign a rectangle to each segment $$I(a,b)\approx\sum_{i=k}^N f(x_k)h.$$ To discretize the continuous variable x, we take the interval $[a,b]$ and divide it into $N$ equally sized portions $h=(b-a)/N$. The sequence of $x$-points then becomes $x_k=a+(k-\frac{1}{2})h$. 

Discretization is a core concept in computational science. An example of the simple rectangular integration is shown in panel a) of the following figure. 

We can now write a little python program that evaluates the sum for the integral. We do this for the function $f(x)=x^4-2x+1$ from $a=0$ to $b=2$. We know the correct answer for this integral and compare the result of our program to it $$\int_0^2 (x^4-2x+1)dx = \left[\frac{1}{5}x^5-x^2+x\right]_0^2=4.4 .$$


```python
def f(x):
    return x**4 - 2*x + 1

N = 10
a = 0.0
b = 2.0
h = (b-a)/N

s = 0.0
for k in range(1,N):
    s += f(a+(k-0.5)*h)

print(h*s)
```

    2.3003400000000007


With only 10 points, the value of the integral is not very good, because the constant value of the function that we assume in each interval does not reflect the curvature of the function $f(x)$ very well. We could improve the accuracy of our integration by using more points (i.e. by making the interval smaller and smaller). This strategy will eventually lead to success, but it might require a large number of points. Maybe we can do something smarter.


### Trapezoidal rule

Instead of replacing the function in each interval with a horizontal line, we could represent the function by a line that runs through the function values at the end point of each interval. Such sloped lines, would capture more rapidly varying functions better. The idea is shown in panel b) of the previous figure.

The two end-points of each segment are $a+(k-1)h$ and $a+kh$. Following the $\textit{trapezoidal rule}$ this gives the simple formula for the area of the segment $$A_k=\frac{1}{2}h\left[ f(a+(k-1)h)+f(a+k h)\right].$$ All we have to do now is sum up the segments again to obtain the integral $$\begin{align} I(a,b)\approx \sum_{k=1}^N A_k &=\frac{1}{2}h\sum_{k=1}^N \left[ f(a+(k-1)h)+f(a+kh)\right]  \\ &=h \left[ \frac{1}{2}f(a)+\frac{1}{2}f(b)+\sum_{k=1}^{N-1} f(a+kh) \right]\end{align}$$ 

We can easily modify the code for the rectangular integration to incorporate the trapezoidal rule:


```python
def f(x):
    return x**4 - 2*x + 1

N = 10
a = 0.0
b = 2.0
h = (b-a)/N

s = 0.5*f(a) + 0.5*f(b)
for k in range(1,N):
    s += f(a+k*h)

print(h*s)
```

    4.50656


We see that with still only 10 discretization points, we immediately achieve a better integration accuracy. With a value of 4.507, we are now $\sim$2\% off.

### Simpson's rule

The trapezoidal rule approximated the integrand with linear segments. To improve on the trapezoidal rule, we make an attempt to represent the function better. We do this by a powerlaw expansion. The simplest expansion beyond linear is quadratic. To fit a quadratic piece to our function requires 3 points, which implies that we now need two discretized segments for each fit. This is schematically shown in the figure. 



Denoting our 3 points $-h$, $0$ and $h$, we fit a quadratic function $Ax^2+Bx+C$ through these points. At the 3 points the function assumes the following values $$\begin{gather} f(-h)=Ah^2-Bh+C, & f(0)=C, & f(h)=Ah^2+Bh+C\end{gather}.$$ Solving this system of equations provides the desired expressions for the 3 coefficients $$\begin{gather} A=\frac{1}{h^2}\left[\frac{1}{2}f(-h)-f(0)+\frac{1}{2}f(h)\right], & B=\frac{1}{2h}\left[f(h)-f(-h)\right], & C=f(0)\end{gather}.$$ These expressions we insert into the integral of the quadratic approximation to the integral under the curve of the first two segments: $$\int_{-h}^{h} (Ax^2+Bx+C)dx=\frac{2}{3}Ah^3+2Ch=\frac{1}{3}\left[f(-h)+4f(0)+f(h)\right] .$$

Now we have to generalize this expression to incorporate also the remaining segments. We do this by sliding the procedure along in pairs of segments and summing up the results. In general, the three points of our first pair of segments are at $a$, $a+h$ and $a+2h$. The three points for the next pair of segments lie at $a+2h$, $a+3h$ and $a+4h$, and so forth. The approximate value of the integral then becomes $$\begin{align}I(a,b)\approx & \frac{h}{3}\left[ f(a)+4f(a+h)+f(a+2h)\right] \\ & +\frac{h}{3} \left[f(a+2h)+4f(a+3h)+f(a+4h) \right] + \ldots \\ & +\frac{h}{3} \left[f(a+(N-2)h+4f(a+(N-1)h)+f(b)\right] \end{align}$$ Note that the total number of slices must be even for this to work. Collecting terms, we now have $$\begin{align} I(a,b)\approx & \frac{h}{3} \left[f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+\ldots + f(b)\right] \\ &=\frac{h}{3} \left[f(a)+f(b)+4\!\!\!\!\sum_{\: \: k \: \textrm{odd} \\ 1\ldots N-1} \!\! f(a+kh)+2 \!\!\!\!\sum_{\: \: k \: \textrm{even} \\ 2\ldots N-2} \!\! f(a+kh) \right] \end{align} $$ This formula is called the $\textit{extended Simpson's rule}$.

The sums over odd and even values of $k$ can be conveniently accomplished in Python by using a for loop of the form "$\texttt{for k in the range(1,N,2)}$" for the odd terms and "$\texttt{for k in the range(2,N,2)}$" for the even terms. Alternatively, we could rewrite Simpson's rule as $$ I(a,b)\approx \frac{h}{3} \left[f(a)+f(b)+4\sum_{k=1}^{N/2} f(a+(2k-1)h)+2\sum_{k=1}^{N/2-1}f(a+2kh) \right] $$ and use an ordinary for loop.

With this insight, we modify our integration code once more to adapt it to the Simpson's rule:


```python
def f(x):
    return x**4 - 2*x + 1

N = 10
a = 0.0
b = 2.0
h = (b-a)/N

s = f(a) + f(b)
for k in range(1,N,2):
    s += 4*f(a+k*h)
    
for k in range(2,N,2):
    s += 2*f(a+k*h)
    
print(h*s/3.0)
```

    4.400426666666668


In summary, we have seen three increasingly more sophisticated integration algorithms. Here is an overview of their performance for 10 integration points: $$ \begin{array}{llrr} \text{Method} & & \text{Value} & \text{Error} \\\hline  \text{Flat} & : & 2.3003 & 52.20 \% \\ \text{Trapezoidal} & : & 4.5066 & 2.40\% \\ \text{Simpson} & : & 4.4004 & 0.01\% \\ \text{Exact} & : & 4.4000 \\  \end{array} $$

# Higher-order integration methods

In the previous section, we have increased the sophistication of approximating the integrand in each method. We went from constants (0th order) to lines (1st order (trapezoidal)) to quadratic curves (2nd order (Simpson's rule)). In principle, we could keep going to higher and higher orders. The general expression for the integral is $$\int_a^b f(x)dx\approx \sum_{k=1}^{N}w_kf(x_k),$$ where $x_k$ are the positions at which we evaluate the integrand and $w_k$ are a set of weights given by the integration method. Since we used homogenuous grids (i.e. fixed grid spacing), we can split off the grid spacing $h$ from the weights $w_k=c_k h$. The different integration methods then only differ in the set of coefficients that make up the integration weights. Different coefficient sets are summarized in the table below.  $$ \begin{array}{lll} \text{Order} & \text{Polynomial} & \text{Coefficients} \\\hline \text{0 (naive)} & \text{constant} & 1,1,1,\ldots,1 \\ \text{1 (trapezoidal rule)} & \text{straight line} & \frac{1}{2}, 1, 1, \ldots, 1, \frac{1}{2} \\ \text{2 (Simpson's rule)} & \text{quadratic} & \frac{1}{3},\frac{4}{3},\frac{2}{3},\frac{4}{3},\ldots,\frac{4}{3},\frac{1}{3} \\ \text{3} & \text{cubic} & \frac{3}{8},\frac{9}{8},\frac{9}{8},\frac{3}{4},\frac{9}{8},\frac{9}{8},\frac{3}{4},\ldots, \frac{9}{8},\frac{3}{8}\\ \text{4} & \text{quartic} & \frac{14}{45}, \frac{64}{45},\frac{8}{15},\frac{64}{45},\frac{28}{45},\frac{64}{45}, \frac{8}{15},\frac{64}{45}, \ldots, \frac{64}{45},\frac{14}{45} \end{array}$$

We can now view integration as a simple sum over integration $\textit{weights}$ and the integrand evaluated at certain discretized $\text{grid}$ points. The weights depend on the integration method and give the weighting of the integrand at a corresponding grid point in the sum. With this new perspective we can rethink our integration program and make it more modular. If we create a function that prefills the weights, we can compute integrals for different methods in the same code framework. Here is an example of our Simpson rule integration adapted to a sum over a grid with the corresponding weights.



```python
from numpy import empty,array
def f(x):
    return x**4 - 2*x + 1

# variable definitions
N = 10
a = 0.0
b = 2.0
h = (b-a)/N
s = 0.0
weights = empty(N+1,float)
grid = empty(N+1,float)

# grid specification (this step is in principle redundant, but we use it to illustrate the concept)
for k in range (0,N+1):  
    grid[k]=a+k*h

# weights specification (only this part would need to be changed for a new method)
for k in range(1,N,2):
    weights[k]=h*4.0/3.0
    
for k in range(2,N,2):
    weights[k]=h*2.0/3.0

weights[0]=h/3.0
weights[N]=h/3.0

# actual integration
for k in range(0,N+1):
    s += weights[k]*f(grid[k])
    
print(s)
```

    4.400426666666667


The program gives the same value as its orginal version. We can now, however, easily upgrade it by modifying the section that specifies the weights. In principle, we could even outsource the weight definitions into subroutines. 

# Gaussian quadrature

Equipped with our new perspective of integration as looping over a grid with weighted integration points, we now consider non-uniform grids. We will choose the position of the grid points such that they are optimal for the integration. The integration weights are then derived correspondingly.

For simplicity, we restrict the domain of integration to $[-1,1]$ and will later scale it back to $[a,b]$. Our objective is then to find the grid points $x_k$ and weights $w_k$ for the integral $$\int_{-1}^{1}f(x)dx\approx\sum_{k=1}^Nw_kf(x_k).$$ Let us assume that our function $f(x)$ is a polynomial in $x$ of degree $2N-1$. We can then use the properties of the Legendre polynomials for our integration method: 
 - The legendre polynomial $P_N(x)$ is orthogonal to every polynomial of lower degree, i.e. $ \int_{-1}^{1}x^k P_N(x)dx=0$ for all integers $k$ in the range $0 \leq k \lt N$. 
 - For all $N$, $P_N(x)$ has $N$ real roots that all lie in the interval from $-1$ to $1$.
 
If we now divide $f(x)$ by the Legendre polynomial $P_N(x)$, we get $$f(x)=q(x)P_N(x)+r(x),$$ where $q(x)$ and $r(x)$ are both polynomials of degree $N-1$ or less. This simplifies our original integral $$\int_{-1}^{1} f(x)dx=\int_{-1}^{1} q(x)P_N(x)dx+\int_{-1}^{1}r(x)dx=\int_{-1}^{1}r(x)dx, $$ because $P_N(x)$ is orthogonal to $q(x)$. This transformation has not gained us anything, since we do not know the function $r(x)$. However, we can make the same substitution in the summation over the grid and weights: $$ \sum_{k=1}^N w_kf(x_k)=\sum_{k=1}^Nq(x_k)P_N(x_k)+\sum_{k=1}^Nw_kr(x_k)$$
We now make use of the second property of the Legendre polynomials, namely, that $P_N(x)$ has $N$ many zeros between $-1$ and $1$. We take these zero positions (roots) as the grid points $x_k$ so that $P_N(x_k)=0$ for all $k$. This simplifies our sum to $$\sum_{k=1}^N w_k f(x_k) = \sum_{k=1}^N w_k r(x_k)=\int_{-1}^{1} r(x) dx = \int_{-1}^{1} f(x) dx .$$ The second equality holds only, because we assumed that $f(x)$ is a polynomial of degree $2N-1$, which makes $r(x)$ a polynomial of degree $N-1$ or less.

To proceed, we have to do the following:

 - Find the roots of $P_N(x)$: this is possible, but tedious. The derviation can be found e.g. in (Abramowitz & Stegun 1972, p. 887).
 
 - Given a set of $x_k$ grid points (corresponding to the roots of $P_N(x)$) compute the corresponding weights.
 
For the second point, we assume that we can find a single polynomial of degree $N-1$ to fit the function $f(x)$ or $r(x)$. For this, we use the _method of interpolating polynomials_. Our _interpolating polynomial_ is $$ \begin{align} \phi_k(x) & =\prod_{m=1\ldots N \\ \:\:\: m\neq k}\frac{(x-x_m)}{(x_k-x_m)} \\ & = \frac{(x-x_1)}{(x_k-x_1)} \times \ldots \times \frac{(x-x_{k-1})}{(x_k-x_{k-1})} \frac{(x-x_{k+1})}{(x_k-x_{k+1})} \times \ldots \times  \frac{(x-x_N)}{(x_k-x_N)} \end{align}$$ Note that the numerator contains one factor for each sample point except the point $x_k$. $\phi_k(x)$ is thus a polynomial of degree $N-1$. It has the property $$\phi_k(x_m)=\delta_{km}.$$ We use this property to define a new function $\Phi(x)$: $$\Phi(x)=\sum_{k=1}^Nf(x_k)\phi_k(x),$$ which is also a polynomial of degree $N-1$. By definition $\Phi(x)$ fits $f(x)$ exactly at $x_m$: $$\Phi(x_m)=\sum_{k=1}^Nf(x_k)\phi_k(x_m)=\sum_{k=1}^Nf(x_k)\delta_{km}=f(x_m).$$
 
We now insert $\Phi(x)$ into our integral $$\int_{-1}^{1}f(x)dx\approx\int_{-1}^{1}\Phi(x)dx=\int_{-1}^{1}\sum_{k=1}^Nf(x_k)\phi_k(x)dx=\sum_{k=1}^Nf(x_k)\int_{-1}^{1}\phi_k(x)dx=\sum_{k=1}^Nf(x_k)w_k.$$ In other words, for an arbitrary set of points $x_k$, the integration weights are given by $$w_k=\int_{-1}^{1}\phi_k(x)dx.$$ 

This is a general expression for finding integration weights for a given set of grid points. For the Gauss Legendre quadrature, we insert the roots of $P_N(x)$ as $x_k$ into the interpolating polynomial and then carry out the integral to find the corresponding weights $w_k$. For small $N$, analytic expressions have been derived. In general, however, roots and weights are calculated numerically. We have done this in the accompanying program gaussxw.py. The figure below shows the positions of the grid points and the values of the corresponding weights for $N=10$ and $N=100$. The grid points are not evenly spaced. Their density increases towards the edges of the integration interval. Concomitantly, the weights are lowest at the interval edges and rise in towards the middle of the interval, where the point density is lowest.   



The calculation of the grid points and weights in Gaussian quadrature takes a bit of effort. In return, Gaussian quadrature exactly integrates functions that are polynomials of degree $2N-1$ with only $N$ grid points. This is due to the properties of the Legendre polynomials and the clever grid choice. As an example, we will consider our test function $x^4-2x+1$. This time, we will integrate it with only $N=3$ grid points using Gaussian quadrature.

Before we perform the integration, we need to briefly consider how to change the integration interval from $[-1,1]$ to $[a,b]$. Since the area under a curve does not depend on where that curve is along the $x$ axis, the sample points can be freely slid up and down the $x$ axis _en masse_. If the desired domain is wider or narrower than the interval from $-1$ to $1$ then we also need to spread the points out or squeeze them together. The stretching or squeezing operation that accomplishes this is $$x_k'=\frac{1}{2}(b-a)x_k+\frac{1}{2}(b+a).$$ Similarly, the weights do not change, if we are simply sliding the sample points up and down the $x$ axis, but if the width of the integration domain changes then the value of the integral will increase or decrease by a corresponding factor. The weights have to be rescaled accordingly $$w_k'=\frac{1}{2}(b-a)w_k.$$

Here, finally, is then then the Python program that integrates $x^4-2x+1$ with only 3 grid points using Gaussian quadrature.


```python
from gaussxw import gaussxw

def f(x):
    return x**4 - 2*x + 1

N = 3
a = 0.0
b = 2.0

# Calculate the sample points and weights, then map them
# to the required integration domain
x,w = gaussxw(N)
xp = 0.5*(b-a)*x + 0.5*(b+a)
wp = 0.5*(b-a)*w

# Perform the integration
s = 0.0
for k in range(N):
    s += wp[k]*f(xp[k])

print(s)
```

    4.4000000000000075


# Choosing an integration method

Which integration method should you use in practice? There is no general answer, because the performance of the methods depends on the nature of the integrand. If in doubt, try several methods. Below is a brief table of advantages and disadvantes of each method. $$ \begin{array}{lll} \text{Method} & \text{Advantage} & \text{Disadvantage} \\\hline  \text{Trapezoidal} & \text{simple and versatile; works for pathological integrands and noisy data} & \text{limited accuracy} \\ \text{Simpson} & \text{simple with good accuracy} & \text{less suitable for pathological integrands and noisy data}  \\ \text{Gaussian} & \text{very good accuracy} & \text{less suitable for pathological integrands and noisy data} \\  \end{array} $$

# Integrals over infinite ranges

Often in physics we encounter integrals over infinite ranges, like $\int_0^\infty f(x)dx$. The techniques we have used so far will not work for such integrals, because we would need an infinite number of integration points. The solution to this problem is to change variables. For an integral over the range from $0$ to $\infty$ the standard change of variables is $$z=\frac{x}{1+x} \quad \text{or equivalently} \quad x=\frac{z}{1-z} .$$ With $dx=dz/(1-z)^2$ we obtain $$\int_0^\infty f(x)dx=\int_0^1 \frac{1}{(1-z)^2} f(\frac{z}{1-z})dz,$$ which can be done using any of the techniques covered earlier in the chapter.

For integrals over a range from a nonzero value $a$ to $\infty$, we have to make two changes of variables. First we change $y=x-a$, which shifts the start of the integration range to $0$, and then apply the previous substitution, but this time in $y$: $z=y/(1+y)$. Or we combine both changes into a single one: $$ z=\frac{x-a}{1+x-a} \quad \text{or equivalently} \quad x=\frac{z}{1-z}+a .$$ Again $dx=dz/(1-z)^2$, so that we end up with $$\int_a^\infty f(x)dx=\int_0^1 \frac{1}{(1-z)^2} f(\frac{z}{1-z}+a)dz .$$

Integrals from $-\infty$ to $a$ can be done the same way by substituting $z \rightarrow -z$ and integrals from $-\infty$ to $\infty$ can be broken up into to integrals from $-\infty$ to $0$ and from $0$ to $\infty$. Alternatively, we could use a single change of variables, such as $$x=\frac{z}{1-z^2}, \quad dx=\frac{1+z^2}{(1-z^2)^2}dz,$$ which would give $$\int_{-\infty}^\infty f(x)dx=\int_{-1}^1 \frac{1+z^2}{(1-z)^2} f(\frac{z}{1-z^2})dz.$$

As an example, we will calculate the value of the following integral using Gaussian quadrature: $$I=\int_0^\infty e^{-t^2}dt$$ We make the change of variables $z=t/(1+t)$ and the integral becomes $$I=\int_0^1 \frac{e^{-z^2/(1-z)^2}}{(1-z)^2}dz .$$ We modify our program for the function $x^4-2x+1$ and perform the integration with $N=50$ grid points.


```python
from gaussxw import gaussxwab
from math import exp

def f(z):
    return exp(-z**2/(1-z)**2)/(1-z)**2

N = 50
a = 0.0
b = 1.0
x,w = gaussxwab(N,a,b)
s = 0.0
for k in range(N):
    s += w[k]*f(x[k])
print(s)
```

    0.8862269254528349


The value of this integral is known analytically: $I=\frac{1}{2}\sqrt{\pi}=0.886226925453\ldots$. Again we see the impressive accuracy of the Gaussian quadrature method. With just 50 sample points, we have calculated an estimate of the integral that is correct to the limits of machine precision.

# Multiple integrals

Integrals over more than one variable are common in physics problems and can be tackled using generalizations of the methods we have already seen. Consider for instance the integral $$I=\int_0^1\int_0^1f(x,y) dx dy.$$ We can rewrite this by defining a function $F(y)$ thus $$F(y)=\int_0^1f(x,y)dx.$$ Then our integral becomes $$I=\int_0^1 F(y) dy.$$ We can thus do multiple integrals numerically by first integrating in one variable and then in the others. For instance, if we do the integrals by Gaussian quadrature with the same number $N$ of points for both $x$ and $y$ integrals, we have $$F(y)\approx \sum_{i=1}^N w_i f(x_i,y) \quad \text{and} \quad I\approx \sum_{j=1}^N w_j F(y_j) .$$ Substituting the first into the second sum not surprisingly gives us the following double sum for the integral $$ I\approx\sum_{i=1}^N\sum_{j=1}^N w_i w_j f(x_i,y_j). $$ This expression has a form similar to the standard integration formula for single integrals with a sum over values of the function $f(x,y)$ at a set of grid points, multiplied by appropriate weights. Now the points are distributed on a two dimensional grid. The approach can be generalized to arbitrary dimensions. The figure below shows the points on a two dimensional Gaussian quadradture grid.




# Derivatives

## Forward and backward differences

The standard definition of a derivative, the one you see in calculus books, is $$\frac{df}{dx}=\lim_{h\rightarrow 0}\frac{f(x+h)-f(x)}{h} .$$ The basic method for calculating numerical derivatives is precisely an implementation of this formula. We cannot take the limit $h \rightarrow 0$ in practice, but we can make $h$ very small and then calculate $$\frac{df}{dx} \approx \frac{f(x+h)-f(x)}{h} .$$ This is called the _forward difference_, because it is measured in the forward (i.e. positive) direction from the point of interest $x$. Analogously, we can define the _backward difference_ $$\frac{df}{dx} \approx \frac{f(x)-f(x-h)}{h} .$$ Both are shown in the figure below. They usually give the same value, provided $h$ is small enough and the function is not pathological.




## Errors

Forward and backward difference are usually not the most accurate. To understand this, we Taylor expand $f(x)$ around $x$: $$f(x+h)=f(x)+hf'(x)+\frac{1}{2}h^2f''(x)+\ldots ,$$ where $f'$ and $f''$ denote the first and second derivatives of $f$, respectively. Rearranging gives us $$f'(x)=\frac{f(x+h)-f(x)}{h}-\frac{1}{2}hf''(x)+\dots .$$ The first part is our numeric forward difference. The following terms we omit in our numeric difference. They contribute to the error and are proportional to $h$ in leading order ($\frac{1}{2}h|f''(x)|$). It would seem that if we make $h$ smaller, we would reduce the error. However, we are bound by machine precision and that affects the accuracy of the numeric difference, which means we have a lower bound for $h$. Since $f(x+h)$ and $f(x)$ are typically close in value the total round error on $f(x+h)-f(x)$ will be $2C|f(x)|$, where $C$ is machine precision (typically $10^{-16}$ in Phython). The total error for our derivative is thus $$\epsilon =\frac{2C|f(x)|}{h} + \frac{1}{2}h|f''(x)|.$$

We want to find the value of $h$ that minimizes this error, so we differentiate with respect to $h$ and set the result equal to zero, which gives $$-\frac{2C|f(x)|}{h^2}+\frac{1}{2}|f''(x)|=0, \quad \text{or equivalently} \quad h=\sqrt{4C\left|\frac{f(x)}{f''(x)}\right|}  .$$ Substituting this back into our expression for the error $\epsilon$, we find that the lowest error on our deivative is $$\epsilon=\sqrt{4C|f(x)f''(x)}.$$ If both $f(x)$ and $f''(x)$ are of order $1$, we should choose $h$ to be of order $\sqrt{C}$, which will typically be $10^{-8}$ and our final error will also be of order $\sqrt{C}$, i.e. $10^{-8}$. In many cases, this might be ok, but it is significantly less accurate than what we have seen so far and we can do better.

## Central differences

A simple improvement on the forward and backward difference is the _central difference_: $$\frac{df}{dx}\approx\frac{f((x+h/2)-f(x-h/2)}{h}.$$ The central difference is similar to the forward and the backward difference, approximating the derivative using the difference between two values of $f(x)$ at points a distance $h$ apart. What has changed is that the two points are now placed symmetrically around $x$, one at a distance $\frac{1}{2}h$ in the forward and one at $-\frac{1}{2}h$ in the backward direction. 

To calculate the approximation error on the central difference we write two Taylor expansions: $$\begin {align} & f(x+h/2)=f(x)+\frac{1}{2}hf'(x)+\frac{1}{8}h^2f''(x)+\frac{1}{48}h^3f'''(x)+\ldots \\ &  f(x-h/2)=f(x)-\frac{1}{2}hf'(x)+\frac{1}{8}h^2f''(x)-\frac{1}{48}h^3f'''(x)+\ldots\end{align}$$ Subtracting the second expression from the first and rearranging for $f'(x)$, we get $$f'(x)=\frac{f((x+h/2)-f(x-h/2)}{h}-\frac{1}{24}h^2f'''(x)+\ldots.$$ To leading order the magitude of the error is now $\frac{1}{24}h^2|f'''(x)|$, which is one order in $h$ higher than before. The size of the rounding error remains unchanged, so the total error on our derivative is $$\epsilon =\frac{2C|f(x)|}{h} + \frac{1}{24}h^2|f'''(x)|.$$ Differentiating to find the minimum and rearranging, we find that the optimal value of $h$ is $$h=\left(24C\left|\frac{f(x)}{f'''(x)}\right| \right)^{\frac{1}{3}}.$$ Substituting this back into the error itself, we find the optimal error to be $$\epsilon=(\frac{9}{8}[f(x)]^2|f'''(x)|)^\frac{1}{3}.$$ If we again assume that $f(x)$ and $f'''(x)$ are of order 1, the ideal $h$ is now of order $C^\frac{1}{3}$, which is typically $10^{-5}$, but the error is of order $C^\frac{2}{3}$, which will be $10^{-10}$. 

Thus, the central difference is about a factor $100$ more accurate than the forward or backwards difference. And we achieve this with a larger value of $h$, which is a bonus.

## Second derivatives

We can also derive numerical approximations for the second derivative of a function $f(x)$. The second derivative is, by definition, the derivative of the first derivative, so we can calculate it by applying our first-derivative formulas twice. If we do this for the central difference formula we can write expressions for the first derivative at $x+h/2$ and $x-h/2$: $$ f'(x+h/2)\approx \frac{f(x+h)-f(x)}{h} \quad \text{and} \quad  f'(x-h/2)\approx \frac{f(x)-f(x-h)}{h} .$$  This gives us for the second derivative: $$\begin{align} f''(x) & \approx \frac{f'(x+h/2)-f'(x-h/2)}{h} \\ & = \frac{(f(x+h)-f(x))-(f(x)-f(x-h))}{h^2} \\ &  = \frac{f(x+h)-2f(x)+f(x-h)}{h^2} \end{align}$$ This is the simplest approximation for the second derivative and we will use it extensively for solving second-order differential equations. 

The error of the second derivative can be estimated analogously to before by Taylor expanding $f(x+h)$ and $f(x-h)$. We obtain $$\epsilon=\frac{4Cf(x)}{h^2}+\frac{1}{12}h^2|f''''(x)|.$$ The optimum values of $h$ and $\epsilon$ are $$h=\left(48C\left|\frac{f(x)}{f''''(x)} \right| \right)^\frac{1}{4} \quad \text{and} \quad \epsilon=(\frac{4}{3}C|f(x)f''''(x))^\frac{1}{2} .$$ This means that the error in the 2nd derivative is of order $\sqrt{C}$, which is typically around $10^{-8}$. Our 2nd derivative therefore has about the same error as the forward or backward difference for the first derivative.


```python

```

56277329ae79f723066790835f0c9ae1c1d4c987
ipynb
Jupyter Notebook
Lecture 2 - Integrals and Derivatives/Integrals_and_Derivatives.ipynb
hlappal/comp-phys
8d78a459bc5849ddf5c6c21d484503136bccccbd
MIT
Lecture 2 - Integrals and Derivatives/Integrals_and_Derivatives.ipynb
hlappal/comp-phys
8d78a459bc5849ddf5c6c21d484503136bccccbd
MIT
Lecture 2 - Integrals and Derivatives/Integrals_and_Derivatives.ipynb
hlappal/comp-phys
8d78a459bc5849ddf5c6c21d484503136bccccbd
MIT
Qwen/Qwen-72B
1. YES
2. YES