Math 451: Computational Mathematics
3.5 Integration
3.5.1 Riemann Sums ¶
Methods in numerical integration arise from the definition of the Riemann Integral
$$
\displaystyle\int_a^b\,f(t)\,dt = \displaystyle\lim_{n\rightarrow\infty} \displaystyle\sum_{k=1}^n f(x_k^*)\Delta x
$$
The Riemann Sum
$$
\displaystyle\sum_{k=1}^n f(x_k^*)\Delta x
$$
is an approximation of the
area under the curve
obtained by
-
partitioning the interval $[a,b]$ into $n$ non-overlapping subintervals $[x_k, x_{k+1}]$
$$ a = x_0 < x_1 < x_2 < \dots < x_{n-1} < x_n = b $$
whose union is all of $[a,b]$. -
choosing a value $x_k^*\in[x_k, x_{k+1}]$ for $1\le k \le n$.
-
evaluating the integrand $f(x)$ at each chosen point to obtain a height of $n$ rectangles.
-
computing the areas of each rectangle $f(x_k^*)\Delta x$, where $\Delta x = \left(\frac{b-a}{n}\right)$
-
summing the areas of the rectangles to obtain the approximation.
riemannsum(@(x) cos(x), 0, pi/4, 10, 'l', 'b', 'FaceAlpha', .7)
The
riemannsum
function used for this notebook can choose the value in each interval of the partition
$$
x_k^*\in[x_k, x_{k+1}]
$$
utilizing one of several methods.
- 'l' left end point of each subinterval (Left Sum)
- 'r' right end point of each subinterval (Right Sum)
- 'm' midpoint of each subinterval (Middle Sum)
- 's' a uniformly distributed random value in each subinterval (True Riemann Sum)
- 't' use both end points to create a trapezoid instead of a rectangle
method = [ 'l', 'r', 'm', 's', 't' ];
name = { 'left sum', 'right sum', 'middle sum', 'Riemann sum', 'trapezoid rule'};
for k=1:5
subplot(2,3,k)
riemannsum(@(x) cos(x), 0, pi/4, 10, method(k), 'b');
title(name{k})
end
3.5.2 Accuracy ¶
Computing the cost or accuracy of these methods of computing Riemann sums depends on the shape of the integrand function. However we can compare the convergence of each method for the same function.
fprintf('%14s %18s %19s %18s %19s\n', ...
'Left Sum', 'Right Sum', 'Middle Sum', 'Riemann Sum', 'Trapezoid Rule')
fprintf('%s\n',repmat('-',1,94))
f = @(x) cos(x);
results = zeros(50,5);
for k=1:4:201
results(k,:) = [ riemannsum(f, 0, pi/2, k, 'l'), ...
riemannsum(f, 0, pi/2, k, 'r'), ...
riemannsum(f, 0, pi/2, k, 'm'), ...
riemannsum(f, 0, pi/2, k, 's'), ...
riemannsum(f, 0, pi/2, k, 't') ];
fprintf('%18.16f %18.16f %18.16f %18.16f %18.16f\n', results(k,:))
end
As we can see, the left sum and right sum do not converge to the limit as quickly as the midpoint sum and trapezoid rule converge.
3.5.3 Adaptive Riemann Sum ¶
Approximating the area with rectangles and trapezoids becomes less accurate with functions that vary their curvature a great deal over our interval of integration $[a,b]$.
f = @(x) sin(x)./x;
riemannsum(f, pi/64, 3*pi, 1000, 't', 'b')
The
adaptive
version of
riemannsum
-
first computes $I_1$ the Riemann sum using only one subinterval $[a,b]$
-
this it recursively completes the following for interval $[a,b]$ passing estimate $I_1$
-
next computes the Riemann sum for both $\left[a, \frac{a+b}{2}\right]$, $I_2$ and $\left[\frac{a+b}{2}, b\right]$, $I_3$.
-
the algorithm compares the values of $|I_1 - (I_2 + I_3)| < $ tolerance.
-
if the values are with the specified tolerance of each other the area approximation is $I_2 + I_3$, otherwise it invokes step 3 twice
- once for interval $\left[a, \frac{a+b}{2}\right]$ passing estimate $I_2$, and
- once for interval $\left[\frac{a+b}{2}, b\right]$ passing estimate $I_3$.
-
recursion continues until a maximum recursion depth is achieved and a warning issued, or the desired estimate is computed.
Adaptive algorithms can be accurate, and efficient. Adaptive methods will accept estimates on wide subintervals over which the function is featureless and continuous. It will only require smaller subintervals where the function is steep or irregular. However they should be used cautiously.
Adaptive methods can be very inefficient if the integrand has a discontinuity in the interval of integration because unnecessary subdivision will occur. One should start out by breaking up such regions into subintervals so that discontinuities occur between subintervals.
Exercise 3.5.1 ¶
Re-write
riemannsum.m
so that the function
riemannsum
has arguments
function area = riemannsum(f, a, b, n, method, tolerance, c, varargin)
riemannsum computes and graphs the Riemann sum of an integrable function
f, on an interval [a,b] using n equal rectangular polygons.
Inputs:
f a handle to an integrable function
a left end point of interval of integration
b right end point of interval of integration
n the maximum depth of recursion
method the method for computing the Riemann sum is a single character
string chosen from 'l','r','m','s','t' for
'left sum', 'right sum', 'middle sum', 'random', 'trapezoid'
option 's' chooses a pseudo-random value in each interval of
the partition of [a,b].
tolerance the required acceptable error in the computation of the area
c color must be a single character string chosen from the list
'r','g','b','c','m','y','w','k', or
an RGB row vector triple, [r g b], the rectangles are
filled with the constant specified color.
varargin the variable argument list modifies one or more Name,Value
argument pairs to set the patch properties of the filled
rectangles of the graph. For example
FaceColor - this is the fill color you set with C
FaceAlpha - alpha blending the faces for transparency [0,1]
EdgeColor - an RGB row vector triple line color specification
LineStyle - for plotting edges
Outputs:
area the area between the curve and the x-axis
You can find the source to
riemannsum
here
. You are adding an new argument
tolerance
to use to stop the algorithm instead of the number of intervals. You are changing the meaning the input argument
n
from the number of subintervals to the
maximum recursion depth
.
You will need to make some simple modifications to the main function
riemannsum.m
along with the addition of the
tolerance
argument.
-
Make sure you set default values for the arguments. You are changing the ordering of the arguments, so pay attention to that. Your default values for the
toleranceshould be1e-20. -
Presently,
riemannsumcallsGraphRiemannSumin the main body to set up the computation, you want to change this to instead call your own function that controls the recursion. For you assignment, you are permitted to ignore the graphing features and completing the assignment should be possible without you interacting withGraphRiemannSumat all. -
You will write a
brand new subfunction
that does the adaptive recursion. This function will use
ComputeRiemannSumto find the area, passing it the correct arguments to return area values. The recursion process will be controlled entirely by this function, which will call bothComputeRiemannSumand itself (name it what you would like) recursively. It should follow the algorithm described above.
Your use of this self-initiated mediated course material is subject to our Creative Commons License .
Creative Commons Attribution-NonCommercial-ShareAlike 4.0
Attribution
You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
Noncommercial
You may not use the material for commercial purposes.
Share Alike
You are free to share, copy and redistribute the material in any medium or format. If you adapt, remix, transform, or build upon the material, you must distribute your contributions under the
same license
as the original.