Calculus 2 coding example: Using NumPy for numerical integration

Here we will give a quick example demonstrating the use of Simpson’s rule in NumPy, as implemented by Patrick Walls.

Mostly we have been using SymPy, which is great for symbolic manipulation, but not for numerical computation. That’s where NumPy (or its sometimes more powerful friend, SciPy) comes in.

NumPy also implements common mathematical functions, like trig functions and square roots. If we are using a routine that is coded using NumPy, we should use the NumPy version of these functions.

Note that functions are treated a bit differently as well: we are not just trying to define a symbolic expression. In SymPy, we can (after defining x as a symbol) enter something like f = sy.sqrt(1+3*(sy.sin(x))**2)), and if we ask the computer to display f, we get a nice expression: \(\sqrt{1+3\sin^2(x)}\).

Now, when we define a function, we need to tell the computer what to do with a given input. To define the same function in NumPy, we would do:

def f(x):
    return np.sqrt(1+3*(np.sin(x))**2)

A shortcut for this uses Python’s lambda syntax: we can also write this in one line as

f = lambda x: np.sqrt(1+3*(np.sin(x))**2)

Let’s see how Simpson’s Rule works.

import numpy as np # import numpy. We don't need matplotlib since we're not graphing anything.

# this is a direct copy-paste from Patrick's website
def simps(f,a,b,N=50):
    '''Approximate the integral of f(x) from a to b by Simpson's rule.

    Simpson's rule approximates the integral \int_a^b f(x) dx by the sum:
    (dx/3) \sum_{k=1}^{N/2} (f(x_{2i-2} + 4f(x_{2i-1}) + f(x_{2i}))
    where x_i = a + i*dx and dx = (b - a)/N.

    Parameters
    ----------
    f : function
        Vectorized function of a single variable
    a , b : numbers
        Interval of integration [a,b]
    N : (even) integer
        Number of subintervals of [a,b]

    Returns
    -------
    float
        Approximation of the integral of f(x) from a to b using
        Simpson's rule with N subintervals of equal length.

    Examples
    --------
    >>> simps(lambda x : 3*x**2,0,1,10)
    1.0
    '''
    if N % 2 == 1:
        raise ValueError("N must be an even integer.")
    dx = (b-a)/N
    x = np.linspace(a,b,N+1)
    y = f(x)
    S = dx/3 * np.sum(y[0:-1:2] + 4*y[1::2] + y[2::2])
    return S

def f(x):
    return np.sqrt(1+4*x**2) #note we need the numpy square root

simps(f,0,4,20) # 

import sympy as sy #now do it exactly with sympy
x = sy.symbols('x')
g = sy.sqrt(1+4*x**2) #this is now a sympy square root
g #display the function, which we can do because this is sympy. Use g so we don't overwrite f above.

sy.integrate(g,(x,0,4))

sy.integrate(g,(x,0,4)).evalf()

Another way to do Simpson’s rule is with the Integrate library from SciPy, which already has a built-in function for this. However, to integrate \(f(x)\) over the interval \([a, b]\) this implementation requires us to first create the partition of \([a,b]\) to be used. The linspace function from NumPy does this for us. Note that we first give names to the values we will use, in case we want to change those later.

import scipy.integrate as spi
a=0; b=4; N=20
x = np.linspace(a,b,N+1) # a partition into N subintervals requires N+1 points
y = f(x)
spi.simps(y,x)

And for comparison:

simps(f,a,b,N)

Looks good!

Finally, what if you do want to work in SymPy first? Maybe you want to define an arclength function that takes \(f(x)\) and \([a,b]\) as input, and gives you the length of \(y=f(x)\), \(a\leq x\leq b\) as output.

And maybe you’d rather calculate the derivative symbolically, instead of numerically. (Then we’d have to talk about how to define numerical differentiation!)

To do that, you need to start out working in SymPy, and then convert. To do this, SymPy has a lambdify command. Let’s see how we could use that. First, let’s define our function: \(f(x)=x^2\).

fx = x**2
fx

Whoops!! What went wrong here? We originally defined x as a SymPy symbol, but later, we defined x to be the partition we used for Simpson’s Rule! We’d better redefine it. Or better yet, use a different variable.

z = sy.symbols('z')
fz = z**2
fz

That’s better. Now, we want to compute the function to be integrated:

s = sy.sqrt(1+(fz.diff(z))**2)
s

This is what we want to integrate. But if we put this directly into our simps function, it’s going to fail. First, we “lambdify”:

F = sy.lambdify(z,s,"numpy") # changes expression s, defined using variable z, to a numpy function

F(x) # check that we can apply this to the array x

Hooray! Now let’s feed this through our simps function:

simps(F,0,4,20)

Bingo. Finally, maybe you want to package all this up into a function. We could do that too:

def arclength(f,a,b,N): 
    '''approximate arclength of y=f(x) from a to b, using N segments.
    
    Assumes that:
    - a function `simps` has been defined to implement Simpson's rule,
      taking arguments a, b, N
    - z has been defined as a SymPy symbol
    - f has been defined as a SymPy function using the variable z'''
    
    F = sy.sqrt(1+(f.diff(z))**2)
    G = sy.lambdify(z, F, "numpy")
    return simps(G,a,b,N)    

z = sy.symbols('z')
f = z**2
arclength(f,0,4,20)