# Numerical Integration

Numerical Integration

# 1. Numerical Integration

## 1.2. Objective

• To introduce the concept of numerical integration
• To have an understanding of the trade-off between the number of calculations and accuracy of an estimate of area.

## 1.3. Introduction

The area under a curve can be estimated by drawing rectangles with heights that are determined by the location of the curve. As the width of the rectangles decrease, the accuracy of the estimate increases.

Alternative estimates of the area can be obtained by drawing alternative shapes - including trapezoids and rectangles with a parabolic top. These alternative methods allow for fewer rectangles to be used to obtain the same level of accuracy as when regular rectangles are used.

## 1.4. Rectangle Method

In the rectangle method, one draws vertical rectangles of a given width upward from the x-axis until a location on the top of rectangle passes intersects the curve.

Three possibilities for setting the rectangle heights are

• Left: The rectangle's height is the height of the curve at the left side of the rectangle
• Middle: The rectangle's height is the height of the line at the middle of the rectangle
• Right: The rectangle's height is the height of the line at the right side of the rectangle

An example of each case is shown below.

Note that the heavy black dots indicate the y-value for the x points at which the curve is evaluated. In the rectangle method, the number of evaluation points is equal to the number of rectangles. In the above images, the area, A, of each rectangle is shown along with the error, ΔA, is shown. For this particular curve, the "middle" method happens to have zero error.

In each of the following examples, the code on the left is a long-hand version that computes the total area. The version on the right is a for loop that performs the same computation.

## 1.5. Using "Left" rule

 Long-hand version w = 2; x = [0,2]; y(1) = 1+x(1)/4; a(1) = w*y(1); y(2) = 1+x(2)/4; a(2) = w*y(2); A = a(1) + a(2)  Short-hand version w = 2; x = [0,2]; for i = [1:length(x)] y(i) = 1+x(i)/4; a(i) = w*y(i); end A = a(1) + a(2) % or A = sum(a); 

## 1.6. Using "Middle" rule

 w = 2; x = [1,3]; y(1) = 1 + x(1)/4; a(1) = w*y(1); y(2) = 1 + x(2)/4; a(2) = w*y(2); A = a(1) + a(2)  w = 2; x = [1,3]; for i = [1:length(x)] y(i) = 1 + x(i)/4; a(i) = w*y(i); end A = a(1) + a(2) % or A = sum(a); 

## 1.7. Using "Right" rule

 w = 2; x = [2,4]; y(1) = 1+x(1)/4; a(1) = w*y(1); y(2) = 1+x(2)/4; a(2) = w*y(2); A = a(1) + a(2)  w = 2; x = [2,4]; for i = [1:length(x)] y(i) = 1+x(i)/4; a(i) = w*y(i); end A = a(1) + a(2) % or A = sum(a); 

## 1.8. Rectangle Method Generalization

In the above examples, all of the code had a similar pattern: The width of the rectangle was set along with the evaluation points x. Then a loop was used to compute the area of each rectangle. For example, for the "left" method using two rectangles, the code was

w = 2;
x = [0,2];
for i = [1:length(x)]
y(i) = 1+x(i)/4;
a(i) = w*y(i);
end
A = sum(a);


For three rectangles, the code would be

w = 4/3;
x = [0,4/3,8/3];
for i = [1:length(x)]
y(i) = 1+x(i)/4;
a(i) = w*y(i);
end
A = sum(a);


and for four rectangles

w = 1;
x = [0,1,2,3];
for i = [1:length(x)]
y(i) = 1+x(i)/4;
a(i) = w*y(i);
end
A = sum(a);


When the number of rectangles used is large, a better approach is to use the function linspace to compute the values for x and to compute w based on the values of x. The function linspace takes three arguments, the left boundary, the right boundary, and the number of points to generate between and including the boundary. To create the equivalent of

x = [0,1.33,2.66]


with linspace (corresponding to the three-rectangle case), we first need to request the generation of four points in the range of x=0-4 to get the correct values for x

a = 0;
b = 4;
N = 3; % Number of rectantles
x = linspace(a,b,N+1) % Gives x = [0,1.33,2.66,4]


and then we will ignore the last point in x in the for loop:

N = 3; % Number of rectangles.
b = 4; % Right boundary
a = 0; % Left boundary
x = linspace(a,b,N+1); % Gives x = [0,1.33,2.66,4]
w = x(2)-x(1); % Gives 1
for i = [1:length(x)-1] % Evaluate for i = 1,2,3; ignores last value in array x
y(i) = 1+x(i)/4;
a(i) = w*y(i);
end
A = sum(a);


Note that the number of steps in the loop is 3 (length(x)-1) instead of length(x). The advantage of this approach is to increase the number of rectangles, one only needs to modify the parameter N.

## 1.9. Rectangle Method Vectorization

The above code can be further simplified using vectorization.

N = 4;
x = linspace(0,4,N+1); % Gives x = [0,1,2,3,4]
x = x(1:end-1); % Remove last value in x
w = x(2)-x(1); % Gives 1
y = 1+x/4; % Creates an array y based on array x
a = w*y; % Create an array y based on array y
A = sum(a)


or

N = 4;
x = linspace(0,4,N+1);
x = x(1:end-1); % Remove last value in x
w = x(2)-x(1);
A = sum(w*(1+x/4));


## 1.10. Trapezoid Method

In the previous example, the "middle" method happens to give the exact answer because the overestimate matches the underestimate for each area used in the estimate. One can observe that if we use the "left" rule rectangle plus a triangle to form the areas, we will also get an exact answer, and there is no under or overestimate. An important difference between the rectangle and trapezoid methods is that for the rectangle method N rectangles require N evaluation points whereas for N trapezoids, N+1 evaluation points are needed.

In the above plot, the line y = 1+x/4 is shown along with two trapezoids. The total area of the rectangles + triangles (trapezoids) can be computed as

  x = [0,2,4]; % Note that the number of evaluation points for x and y is equal to the number of trapezoids + 1.
% or
% N = 2;
% x = linspace(0,4,N);
w = x(2) - x(1);
y(1) = 1+x(1)/4;
y(2) = 1+x(2)/4;
a(1) = w*y(1) + 0.5*w*(y(2)-y(1)); %  w*y(1) is rectangle area, 0.5*w*(y(2)-y(1)) is triangle area.
y(3) = 1+x(3)/4;
a(2) = w*y(2) + 0.5*w*(y(3)-y(2));
A = sum(a)


This can be simplified algebraically to

  x = [0,2,4];
w = x(2) - x(1);
y(1) = 1+x(1)/4;
y(2) = 1+x(2)/4;
a(1) = w*0.5*(y(2)+y(1)); % Corresponds to (width)*(average height)
y(3) = 1+x(3)/4;
a(2) = w*0.5*(y(3)+y(2));
A = sum(a)


To better determine how to write this using a for loop, write out the calculation in long-hand form using four trapezoids:

  clear
x = [0,1,2,3,4];
w = x(2) - x(1);
y(1) = 1+x(1)/4;
y(2) = 1+x(2)/4;
a(1) = w*0.5*(y(2)+y(1));
y(3) = 1+x(3)/4;
a(2) = w*0.5*(y(3)+y(2));
y(4) = 1+x(4)/4;
a(3) = w*0.5*(y(4)+y(3));
y(5) = 1+x(5)/4;
a(4) = w*0.5*(y(5)+y(4));
A = sum(a)


Here the repeated pattern, for i=1,2,3,4 is

 y(i+1) = 1+x(i+1)/4;
a(i) = w*0.5*(y(i+1)+y(i));


so the for loop version is

  clear
x = [0,1,2,3,4];
w = x(2) - x(1);
y(1) = 1+x(1)/4;
for i = 1:length(x)-1
y(i+1) = 1+x(i+1)/4;
a(i) = w*0.5*(y(i+1)+y(i));
end
A = sum(a)


## 1.11. Simpson's Method

In the trapezoid method, the top part of the each area was determined by connecting the left and right points with a straight line and the area below the straight line was used to estimate a part of the total area. The next level of approximation is to make the top part of each area a quadratic. To uniquely determine a quadratic, three points are needed. It can be shown that the exact value for the area under a quadratic that passes through the evenly spaced points y(1), y(2), and y(3) is

(h/3)*(y(1) + 4*y(2) + y(3))


where h is the horizontal separation between the equally spaced points and y(1), y(2), and y(3) are height of the curve at the left, middle, and right sides of the area, respectively.

In the case where the area to be approximated is under a quadratic curve, this formula will give an exact answer, independent of the number of areas used. In the following image, two areas are shown.

In Simpson's approximation, the two areas and their sum can be computed with

x = [0:4];
% or
% N = 2;
% x = linspace(0,4,2*N + 1);
% h = x(2) - x(1);
y = 1 + x.^2/16;
h = 1; % Separation between each blue point; is 1/2 width of each area.
a(1) = (h/3)*(y(1) + 4*y(2) + y(3));
a(2) = (h/3)*(y(3) + 4*y(4) + y(5));
A = sum(a) % Gives exact answer for area = 4 + 4^3/(16*3) = 5.3333


# 2. Problems

## 2.1. Rectangle method hand calculation

By hand, compute the area under the curve y = ex between x = 0 and x = 2 using two rectangles and the "left" rule.

Write a program that uses a for loop to compute the area of this curve and verify that your program gives the correct total area and areas of each rectangle.

Repeat the previous calculation using the "middle" rule.

## 2.2. Rectangle Method I

1. Use all three methods of the rectangle method with width of 0.5 to determine the area under the curve y=x in for x=0-4. First, write your code in long-hand form and then translate the long-hand form to a program that uses a for loop. In all cases, draw a sketch on graph paper and manually compute the expected areas. Use this calculation to verify that your program is producing correct results.
2. For the "left" method, what rectangle width is required for the estimated area to match the actual area to three decimal places (That is, so that the difference between the estimate and the exact area is 0.001 0.01 or less.)

## 2.3. Rectangle Method II

The following program is used to estimate the area under the curve y = x2/4 between x = 0 and x = 5 using the rectangle method:

A = 0;
for i = [1:5]
A = A + (i)^2/4;
end


On a piece of paper, sketch the curve and the rectangles that correspond to the program above. Indicate the numeric value for the height of each rectangle.

## 2.4. Rectangle Method III

The following program is used to estimate the area under the curve y = x2/4 between x = 0 and x = 5 using the rectangle method:

A = 0;
for i = [1:5]
A = A + (i-1/2)^2/4;
end


On a piece of paper, sketch the curve and the rectangles that correspond to the program above. Indicate the numeric value for the height of each rectangle.

## 2.5. Rectangle Method IV

The following program computes the area under a curve using four rectangles between x = 0-4.

xf = 4;
x = [0:xf];
w = x(2) - x(1);
y = 1 + x/4;
for i = 1:length(x)-1
a(i) = w*y(i);
end


Modify this program so that it computes the area under a curve using for rectangles for arbitrary values of xf > 0. To compute an estimate of the area for a different xf but using four rectangles, I should only have to modify the variable xf.

## 2.6. Rectangle Method V

For the curve y = 1 + x + x2/16 in the range of x = 0-4,

• Compute the exact value of the total area under the curve in this range
• Use the "left" rectangle method with two areas to estimate the total area under the curve in this range
• Use the "middle" rectangle method with two areas to estimate the total area under the curve in this range
• Use the "right" rectangle method with two areas to estimate the total area under the curve in this range
• Compute the minimal number of areas that are needed to match the exact value to three decimal places for each method

Sketch the rectangles on a piece of graph paper for each method. Label the heights of each rectangle and the x-values at which the heights were computed.

## 2.7. Rectangle Method VI

The following is the header of a function that specifies how it is called along with a description of what it does. Write the code for this function below the header information.

rectint.m

function A = rectint(xo,xf,N)
% RECTINT Estimate area under y = x^2 using rectangles.
%
%    A = RECTINT(xo,xf,N)
%
%    Estimates the area A under the curve y = x^2 between xo and xf using the midpoint
%    rectangle rule using N rectangles.


Create a file named rectint_test.m that has 5 tests that verify your function will work for arbitrary values of inputs. For example, one test could be

rectint_test.m

xo = 0;
xf = 4;
N  = 4;
A = rectint(xo,xf,N);
Ae = 0.5^2 + 1.5^2 + 2.5^2 + 3.5^2; % Expected value based on hand calculation with xo=0, xf=4, and N=4.
fprintf('xo = %.1f, xf = %.1f, Computed - Expected = %0.8f\n', xo,xf,A-Ae);


When a user executes rectint_test.m and inspects the results displayed, they should be confident that the function has been thoroughly tested and works as indicated by the documentation.

## 2.8. Rectangle Method VII

Write a function named rect that takes an input of the left and right boundaries and the total number of areas to use and returns the estimate of the area under the curve y = ex using the middle rectangle method. To estimate the area in the range of x = 0-1 using two rectangles, I could enter

A = rect(0,1,2)


## 2.9. Rectangle Method VIII

Compute the area under the curve y = ex between x = 0 and x = 2 using two rectangles and the "left" rule.

## 2.10. Rectangle Method IX

Compute the area under the curve y = ex using 22, 24, 26, and 28 areas and compare your answer with the exact value for the area.

## 2.11. Rectangle Method Generalization

Write a program that estimates the area under the curve y = 1 + x + x2/16 in the range of x = a-b using N rectangles. The first four lines of your program should be

clear
N = 2;
a = 0;
b = 4;


and your program should give the correct answer for the estimate for any value of N > 1 and arbitrary values of a and b (assume b > a). To test your program, I will select arbitrary values of N, a, and b. I should not need to modify any other part of your program to get an estimate.

## 2.12. Rectangle Method Generalization

Write a program that will estimate the area under the curve y = x2 in the range 0-1 using the middle rectangle method and N rectangles. Your program should work for any integer value of N greater than zero.

Compute the exact area and save the value in a variable named Aexact

Use this code to create an array A of area estimates associated with N values from 1 through 100.

Plot (A-Aexact)/Aexact versus N with appropriate legend and axis labels.

## 2.13. Trapezoid Method I

Write a program that estimates the area under the curve y = 1 + x2/16 in the range of x = 0-4 using the trapezoid method with four trapezoids:

1. Without using a for loop
2. Using a for loop

then draw a sketch of the four trapezoids on a piece of paper and compute the areas manually to verify your computation. Then,

• Determine the minimum number of trapezoids that are needed for the estimate to match the exact answer within +/- 0.01 or less.

## 2.14. Trapezoid Method II

One form of the Trapezoid method for two areas

 A = w*(0.5*y(1) + y(2) + 0.5*y(3))


For three areas, the formula is

 A = w*(0.5*y(1) + y(2) + y(3) + 0.5*y(4))


Write a for loop that computes the area under the curve y = 1 + x + x2/16 in the range of x = 0-4 using both formulas.

Generalize your program so that it works for 10 areas.

## 2.15. Trapezoid Method IIIa

By hand and on a sheet of paper, use the trapezoid method with two trapezoids to estimate the area under the curve y = 1 + x + x2 in the range x = [0,1].

## 2.16. Trapezoid Method IIIb

For the curve y = 1 + x + x2/16 in the range of x = 0-4, use the trapezoid method with two areas to estimate the total area under the curve.

(To verify your answer, draw out the trapezoids on a piece of paper and compute their areas by hand and compare with the values your code produces.)

## 2.17. Trapezoid Method IV

Use the trapezoid method to estimate the area under the curve ex from x = 0 to x = 3. State how many trapezoids are needed to match the exact answer to three decimal places.

## 2.18. Simpson's Method I

For the curve y = 1 + x + x2/16 in the range of x = 0-4,

• Compute the exact value of the total area under the curve in this range
• Use Simpson's method with two areas to estimate the total area under the curve in this range
• Compute the minimal number of areas that are needed to match the exact value to three decimal places

## 2.19. Simpson's Method II

For the curve y = 1 + x3/81 in the range of x = 0-4,

• Compute the exact value of the total area
• Use Simpson's method with two areas to estimate the area
• Use the Trapezoid method with two areas to estimate the area
• Compute the number of trapezoids that are needed to match the two-area Simpson's method estimate to two decimal places

## 2.20. Simpson's Method III

For the curve y = ex in the range of x = 0-4,

1. Compute the exact value of the integral
2. Estimate the integral using three areas and Simpson's method by hand (using only a pencil, paper, and calculator). Turn in your solution on a sheet of paper.
3. Write a program that does the calculation for 2. using MATLAB.
4. Write a program that computes the area using 100 areas. When your program is executed, the last number displayed should be this estimate.

## 2.21. Simpson's Method IV

For the curve y = 2ex/e4 in the range of x = 0-4,

• Compute the exact value of the total area
• Use the Trapezoid method with n = 21, 22, 23, ..., 28 areas to estimate the area
• Use Simpson's method with n = 21, 22, 23, ..., 28 areas to estimate the area

Hint: For both methods, write out the calculation that you need to make in long-hand form for the first three values of n. For example, for n = 2 the Trapezoid calculation is

x = [0:2:4];
y = 2*exp(x)/exp(4);
A(1) = 2.0*(0.5*y(1) + y(2) + 0.5*y(3))


and for Simpson's method it is

x = [0:1:4];
y = 2*exp(x)/exp(4);
A(1) = (h/3)*(y(1) + 4*y(2) + y(3)   +   y(3) + 4*y(4) + y(5) )


Once you have done this for n = 22 and n = 23, look for a pattern that will allow you to do the calculation for all n using a for loop.

Create a plot showing the errors for the two methods as a function of the number of areas. For example, if errT is an array containing the difference between the exact and Trapezoid-estimated error and errS is an array containing the difference between the exact and Simpson's-estimated error, and E is the exact area,

n = [2^1, 2^2, 2^3, 2^4, 2^5, 2^6, 2^7, 2^8]; % Or n = 2.^[1:8];
figure(1)
loglog(n, errT*100/E);
hold on; grid on;
loglog(n, errS*100/E);
legend('Trapezoid % Error', 'Simpsons % Error');
xlabel('# of areas');

Summary of results for this problem. Note that the above plotting code will result in a different scale on the x-axis, but your results should be the same.

## 2.22. Trapezoid and Simpson's Method Comparison

For the curve y = 1 + x3/81 in the range of x = 0-4,

• Compute the exact value of the total area
• Use Simpson's method with two areas to estimate the area
• Use the Trapezoid method with two areas to estimate the area
• Compute the number of trapezoids that are needed to match the two-area Simpson's method estimate to two decimal places

## 2.23. Simpson's Method Derivation

Consider a polynomial of the form

$\frac{}{}y(x) = a + bx + cx^2$

In terms of a, b, and c, compute

$A = \int_0^2 y(x)dx$

Compare this with a re-write the following formula in terms of a, b, and c

$A = \frac{w}{3}\left(y(0) + 4y(1) + y(2)\right)$

where w = 1.

## 2.24. Integration Functions

Use the MATLAB function integrate integral to estimate the area under the curve y = 10 + x2 from 0-1.

## 2.25. Integration functions

MATLAB has a numerical integration function, integrate.

The function has two input parameters that can be used to control the accuracy of the solution: AbsTol and RelTol. In contrast, the rectangle, trapezoid, and simpson methods have only one parameter to control the accuracy of the solution: the width of each area (or equivalently, the number of areas).

Use integrate function to estimate the integral of various curves for which you can compute the exact answer. Using trial-and-error, develop an understanding of how AbsTol and RelTol affect the accuracy of the solution.

Provide a summary of your results using any combination of words, tables, and figures.

You may also want to include information about how long MATLAB's function took to execute in comparison to the methods we have covered previously. To time an operation, use the functions tic() and toc(). For example, the following program can be used to time various ways of populating an array.

clear
N = 100000;
tic();
for i = 1:N
A(i) = i;
end
fprintf('Time to populate %d element array without pre-allocation: %.5f\n',N,toc());

clear
N = 100000;
tic();
A = zeros(N,1);
for i = 1:N
A(i) = i;
end
fprintf('Time to populate %d element array with pre-allocation: %.5f\n',N,toc());

clear
N = 100000;
tic();
A = [1:N];
fprintf('Time to populate %d element array using colon operator: %.5f\n',N,toc());


## 2.26. Monte Carlo integration

Write a program that uses the Monte Carlo Method to estimate the area under the curve y = 10 + x2 from 0-1. How many sample points are needed in order to match the exact solution to two decimal places?

## 2.27. Physics application: non-ideal spring

A spring has a restoring force F = kx, where k = e − | x | / L and L = 0.1 meters.

Using numerical integration,

1. estimate the work required to compress the spring 0.01 meters.
2. estimate the work required to compress the spring 0.1 meters.

Create a plot that shows the work required to compress the spring as a function of x in the range of 0 to 0.5 and includes a line that shows the work required to compress the spring assuming k = 1. Your plot should have legend labels and properly labeled axes (with units).

## 2.28. Non-ideal spring

A spring has a restoring force F = kx, where k = ko(1 + x2 / L2), ko = 1N / m. and L = 10 meters.

A horizontally traveling 1 kg particle with a speed of 2 m/s strikes this spring.

1. Use a graphical method for estimating the amount of spring compression.
2. Use numerical integration to estimate the amount of compression.

# 3. Activity

## 3.1. Numerical Integration

Use vertical rectangles to estimate the areas of the objects shown on the following hand-outs.

## 3.2. Numerical Integration

 Estimate the area of the blue shape using 4 vertical rectangles. Draw the rectangles and show your calculations. State the algorithm that you used to determine the height of each rectangle. Did your algorithm lead to an over- or under-estimate of the area? If twice as many rectangles are used How many more calculations are required to estimate the area? Will the answer be closer to the analytic value of the area?

## 3.3. Numerical Integration

 Estimate the area of the blue shape using 4 vertical rectangles. Draw the rectangles and show your calculations. State the algorithm that you used to determine the height of each rectangle. Did your algorithm lead to an over- or under-estimate of the area If twice as many rectangles are used How many more calculations are required to estimate the area? Will the answer be closer to the analytic value of the area?

## 3.4. Numerical Integration

In this problem, we are given an equation that represents the height so we can compute the area by hand calculations and by iteration.

Compute the area between the lines y=0, x=3, and y=x2 using rectangles of width 1. Use this graph paper and let 1 unit be the width of a square.

Repeat the above using rectangles of width 1/2 of a square.

Write down your rule for selecting the heights of the rectangles in 1. and 2. What would an advantage be of having a simple rule or a rule that is the same for all rectangles as opposed to a rule that depends on x?

Write a program without a for to compute the area using w=1.

Re-write the program with a for loop.

Answer the last two questions using w = 0.5.

## 3.5. Numerical Integration

In this problem, we are given an equation that represents the height so we can compute the area by hand calculations and by iteration.

• Compute the area between the lines y=0, x=3, and the curve y=0.3333·x3 using rectangles of width 1. Use this graph paper and let 1 unit be the width of a square. Turn in your diagram at the start of class. Include your calculations on your paper. You do not need to put your answers to this question on your wiki page.
• Write a program without a for to compute the area using w=1.
• Write a program with a for to compute the area using w=1.
• Repeat all of the above using using w = 0.5.