import numpy as np
import matplotlib.pyplot as plt
Numerical Differentiation¶
Given a set of points $x$ and $f(x)$, how do I compute $f'(x)$?
The best you can do is compute an estimate for $f'(x)$ at each of the points. Recall the definition of a derivative:
$$f'(x) = \lim_{h\rightarrow 0} \frac{f(x + h) - f(x)}{h}$$An easy way to estimate at a particular point, say the 4th data point, $f'(x_4)$ is to use:
$$f'(x_4) \approx \frac{f(x_5) - f(x_4)}{x_5 - x_4}$$[Show secant picture]
The Forward Numerical Derivative¶
$$f'(x_i) \approx \frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i}$$x = np.arange(0,5,0.5) #Create some points
y = x ** 2 #Create the square of those points, now we have f(x_i)
derivative = [] #This is where I will store the answer
#Since this if the forward derivative, we can only go up to the second to last point
for i in range(len(x) - 1):
dx = x[i + 1] - x[i]
dy = y[i + 1] - y[i]
derivative.append(dy / dx)
#We have to "repair" our problem of skipping the last point
#We just repeat the last element
derivative.append(derivative[-1])
plt.plot(x,derivative)
plt.show()
Forward Rule in 1 line with Numpy¶
#Using numpy
#Note x[1:] is all but the first element
#x[:-1] is all but the last element
#So x[1:] - x[:-1] is an idiom for x[1] - x[0], x[2] - x[1], x[3] - x[2], .....
derivative = (y[1:] - y[:-1]) / (x[1:] - x[:-1])
#Here I just skip the last point instead of adding it to the derivative
plt.plot(x[:-1], derivative)
plt.show()
Central Difference Rule¶
Obviously going only forward seems unintuitive. That leads to the central difference rule. You approximate the derivative with the average of the forward and backward approximations. Be careful with the indexes!
$$f'(x_i) = \frac{1}{2}\left[\frac{f(x_{i+1})- f(x_i)}{x_{i + 1} - x_{i}} + \frac{f(x_{i})- f(x_{i-1})}{x_{i} - x_{i-1}}\right]$$forward = (y[1:] - y[:-1]) / (x[1:] - x[:-1])
backward = forward #They're the same numbers, but in forwards we start at 0 go to N-1 and
#in backwards we start at 1 and go to N
central = 0.5 * (forward[:-1] + backward[1:])# Combine the two
plt.plot(x[1:-1], central)
plt.show()
#To fill in, we can use the backward and forward respectively
#We have to declare it ahead of time, because we can't create it all at once. Only in pieces
central = np.zeros(len(x))
#this code is the same as above
forward = (y[1:] - y[:-1]) / (x[1:] - x[:-1])
backward = forward #All that's different is which x-value the points are attached to.
central[0] = forward[0] #Forwards exists at the left point
central[-1] = backward[-1] #Backwards exists at the right point
central[1:-1] = 0.5 * (forward[1:] + backward[:-1]) #Set the rest of the points
plt.plot(x[1:], backward, label="Backward Difference")
plt.plot(x[:-1], forward, label="Forward Difference")
plt.plot(x, central, label="Central Difference")
plt.legend(loc="upper left")
plt.show()
Numerical Integration¶
One strategy is the Riemann integration strategy, which looks like:
$$\int f(x)\,dx \approx \sum^N f(x_i) \Delta x$$[draw picture]
We can do better than that. [Draw Picture]. The trapezoidal rule allows us to better approximate integrals. It's definition is (for non-uniform grid):
$$\int_a^b f(x)\,dx \approx \frac{1}{2} \sum_{i=1}^N (x_{i + 1} - x_i)(f(x_{i + 1}) + f(x_i))$$from math import pi
x = np.linspace(0,pi,10)
y = np.sin(x)
little_traps = 0.5 * (x[1:] - x[:-1]) * (y[1:] + y[:-1])
trap = np.sum(little_traps)
print(trap)
Comparison of Riemann vs Trapezoidal¶
def trap(x,y):
return 0.5 * np.sum( (x[1:] - x[:-1]) * (y[1:] + y[:-1]) )
def rie(x,y):
return np.sum(y[:-1] * (x[1:] - x[:-1]) )
rie_error = []
trap_error = []
points = range(2, 100,1)
for i in points:
x = np.linspace(0,1,i)
y = np.exp(x)
rie_error.append(rie(x,y))
trap_error.append(trap(x,y))
plt.hlines(np.exp(1) - 1, 0,100, label="exact")
plt.plot(points, rie_error, label="Riemann")
plt.plot(points, trap_error, label="Trapezoidal")
plt.xlabel('Number of intervals (trapezoids)')
plt.ylabel('Value of Integral')
plt.xlim(0,100)
plt.legend()
plt.title("Error in $e^x$")
plt.show()
rie_error = []
trap_error = []
points = range(2, 100,1)
for i in points:
x = np.linspace(0,pi,i)
y = np.cos(x)
rie_error.append(rie(x,y))
trap_error.append(trap(x,y))
plt.hlines(np.sin(pi), 0,100, label="exact")
plt.plot(points, rie_error, label="Riemann")
plt.plot(points, trap_error, label="Trapezoidal")
plt.xlabel('Number of intervals (trapezoids)')
plt.ylabel('Value of Integral')
plt.xlim(0,100)
plt.legend()
plt.title("Error in $cos(x)$")
plt.show()
Passing Functions around¶
You can pass functions to functions. Example:
def execute(f):
print(f(5)) #<--- I just assume f is a function!
execute(np.cos)
execute(5)
Why the heck would you want to do this?!¶
Create a function that:
- Plots the given function
- Evaluates an integral on a function
- Samples from the given function
- Calls the function and histograms its output
def plot_function(fxn, xlo, xhi):
x = np.linspace(xlo, xhi, 1000)
plt.plot(x, fxn(x))
plt.show()
plot_function(np.exp, -1, 1)
def trap_integrate(fxn, a, b, points=1000):
x = np.linspace(a,b,points)
y = fxn(x) # <----- Just execute function. Note this only works if fxn is a numpy function!!!
return 0.5 * np.sum( (x[1:] - x[:-1]) * (y[1:] + y[:-1]) )
print(trap_integrate(np.exp, 0,1))
print(np.exp(1) - 1)
Integrating Functions¶
Lots of people know about this hidden feature and some great integration functions exist. Such as the adaptive Gaussian-quadrature called scipy.integrate.quad. This works for functions which can be evaluated anywhere
If you are working with a function (not data) you use scipy.integrate.quad
If you are working with data (not a function) you use your own trapezoidal code
def my_special_function(x):
return x ** 2
from scipy.integrate import quad
#THE FUNCTION MUST BE A NUMPY COMPATIBLE FUNCTION
#So, use vectorize or if it only uses numpy functions (+,-,*,**,np.sin,np.cos...) you're OK
#If statements ARE NOT NUMPY FUNCTIONS
ans, error = quad(my_special_function,0,2) #<-- Notice that this is just like the np.histogram function, which gives two values
print('The answer is {} with an estimated error of {}'.format(ans,error))
def my_special_function(x):
return np.exp(-x**2)
ans,err = quad(my_special_function, 0, np.inf)
print(ans, err)
Let's try a more challlenging integral¶
def confusing_function(x):
if x < 0:
return np.exp(x)
if x < 3:
return x**2
if x == 3:
return 0
if x > 3:
return np.sqrt(x)
if x >= 5:
return np.exp(-x)
#It converts my regular function into a numpy function
np_confusing_function = np.vectorize(confusing_function)
x = np.arange(-4,6,0.01)
plt.plot(x, np_confusing_function(x))
plt.show()
ans,err = quad(np_confusing_function, -4, 6) #Integrate it from -4 to 6 using the numpy version
print(ans,err)
Giving Hints to Quad¶
You can tell quad where it will run into problems, such as dicontinuities.
ans,err = quad(np_confusing_function, -4, 6, points=[0,3,5])
print(ans,err)
Lambda Functions¶
Sometimes it's tedious to define a function just to integrate. For short, you can use the special lambda word:
my_special_function = lambda x: x**2
print(my_special_function(5))
ans, error = quad(lambda x: np.sin(x), 0, 2)
print('The answer is {} with an estimated error of {}'.format(ans,error))
Plotting Integrals¶
To plot an integral, you can loop over multiple integrations
x = np.arange(0,5,0.1)
y = []
for xi in x:
ans,err = quad(lambda x: x ** 2, 0, xi)
y.append(ans)
plt.plot(x,y, label='$\int_0^x\, s^2\,ds$')
plt.plot(x, x**2, label='$x^2$')
plt.legend(loc='upper left')
plt.show()
Summary of Integration and Differentiation¶
If you are integrating a function (not data) you use scipy.integrate.quad
If you are integrating data you use write your trapezoidal code
If you are differentiating a function (not data) just take the derivative using calculus
If you are differentiating data you use the central difference