Contact Us!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download
Views: 240
License: MIT
Image: ubuntu2004
Tweet
Share
Share
Kernel: Python 3 (system-wide)

Numerical Integration Python Lab

In this lab, we look at using computers to approximate the values of definite integrals. Recall that this corresponds to the area under a curve.

Remember that a line that starts with # is a comment line, not executable code. Also, if multiple lines are enclosed in three quote marks, those lines are also comments.

To execute a line of code, click anywhere inside the box, hold down the Shift key and press Enter.

Enter your name and your group members' names in the line below. The \n starts a new line. Then execute the line.

print('\r Trapezoidal Rule and Simpson\'s Rule\n Methods of Approximating Area\n Your Name\n My group members are Name1 and Name2.')
Trapezoidal Rule and Simpson's Rule Methods of Approximating Area Your Name My group members are Name1 and Name2.

Now we need to call up the numeric Python library, numpy, and the plotting or graphing library, pyplot. Execute the following line of code.

from matplotlib import pyplot as plt
import numpy as np

We will use the function x3ex\frac{x^3}{e^x} over the interval [0, 8]. That is, we are going to approximate ∫08x3exdx\large\int_0^8{\frac{x^3}{e^x}}dx. While this function does, in fact, integrate, we will use it to demonstrate the techniques.

The true, exact value of the area is 6−758e−8=6−758e8≈5.745719328049896.6-758e^{-8} = 6 - \frac{758}{e^8} \approx 5.745719328049896. Let us save this value, so we can judge how close our approximations are to the true value.

This line of code is the first that you will copy into Practice 1.

"""Enter the next line in one input line, then hold the Shift key while hitting Enter. """
exact = 6-758*np.exp(-8)

# Print out the value for exact, so you can verify it is correct.
exact
5.745719328049896

Part I: The Trapezoidal Rule

In the following input line, we define a program called trapezoidal, which computes the approximate area under a curve over a given interval using the Trapezoidal Rule.

def trapezoidal(f, a, b, n):
    N = 10 # Use N*n+1 points to plot the function smoothly
    dx = float(b-a)/n
    result = 0.5*f(a) + 0.5*f(b)
    for i in range(1, n):
        result += f(a + i*dx)
    result *= dx
        
    x = np.linspace(a,b,n+1)
    y = f(x)
    
    rightx = np.linspace(a, b, N*n+1)
    righty = f(rightx)
    
    plt.plot(rightx,righty)
    
    for i in range(n):
        xs = [x[i], x[i], x[i + 1], x[i + 1]]
        ys = [0,f(x[i]),f(x[i+1]),0]
        plt.fill(xs,ys,'b',edgecolor='b',alpha=0.2)
    
    plt.title('Trapezoid Rule, n = {}'.format(n))
    plt.show()
    
    return result

In the next input line, we use the lambda command to define a funtion quickly. We set the number of sub-intervals, n, to be 4.

Next we set the variable trap_4 to be the output of the trapezoidal program. Then we print out the value of it.

Finally, we compute the error in our approximation by subtracting the value from the exact value.

For Practice 1, you will copy the following input line and paste it into the second empty slot after Practice 1.

v = lambda t: t**3/np.exp(t) # The lambda command is a quick way to define a function without def.
n = 4
trap_4 = trapezoidal(v, 0, 8, 4)

print('The number of subintervals is n = ',n,'.\n The approximate area is trap_4 = ',trap_4)

print ('Trapezoidal Rule, n = 4, error is ', exact - trap_4)
Image in a Jupyter notebook
The number of subintervals is n = 4 . The approximate area is trap_4 = 5.7523441153497314 Trapezoidal Rule, n = 4, error is -0.006624787299835511

Now run the program for n = 8. Change the name of the variable to trap_8.

For Practice 1, you will copy the following input line and paste it into the third empty slot after Practice 1.

n = 8
trap_8 = trapezoidal(v, 0, 8, 8)

print('The number of subintervals is n = ',n,'.\n The approximate area is trap_8 = ',trap_8)

print ('Trapezoidal Rule, n = 8, error is ', exact - trap_8)
Image in a Jupyter notebook
The number of subintervals is n = 8 . The approximate area is trap_8 = 5.743321233849517 Trapezoidal Rule, n = 8, error is 0.002398094200378509

Now run the program for n = 16. Change the name of the variable to trap_16.

For Practice 1, you will copy the following input line and paste it into the fourth empty slot after Practice 1.

n = 16
trap_16 = trapezoidal(v, 0, 8, 16)

print('The number of subintervals is n = ',n,'.\n The approximate area is trap_16 = ',trap_16)

print ('Trapezoidal Rule, n = 4, error is ', exact - trap_16)
Image in a Jupyter notebook
The number of subintervals is n = 16 . The approximate area is trap_16 = 5.743975569597121 Trapezoidal Rule, n = 4, error is 0.0017437584527746353

For Practice 1 below,

  1. In the first empty slot below, copy the third input line. Redefine the value of exact and name it exact1. Let exact1 = 1/2−1/2cos⁡(25)≈0.004398594068263317.1/2 -1/2\cos(25) \approx 0.004398594068263317. Remember to use the command np.cos(25).

  2. Copy the three previous input lines into the next three empty slots.

    1. In the second slot,

      1. change the name of the function to w and define it using the lambda command. Don't forget to use t*np.sin(t**2). The * marks explicit multiplication and the ** marks an exponent.

      2. Leave n as 4.

      3. Remember to change the name of the function to w when you call Trapezoidal.

      4. Remember to change the interval to 0 to 5 when you call Trapezoidal.

      5. Change the error to exact1 - trap_4.

    2. In the third slot,

      1. Leave n as 8.

      2. Remember to change the name of the function to w when you call Trapezoidal.

      3. Remember to change the interval to 0 to 5 when you call Trapezoidal.

      4. Change the error to exact1 - trap_8.

    3. In the fourth slot,

      1. Leave n as 16.

      2. Remember to change the name of the function to w when you call Trapezoidal.

      3. Remember to change the interval to 0 to 5 when you call Trapezoidal.

      4. Change the error to exact1 - trap_16.

