Some material comes from the Computational Physics book by Mark Newman at University of Michigan, materials developed by Matt Moelter and Jodi Christiansen for PHYS 202 at Cal Poly, as well as some NumPy tutorials.
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
Derivatives are everywhere in physics. Velocity is the time derivative of position. Acceleration is the time derivative of velocity. Derivatives tell us about the slope of a function and the minimum or maximum of a function occurs where the derivative is zero. We use these facts all the time so let’s investigate some numerical derivatives and their applications.
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} $$This is the 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} $$The forward and backward differences typically give about the same answer and in many cases you can use either. Most often one uses the forward difference. There are a few special cases where one is preferred over the other, particularly when there is a discontinuity in the derivative of the function at the point $x$ or when the domain of the function is bounded and you want the value of the derivative on the boundary, in which case only one or other of the two difference formulas will work. The rest of the time, you just pick whichever one is most convenient. A third option is the centered difference, which we'll learn about in a minute.
Consider a function like, $y = x^5$. Numerically in a computer, this function is represented as an array of $y$ values evaluated at $x$ values.
In [ ]:
x = np.arange(0.,10.01,0.1)
y = x**5
The numerical derivative is calculated using the pre-calculus method of finding $\Delta y$ and dividing by $\Delta x$. This is only a good approximation to the derivative if $\Delta$ is small, but it’s the best we’ve got. We are just going to calculate the slope, $\Delta y/\Delta x$ at every point.
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)
The vector 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()
They are 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()
Since the center difference is the more accurate of the two, we'll generally use that one. The overall accuracy of the derivative depends on the spacing of points. You can explore this phenomenon further with the accompanying exercises.