**Instructions:** Create a new directory called `Differentiation`

with a notebook called `DifferentiationTour`

. Give it a heading 1 cell title **Numerical Differentiation**. Read this page, typing in the code in the code cells and executing them as you go.

**Do not copy/paste**.

Type the commands yourself to get the practice doing it. This will also slow you down so you can think about the commands and what they are doing as you type them.</font>

Save your notebook when you are done, then try the accompanying exercises.

```
In [ ]:
```%pylab inline
import numpy as np
import matplotlib.pyplot as plt

You hear a lot less about numerical derivatives than integrals, for a number of reasons:

The basic techniques for numerical derivatives are quite simple, so they don’t take long to explain.

Derivatives of known functions can always be calculated analytically, so there’s less need to calculate them numerically.

There are some significant practical problems with numerical derivatives, which means they are used less often then numerical integrals.

The standard definition of a derivative, the one you see in the calculus books, is

$$ \frac{df}{dx} = \lim\limits_{h\rightarrow 0} \frac{f(x+h)−f(x)}{h} $$*forward difference* method - it’s simply the slope of the curve $f(x)$ measured over a small interval of width $h$ in the forward direction from $x$.

The backward difference, which has the mirror image definition, is given by

$$ \frac{df}{dx} \simeq \frac{f(x) - f(x-h)}{h} $$```
In [ ]:
```x = np.arange(0.,10.01,0.1)
y = x**5

```
In [ ]:
```dy = np.diff(y) #y(n+1)-y(n)
#Check out the help on the diff method if you want to learn more:
#help(np.diff)
print len(y), len(dy)
dx = np.diff(x)
dydx = (dy/dx)
print len(dydx)

The following won’t work. Can you see why?

```
In [ ]:
```plt.plot(x,dydx)

`x`

is longer by one element than the derivative vector `dydx`

. If we want to have the `dydx`

with the same length as our original `x,y`

, we can do:

```
In [ ]:
```#specify the size of dy ahead because diff returns an array of n-1 elements
dydx1 = np.zeros(y.shape,float) #we know it will be this size
dydx1[0:-1] = np.diff(y)/np.diff(x)
dydx1[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

which will use the forward finite difference for all elements but the last one, where we use a backward difference.

Recall the slicing syntax:

```
array[:] means take the whole array
array[2:7] means take the third through 6th elements (Python numbering)
array[-1] means return the last element in the array (negative numbers count back from the end)
array[1:-1] means take everything but the first and last elements in the array
```

To appreciate how array slicing like `array[2:7]`

works, think of the left-hand-side of the colon as one index and the right-hand-side as another and implement any logic for each side separately. So for example, to implement the two-point forward difference using any arbitrary start and end points, `i`

and `j`

you would use:

`dydx[i:j] = (y[i+1:j+1] - y[i:j])/(x[i+1:j+1] - x[i:j])`

As an alternative to the two-point forward difference, we can calculate the slope at the center of each x-bin. In this case, we have to treat the first and last elements specially, but we can use array slices for everything in between:

```
In [ ]:
```#calculate dydx by center differencing using array slices
dydx2 = np.zeros(y.shape,float) #we know it will be this size
dydx2[1:-1] = (y[2:] - y[:-2])/(x[2:] - x[:-2]) #center difference
dydx2[0] = (y[1]-y[0])/(x[1]-x[0]) #forward difference
dydx2[-1] = (y[-1] - y[-2])/(x[-1] - x[-2]) #backward difference

Let’s compare our numerical derivatives to the exact value:

$$\frac{dy}{dx} = 5x^4$$```
In [ ]:
```dydxExact = 5*x**4

```
In [ ]:
```plt.plot(x,dydx1,label='forward difference')
plt.plot(x,dydx2,label='center difference')
plt.plot(x,dydxExact,label='exact')
plt.legend(loc='upper left')
plt.show()

*slightly* different but it is hard to see. Let’s plot the error introduced by our numerical derivatives. Start from the second value in the array to avoid dividing by zero.

```
In [ ]:
```percentError1 = 100*abs( dydx1[1:] - dydxExact[1:] )/dydxExact[1:]
percentError2 = 100*abs( dydx2[1:] - dydxExact[1:] )/dydxExact[1:]
#Use a semilog y-axis to see the variation better
plt.semilogy(x[1:], percentError1, label='forward difference')
plt.semilogy(x[1:], percentError2, label='center difference')
plt.ylabel("percent error with dydxExact")
plt.legend(loc="upper right")
plt.show()