# 10 Differentiation – Student Activity
See *Computational Physics* (Landau, Páez, Bordeianu), Chapter 5.5

## Implementation of Finite Difference Algorithms in Python

*Forward difference* and *central difference* (discussed in class):

\begin{align}
D_\text{fd} y(t) &\equiv \frac{y(t+h) - y(t)}{h} \\
D_\text{cd} y(t) &\equiv \frac{y\Big(t + \frac{h}{2}\Big) - y\Big(t - \frac{h}{2}\Big)}{h}
\end{align}

and also the *extended difference algorithm*

\begin{align}
D_\text{ep} y(t) &\equiv \frac{4 D_\text{cd}y(t, h/2) - D_\text{cd}y(t, h)}{3} \\
  &= \frac{8\big(y(t+h/4) - y(t-h/4)\big) - \big(y(t+h/2) - y(t-h/2)\big)}{3h}
\end{align}



```python
def D_fd(y, t, h):
    """Forward difference"""
    return (y(t + h) - y(t))/h

def D_cd(y, t, h):
    """Central difference"""
    # implement

def D_ep(y, t, h):
    """Extended difference"""
    # implement
```

### Test your implementations
Test function: $y(t) = \cos t$
1. What is the analytical derivative $\frac{d\cos(t)}{dt}$?
1. Calculate the derivative of $y(t) = \cos t$ at $t=0.1, 1, 100$.
1. Print derivative and relative error $E = \frac{D y(t) - y^{(1)}(t)}{y^{(1)}(t)}$ (finite difference value $D y(t)$ compared to the analystical value $y^{(1)}(t)$– use numpy functions for "exact" values) as function of $h$.
1. Reduce $h$ until you reach machine precision, $h \approx \epsilon_m$
1. Plot $\log_{10} |E(h)|$ against $\log_{10} h$.

Try to do the above for all three algorithms

#### Function definitions 


```python
import numpy as np
# test function
y = np.cos

# Analytical derivative (use np.cos, np.sin, np.exp or whatever else you need)
def y1(t):
    # IMPLEMENT ME (remove `pass`)
    pass

t_values = np.array([0.1, 1, 100], dtype=np.float64)
```

Use numpy functions for everything because then you can operate on all `t_values` at once.

#### Evaluate the finite difference derivatives
Note that we pass *a function* `y` to the forward difference function `D_fd` and we can also pass a whole array of `t_values`!


```python
D_fd(y, t_values, 0.1)
```

#### Evaluate the exact derivatives
Compute the exact derivatives (again, operate on all $t$ together... start thinking in numpy arrays!)


```python
y1(t_values)
```

Calculation of the **absolute error**: subtract the two arrays that you got previously:


```python
# IMPLEMENT ME
```

#### Calculate the relative error $E$


```python
def error(Dxx, y, y1, t, h):
    """Relative error."""
    # INCOMPLETE ... replace None with appropriate terms
    return (Dxx(y, t, h) - None)/None
```

Note that we pass again a general function for the difference operator so that we can use `error()` with `D_fd()`, `D_cd()` and `D_ep()`.


```python
error(D_fd, y, y1, t_values, 0.1)
```

#### Plot $|E|$
Plot $\log_{10} |E(h)|$ against $\log_{10} h$.


```python
import matplotlib.pyplot as plt
%matplotlib inline
```


```python
h_values = 10**(np.arange(-15, -1, 0.1))
abs_errors = np.abs(error(D_fd, y, y1, 0.1, h_values))
```


```python
# INCOMPLETE, replace None with appropriate terms
plt.loglog(None, None, label="t=0.1")
```

Plot the three different $t$ values together in one plot:


```python
# INCOMPLETE: replace None values by appropriate terms
for t in t_values:
    abs_errors = None
    plt.loglog(None, None, label=r"$t={}$".format(t))
ax = plt.gca()
ax.legend(loc="best")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$|E|$")
```

* error behavior depends on $t$ and on cancellation of errors (e.g. for $t=1$
* algorithmic error decreases for decreasing $h$ until the round of error starts dominating


```python

```

1da979e136622ec2482a39a9015d9ed6d009cc77
ipynb
Jupyter Notebook
10_differentiation/10-differentiation-students.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
10_differentiation/10-differentiation-students.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
10_differentiation/10-differentiation-students.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES