Lab 05 - Numerical Integration - The Trapezoidal Rule and Simpson's Rule.ipynb

Kernel: SageMath 9.2

Lab 05 - Numerical Integration: The Trapezoidal Rule and Simpson's Rule

Overview

As we learned in Calculus I, there are two ways to evaluate a definite integral: using the Fundamental Theorem of Calculus or using numerical approximations such as Reimann Sums. While the FTC is nice in theory, it cannnot be applied in many cases, as antiderivatives are often difficult or even impossible to find. With the help of computers, using numerical integration is the most practical way to evaluate definite integrals. In this lab, we will use the Trapezoidal Rule and Simpson's Rule to approximate definite integrals.

Important SageMath Commands Introduced in this Lab

ParseError: KaTeX parse error: Undefined control sequence: \hfill at position 32: …|l|l|} \hline \̲h̲f̲i̲l̲l̲ ̲\textbf{Command…

Related Course Material

Section 8.7

In Calculus I, we learned how to use rectangles to approximate the definite integral. Another geometic figure which can be used to approximate the definite integeral is a trapezoid. Recall that the area of a trapzoid is $A = \frac{1}{2}h(b_1 + b_2)$. In order to use trapezoids to approximate the definite integral, we subdivide the interval $[a,b]$ into $n$ subintervals of equal length $\Delta x = \frac{b-a}{n}$. Then, on each subinterval we let $h = \Delta x$ and $b_1 = f(x_{i-1})$ and $b_2 = f(x_i)$. Thus, the area of the trapzoid in the $i^\text{th}$ subinterval is $A = \frac{\Delta x}{2}(f(x_{i-1}) + f(x_i))$. Therefore, we have the following approximation of the definite integral using the Trapezoid Rule: $$\displaystyle \int_a^b f(x) \ dx \approx \dfrac{\Delta x}{2}\sum_{i = 1}^{n} \big(f(x_{i-1}) + f(x_i)\big),$$ where $n$ is the number of trapezoids we are using in our approximation and $x_i = a + i \Delta x$.

By using trapezoids, we are essentially approximating the curve $y=f(x)$ using piecewise linear functions. In Simpson's Rule, we will instead approximate $y = f(x)$ using piecewise quadratic functions. This time, we start by subdividing the interval $[a,b]$ into $2n$ subintervals of equal length $\Delta x$ = $\frac{b - a}{2n}$. On each pair of subintevals $[i-1,i]$ and $[i, i+1]$ with $i$ odd, we approximate the area spanned by the two subintervals with $A = \frac{\Delta x}{3}(f(x_{i-1}) + 4f(x_i) + f(x_{i+1})$. Therefore, we have the following approximation of the definite integral using Simpson's Rule: $$\displaystyle \int_a^b f(x) \ dx \approx \dfrac{\Delta x}{3} \sum_{i=0}^{n-1} \big( f(x_{2i}) + 4f(x_{2i+1}) + f(x_{2i + 2})\big)$$ where $n$ is half of the number of subintervals we are using in our approximation and $x_i = a + i \Delta x$.

Example 1

Let us begin by practicing using lists (or arrays in some programming languages) and $\textbf{for}$ loops. In order to create a list of elements, simply surround the elements in [ $\$ ].

fruits = ['apple', 'banana', 'cherry']

We can return specific values of the list by adjoining ParseError: KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at end of input: \textbf{[iParseError: KaTeX parse error: Expected 'EOF', got '}' at position 2: ]}̲ immediately after the name of the name of the list. The number $i$ is refered to as the index which is the position of the element in the list.

$\textbf{Caution:}$ The starting index of a list is $0$, not $1$ as you might expect. Therefore, a list which has five elements would use the indices $0$ through $4$ and not $1$ through $5$. This is common to most programming languages.

print(fruits[0])
print(fruits[1])
print(fruits[2])

We can adjoin two lists together by using the $\textbf{append}$ command.

fruits.append('durian')
fruits

$\textbf{Caution:}$ Be careful when appending elements to a list since if you reexecute the above code, $\textbf{fruits}$ will contain two elements of 'durian' .

Note that all of the elements in the list $\textbf{fruit}$ are strings. We can add non-strings to the list as well if we prefer. Also, we can adjoin multiple elements at once by using the $\textbf{extend}$ command.

fruits.extend([1, 3/2, 2.5])
fruits

Suppose we wanted to print all of the elements of $\textbf{fruit}$ on individual lines. One way to do this is to have multiple lines containing print statements as we did previously. This method is not ideal, however, if our list contains many elements. A better way to accomplish this is to use a $\textbf{for}$ loop. If we type the code ParseError: KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at end of input: \textbf{"for iParseError: KaTeX parse error: Expected 'EOF', got '}' at position 13: in fruits:"}̲, then SageMath will iterate the variable $i$ through each element of the list $\textbf{fruits}$. We can then follow the $\textbf{for}$ loop with the line $\textbf{"print(i)"}$ which will print the element $i$ as $i$ iterates through all of the elements of $\textbf{fruits}.$

Note: the variable $i$ is arbitrary. You can use any letter or word you want as the iterative variable in the $\textbf{for}$ loop.

for i in fruits:
    print(i)

Example 2

Create a list called "lst" which contains the numbers $1$ through $10$.

One way to do the above is to list all ten numbers from 1 through 10. This is doable, but not if we wanted to list the first one thousand numbers. The command ParseError: KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at end of input: \textbf{range(iParseError: KaTeX parse error: Expected 'EOF', got '}' at position 2: )}̲ returns the list containing the numbers $0$ through $i-1$. This does not give us exactly the list we want, but we can use this along with a $\textbf{for}$ loop to create the desired list.

lst = []
for i in range(10):
    lst.append(i+1)
lst

The above code creates a blank list, and then inserts the the numbers which are one more than the elements of $\textbf{range($10)$}$into the list. Since$\textbf{range($10$)}$ was the list containing the numbers 0 through 9, the new list contains the numbers $1$ through $10$.

A third way we could do this is to use the $\textbf{range(m,n)}$ command.

lst = range(1,11)
lst

If you are using SageMath 9.2 or other versions of SageMath which use Python3, then $lst$ did return a list, but instead it just returned the command $\textit{range(1,11)}.$ In order to convert this command into a list, you need to use the ParseError: KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at end of input: \textbf{list(\dotsParseError: KaTeX parse error: Expected 'EOF', got '}' at position 2: )}̲ command.

lst = list(range(1,11))
lst

Note that we used $\textbf{range(1,11)}$ to obtain the list of numbers 1 through 10. This is because the range command starts at the first number and stops at the last number minus 1.

Now, add the numbers 11 through 13 to $\textit{lst}$.

Return the sum of the fourth and thirteenth elements of $\textit{lst}$.

Use a $\textbf{for}$ loop to print the square of each number in $\textit{lst}$ on its own line.

If all of the elements in a list are numbers, then you can use the $\textbf{sum}$ command to return the sum of all of the elements in the list.

sum(lst)

Example 3

We will use SageMath to approximate $\displaystyle \int_0^1 e^{-x^2} \ dx$ using the Trapezoid Rule with $n = 10$ trapezoids. First, we define our function and assign $a, b, n,$ and $\Delta x$ to their appropriate values.

def f(x):
    return e^(-x^2)
a = 0
b = 1
n = 10
DeltaX = (b-a)/n

Now, we will create a list which contains the summands of our Trapezoid Rule formula. Recall that the Trapezoid Rule is $$\displaystyle \int_a^b f(x) \ dx \approx\dfrac{\Delta x}{2}\sum_{i = 1}^{n} \big(f(x_{i-1}) + f(x_i)\big).$$

We need to add the terms in the summation to our list. Note that the summation is a sum containing $n$ terms. Therefore, we will use a $\textbf{for}$ loop which goes through $n$ terms. Remember that if we use the ParseError: KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at end of input: \textbf{range(nParseError: KaTeX parse error: Expected 'EOF', got '}' at position 2: )}̲ command, then $i$ will go from $0$ to $n-1$, so we will instead use $\textbf{range(1,n+1)}$ to go from $1$ to $n$.

summands = []     ## Creating a blank list
for i in range(1,n+1):
    summands.append(f(a + (i-1)*DeltaX) + f(a + i*DeltaX))     ## Adding f(x_(i-1)) + f(x_i) to the list

Our list now contains all $10$ of our summands. We can check by calling the name of the list to see what it contains.

summands

We can also have SageMath tell us exactly how many elements a list contains by using the $\textbf{len}(\textit{lst})$ command.

len(summands)

Lastly, we can complete our approximation by adding all of the summands together using the $\textbf{sum}$ command and multiplying by $\dfrac{\Delta x}{2}$.

TrapApprox = (DeltaX / 2) * sum(summands)
TrapApprox

In order to get a comprehensable answer, we can use the $\textbf{round}(\textit{num, n})$ to round the answer to $n$ decimal places.

round(TrapApprox, 10)

Therefore, by using the Trapezoid Rule, we get $$\displaystyle \int_0^1 e^{-x^2} \ dx \approx 0.7462107961.$$

We can use SageMath to see how close our approximation is to the actual answer.

Answer = integrate(f(x), x, 0,1)
Error = abs(Answer - TrapApprox)
round(Error, 10)

Therefore, our approximation is accurate to within three decimal places.

Note that, we can combine all of the steps we used to compute the Trapezoid Rule approximation into a single function which will allow us to approximate the definite integral quicker.

def TrapezoidRule(f, a, b, n): 
    DeltaX = (b-a)/n
    summands = []     ## Creating a blank list
    for i in range(1,n+1):
            summands.append(f(a + (i-1)*DeltaX) + f(a + i*DeltaX))     ## Adding f(x_(i-1)) + f(x_i) to the list
    return (DeltaX / 2.0) * sum(summands)

We have just created a function which accepts arguments $f$ for our function, $a$ for our lower bound, $b$ for our upper bound, and $n$ for the number of trapezoids. We can check that this function does return what we expect by using it to calculate the Trapezoid Rule approximation we just did in detail.

def f(x):
    return e^(-x^2)
TrapezoidRule(f, 0, 1, 10)

Again, to make sense of this answer, we can use the $\textbf{round}$ command.

round(TrapezoidRule(f, 0, 1, 10), 10)

Example 4

Now, we will use SageMath to approximate $\displaystyle \int_0^1 e^{-x^2} \ dx$ using Simpson's Rule with $n = 5$.

def f(x):
    return e^(-x^2)
a = 0
b = 1
n = 5
DeltaX = (b-a) / (2*n)

Again, we will create a list which contains all of our summands of our Simpson's Rule formula. Recall that Simpson's Rule is $$\displaystyle \int_a^b f(x) \ dx \approx \dfrac{\Delta x}{3} \sum_{i=0}^{n-1} \big( f(x_{2i}) + 4f(x_{2i+1}) + f(x_{2i + 2})\big).$$ Note that this time our summation starts at $i = 0$ and stops at $i = n-1$, so we can use $\textbf{range(n)}.$

summands = []
for i in range(n):
    summands.append(f(a + 2*i*DeltaX))     ## Adding f(x_2i) to summands
    summands.append(4 * f(a + (2*i + 1)*DeltaX))     ## Adding 4f(x_(2i+1)) to summands
    summands.append(f(a + (2*i + 2)*DeltaX))     ## Adding f(x_(2i+2)) to summands

Now, we can get our approximation by adding all of the elements in summand and multiplying by $\dfrac{\Delta x}{3}.$

SimpApprox = (DeltaX / 3.0) * sum(summands)
SimpApprox

round(SimpApprox,10)

Therefore, Simpson's Rule gives the approximation $$\displaystyle \int_0^1 e^{-x^2} \ dx \approx .7468249483.$$ Now, check how close the approximation is to the actual answer.

It follows that our approximation using Simpson's Rule is accurate to 6 decimal places.

Create the function $\textbf{SimpsonsRule}$ which combines all the steps used to estimate the definite integral using Simpson's Rule. Then check that your function gives the same result as above.

def SimpsonsRule(f, a, b, n):
    ## Insert your code here

Example 5

Use the functions $\textbf{TrapezoidRule}$ and $\textbf{SimpsonsRule}$ to approximate the following definite integrals with $n = 4$. Also, determine the abosolute error in each approximation. Remember that in order to use a variable other than $x$, you have to define it as a variable using the $\textbf{var}$ command.

1. $\displaystyle \int_0^2 (t^3 + t) \ dt$

2. $\displaystyle \int_1^2 \dfrac{1}{s^2} \ ds$

3. $\displaystyle \int_0^\pi \sin(t) \ dt$