Don't forget to Shift + Enter to execute the code.

There are several ways to copy a line of code. The easiest here is to

  1. click anywhere in the input line you want to copy.

  2. Use CTRL + A (pc) or CMD + A (Mac) to select all in the input line. (Recall that an input "line" will probably be multiple lines of code. A single input line has In followed by square brackets (with a number in the middle after execution).)

  3. Then use CTRL + c or CMD + c to copy the code.

  4. Next, click inside the empty input line where you want to copy the code.

  5. Use CTRL + v or CMD + v to paste.

Practice 1: Estimate the area under the curve of the function xsin⁡(x2)x\sin(x^2) over the interval [0,5][0, 5]. That is, approximate ∫05xsin⁡x2dx\int_0^5{x \sin{x^2}}dx using the Trapezoidal Rule.

Part II: Simpson's Rule

In the following input line, we define a program called simps_with_plot, which computes the approximate area under a curve over a given interval using Simpson's Rule.

def simps_with_plot(f, a, b, n): #Changed code
    if n % 2 == 1:
        raise ValueError("n must be an even integer.")
    if a > b:
        print('Incorrect bounds. a should be smaller than b, or you will get negative area.')
    dx = float(b - a)/n
    x = np.linspace(a, b, n+1)
    y = f(x)
    SimpSum = dx/3 * np.sum(y[0:-1:2] + 4*y[1::2] + y[2::2]) #This is all but the last line of the program for Simpson's Rule.
    
    # The lines below build the plot.
    plt.figure() 
    for i in range(0, n, 2): # start with x_0, go to x_n by 2s. Two boxes per parabolic top.
        
        """The following two lines produce a matrix representing the equations with three unknowns. """
        coeff_matrix = np.array([[x[i]**2, x[i], 1],[x[i + 1]**2, x[i + 1], 1], [x[i + 2]**2, x[i + 2], 1]])
        const_matrix = np.array([f(x[i]), f(x[i + 1]), f(x[i + 2])])
        """ this is an array representing the solution for the parabola, ax^2 + bx + c"""
        soln = np.linalg.solve(coeff_matrix, const_matrix) 
        """ this is the actual equation for the parabola, ax^2 + bx + c. Don't confuse the a and b here with the a and b for our definite integral. They are not the same."""
        def parab_eqn(t):
            return soln[0]*t**2 + soln[1]*t + soln[2]
        
        # Now to graph the parabola top in the correct space
        x_para = np.arange(x[i], x[i + 2], 0.01)
        
        #plt.figure()
        plt.plot(x_para, parab_eqn(x_para), 'b--', x_para, f(x_para), 'k')
        plt.fill_between(x_para, parab_eqn(x_para), color = 'g')
        
        """The code below would graph the Trapezoidal Rule with Simpson's Rule. Right now it’s turned off.
        for i in range(n):"""
        """xs = [x[i], x[i], x[i + 1], x[i + 1]]
        ys = [0,f(x[i]),f(x[i+1]),0]
        plt.fill(xs,ys,'b',edgecolor='b',alpha=0.2)"""
        
    plt.vlines(x, 0, y, colors = 'r', linestyles = 'solid')
    plt.grid()
    plt.title('Simpsons Rule, n = {}'.format(n))
    plt.show() # This is the end of the plotting part of the code.
    
    return SimpSum # This is the final line of the program.

In this example, we are going to call up the program three times, with different values for n, the number of subintervals. We will return to the function v with exact value given by exact, over the interval [0,8].[0, 8].

In Practice 2, you will copy the following input line to use with the function w.

simps_4 = simps_with_plot(v, 0, 8, 4)
print('Simpson’s Rule with n = 4 is', simps_4)
print ('Simpson’s Rule, n = 4, error is ', exact - simps_4)

simps_8 = simps_with_plot(v, 0, 8, 8)
print('Simpson’s Rule with n = 8 is', simps_8)
print ('Simpson’s Rule, n = 8, error is ', exact - simps_8)

simps_16 = simps_with_plot(v, 0, 8, 16)
print('Simpson’s Rule with n = 16 is', simps_16)
print ('Simpson’s Rule, n = 16, error is ', exact - simps_16)

For Practice 2 below,

  1. Copy the previous input line into the empty slot below Practice 2.

    1. Remember to change the name of the function to w when you call simps_with_plot. It is already defined from before.

    2. Remember to change the interval to 0 to 5 when you call simps_with_plot.

    3. Change the error to exact1 - trap_#, with the appropriate number for all three calls.

Don't forget to Shift + Enter to execute the code.

There are several ways to copy a line of code. The easiest here is to

  1. click anywhere in the input line you want to copy.

  2. Use CTRL + A (pc) or CMD + A (Mac) to select all in the input line. (Recall that an input "line" will probably be multiple lines of code. A single input line has In followed by square brackets (with a number in the middle after execution).)

  3. Then use CTRL + c or CMD + c to copy the code.

  4. Next, click inside the empty input line where you want to copy the code.

  5. Use CTRL + v or CMD + v to paste.

Practice 2: Estimate the area under the curve of the function xsin⁡(x2)x\sin(x^2) over the interval [0,5][0, 5]. That is, approximate ∫05xsin⁡x2dx\int_0^5{x \sin{x^2}}dx using Simpson's Rule.

End of Lab