<a href="https://colab.research.google.com/github/cohmathonc/biosci670coh/blob/master/IntroductionComputationalMethods/03_IntroCompMethods_NumericalDifferentiation.ipynb" target="_parent"></a>


```python
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
```

## Numerical Differentiation

The continuous derivative of an analytic function $f(x)$ is defined as the limit

$$f'(x)=\lim_{h\rightarrow 0}\frac{f(x+h)-f(x)}{h} \; , \tag{1}$$

where $h$ is in infinitesimal small quantity.

### Discretization

In the discussion of numeric precision (see [this notebook](https://github.com/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/02_IntroPythonForScientificComputing.ipynb)), we saw that *any*  digital number representation is discrete.
This implies that $h$ is finite and we cannot compute the limit $h\rightarrow 0$. 

Often we have to resort to numerical differentiation in cases when the functional expression is unknown, which means that $f(x)$ cannot be evaluated at any arbitrary point $x$ or $x+h$.

The goal of numerical differentiation is to approximate the n-th derivative, given a sequence of values $f(x_0), f(x_1), \ldots f(x_N)$ of an unkown function $f(x)$ at known evaluation points $x_0, x_1, \ldots, x_N$. 
Often the points at which values of a function are known are equally spaced. In that case, the spacing $h$ corresponds to the spacing between subsequent points, i.e. $h=x_{i+1} - x_i = \Delta x$.



```python
## domain
x_min = 0
x_max = 4*np.pi

# we can either 

# (1) fix the number of points, or
# n_points = 100
# x = np.linspace(x_min, x_max, n_points)
# h = (x_max-x_min)/(n_points-1)

# (2) fix the grid spacing
delta_x = 0.1
x = np.arange(x_min, x_max, delta_x)
n_points = len(x)

## function values
y = np.cos(x)

plt.plot(x, y, label='cos(x)')
plt.legend()
```

**Discretization**



### Forward Difference

Identifying $h$ in Eq. (1) with $\Delta x$, we can approximate the derivative of the function at a specific point $x_i$ by **finite differences**:

$$f'(x_i)\approx \frac{f(x_i+\Delta x)-f(x_i)}{\Delta x} = \frac{f(x_{i+1})-f(x_i)}{\Delta x}\;. \tag{2}$$

The formula in Eq. (2) is called the **forward difference formula** because it estimates the function's derivative at $x_i$ using information about the function's value at the point $x_i$ and the 'next' point $x_{i+1}$. 

---

**Exercise:**

- Approximate the derivative of $\cos(x)$ using the forward difference formula (2). 
- Plot $\cos(x)$, $\cos'(x)=-\sin(x)$ and the numerical approximation.
- Examine the behavior of the numerical approximation for varying size of the grid spacing $\Delta x$. 

---


```python
def derivative_forward(y, delta_x):
    """
    Computes forward difference on elements in 'y' array, assuming constant grid spacing 'delta_x'.

    Args:
    - y: an array of data points / function values
    - x: scalar, indicating the spacing between subsequent observation/evaluation points
    """
    n_eval_points = y.shape[0]
    derivative = np.ones(n_eval_points) * np.nan
    for i in range(n_eval_points-1):
        derivative[i] = (y[i+1]-y[i])/delta_x
    return derivative
```

This implementation 
- creates a new array $y'$ (`derivative`)
- scans the input array `y`, and for each index $i$ of `y` computes the finite difference by accessing the $i$-th and $(i+1)$-th element of `y`.
- assigns the finite difference to the $i$-th element of vector $y'$ (`derivative`)



Note that the last element ($i=N$) of the array $y'$ (`derivative`) remains empty because the operation `y[i+1]-y[i]` can only be performed for $i\leq N-1$.

We have multiple options to handle this in our implementation. We could:
1. Create an array `derivative` that has one element less than the array `y`.
2. Create an array `derivative` of same size as array `y` and explicitly indicate that no specific value is assigned to the last element.

Above implementation uses approach (2) by initializing the array `derivative` with elements of the `np.nan` data type ([Not a Number](https://en.wikipedia.org/wiki/NaN)) instead of a specific value, such as ones or zeros.


```python
## THIS IS EQUIVALENT TO ABOVE BUT AVOIDS FOR LOOP BY VECTORIZATION

def derivative_forward_array(y, delta_x):
    """
    Computes forward difference on elements in 'y' array, assuming constant grid spacing 'delta_x'.

    Args:
    - y: an array of data points / function values
    - x: scalar, indicating the spacing between subsequent observation/evaluation points
    """
    # instead of iterating over each pair i, i+1, we create a new shifted array
    y_forward = np.roll(y, -1)     # i-th item in this array corresponds to value at i+1 in y array
    derivative = (y_forward - y)/delta_x
    derivative[-1] = np.nan
    return derivative
```


```python
## Domain
x_min = 0
x_max = 4*np.pi

## Evaluation Grid
delta_x = 0.5
x = np.arange(x_min, x_max, delta_x)

## function values
y = np.cos(x)
y_p = -np.sin(x)                                  # analytic derivative
y_p_approx = derivative_forward(y, delta_x)       # numeric approximation of derivative


# plot f, f', f'_approximated
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 
ax.plot(x, y, 
        label='f(x) = cos(x)' )
ax.plot(x, y_p, 
        label="exact f'(x) = -sin(x)")
ax.plot(x, y_p_approx, 
        label="approx f'(x)", color='red', linestyle=':')

# change x axis ticks to multiples of pi
x_ticks_values = np.arange(0, 4*np.pi, np.pi/2 )
x_ticks_labels = np.arange(0, 4, 0.5)
ax.set_xticks(x_ticks_values)
ax.set_xticklabels(x_ticks_labels)
ax.set_xlabel("x in multiples of $\pi$")

ax.legend()
plt.show()

```

### Sources of Errors

Above example shows that the numerical values of the derivative computed by approximation (2) deviate from the analytic derivative with increasing $\Delta x$.

To estimate the error made by this approximation we focus on a specific point $x_0$ and expand the function $f(x)$ as a [*Taylor Series*](https://en.wikipedia.org/wiki/Taylor_series) around that point.
A Taylor series is a representation of a function at a single point $x_0$ as an infinite sum of terms that are calculated from the values of the function's derivatives at this point:

$$f(x)\big|_{x=x_0} = \sum_{n=0}^\infty\frac{f^{(n)}(x_0)}{n\,!}(x-x_0)^n \; . \tag{3}$$

Only taking into account the 0th ($n=0$), 1st ($n=1$) and 2nd ($n=2$) derivatives, we obtain the following approximation:
$$f(x)\big|_{x=x_0} \approx f(x_0) +(x-x_0)\, f'(x_0) + \frac{(x-x_0)^2}{2}f''(x_0) \;. \tag{4}$$

Where we neglected terms of third order and higher in $(x-x_0)$.


How does this help us in estimating the approximation error in the *forward difference formula* (2) introduced above?

Substituting $h=x-x_0$ in the truncated Taylor expansion (4) gives
$$f(x_0+h) \approx  f(x_0) + h\,f'(x_0) + \frac{h^2}{2}f''(x_0) ;.$$

We can now bring $f'(x_0)$ to the left-hand side and recover the *forward difference formula* from before (2):
    $$f'(x_0)\approx \frac{f(x_0+h)-f(x_0)}{h} -  \frac{h}{2}f''(x_0)\;. $$
    
Thus, the error made by this approximation is:
$$\underbrace{\frac{f(x_0+h)-f(x_0)}{h}}_{\text{forward difference approximation} } - f'(x_0)\approx   \underbrace{\frac{h}{2}f''(x_0)}_{\text{truncation error} } \;. \tag{5}$$


This error is called **truncation error** because the *forward difference formula* is obtained by truncating this term from the *actual* derivative, see Eq. (3). 
For the *forward difference formula*, the truncation error  is proportional to $h=h^1$ ($\mathcal{O}(h)$)
The *forward difference formula* (2) is therefore called a **first order method**.

*Note 1*:
In step from equation (3) to (4), we have truncated multiple terms of higher order in $h$ (all terms of order $n\geq3$). Since we assume that $h$ is very small ($h<<1$), it follows that $h^n>h^{n+1}$. The most important error term is therefore the one of lowest order in $h$ and we can neglect all other terms.

*Note 2*:
The truncation error is proportional to the second derivative (and higher derivatives) of the function whose first derivative we aim to approximate. 
This means that the forward difference formula is the *exact* derivative for linear functions.


---
**Exercise:**

Let's estimate the  error made by the  *forward difference* approximation (2) and compare it with the theoretical truncation error.

For this example we consider the function $f(x)=\cos(x)$ with first derivative $f'(x)=-\sin(x)$ and second derivative $f''(x)=-\cos(x)$. 

We approximate its derivative $f'_{\mathrm{num}}(x)$ at the point $x_0=\pi/6$ and compute the following errors:

\begin{align}
\text{approximation error:}  \qquad\epsilon_{\mathrm{num}} &= \bigl| f'(x_0) - f'_{\mathrm{num}}(x_0) \bigr| \\
\text{truncation error:}  \qquad\epsilon_{\mathrm{trun}} &= -\frac{h}{2}\, f''(x_0) = \frac{h}{2}\cos(x_0)
\end{align}

We repeat this computation for multiple values of $h$, starting from $h=0.1$ and halving its value in each subsequent evaluation, so that $h_i =\frac{0.1}{2^i}$ for $i\in [0, 1 ...,40]$.

---



```python
## Functions and evaluation points
# f=cos(x_0), f'=sin(x_0); x_0 = pi/6
x_0 = np.pi/6
y   = np.cos(x_0)
y_p = -np.sin(x_0)

## Initalization
n_steps = 40
inital_h = 0.1

h_values = np.zeros(n_steps)                  # for h values
y_p_approx_num = np.zeros(n_steps)            # for numeric approximation of f'
eps_num = np.zeros(n_steps)                   # for approximation error
eps_trun = np.zeros(n_steps)                  # for truncation error

## Computations in loop
for i in range(n_steps):
    # compute current h
    h_values[i] = inital_h/(2**i)
    # compute discretization (we only consider x_0 and its neighbouring points)
    x = np.array([x_0-h_values[i], x_0, x_0+h_values[i]])           # x_0 at index 1
    y = np.cos(x)
    # compute numerical approximation of cos(x_0)' with current h
    y_p_approx_num_x0 = derivative_forward(y, h_values[i])[1]       # get approximated derivative at x_0 
    y_p_approx_num[i] = y_p_approx_num_x0                           
    # compute error between analytic and numeric derivative
    eps_num[i] = np.abs(y_p - y_p_approx_num_x0)
    # compute the theoretical (truncation) error from taylor expansion 
    eps_trun[i] = np.cos(x_0)*h_values[i]/2 

## Print results
print("eps_num: ",eps_num)
print("eps_trun: ",eps_trun)


```

Now, we compare those two errors graphically by plotting the errors $\epsilon_{\mathrm{num}}$, $\epsilon_{\mathrm{trun}}$ against $h$. Since both quantities, errors and $h$, vary by multiple orders of magnitude, we use logarithmic axes.


```python
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 

ax.loglog(h_values, eps_num, marker='o', 
          label='error of numerical approximation: $\epsilon_{\mathrm{num}}$' )
ax.loglog(h_values, eps_trun, marker='x', 
          label='theoretical truncation error: $\epsilon_{\mathrm{trun}}$' )

ax.set_xlabel('$h$')
ax.set_ylabel("absolute error")
ax.set_title("Errors in forward-difference approximation of 1st " \
             "derivative of $\cos(\pi/6)$")
ax.legend()
plt.show()
```

For large values $h$, the error of our numerical approximation $\epsilon_{\mathrm{num}}$ agrees with the truncation error $\epsilon_{\mathrm{trun}}$ that is theoretically expected for this *first order* differentiation method.
However, with decreasing $h$, the numerical error *increases* and begins to deviate from the expected truncation error.

This increase in $\epsilon_{\mathrm{num}}$ is not explained by truncation but is caused by *rounding* due to the finite precision of standard numeric data types. This error is called **rounding error**.
**Rounding** and **truncation** are the *two primary sources of errors* when using numerical methods. 

Such errors can have serious consequences. A few extreme cases where bad numerical computing practices lead to disasters are collected [here](http://www-users.math.umn.edu/~arnold/disasters/) and [here](https://www5.in.tum.de/persons/huckle/bugse.html).


### Convergence

We are concerned with numerical methods that approximate the solution $u$ of a continuous problem by computing the solution $\tilde{u}$ of a discretized problem. 
In the previous section we have shown that the accuracy of this approximation depends on the granularity of the discretization, i.e. the spacing $h$. Therefore, we denote the numerical approximation by $\tilde{u}_h$.

An important characteristic of a numerical method is its convergence behavior, that is *whether* and *how quickly* the numerical estimate $\tilde{u}_h$ converges to the exact value $\tilde{u}$ with decreasing $h$.

In Eq. (5), we showed that the truncation error of the *forward difference method*, Eq  (2), is expected to be proportional to $h$ and therefore the method is said to be *of order 1*.
Decreasing $h$ by a factor of e.g. $2$ will decrease the error in the estimation $\tilde{u}_{h/2}$ by the same factor $2$.

In general, when a numerical method is said to be *of order $p$*,  this means that for sufficiently small $h$:
$$\left| \tilde{u}_h-u\right|\leq C\, h^p \tag{5}\; ,$$
where $C$ is a number independent of $h$ that typically depends on the exact solution.
The speed at which the approximation by this numerical method approaches the true solution is called the *convergence rate* of the method and is $h^p$.

If the order $p$ of a given numerical method is known, estimating $p$ from a sequence of approximations is a good way to check the implementation of this method. 

If the exact value $u$ is known, we can compute
$$\log{\left| \tilde{u}_h-u\right|} \leq \log{\left|C\right|} + p\log{h}$$
for multiple different values $h$ and fit a linear function of $\log{h}$ to $\log{\left| \tilde{u}_h-u\right|}$ to approximate $p$.

The standard way to obtain precise estimates for $p$ is to halve the spacing $h$ and compare the ratios of the errors resulting from subsequent spacings:

$$ \frac{\left|\tilde{u}_{h}-u\right|}{\left|\tilde{u}_{h/2}-u\right|} \leq \frac{C\, h^p }{C\, (h/2)^p }=2^p\tag{6}$$
and thus
$$\log_2{\left|\frac{\tilde{u}_{h}-u}{\tilde{u}_{h/2}-u}\right| \leq p } \tag{7}\; .$$



Alternatively, if the the exact value $u$ is not known, comparing the ratios of differences between $\tilde{u}_{h}$ for different $h$ yields a similar result:

$$ \frac{\left|\tilde{u}_{h}-\tilde{u}_{h/2}\right|}{\left|\tilde{u}_{h/2}-\tilde{u}_{h/4}\right|} \leq \frac{1-2^{-p}}{2^{-p}-2^{-2p}}=2^p\tag{8}$$

---
**Exercise (2):**

In the previous exercise we chose $h_{i+1} = \frac{1}{2}h_i$, so that we can now apply formula (8) to verify whether our implementation of the *forward difference method* shows the expected convergence behavior.

---




```python
# compute ratio of errors between steps: 
# ratio(n) = error(n-1)/error(n)
#          = error(h_n)/error(h_n/2)

eps_ratio = np.ones_like(eps_num)*np.nan

for i in range(1,eps_num.shape[0]):         # start at 1 to access index i-1
    eps_ratio[i] = eps_num[i-1]/eps_num[i]

convergence_order = np.log2(eps_ratio)

print("ratio of subsequent error: ", eps_ratio)
print("convergence order: ",convergence_order)
```

### Other Finite Difference Schemes


We will now apply the analysis methods developed above to various finite difference schemas.
To facilitate this analysis we define functions for recurring tasks.


```python
   
def compute_approx_different_h(function_num_deriv, function, 
                               inital_h=0.1, x_0=0, n_steps=40):
    """
    Computes 'n_steps' approximations of derivative of 'function' around 'x_0'
    using 'function_num_deriv', halving spacing delta x in each iteration.

    Args:
    - function_num_deriv: function object that computes numerical derivative from 
                          - array of function evaluations / data points and 
                          - spacing of evaluation points delta_x
    - function: function object for function whose derivative is to be estimated
    - initial_h: float, initial spacing
    - x_0: float, point where derivative of 'function' is to be estimated
    - n_steps: integer number of different spacings steps
    - n: degree of derivative; only defined for first, second and third derivative.

    Returns:
    - array of spacing values, array of estimated derivatives at x_0 for each choosen spacing
    """
    # initialize arrays
    h_values = np.zeros(n_steps)                  # for h values
    y_p_approx_num = np.zeros(n_steps)            # for numeric approximation of f'
    # compute approximations for different h
    for i in range(n_steps):
        # compute current h
        h_values[i] = inital_h/(2**i)
        # compute discretization (we only consider x_0 and its neighbouring points)
        x = np.array([x_0-h_values[i], x_0, x_0+h_values[i]])      # x_0 at index 1
        y = function(x)
        # compute numerical approximation of cos(x_0)' with current h
        y_p_approx_num[i] = function_num_deriv(y, h_values[i])[1]  # get approximated derivative at x_0                         
    return h_values, y_p_approx_num

  
def compute_convergence_order(eps):
    """
    Estimates convergence order from ratio of errors between 
    finite difference estimations with different spacings.
    Assumes that step size between subsequent estimations has been halved.
    """
    eps_shifted = np.roll(eps, shift=1)
    eps_shifted[0] = np.nan
    eps_ratio = eps_shifted / eps
    convergence_order = np.log2(eps_ratio)
    return convergence_order
  
   
def plot_errors(h_vals, eps_num, eps_trun=None, title=None):
    """
    Plots errors vs stepsize on LogLog axes.

    Args:
    - eps_num: array of errors of the numeric approximation
    - eps_trun: array of theoretic truncation errors
    - h_vals: array of h values

    All arrays are expeceted to be of equal length
    """
    # plot
    fig = plt.figure(figsize=plt.figaspect(0.5))
    ax = fig.add_subplot(111) 
    ax.loglog(h_vals, eps_num, marker='o', 
            label='error of numerical approximation: $\epsilon_{\mathrm{num}}$')
    if eps_trun is not None:
        ax.loglog(h_vals, eps_trun, marker='x', 
              label='theoretical truncation error: $\epsilon_{\mathrm{trun}}$')
    ax.set_xlabel('spacing $h$')
    ax.set_ylabel("absolute error")
    if title:
        ax.set_title(title)
    ax.legend()
    plt.show()

```

#### Backward Difference

Similary to the *forward difference* (2), we can define the *backward difference*:

$$f'(x)\approx \frac{f(x)-f(x-h)}{h}\, \qquad \text{with}\qquad h>0 \;. \tag{9}$$

You can easily derive this formula from (2) by assuming a negative $h$, say $h=-\tilde{h}$ where $\tilde{h}>0$.

Like the *forward difference*, this is a *one-sided difference* first order method.

---
**Exercise (3):**

Implement the backward difference method.

---





```python
def derivative_backward(y, delta_x):
    """
    Computes backward difference on elements in 'y' array, assuming constant grid spacing 'delta_x'.

    Args:
    - y: an array of data points / function values
    - x: scalar, indicating the spacing between subsequent observation/evaluation points
    """
    n_eval_points = y.shape[0]
    derivative = np.ones(n_eval_points) * np.nan
    for i in range(1,n_eval_points):
        derivative[i] = (y[i]-y[i-1])/delta_x
    return derivative
```


```python
## THIS IS EQUIVALENT TO ABOVE BUT AVOIDS FOR LOOP BY VECTORIZATION

def derivative_backward_array(y, delta_x):
    """
    Computes backward difference on elements in 'y' array, assuming constant grid spacing 'delta_x'.

    Args:
    - y: an array of data points / function values
    - x: scalar, indicating the spacing between subsequent observation/evaluation points
    """
    # instead of iterating over each pair i, i+1, we create a new shifted array
    y_backward = np.roll(y, -1)     # i-th item in this array corresponds to value at i+1 in y array
    derivative = (y_backward - y)/delta_x
    derivative[0] = np.nan
    return derivative
```


```python
## Functions and evaluation points
x_0 = np.pi/6
function = np.cos
n_steps=40
inital_h=0.1

# Compute approximations
h_values, y_p_approx_num = compute_approx_different_h(function_num_deriv=derivative_backward, 
                                                        function=function, 
                                                        x_0=x_0, n_steps=n_steps, inital_h=inital_h)
# Compute numerical error
y_p = -np.sin(x_0)
eps_num = np.abs(y_p - y_p_approx_num)
# Compute convergence order
convergence_order  = compute_convergence_order(eps_num) 
# compute theoretical (truncation) error from taylor expansion 
eps_trun = np.cos(x_0)*h_values/2 
# plot
plot_errors(h_values, eps_num, eps_trun, 
            title="Errors of backward difference approximation of 1st derivative")

print("Convergence order: ",convergence_order)

```

#### Centered Difference

Forward and backward difference methods are both based on *one-sided differences*.
Could we use information from both sides of the current evaluation point $x_0$ to inform the numerical derivative?


From the Taylor expansion, equation (3), we can derive truncated approximations for the forward and backward difference:
\begin{align}
 f(x+h) &= f(x) + h\, f'(x) + \frac{1}{2}h^2\, f''(x) + \frac{1}{6}h^3\, f'''(x) \tag{10a}\\
 f(x-h) &= f(x) - h\, f'(x) + \frac{1}{2}h^2\, f''(x) - \frac{1}{6}h^3\, f'''(x) \tag{10b}
\end{align}
Substracting (10b) from (10a) gives:
$$f(x+h) - f(x-h) = 2h\, f'(x) + \frac{2}{6}h^3\, f'''(x)\,$$
and after rearranging:
$$f'(x) = \frac{f(x+h)-f(x-h)}{2h} - \frac{1}{6}h^2\, f'''(x)\,. \tag{11}$$

This approximation of the derivative is called **centered difference** formula. The highest order truncation error term is of *second order* in $h$ ($\mathcal{O}(h^2)$), which means that the *centered difference* approximation is a second order method.

Higher-order approximations can be derived in a similar manner, using Taylor expansions with more terms.

---
**Exercise (4):**

Implement the centered difference method for the first derivative (10) and verify that your implementation shows the expected convergence behavior.

---


```python
       
def derivative_centered(y, delta_x):
    """
    Computes centered difference on elements in 'y' array, assuming constant grid spacing 'delta_x'.
    Args:
    - y: an array of data points / function values
    - x: scalar, indicating the spacing between subsequent observation/evaluation points
    """
    n_eval_points = y.shape[0]
    derivative = np.ones(n_eval_points) * np.nan
    for i in range(1, n_eval_points-1):
      derivative[i] = (y[i+1]-y[i-1])/(2*delta_x)
    return derivative

```


```python
## THIS IS EQUIVALENT TO ABOVE BUT USES NUMPY'S VECTORIZATION INSTEAD OF FOR LOOP

def derivative_centered_array(y, delta_x):
    """
    Computes centered difference on elements in 'y' array, assuming constant grid spacing 'delta_x'.

    Args:
    - y: an array of data points / function values
    - x: scalar, indicating the spacing between subsequent observation/evaluation points
    """
    y_backward = np.roll(y, 1)     # i-th item corresponds to value at i-1 
    y_forward = np.roll(y, -1)     # i-th item corresponds to value at i+1
    derivative = (y_forward-y_backward)/(2*delta_x)
    derivative[0] = np.nan          
    derivative[-1] = np.nan
    return derivative
```


```python
## Implementation where you can provide either
#  - scalar value for constant grid spacing, or
#  - array of same size as y for actual evaluation points

def derivative_centered_alternative(y, delta_x):
    """
    Computes centered derivative of data at point `a` with.
    Warning: Derivative at boundary points is not computed

    Args:
    - data: an array of data points / function values
    - x:    can be a scalar, indicating the spacing between subsequent 
            observation/evaluation points, or
            an array of same length as data, indicating the observation points themselves
    """
    n_eval_points = data.shape[0]
    derivative = np.ones(n_eval_points) * np.nan
    if isinstance(x, np.ndarray):      # x is array of coordinates
        for i in range(1, n_eval_points-1):
            derivative[i] = (data[i+1]-data[i-1])/(x[i+1]-x[i-1])
    elif (isinstance(x, float) or isinstance(x, int)):
        for i in range(1, n_eval_points-1): # x is spacing between points
            derivative[i] = (data[i+1]-data[i-1])/(2*x)
    return derivative
```


```python
## THIS IS EQUIVALENT TO ABOVE BUT USES NUMPY'S VECTORIZATION INSTEAD OF FOR LOOP

def derivative_centered_alternative_array(data, x):
    """
    Returns centered derivative of function `function` at point `a` with.
    Warning: Derivative at boundary points is not computed

    Args:
    - data: an array of data points
    - x:    can be a scalar, indicating the spacing between subsequent 
            observation/evaluation points, or
            an array of same length as data, indicating the observation points
    """
    data_backward = np.roll(data, 1)     # i-th item corresponds to value at i-1 
    data_forward = np.roll(data, -1)     # i-th item corresponds to value at i+1
    if isinstance(x, np.ndarray):
        x_backward = np.roll(x, 1)       # i-th item corresponds to value at i-1 
        x_forward = np.roll(x, -1)       # i-th item corresponds to value at i+1
        derivative = (data_forward-data_backward)/(x_forward - x_backward)
    elif (isinstance(x, float) or isinstance(x, int)):
        derivative = (data_forward-data_backward)/(2*x)
    # above computations use 'periodic' boundary conditions, i.e. x_0 = x_{N+1}
    # but in general, we do not have sufficient information to approximate 
    # f' at the boundary: x_0 and x_N
    derivative[:1] = np.nan          
    derivative[-1:] = np.nan
    return derivative
```


```python
## Functions and evaluation points
x_0 = np.pi/6
function = np.cos
n_steps=40
inital_h=0.1

# Compute approximations
h_values, y_p_approx_num = compute_approx_different_h(function_num_deriv=derivative_centered_array, 
                                                    function=function, 
                                                    x_0=x_0, n_steps=n_steps, inital_h=inital_h)
# Compute numerical error
y_p = -np.sin(x_0)
eps_num = np.abs(y_p - y_p_approx_num)
# Compute convergence order
convergence_order  = compute_convergence_order(eps_num) 
# compute theoretical (truncation) error from taylor expansion 
eps_trun = np.sin(x_0)*h_values**2/6 
# plot
plot_errors(h_values, eps_num, eps_trun, 
            title="Errors of centered difference approximation of 1st derivative")

print("Convergence order: ",convergence_order)

```

#### Higher order Derivatives

A similar approach as for the *centered difference* scheme can be used to obtain finite approximations of higher derivatives, such as $f''(x)$, $f'''(x)$, etc.

Let's derive a second order scheme for the second derivative:

We add the fourth derivative term ($n=4$) to equations (10) and consider the sum of these extended (10a) and (10b):
$$f(x+h) + f(x-h) = 2\, f(x) + h^2\, f''(x) + \frac{2}{24}h^4\, f''''(x)\;.$$

Rearrangement yields a second order finite difference approximation of the second derivative:
$$f''(x) = \frac{f(x+h) - 2\,f(x) + f(x-h)}{h^2}-\frac{2}{24}h^2\, f''''(x) \tag{12}$$

---
**Exercise (5):**

Implement Eq. (12) and verify that your implementation shows the expected convergence behavior.

---

### Derivatives at Domain Boundaries 

Note that our ability to approximate the derivative of a function at the domain boundary with a given finite difference scheme may be limited by the available information.
Derivatives at these points can typically be approximated by a combination of one-sided and centered finite difference schemes.

For example, the centered difference scheme (11) requires function values from two neighboring points $f(x_{i-1})$, $f(x_{i+1})$ to approximate $f'(x_i)$, and therefore cannot provide estimates for the derivative at points $x_0$ and $x_N$.

On the other hand, forward and backward difference formula only require information from the current point $x_i$ and the 'next' $x_{i+1}$ (forward) or  'previous' point $x_{i-1}$ (backward), respectively.
Hence, the first derivative on the entire grid $x_0, x_1, \ldots, x_N$ could be estimated using a combination of forward, centered and backward difference schemes.

### 'Exact' Numerical Derivative

When the truncation error of a finite difference method vanishes for a given function, the method computes the exact derivative of that function.

Consider the following polynomials:

- Linear: $f(x)=1+x$
- Quadratic: $f(x)=1+x+x^2$
- Cubic: $f(x)=1+x+x^2+x^3$

We now approximate the first derivative $f'(x)$ at $x=10$ for each of these Polynomials using the centered difference scheme and compute the approximation error in function of spacing $h$:


```python
# Functions and derivatives
def linear_function(x):
    return 1+x

def linear_function_deriv(x):
    return 1

def quad_function(x):
    return 1+x+x**2

def quad_function_deriv(x):
    return 1+2*x

def cube_function(x):
    return 1+x+x**2+x**3

def cube_function_deriv(x):
    return 1+2*x+3*x**2


## Evaluation points
x_0 = 10
n_steps=40
inital_h=0.1

# Linear
h_values_1, y_p_approx_num_1 = compute_approx_different_h(function_num_deriv=derivative_centered_array, 
                                                        function=linear_function, 
                                                        x_0=x_0, n_steps=n_steps, inital_h=inital_h)
y_p_1 = linear_function_deriv(x_0)
eps_num_1 = np.abs(y_p_1 - y_p_approx_num_1)

# Quadratic
h_values_2, y_p_approx_num_2 = compute_approx_different_h(function_num_deriv=derivative_centered_array, 
                                                        function=quad_function, 
                                                        x_0=x_0, n_steps=n_steps, inital_h=inital_h)
y_p_2 = quad_function_deriv(x_0)
eps_num_2 = np.abs(y_p_2 - y_p_approx_num_2)

# Cubic
h_values_3, y_p_approx_num_3 = compute_approx_different_h(function_num_deriv=derivative_centered_array, 
                                                        function=cube_function, 
                                                        x_0=x_0, n_steps=n_steps, inital_h=inital_h)
y_p_3 = cube_function_deriv(x_0)
eps_num_3 = np.abs(y_p_3 - y_p_approx_num_3)


# Plot
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 

ax.loglog(h_values_1, eps_num_1, marker='o', 
        label='error of numerical approximation: Linear Function')
ax.loglog(h_values_2, eps_num_2, marker='o', 
        label='error of numerical approximation: Quadratic Function')
ax.loglog(h_values_3, eps_num_3, marker='o', 
        label='error of numerical approximation: Cubic Function')

ax.set_xlabel('spacing $h$')
ax.set_ylabel("absolute error")
ax.legend()
plt.show()
```

The truncation error of the centered difference scheme is proportional to the third derivative of the function and therefore vanishes for linear and quadratic functions.
The approximation error for those cases is dominated by numeric rounding errors for all choices of $h$ and increases monotonically with decreasing spacing $h$.
This is not the case for the cubic function where the truncation error does not vanish. 

## Useful Python / Numpy Functions




```python
# np.gradient(); 1st derivative, 2nd order
# check help(np.gradient) for info

x_min = 0
x_max = 4*np.pi
n_points = 100

x = np.linspace(x_min, x_max, n_points)
y = np.sin(x)
y_p = np.gradient(y, x)

# plot f, f'_approximated
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 
ax.plot(x, y, 
        label='f(x) = cos(x)' )
ax.plot(x, y_p, 
        label="approx f'(x)", color='red', linestyle=':')

# change x axis ticks to multiples of pi
x_ticks_values = np.arange(0, 4*np.pi, np.pi/2 )
x_ticks_labels = np.arange(0, 4, 0.5)
ax.set_xticks(x_ticks_values)
ax.set_xticklabels(x_ticks_labels)
ax.set_xlabel("x in multiples of $\pi$")

ax.legend()
plt.show()
```


```python
# np.diff(x, n) computes the difference x_i+n - x_i and shortens the array by n

x = np.linspace(0, 4*np.pi, 40)
y = np.sin(x)

# using np.diff, the forward difference scheme becomes
y_p = np.diff(y, n=1) / np.diff(x, n=1)
x_shortened = x[1:]    

# plot f, f'_approximated
fig = plt.figure(figsize=plt.figaspect(0.5))
ax = fig.add_subplot(111) 
ax.plot(x, y, 
        label='f(x) = cos(x)' )

# note that we removed x[0]
ax.plot(x_shortened, y_p, 
        label="approx f'(x)", color='red', linestyle=':')

# change x axis ticks to multiples of pi
x_ticks_values = np.arange(0, 4*np.pi, np.pi/2 )
x_ticks_labels = np.arange(0, 4, 0.5)
ax.set_xticks(x_ticks_values)
ax.set_xticklabels(x_ticks_labels)
ax.set_xlabel("x in multiples of $\pi$")

ax.legend()
plt.show()
```

## Exercises

- In [this](https://github.com/cohmathonc/biosci670/blob/master/IntroductionComputationalMethods/exercises/04_NumericalDifferentiation.ipynb) exercise you extend the centered difference scheme to the domain boundary and approximate the first derivative of a sequence of given data points.

###### About 
This notebook is part of the *biosci670* course on *Mathematical Modeling and Methods for Biomedical Science*.
See https://github.com/cohmathonc/biosci670coh for more information and material.

5291640e7871454771eb23d330d5e5673c7ccc07
ipynb
Jupyter Notebook
IntroductionComputationalMethods/03_IntroCompMethods_NumericalDifferentiation.ipynb
sbranciamore/biosci670coh
92e91aa1e0266b54b0831537b857796eaf0b949c
CC-BY-4.0
IntroductionComputationalMethods/03_IntroCompMethods_NumericalDifferentiation.ipynb
sbranciamore/biosci670coh
92e91aa1e0266b54b0831537b857796eaf0b949c
CC-BY-4.0
IntroductionComputationalMethods/03_IntroCompMethods_NumericalDifferentiation.ipynb
sbranciamore/biosci670coh
92e91aa1e0266b54b0831537b857796eaf0b949c
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES