Numerical Integration
From ComputingForScientists
1. Numerical Integration
1.1. References
 Chapter 7 of An Introduction to Programming and Numerical Methods in MATLAB
 Lecture 21 and 22 of Introduction to Numerical Methods and Matlab Programming for Engineers
1.2. Objective
 To introduce the concept of numerical integration
 To have an understanding of the tradeoff 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 xaxis 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 yvalue 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 longhand 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
Longhand 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) 
Shorthand 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 threerectangle case), we first need to request the generation of four points in the range of x=04 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:end1); % 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:end1); % 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 longhand 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
1.12. Monte Carlo Method
2. Problems
2.1. Rectangle method hand calculation
By hand, compute the area under the curve y = e^{x} 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.
Answer 

For the hand calculation: y = e^x initial x = 0, final x = 2 number of rectangles = 2 w = (final x  initial x) / number of rectangles = (2  0) / 2 = 1 reference point for first rectangle: x1 = 0 first y = e^x = e^(0) = 1 reference point for second rectangle: x2 = x1 + w = 0 + 1 = 1 second y = e^x = e^(1) = e Turned into code: x = [0,1]; % Array of the reference points w = 1; % Width of each rectangle for i = [1:length(x)] y(i) = exp(x(i)); % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles Alternatively, using linspace: x = linspace(0,2,3); % Splits the interval between x = 0 and x = 2 into 3 points, the first two of which are our reference points w = x(2)  x(1); for i = [1:length(x)1] y(i) = exp(x(i)); % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles 
Repeat the previous calculation using the "middle" rule.
Answer 

For the hand calculation: y = e^x initial x = 0, final x = 2 number of rectangles = 2 w = (final x  initial x) / number of rectangles = (2  0) / 2 = 1 reference point for first rectangle: x1 = w/2 = 0.5 first y = e^x = e^(0.5) reference point for second rectangle: x2 = x1 + w = 0.5 + 1 = 1.5 second y = e^x = e^(1.5) Turned into code: x = [0.5,1.5]; % Array of the reference points w = 1; % Width of each rectangle for i = [1:length(x)] y(i) = exp(x(i)); % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles Alternatively, using linspace: x = linspace(0,2,3); % Splits the interval between x = 0 and x = 2 into 3 points, the first two of which are our reference points w = x(2)  x(1); x = x + w/2; % Moves the reference points to the middles of the rectangles by adding half the width for i = [1:length(x)1] y(i) = exp(x(i)); % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles 
2.2. Rectangle Method I
 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=04. First, write your code in longhand form and then translate the longhand 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.  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.0010.01 or less.)
Answer 

% Equation of line y = x % For four rectangles, w = 1; format long; clear; w = 1; a(1) = w*0; a(2) = w*1; a(3) = w*2; a(4) = w*3; sum(a) clear; w = 1.0; for i = 1:4 a(i) = w*(i1); end sum(a) % For eight rectangles, w = 0.5; clear; w = 0.5; a(1) = w*0; a(2) = w*0.5; a(3) = w*1.0; a(4) = w*1.5; a(5) = w*2.0; a(6) = w*2.5; a(7) = w*3.0; a(8) = w*3.5; sum(a) % In above, area number (i) is related to height according to % h = (i1)/2.0. So we can write an equivalent for loop: clear; w = 0.5; for i = 1:8 a(i) = w*(i1)/2.0; end sum(a) % Now, compare the two for loops clear; w = 1.0; for i = 1:4 a(i) = w*(i1); end sum(a) clear; w = 0.5; for i = 1:8 a(i) = w*(i1)/2.0; end sum(a) % It appears that w is halved, the number of iterations doubles and the % denominator in the area equation is halved. We can generalize this to clear; n = 100; % or n = 1; w = 1/n; for i = 1:4*n a(i) = w*(i1)/(n); end sum(a) % And verify that when n = 1 and n = 2, we get the same result as the two % for loops for w = 1.0 and w = 0.5. % The manual way of determining the number of iterations required is to % guess value of n until the estimated is greater than or equal to 7.99 % (the exact area is 8.0 and this method always underestimates the exact % answer for this particular curve). Doing so gives % n = 199, w = 1/199, sum(a) = 7.98994975 % n = 200, w = 1/200, sum(a) = 7.99000000 % n = 201, w = 1/201, sum(a) = 7.99004975 % so the required width is w = 1/200 % The automated way of doing this is to use a for loop: clear for n = 1:1000 % Guess the maximum value w = 1/n; for i = 1:4*n a(i) = w*(i1)/n; end fprintf('n = %4d, Estimated area = %.8f\n',n,sum(a)); if abs(sum(a)  8) <= 0.01 break; % Stop execution of outer for loop (loop over n) end end 
2.3. Rectangle Method II
The following program is used to estimate the area under the curve y = x^{2}/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.
Answer 

2.4. Rectangle Method III
The following program is used to estimate the area under the curve y = x^{2}/4 between x = 0 and x = 5 using the rectangle method:
A = 0; for i = [1:5] A = A + (i1/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.
Answer 

2.5. Rectangle Method IV
The following program computes the area under a curve using four rectangles between x = 04.
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 + x^{2}/16 in the range of x = 04,
 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 xvalues at which the heights were computed.
Answer 

The integral of y = 1 + x + x^{2}/16 is y = x + x^{2}/2 + x^{3}/48 Exact value: x = [0,4]; % Lower and upper limits of the x range y = x + (1/2)*x.^2 + (1/48)*x.^3; % Computes the y values of the integral at the lower and upper limit of x a = y(2)  y(1) % Calculates the exact area by taking the difference between the values of the integral at the limits Left rectangles: clear x = linspace(0,4,3); % Splits the interval between x = 0 and x = 2 into 3 points % Or, % x = [0,2,4]; w = x(2)x(1); % Width of each rectangle for i = [1:length(x)1] y(i) = 1 + x(i) + (1/16)*x(i)^2; % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles Middle rectangles: clear x = linspace(0,4,3); % Splits the interval between x = 0 and x = 2 into 3 points % Or, % x = [0,2,4]; w = x(2)x(1); % Width of each rectangle for i = [1:length(x)1] x(i) = x(i) + w/2; % Moves the x values to the middle of each area y(i) = 1 + x(i) + (1/16)*x(i)^2; % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles Right rectangles: clear x = linspace(0,4,3); % Splits the interval between x = 0 and x = 2 into 3 points w = x(2)x(1); % Width of each rectangle % Or, % x = [0,2,4]; for i = [1:length(x)1] y(i) = 1 + x(i+1) + (1/16)*x(i+1)^2; % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles 
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,AAe);
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 = e^{x} using the middle rectangle method. To estimate the area in the range of x = 01 using two rectangles, I could enter
A = rect(0,1,2)
2.9. Rectangle Method VIII
Compute the area under the curve y = e^{x} between x = 0 and x = 2 using two rectangles and the "left" rule.
Answer 

x = [0:2/100:2]; % Array of the reference points w = x(2)x(1); % Width of each rectangle for i = [1:length(x)1] y(i) = exp(x(i)); % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles Using linspace: x = linspace(0,2,101); % Splits the interval between x = 0 and x = 2 into 100 points w = x(2)x(1); % Width of each rectangle for i = [1:length(x)1] y(i) = exp(x(i)); % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles 
2.10. Rectangle Method IX
Compute the area under the curve y = e^{x} using 2^{2}, 2^{4}, 2^{6}, and 2^{8} areas and compare your answer with the exact value for the area.
Answer 

for areas = [2^2,2^4,2^6,2^8] x = [0:2/areas:2]; % Array of the reference points w = x(2)x(1); % Width of each rectangle for i = [1:length(x)1] y(i) = exp(x(i)); % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a) % Sums the total area of all rectangles end Using linspace: for areas = [2^2,2^4,2^6,2^8] x = linspace(0,2,areas+1); % Splits the interval between x = 0 and x = 2 into a number of points depending on the current area w = x(2)x(1); % Width of each rectangle for i = [1:length(x)1] y(i) = exp(x(i)); % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a) % Sums the total area of all rectangles end 
2.11. Rectangle Method Generalization
Write a program that estimates the area under the curve y = 1 + x + x^{2}/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.
Answer 

Using middle rectangle method: clear N = 2; a = 0; b = 4; x = linspace(a,b,N+1); % Splits the interval between x = 0 and x = 2 into 100 points w = x(2)x(1); % Width of each rectangle for i = [1:length(x)1] x(i) = x(i) + w/2; % Moves the x values to the middle of each area y(i) = 1 + x(i) + (1/16)*x(i)^2; % Calculates the y values at the reference points a(i) = w*y(i); % Calculates the area of the given rectangles end A = sum(a); % Sums the total area of all rectangles 
2.12. Rectangle Method Generalization
Write a program that will estimate the area under the curve y = x^{2} in the range 01 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 (AAexact)/Aexact
versus N
with appropriate legend and axis labels.
Answer 

clear % The following is based on the example in the notes % for the "left" rectangle method. The lines % changes are the line after % w = x(2)  x(1) % the line 1+x(i)/4 was changed to x(i).^2 and % the b = 4 was changed to b = 1. 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 x = x + w/2; % Shift evaluation points by w/2. for i = [1:length(x)1] % Evaluate for i = 1,2,3; ignores last value in array x y(i) = 1+x(i).^2; a(i) = w*y(i); end A = sum(a); % To compute areas for values of N = 1:100, wrap the above % in a for loop; change N to k and make A an array. % Note that b = 1 and a = 0 were kept inside the for % loop for presentation clarity. In practice, they % should be placed before the k loop because they % do not depend on the number of rectangles. clear; N = 1:100; for k = 1:length(N) b = 1; a = 0; x = linspace(a,b,k+1); w = x(2)x(1); x = x + w/2; for i = [1:length(x)1] y(i) = x(i).^2; a(i) = w*y(i); end A(k) = sum(a); end % Don't use function integrate for this. % That function gives an estimate, not the exact area. Aexact = 1/3; figure(1);clf;hold on plot(N,(AAexact)/Aexact,'k.','MarkerSize',20); % Plot dots. plot(N,(AAexact)/Aexact,'k'); % Plot lines xlabel('Number of rectangles'); ylabel('Fractional error'); grid on; % To better show the features, use a loglog plot. figure(2);clf; loglog(N,abs(AAexact)/Aexact,'k.','MarkerSize',20); % Plot dots. hold on; % Bug in MATLAB  can't use hold on before first loglog plot command (as was done above). loglog(N,abs(AAexact)/Aexact,'k'); % Plot lines xlabel('Number of rectangles'); ylabel('Fractional error'); title('Integration error using middle rectangle method for y(x) = x^2'); grid on; 
2.13. Trapezoid Method I
Write a program that estimates the area under the curve y = 1 + x^{2}/16 in the range of x = 04 using the trapezoid method with four trapezoids:
 Without using a for loop
 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,
 Compare your estimate with the exact answer for the area.
 Determine the minimum number of trapezoids that are needed for the estimate to match the exact answer within +/ 0.01 or less.
Answer 

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 + x^{2}/16 in the range of x = 04 using both formulas.
Generalize your program so that it works for 10 areas.
Answer 

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 + x^{2} in the range x = [0,1].
Answer 

y1 = 1 y2 = 1+0.5+0.5*0.5 y3 = 1+1.0+1.0*1.0 w = 0.5; a1 = w*0.5*(y2+y1); a2 = w*0.5*(y3+y2); A = a1 + a2 which simplifies to 15/8 
2.16. Trapezoid Method IIIb
For the curve y = 1 + x + x^{2}/16 in the range of x = 04, 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.)
Answer 

clear x = linspace(0,4,3); % Array of the reference points w = x(2)x(1); % Width of each area y(1) = 1 + x(1) + (1/16)*x(1)^2;; % Calculates the y value of the first reference point for i = 1:length(x)1 y(i+1) = 1 + x(i+1) + (1/16)*x(i+1)^2;; % Calculates the y values at the next reference point, which is needed for calculating the triangle's area a(i) = w*y(i) + w*0.5*(y(i+1)y(i)); % Calculates the area of the given rectangles and triangles, adding them together end A = sum(a); % Sums the total area of all rectangles and triangles 
2.17. Trapezoid Method IV
Use the trapezoid method to estimate the area under the curve e^{ − x} from x = 0 to x = 3. State how many trapezoids are needed to match the exact answer to three decimal places.
Answer 

clear; % Exact answer Ae = exp(0)  exp(3); for n = 1:100 x = linspace(0,3,n+1); % Evaluation points. w = x(2)x(1); y(1) = exp(x(1)); for i = 1:length(x)1 y(i+1) = exp(x(i+1)); a(i) = w*0.5*(y(i+1)+y(i)); end A(n) = sum(a); fprintf('N = %d\n',n); fprintf('Exact Area = %.5f\n',Ae); fprintf('Estimated Area = %.5f\n',A(n)); end % First three decimal places are both 0.950 at N = 31. loglog([1:100],(AAe)/Ae); ylabel('(AA_{exact})/A_{exact}'); xlabel('Number of trapezoids'); title('Trapezoid method to estimate e^{x} from 03') grid on; 
2.18. Simpson's Method I
For the curve y = 1 + x + x^{2}/16 in the range of x = 04,
 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
Answer 

Code: [5]. 
2.19. Simpson's Method II
For the curve y = 1 + x^{3}/81 in the range of x = 04,
 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 twoarea Simpson's method estimate to two decimal places
2.20. Simpson's Method III
For the curve y = e^{x} in the range of x = 04,
 Compute the exact value of the integral
 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.
 Write a program that does the calculation for 2. using MATLAB.
 Write a program that computes the area using 100 areas. When your program is executed, the last number displayed should be this estimate.
Answer 

clear Aexact = exp(4)  exp(0); fprintf('Exact area = %.4f\n',Aexact); % For two areas, the code would be x = [0:4]; % 2 areas requires 5 evaluation points. % or % N = 2; % x = linspace(0,4,2*N + 1); % h = x(2)  x(1); y = exp(x); h = 1; 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); fprintf('Simpson'' (2 area) = %.4f\n',A); % For three areas, the code would be x = [0:2/3:4]; % 3 areas requires 10 evaluation points. % or % N = 3; % x = linspace(0,4,2*N + 1); % h = x(2)  x(1); y = exp(x); h = 2/3; a(1) = (h/3)*(y(1) + 4*y(2) + y(3)); a(2) = (h/3)*(y(3) + 4*y(4) + y(5)); a(3) = (h/3)*(y(5) + 4*y(6) + y(7)); A = sum(a); fprintf('Simpson'' (3 area) = %.4f\n',A); % Notice the pattern for the indices for elements multiplied by 4. % They are alway 2x the index value in array a. So the following is % a for loop version of the above. Verify that it % produces the same result by writing the following loop % in longhand when N = 3 and checking that it is the same code as % the above. N = 100; x = linspace(0,4,2*N + 1); h = x(2)  x(1); y = exp(x); for i = 1:N a(i) = (h/3)*(y(2*i1)+4*y(2*i)+y(2*i+1)); end A = sum(a); fprintf('Simpson'' (100 area) = %.4f\n',A); 
2.21. Simpson's Method IV
For the curve y = 2e^{x}/e^{4} in the range of x = 04,
 Compute the exact value of the total area
 Use the Trapezoid method with
n
= 2^{1}, 2^{2}, 2^{3}, ..., 2^{8} areas to estimate the area  Use Simpson's method with
n
= 2^{1}, 2^{2}, 2^{3}, ..., 2^{8} areas to estimate the area
Hint: For both methods, write out the calculation that you need to make in longhand 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
= 2^{2} and n
= 2^{3}, 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 Trapezoidestimated error and errS
is an array containing the difference between the exact and Simpson'sestimated 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');
2.22. Trapezoid and Simpson's Method Comparison
For the curve y = 1 + x^{3}/81 in the range of x = 04,
 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 twoarea Simpson's method estimate to two decimal places
Answer 

Aexact = 4 + 4^4/(81*4); fprintf('Exact area = %.4f\n',Aexact); % The following is from the notes on Simpson's method. % the line y = 1 + x^2/16 was changed to % y = 1 + x^3/16. clear x = [0:4]; % or % N = 2; % x = linspace(0,4,2*N + 1); % h = x(2)  x(1); y = 1 + x.^3/81; h = 1; 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); fprintf('Simpson''s (2 areas) = %.4f\n',A); clear % The following is based on the notes from the trapezoid method using % for trapezoids. To make it apply to two trapzoids, % the line x = [0,1,2,3,4] was changed to % x = [0,2,4]. Also the line % 1+x(1)/4 was changed to 1 + x(i).^3/81 and the line % 1+x(i+1)/4 was changed to 1 + x(i+1).^3/81 x = [0,2,4]; w = x(2)  x(1); y(1) = 1 + x(1).^3/81; for i = 1:length(x)1 y(i+1) = 1 + x(i+1).^3/81; a(i) = w*0.5*(y(i+1)+y(i)); end A = sum(a); fprintf('Trapezoid (2 areas) = %.4f\n',A); % To generalize the above code, we need to create an % array of x values that depends on the number of trapezoids. N = 2; x = linspace(0,4,N+1); % will create the same array x as above. % The following is the same as the previous calculation % except the line x = [0,2,4] has been replaced with % the following two lines. When N = 2, the area result % is the same as when x = [0,2,4] was used. Now change % N until the area matches that for the twoarea Simpson's method, % which gave A = 4.7901. I find that N = 13 gives % A = 4.7948, which rounds to 4.79. N = 13; x = linspace(0,4,N+1); w = x(2)  x(1); y(1) = 1 + x(1).^3/81; for i = 1:length(x)1 y(i+1) = 1 + x(i+1).^3/81; a(i) = w*0.5*(y(i+1)+y(i)); end A = sum(a); fprintf('Trapezoid (%d areas) = %.4f\n',N,A); 
2.23. Simpson's Method Derivation
Consider a polynomial of the form
In terms of a, b, and c, compute
Compare this with a rewrite the following formula in terms of a, b, and c
where w = 1.
2.24. Integration Functions
Use the MATLAB function integrateintegral
to estimate the area under the curve y = 10 + x^{2} from 01.
Answer 

f = @(x) 10 + x.^2 integral(f,0,1) % Gives 10.333333333333332 
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 trialanderror, 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 preallocation: %.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 preallocation: %.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 + x^{2} from 01. How many sample points are needed in order to match the exact solution to two decimal places?
Answer 

Typically 10^{6} samples are needed to get area of 10.33. Note that this is much larger than the number of rectangles or trapezoids that are needed to get this level of accuracy. Monte Carlo integration is only more efficient when the number of dimensions integrated over is large. clear; % Create random points in rectangle % x = 01 % y = 011 % Integral is % (Area of rectangle)*(fraction of points below 10+x^2) % Area of rectangle = 11; N = 100000; x = rand(1,N); % Random values in range 01 y = 11*rand(1,N); % Random values in range 011 ycurve = 10 + x.^2; Ibelow = find(y < ycurve); fbelow = length(Ibelow)/N; A = fbelow*11 
2.27. Physics application: nonideal spring
A spring has a restoring force F = kx, where k = e^{ −  x  / L} and L = 0.1 meters.
Using numerical integration,
 estimate the work required to compress the spring 0.01 meters.
 estimate the work required to compress the spring 0.1 meters.
Compare your answers to the case where k = 1.
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).
Answer 

clear; % Estimate the work required to compress the spring 0.01 m L = 0.1; xf = 0.01; % Final compression. xo = 0.00; % Initial compression. % Use two rectangles x = [xo, xf/2]; dx = x(2)  x(1); % Area of first rectangle w(1) = dx*exp(abs(x(1))/L)*x(1); % Area of second rectangle w(2) = dx*exp(abs(x(2))/L)*x(2); W = sum(w) % Total work. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % To generalize the above (make sure to verify that % when N = 2, you get the same answer as above). % Also add a calculation for the work for k=1. N = 2; x = linspace(xo,xf,N+1); dx = x(2)  x(1); for i = 1:length(x)1 w(i) = dx*exp(abs(x(i))/L)*x(i); wo(i) = dx*x(i); end % With variable k W = sum(w) % With constant k Wo = sum(wo) % Exact answer for k = 1. As N increases (more rectangles), % Wo should approach We. We = (1/2)*(1)*xf^2 % Next we generalize the above code to allow xf to vary. xf = [0.0:0.01:0.5]; N = 2; for j = 1:length(xf) x = linspace(xo,xf(j),N+1); % Note xf now depends on j dx = x(2)  x(1); for i = 1:length(x)1 w(i) = dx*exp(abs(x(i))/L)*x(i); wo(i) = dx*x(i); end % With variable k W(j) = sum(w); % Note W depends on j % With constant k Wo(j) = sum(wo); % Note Wo depends on j end plot(xf,W,'r','LineWidth',2); hold on; grid on; plot(xf,Wo,'b','LineWidth',2); xlabel('x_f [m]'); ylabel('Work [J]'); legend('k = e^{x/L}','k = 1','Location','NorthWest'); 
2.28. Nonideal spring
A spring has a restoring force F = kx, where k = k_{o}(1 + x^{2} / L^{2}), k_{o} = 1N / m. and L = 10 meters.
A horizontally traveling 1 kg particle with a speed of 2 m/s strikes this spring.
 Use a graphical method for estimating the amount of spring compression.
 Use numerical integration to estimate the amount of compression.
Answer 

The energy equation is
PE_{o} = 0 because spring is initially uncompressed. At maximum compression, v_{f} = 0, so KE_{f} = 0. This leaves
When k is constant, the analytical solution is
If k is not constant, note that the energy equation comes from
where W is the work done by any forces on the object as it moves along its path.
Using gives
Doing the integral and assuming x_{o} = 0 gives
So the energy equation
can be written as
This is a 4thorder polynomial. One could use the function
and plot y vs. x_{f}. The x_{f} value where the curve crosses zero is the solution. To solve this, first use the bruteforce method using a constant k and compare with the analytical solution. Then ask what is expected to happen physically when k increases with x. The spring gets stiffer, so we expect the nonconstant k case to result in a smaller x_{f}. The following creates the plot and the zoom tool was used to find . (The zoom tool is available when the figure window is undocked.) This is smaller than that computed for k = constant, which was x_{f} = 2.0. Additional code is given for finding the solution using other methods. L = 10; ko = 1; vo = 2; m = 1; xf = [0:0.0001:3]; y1 = ko*xf.^2  m*vo^2; y2 = ko*(1 + xf.^2/(2*L^2)).*xf.^2  m*vo^2; figure(1);clf plot(xf,y1,'r','LineWidth',2) hold on;grid on; plot(xf,y2,'b','LineWidth',2) xlabel('x_f') ylabel('\Delta KE  W') legend('k=k_o','k=k_o(1+x^2/L^2)') % Additional code below. % To solve this using <code>roots</code>, simplify equation to <math>y=x_f^4/200 + x_f^2  4</math>. r = roots([1/(2*L^2),0,1,0,4]); fprintf('Estimate using roots: %.5f\n',real(r(end))); %Another method of obtaining <math>x_f</math> is to search for the location where the <code>y</code> array transitions from below to above zero: % Find index where crossing. I = find(y2(1:end1) < 0 & y2(2:end) >= 0); %y2(I) % Should be below zero %y2(I+1) % Should be above or equal to zero %xf(I) % Value when y2 below zero %xf(I+1) % Value when y2 above or equal to zero % Estimate as average value before and after crossing zero xf_estimate = (xf(I) + xf(I+1))/2; fprintf('Estimate using energy equation: %.5f\n',xf_estimate); % TODO: Use numerical integration of W with trialanderror guesses for x_f until it matches change in KE. 
3. Activity
3.1. Numerical Integration
Use vertical rectangles to estimate the areas of the objects shown on the following handouts.
 Handout #1: pdf version  svg version
 Handout #2: pdf version  svg version
3.2. Numerical Integration

3.3. Numerical Integration

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=x^{2} using rectangles of width 1. Use this graph paper and let 1 unit be the width of a square.
Answer 

In the image below, three dots were drawn at xvalues of 1, 2, and 3. The yamplitude was computed using the formula y=x^{2} to give yvalues for the dots of 1^{2}=1, 2^{2}=4, and 3^{2}=9. The width of each rectangle is 1, so the total area is the sum of the areas of the three rectangles, or 1*1 + 1*4 + 1*9 = 14. 
Repeat the above using rectangles of width 1/2 of a square.
Answer 

In the image below, three dots were drawn at xvalues of 0.5, 1, 1.5, 2, 2.5, and 3. The yamplitude was computed using the formula y=x^{2} to give yvalues for the dots of 0.5^{2}=0.25, 1^{2}=1, ..., 2.5^{2}=6.25, and 3^{2}=9. The width of each rectangle is 0.5, so the total area is the sum of the areas of the six rectangles, or 0.5*0.25 + 0.5*1 + 0.5*2.25 + 0.5*4+ 0.5*6.25 + 0.5*9 = 11.375. 
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?
Answer 

For the first problem, the rule was to draw a horizontal line from x=0 to x=1 on the xaxis. Then draw a vertical line up to x=1^{2} and draw a dot. The two lines represent two sides of the first rectangle. Next draw a line from x=1 to x=2 on the xaxis. Then draw a vertical line up to x=2^{2} and draw a dot. The two lines represent two sides of the second rectangle. Etc. For second problem, the rule was to draw a horizontal line from x=0 to x=0.5 on the xaxis. Then draw a vertical line up to x=0.5^{2} and draw a dot. The two lines represent two sides of the first rectangle. Next draw a line from x=0.5 to x=1.0 on the xaxis. Then draw a vertical line up to x=1^{2} and draw a dot. The two lines represent two sides of the second rectangle. Etc. 
Write a program without a for
to compute the area using w=1.
Answer 

To figure out how to write a program to do the computation with a for loop, write out the steps in longhand form. clear; w = 1; A(1) = w*1*1; A(2) = w*2*2; A(3) = w*3*3; Atotal = A(1)+A(2)+A(3) now, rewrite the above using an iteration variable. clear; w = 1 i = 1; A(i) = w*i*i; i = 2; A(i) = w*i*i; i = 3; A(i) = w*i*i; Atotal = A(1)+A(2)+A(3) 
Rewrite the program with a for
loop.
Answer 

Both programs above correctly compute the area. The last program can be rewritten as clear; w = 1; for i = [1:3] A(i) = w*i*i; end Atotal = sum(A) Note that the last line was rewritten using the z = sum([1,2,4]) will assign the value 
Answer the last two questions using w = 0.5.
Answer 

Here is a program that almost does what we want. The problem is that the array clear; w = 0.5 i = 0.5; A(i) = w*i*i; i = 1.0; A(i) = w*i*i; i = 1.5; A(i) = w*i*i; i = 2.0; A(i) = w*i*i; i = 2.5; A(i) = w*i*i; i = 3.0; A(i) = w*i*i; Atotal = sum(A) the above (nonvalid) program can be rewritten as clear; w = 0.5 for i = [0.5:0.5:3.0] A(i) = w*i*i end Atotal = sum(A) To make this program work, we need the array clear; w = 0.5 k = 1; for i = [0.5:0.5:3.0] A(k) = w*i*i; k = k+1; end Atotal = sum(A) If you rewrite the above in longhand notation, you will see that it does the correct calculation. 
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·x^{3} 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.
Answer clear; w = 1; for i = 1:3 A(i) = w*0.3333*i*i*i; end Atotal = sum(A) % Gives Atotal = 11.998799999999999
 Write a program with a
for
to compute the area using w=1.
Answer clear; w = 1; A(1) = w*0.3333*1*1*1; A(2) = w*0.3333*2*2*2; A(3) = w*0.3333*3*3*3; Atotal = A(1) + A(2) + A(3) % Gives Atotal = 11.998799999999999
or
clear; w = 1; i = 1; A(i) = w*0.3333*i*i*i; i = 2; A(i) = w*0.3333*i*i*i; i = 3; A(i) = w*0.3333*i*i*i; Atotal = sum(A) % Gives Atotal = 11.998799999999999
 Repeat all of the above using using w = 0.5.
Answer clear; w = 0.5; A(1) = w*0.3333*0.5*0.5*0.5; A(2) = w*0.3333*1*1*1; A(3) = w*0.3333*1.5*1.5*1.5; A(4) = w*0.3333*2*2*2; A(5) = w*0.3333*2.5*2.5*2.5; A(6) = w*0.3333*3*3*3; Atotal = A(1) + A(2) + A(3) + A(4) + A(5) + A(6) % Gives Atotal = 9.186581250000000
or
clear; w = 0.5; i = 1; A(i) = w*0.3333*(i/2)*(i/2)*(i/2); i = 2; A(i) = w*0.3333*(i/2)*(i/2)*(i/2); i = 3; A(i) = w*0.3333*(i/2)*(i/2)*(i/2); i = 4; A(i) = w*0.3333*(i/2)*(i/2)*(i/2); i = 5; A(i) = w*0.3333*(i/2)*(i/2)*(i/2); i = 6; A(i) = w*0.3333*(i/2)*(i/2)*(i/2); Atotal = sum(A) % Gives Atotal = 9.186581250000000
Answer clear w = 0.5; k = 1; for i = [0.5:0.5:3] A(k) = w*0.3333*i*i*i; k = k+1; end Atotal = sum(A) % Gives Atotal = 9.186581250000000
or
clear; w = 0.5; for i = 1:6 A(i) = w*0.3333*(i/2)*(i/2)*(i/2); end Atotal = sum(A) % Gives Atotal = 9.186581250000000
Note Note that the exact answer is
The following program computes the area using 100 rectangles with widths of 1, 1/2, 1/3, ..., 1/100; the first element in the array
Atotal
is for width=1 and the last is for width = 100.clear; for j = 1:100 w(j) = 1/j; k = 1; for i = [w(j):w(j):3.0] A(k) = w(j)*0.3333*i*i*i; k = k+1; end Atotal(j) = sum(A); end Atotal(end) % Atotal(end) is 6.794395492500002 (to be compared with the % exact value of 6.749325000000000)