The rectangle rule - SOLVED

The first step in this process is to place our rectangles. We will begin by deciding where the edges of the rectangles are. To do so we will write a function rectangle_edges which takes 3 inputs N the number of edges that we want, lower_val the value of the smallest edge, upper_val the value of the highest edge. In this function you will need to

1. Define an array of N linearly spaced x values between lower_val and upper_val.

2. Return that array

import numpy as np
import matplotlib.pyplot as plt

def rectangle_edges(N, lower_val, upper_val):
    # COMPLETE THIS FUNCTION AS DESCRIBED
    x_values = np.linspace(lower_val, upper_val, N)
    return x_values

To test your function, run the following code.

print("This should read [0.5, 0.8, 1.1, 1.4, 1.7, 2]")
print(rectangle_edges(6, 0.5, 2.0))

This should read [0.5, 0.8, 1.1, 1.4, 1.7, 2]
[0.5 0.8 1.1 1.4 1.7 2. ]

Guided exercise part 2

Write a function compute_stepsize that computes the gap between the rectangles. HINT the rectangles are equally spaced so you only have to compute the gap between the first two points.

def compute_step_size(rect_edges):
    # COMPLETE CODE
    return rect_edges[1] - rect_edges[0]

rect_edges = rectangle_edges(6, 0.5, 2.0)
stepsize = compute_step_size(rect_edges)

print("Stepsize should be 0.3")
print(f"Stepsize is {stepsize}")

Stepsize should be 0.3
Stepsize is 0.30000000000000004

Guided exercise part 3

We now need to identify the midpoint (on the x-axis) of each of the rectangles. Again, let’s describe what is needed:

1. Write a function f_midpoint. This function should take 2 input values lower_edge and upper_edge and return the midpoint. For example if the lower_edge is 0.5 and the upper_edge is 1, the midpoint is 0.75 (which is halfway between 0.5 and 1).

2. Use the f_midpoint function to find the midpoint of all your rectangles. HINT our example above used 5 rectangles. But 5 rectangles have 6 edges. Therefore if we use N=6 when calling rectangle_edges we will expect only 5 rectangles, and so 5 midpoints. You could do this in a for loop, but is it possible to do it using slicing of numpy arrays? (Do whatever you are happiest with, and talk to us if it is unclear.)

# Write your function f_midpoint here
def f_midpoints(lower_edge, upper_edge):
    return (upper_edge + lower_edge)/2

# Create our set of 5 rectangles from before
rect_edges = rectangle_edges(6, 0.5, 2.0)
stepsize = compute_step_size(rect_edges)
# Now write code to compute the midpoints of these rectangles
# INSERT CODE HERE
#rect_midpoints = []
#for i in range(5):
#    rect_midpoints.append(f_midpoint(rect_edges[i], rect_edges[i+1]))

rect_midpoints = f_midpoints(rect_edges[:-1], rect_edges[1:])

print("The midpoints should be [0.65, 0.95, 1.25, 1.55, 1.85")
print(rect_midpoints) # Change the variable name if you called it something else

The midpoints should be [0.65, 0.95, 1.25, 1.55, 1.85
[0.65 0.95 1.25 1.55 1.85]

Guided exercise part 4

Write a function to compute \(e^{-x^2}\). It should take as input a numpy array of x values, and return an array of \(e^{-x^2}\) for all inputs.

def compute_ex2(x_values):
    # Complete function below
    return np.exp(- x_values**2)


rect_edges = rectangle_edges(6, 0.5, 2.0)
stepsize = compute_step_size(rect_edges)
rect_midpoints = f_midpoints(rect_edges[:-1], rect_edges[1:])

y_values = compute_ex2(rect_midpoints)

print("The y_values should be [0.65540625 0.40555451 0.20961139 0.09049144 0.03263076]")
print(f"The y_values are {y_values}")

The y_values should be [0.65540625 0.40555451 0.20961139 0.09049144 0.03263076]
The y_values are [0.65540625 0.40555451 0.20961139 0.09049144 0.03263076]

Guided exercise part 5

Now we can make a plot to illustrate this. If all your functions above work, you should just be able to run the code below and produce the plots we used for illustration earlier!

# These 4 values are the things that are specific to this example. If we changed these we could compute any integral in this way
N_rectangles = 5
lower_bound = 0.5
upper_bound = 2.0
function = compute_ex2

rect_edges = rectangle_edges(N_rectangles+1, lower_bound, upper_bound)
stepsize = compute_step_size(rect_edges)
rect_midpoints = f_midpoints(rect_edges[:-1], rect_edges[1:])
y_values = function(rect_midpoints)

# Plot the red crosses
plt.plot(rect_midpoints, y_values, 'rx') # Plot the line

# Plot the rectangles, this is a bit more involved
for index in range(N_rectangles):
    plt.fill([rect_edges[index], rect_edges[index], rect_edges[index+1], rect_edges[index+1]],
             [0, y_values[index], y_values[index], 0], alpha=0.5, color='green', linestyle='-')
plt.xlabel('x')
plt.ylabel('$e^{-x^2}$')
plt.xlim(0.5,2)
plt.show()

Guided exercise part 6

Now we can compute the integral. To do this we write another function compute_rect_area. This function takes as input y_values the height of the rectangles and stepsize the width of the rectangles. Note the y_values should be a numpy array (as it is different for each rectangle), but stepsize is just a single value (as all rectangles have the same width). This function should then:

1. Compute the area of each rectangle. Remember area = width * height, so you can quickly compute the area for each rectangle.

2. Sum the area of all the rectangles up

3. Return the sum of the areas of all the rectangles

def compute_rect_area(y_values, stepsize):
    # COMPLETE CODE HERE
    return np.sum(y_values * stepsize)

rect_edges = rectangle_edges(6, 0.5, 2.0)
stepsize = compute_step_size(rect_edges)
rect_midpoints = f_midpoints(rect_edges[:-1], rect_edges[1:])
y_values = compute_ex2(rect_midpoints)

# Now the new bit
integral = compute_rect_area(y_values, stepsize)

print("The integral should be 0.4181083032593553")
print(f"The integral was computed to be {integral}")

The integral should be 0.4181083032593553
The integral was computed to be 0.4181083032593553

Guided exercise part 7

Well done. You’ve just computed the numerical integral of \(e^{-x^2}\). Now let’s tie all this together so we can more easily compute other integrals, and change the ranges. To do this:

Write a function compute_rectangle_integral(function, lower_val, upper_val, num_rectangles). The function should be the function to compute (ie. compute_ex2) in the example above, lower_val is the lower value of integration, upper_val is the upper value and num_rectangles is the number of rectangles to use.

Within this function you should call all the things you’ve already written, so:

1. Call rectangle_edges to get the edges of your rectangles

2. Call compute_step_size to get the step size.

3. Call f_midpoints to get the midpoints of your rectangles

4. Call function to get the y_values (the heights) of your rectangles.

5. Call compute_rect_area to get the integral.

6. Return the integral.

def compute_rectangle_integral(function, lower_val, upper_val, num_rectangles):
    # COMPLETE THIS CODE BELOW
    rect_edges = rectangle_edges(num_rectangles+1, lower_val, upper_val)
    stepsize = compute_step_size(rect_edges)
    rect_midpoints = f_midpoints(rect_edges[:-1], rect_edges[1:])
    y_values = function(rect_midpoints)
    integral = compute_rect_area(y_values, stepsize)
    return integral


integral = compute_rectangle_integral(compute_ex2, 0.5, 2, 100)
print("The integral should be 0.42079376964408877")
print(f"The integral is {integral}")

The integral should be 0.42079376964408877
The integral is 0.42079376964408877

Exercise

Now compute the integral of \(e^{-x^2}\) between \(x=0.5\) and \(x=2\) computed using the “rectangle rule using 100 rectangles.

# ADD CODE BELOW
integral = compute_rectangle_integral(compute_ex2, 0.5, 2, 100)
print(integral)

0.42079376964408877

Exercise

Make a plot of the accuracy of the integral of \(e^{-x^2}dx\) between \(x=0.5\) and \(x=2\) computed using the “rectangle rule” as a function of the number of rectangles used.

The x-axis of the plot should be the number of rectangles used in the integral. Vary this between 1 rectangle and 100 rectangles. The y-axis should show the value of the integral.

HINT You don’t need to write any integration code here, it’s all provided above. Just remove the plotting calls, then call the existing code for many values of the number of rectangles, store the outputs in a list, and then plot.

# ADD CODE BELOW
integral_list = []
for num_rectangles in range(1, 101):
    integral_list.append(compute_rectangle_integral(compute_ex2, 0.5, 2, num_rectangles))

plt.plot(range(1,101), integral_list)
plt.xlabel("Num rectangles")
plt.ylabel(r"Integral of $e^{-x^2}$ between 0.5 and 2")
plt.show()

print(integral_list)

[0.3144170807266467, 0.40226990312386074, 0.41308241713298993, 0.4165522085393729, 0.4181083032593553, 0.41894082588235804, 0.41943854615366943, 0.41975988647889984, 0.4199794284472106, 0.420136083977214, 0.4202517871569165, 0.4203396728984657, 0.42040799942935114, 0.42046217137315095, 0.4205058467893739, 0.4205415734959833, 0.4205711703587752, 0.42059596405710226, 0.4206169407248601, 0.4206348457397349, 0.4206502509823336, 0.4206636011465935, 0.42067524624143, 0.42068546479747904, 0.4206944806978831, 0.42070247555968837, 0.42070959796131263, 0.4207159704020785, 0.4207216946092978, 0.42072685562662393, 0.42073152499337935, 0.42073576323876505, 0.42073962185462843, 0.4207431448677929, 0.420746370102209, 0.4207493301989571, 0.4207520534457769, 0.42075456445571785, 0.4207568847254853, 0.4207590330972486, 0.420761026142529, 0.42076287848284183, 0.4207646030587309, 0.4207662113564653, 0.42076771359985243, 0.42076911891317265, 0.42077043546009235, 0.42077167056252907, 0.4207728308026966, 0.42077392211100656, 0.4207749498420215, 0.4207759188402755, 0.4207768334974816, 0.42077769780239327, 0.4207785153843706, 0.42077928955155797, 0.4207800233243956, 0.4207807194651442, 0.42078138050391395, 0.42078200876169225, 0.42078260637075277, 0.42078317529276826, 0.4207837173349381, 0.420784234164351, 0.42078472732083083, 0.42078519822841076, 0.42078564820563563, 0.42078607847479804, 0.4207864901702437, 0.4207868843458698, 0.4207872619818617, 0.4207876239908093, 0.4207879712232073, 0.42078830447247667, 0.42078862447948256, 0.4207889319366677, 0.42078922749179765, 0.4207895117513808, 0.42078978528376815, 0.42079004862201885, 0.4207903022664767, 0.4207905466871705, 0.42079078232597256, 0.4207910095986147, 0.4207912288965101, 0.42079144058844903, 0.4207916450221471, 0.42079184252567303, 0.4207920334087667, 0.42079221796405536, 0.42079239646817407, 0.4207925691828022, 0.4207927363556325, 0.4207928982212434, 0.4207930550019355, 0.420793206908486, 0.4207933541408684, 0.42079349688890494, 0.42079363533287456, 0.42079376964408877]

Hopefully you can see that the integral very quickly approaches one value. Surprisingly (at least to me) we don’t need many rectangles at all to accurately compute the integral of this function!

Exercise

Now repeat the process for the other two integrals:

Make a plot of the accuracy of \(\int_1^3\frac{\sin{x}}{x} dx\) computed using the “rectangle rule” as described above.

The x-axis of the plot should be the number of rectangles used in the integral vary this between 1 rectangle and 100 rectangles. The y-axis should show the value of the integral.

HINT All you have to do is copy the previous solution and change the function used, and the range used. You will need to write a function to compute \(\frac{\sin{x}}{x}\) to replace the compute_ex2 function.

HINT I compute this integral to be roughly 0.90257

def sinc(x_values):
    return np.sin(x_values) / x_values

integral_list = []
for num_rectangles in range(1, 101):
    integral_list.append(compute_rectangle_integral(sinc, 1, 3, num_rectangles))

plt.plot(range(1,101), integral_list)
plt.xlabel("Num rectangles")
plt.ylabel(r"Integral of $\sin(x) / x$ between 1 and 3")
plt.show()

print(integral_list)

[0.9092974268256817, 0.9043855153776189, 0.9033862496158365, 0.9030307522160774, 0.9028652282287184, 0.9027750573223239, 0.9027206003315867, 0.9026852208380283, 0.9026609489754476, 0.9026435795829132, 0.9026307239754605, 0.9026209438298128, 0.9026133311416681, 0.9026072898133084, 0.902602415401979, 0.9025984256708982, 0.902595118818477, 0.9025923474618903, 0.9025900019332169, 0.9025879992647904, 0.9025862757524618, 0.902584781827737, 0.9025834784544708, 0.9025823345532934, 0.9025813251319593, 0.9025804299088317, 0.9025796322863043, 0.902578918576032, 0.9025782774077981, 0.9025776992739187, 0.9025771761747701, 0.9025767013406085, 0.902576269011454, 0.9025758742616111, 0.9025755128587482, 0.9025751811499988, 0.902574875969298, 0.9025745945615691, 0.9025743345203145, 0.902574093736008, 0.902573870353182, 0.9025736627345743, 0.9025734694310431, 0.9025732891562112, 0.9025731207650239, 0.9025729632354961, 0.9025728156532018, 0.9025726771979454, 0.902572547132362, 0.9025724247920736, 0.9025723095771874, 0.902572200944926, 0.9025720984032277, 0.9025720015051675, 0.9025719098440743, 0.9025718230492532, 0.9025717407822371, 0.9025716627334843, 0.9025715886194579, 0.9025715181800632, 0.9025714511763202, 0.9025713873883856, 0.9025713266137164, 0.9025712686654759, 0.9025712133710683, 0.9025711605708968, 0.902571110117169, 0.9025710618728725, 0.9025710157108756, 0.9025709715130408, 0.9025709291694998, 0.9025708885779519, 0.9025708496430725, 0.902570812275902, 0.9025707763934079, 0.9025707419179673, 0.9025707087769529, 0.9025706769023812, 0.9025706462305413, 0.9025706167016724, 0.9025705882596895, 0.9025705608518823, 0.9025705344287286, 0.9025705089436082, 0.9025704843526534, 0.9025704606145181, 0.9025704376902068, 0.9025704155429491, 0.902570394138019, 0.9025703734425974, 0.9025703534256834, 0.9025703340579398, 0.9025703153115772, 0.9025702971603098, 0.9025702795792201, 0.9025702625446546, 0.9025702460342105, 0.9025702300266039, 0.9025702145016075, 0.9025701994400331]

Make a plot of the accuracy of \(\int_1^3 \sqrt{\sin{x}} dx\) computed using the “rectangle rule” as described above.

The x-axis of the plot should be the number of rectangles used in the integral vary this between 1 rectangle and 100 rectangles. The y-axis should show the value of the integral.

I computed the integral to be 1.717835 here.

def root_sin(x_values):
    return (np.sin(x_values))**(0.5)

integral_list = []
for num_rectangles in range(1, 101):
    integral_list.append(compute_rectangle_integral(root_sin, 1, 3, num_rectangles))

plt.plot(range(1,101), integral_list)
plt.xlabel("Num rectangles")
plt.ylabel(r"Integral of $\sqrt{\sin(x)}$ between 1 and 3")
plt.show()

print(integral_list)

[1.9071417638190211, 1.7723565217937542, 1.7437729029117706, 1.7329954111513168, 1.727770661488054, 1.7248426176024818, 1.7230375411639167, 1.7218466581941605, 1.7210199816151037, 1.7204229203517831, 1.7199777632075626, 1.7196370867075335, 1.719370618598097, 1.7191582997051678, 1.718986412748108, 1.7188453208416572, 1.7187280937317764, 1.7186296447344398, 1.7185461723945057, 1.7184747898477504, 1.7184132730979522, 1.7183598865329102, 1.7183132597415105, 1.7182722990956907, 1.71823612332364, 1.718204015915369, 1.7181753895178602, 1.718149758989429, 1.7181267207884428, 1.7181059370506904, 1.718087123175326, 1.7180700380632599, 1.7180544763797039, 1.7180402623752742, 1.7180272449169796, 1.7180152934659012, 1.7180042948009258, 1.7179941503345941, 1.7179847739018166, 1.7179760899287055, 1.717968031908631, 1.717960541127982, 1.7179535655959823, 1.71794705914205, 1.7179409806514272, 1.717935293415304, 1.7179299645764123, 1.7179249646542096, 1.7179202671370297, 1.7179158481305128, 1.7179116860536578, 1.7179077613752525, 1.7179040563846877, 1.7179005549921107, 1.7178972425536898, 1.7178941057184627, 1.7178911322937904, 1.7178883111268595, 1.7178856320000748, 1.7178830855385634, 1.717880663128065, 1.7178783568421727, 1.7178761593774239, 1.7178740639955004, 1.717872064471519, 1.7178701550478772, 1.7178683303927191, 1.7178665855627375, 1.7178649159697463, 1.7178633173504256, 1.7178617857391631, 1.7178603174434266, 1.7178589090215548, 1.7178575572624928, 1.7178562591676247, 1.7178550119340337, 1.7178538129393859, 1.7178526597282013, 1.7178515499992437, 1.7178504815940254, 1.7178494524863517, 1.7178484607726412, 1.7178475046632373, 1.7178465824742115, 1.7178456926200598, 1.7178448336068037, 1.717844004025728, 1.71784320254767, 1.7178424279176323, 1.7178416789498685, 1.7178409545234132, 1.7178402535778188, 1.7178395751092563, 1.7178389181670255, 1.7178382818501456, 1.7178376653042287, 1.7178370677187385, 1.7178364883242145, 1.7178359263898093, 1.71783538122107]

Rectangle rule - summarized

The rectangle rule is a simple approximation for numerically integrating a function. It can also be done on pen and paper if you want to integrate a function you have plotted …. But we live in a world with computers, so do we really want to do that in the 2020s?

We have shown for the example functions shown that we only need a relatively small number of rectangles to compute this, with a very low computational cost. We have demonstrated how we can show that the integral is “converging” to the right answer by plotting the integral as a function of the number of rectangles used. We will not get a different answer if we use 100000000 rectangles, the integral has already converged. However, there may be functions that require many more points than we used here … But you now know how to check if the integral is accurate enough!

It is possible to compute this using fewer points if we change our approximation to not use rectangles, but to use something else. Let’s explore that now.