Figure 13: Illustration of the rectangle method with evaluating the rectangle height by either the left or right point.
Compute by hand the area composed of two trapezoids (of equal width) that
approximates the integral \( \int_1^3 2x^3dx \). Make a test function
that calls the trapezoidal function in trapezoidal.m
and compares the return value with the hand-calculated value.
Filename: trapezoidal_test_func.m.
Compute by hand the area composed of two rectangles (of equal width) that
approximates the integral \( \int_1^3 2x^3dx \). Make a test function
that calls the midpoint function in midpoint.m
and compares the return value with the hand-calculated value.
Filename: midpoint_test_func.m.
Apply the trapezoidal and midpoint functions to compute
the integral \( \int_2^6 x(x -1)dx \) with 2 and 100 subintervals.
Compute the error too.
Filename: integrate_parabola.m.
We consider integrating the sine function: \( \int_0^b \sin (x)dx \).
a) Let \( b=\pi \) and use two intervals in the trapezoidal and midpoint method. Compute the integral by hand and illustrate how the two numerical methods approximates the integral. Compare with the exact value.
b) Do a) when \( b=2\pi \).
Filename: integrate_sine.pdf.
Modify the file test_trapezoidal.m such that the
three tests are applied to the function midpoint implementing
the midpoint method for integration.
Filename: test_midpoint.m.
The trapezoidal method integrates linear functions exactly, and this
property was used in the test function test_trapezoidal_linear in the file
test_trapezoidal.m. Change the function used in
the section Proper test procedures to \( f(x)=6\cdot 10^8 x - 4\cdot 10^6 \)
and rerun the test. What happens? How must you change the test
to make it useful? How does the convergence rate test behave? Any need
for adjustment?
Filename: test_trapezoidal2.m.
We want to test how the trapezoidal function works for the integral
\( \int_0^4\sqrt{x}dx \). Two of the tests in test_trapezoidal.m
are meaningful for this integral. Compute by hand the result of
using 2 or 3 trapezoids and modify the test_trapezoidal_one_exact_result
function accordingly. Then modify test_trapezoidal_conv_rate
to handle the square root integral.
Filename: test_trapezoidal3.m.
The convergence rate test fails. Printing out r shows that the
actual convergence rate for this integral is \( -1.5 \) and not \( -2 \).
The reason is that the error in the trapezoidal method
is \( -(b-a)^3n^{-2}f''(\xi) \) for some (unknown) \( \xi\in [a,b] \).
With \( f(x)=\sqrt{x} \), \( f''(\xi)\rightarrow -\infty \) as \( \xi\rightarrow 0 \),
pointing to a potential problem in the size of the error.
Running a test with \( a>0 \), say \( \int_{0.1}^4\sqrt{x}dx \) shows that
the convergence rate is indeed restored to -2.
The midpoint method divides the interval of integration into equal-sized subintervals and approximates the integral in each subinterval by a rectangle whose height equals the function value at the midpoint of the subinterval. Instead, one might use either the left or right end of the subinterval as illustrated in Figure 13. This defines a rectangle method of integration. The height of the rectangle can be based on the left or right end or the midpoint.
Figure 13: Illustration of the rectangle method with evaluating the rectangle height by either the left or right point.
a)
Write a function rectangle(f, a, b, n, height='left') for
computing an integral \( \int_a^bf(x)dx \) by the rectangle method
with height computed based on the value of height, which is either
left, right, or mid.
b)
Write three test functions for the three unit test procedures
described in the section Proper test procedures. Make sure you
test for height equal to left, right, and mid. You may
call the midpoint function for checking the result when height=mid.
Hint.
Edit test_trapezoidal.m.
Filename: rectangle_methods.m.
Suppose we want to use the trapezoidal or midpoint method to compute an integral \( \int_a^b f(x)dx \) with an error less than a prescribed tolerance \( \epsilon \). What is the appropriate size of \( n \)?
To answer this question, we may enter an iterative procedure where we compare the results produced by \( n \) and \( 2n \) intervals, and if the difference is smaller than \( \epsilon \), the value corresponding to \( 2n \) is returned. Otherwise, we halve \( n \) and repeat the procedure.
Hint.
It may be a good idea to organize your code so that the function adaptive_integration can be used
easily in future programs you write.
a) Write a function
adaptive_integration(f, a, b, eps, method=midpoint)
eps corresponds to the tolerance
\( \epsilon \), and method can be midpoint or trapezoidal).
b) Test the method on \( \int_0^2x^2dx \) and \( \int_0^2\sqrt{x}dx \) for \( \epsilon = 10^{-1}, 10^{-10} \) and write out the exact error.
c) Make a plot of \( n \) versus \( \epsilon \in [10^{-1}, 10^{-10}] \) for \( \int_0^2\sqrt{x}dx \). Use logarithmic scale for \( \epsilon \).
Filename: adaptive_integration.m.
The type of method explored in this exercise is called adaptive, because it tries to adapt the value of \( n \) to meet a given error criterion. The true error can very seldom be computed (since we do not know the exact answer to the computational problem), so one has to find other indicators of the error, such as the one here where the changes in the integral value, as the number of intervals is doubled, is taken to reflect the error.
Consider the integral
$$ I = \int_0^4 x^x\,dx\thinspace .$$
The integrand \( x^x \) does not have an anti-derivative that can be expressed in terms of standard functions (visit http://wolframalpha.com and type
integral(x^x,x) to convince yourself that
our claim is right. Note that Wolfram alpha does give you an answer, but that answer is an approximation, it is not exact. This is because
Wolfram alpha too uses numerical methods to arrive at the answer, just as you will in this exercise). Therefore, we are forced to compute the integral by numerical methods. Compute a result that is right to four digits.
Hint. Use ideas from Exercise 3.9: Adaptive integration.
Filename: integrate_x2x.m.
In this exercise we shall integrate $$ I_{j,k} = \int_{-\pi}^{\pi} \sin(jx)\sin(kx)dx,$$ where \( j \) and \( k \) are integers.
a) Plot \( \sin(x)\sin(2x) \) and \( \sin(2x)\sin(3x) \) for \( x\in ]-\pi,\pi] \) in separate plots. Explain why you expect \( \int_{-\pi}^{\pi}\sin x\sin 2x\,dx=0 \) and \( \int_{-\pi}^{\pi}\sin 2x\sin 3x\,dx=0 \).
b) Use the trapezoidal rule to compute \( I_{j,k} \) for \( j=1,\ldots,10 \) and \( k=1,\ldots,10 \).
Filename: products_sines.m.
This is a continuation of Exercise 2.18: Fit sines to straight line. The task is to approximate a given function \( f(t) \) on \( [-\pi,\pi] \) by a sum of sines, $$ \begin{equation} S_N(t) = \sum_{n=1}^{N} b_n \sin(nt)\thinspace . \tag{28} \end{equation} $$ We are now interested in computing the unknown coefficients \( b_n \) such that \( S_N(t) \) is in some sense the best approximation to \( f(t) \). One common way of doing this is to first set up a general expression for the approximation error, measured by "summing up" the squared deviation of \( S_N \) from \( f \): $$ E = \int_{-\pi}^{\pi}(S_N(t)-f(t))^2dt\thinspace .$$ We may view \( E \) as a function of \( b_1,\ldots,b_N \). Minimizing \( E \) with respect to \( b_1,\ldots,b_N \) will give us a best approximation, in the sense that we adjust \( b_1,\ldots,b_N \) such that \( S_N \) deviates as little as possible from \( f \).
Minimization of a function of \( N \) variables, \( E(b_1,\ldots,b_N) \) is mathematically performed by requiring all the partial derivatives to be zero: $$ \begin{align*} \frac{\partial E}{\partial b_1} & = 0,\\ \frac{\partial E}{\partial b_2} & = 0,\\ &\vdots\\ \frac{\partial E}{\partial b_N} & = 0\thinspace . \end{align*} $$
a) Compute the partial derivative \( \partial E/\partial b_1 \) and generalize to the arbitrary case \( \partial E/\partial b_n \), \( 1\leq n\leq N \).
b) Show that $$ b_n = \frac{1}{\pi}\int_{-\pi}^{\pi}f(t)\sin(nt)\,dt\thinspace .$$
c)
Write a function integrate_coeffs(f, N, M) that computes \( b_1,\ldots,b_N \)
by numerical integration, using \( M \) intervals in the trapezoidal rule.
d)
A remarkable property of the trapezoidal rule is that it is exact for integrals
\( \int_{-\pi}^{\pi}\sin nt\,dt \) (when subintervals are of equal size). Use this
property to create a function test_integrate_coeff to verify the
implementation of integrate_coeffs.
e)
Implement the choice \( f(t) = \frac{1}{\pi}t \) as a Matlab function
f(t) and call integrate_coeffs(f, 3, 100) to see what the
optimal choice of \( b_1, b_2, b_3 \) is.
f)
Make a function plot_approx(f, N, M, filename) where you plot f(t)
together with the best approximation \( S_N \) as computed above,
using \( M \) intervals for numerical integration. Save the plot to a file
with name filename.
g)
Run plot_approx(f, N, M, filename) for \( f(t) = \frac{1}{\pi}t \)
for \( N=3,6,12,24 \). Observe how the approximation improves.
h)
Run plot_approx for \( f(t) = e^{-(t-\pi)} \) and \( N=100 \).
Observe a fundamental problem: regardless of \( N \), \( S_N(-\pi)=0 \), not
\( e^{2\pi}\approx 535 \). (There are ways to fix this issue.)
Filename: autofit_sines.m.
Use ideas in the section The midpoint rule for a double integral to derive a formula for computing a double integral \( \int_a^b\int_c^d f(x,y)dydx \) by the trapezoidal rule. Implement and test this rule.
Filename: trapezoidal_double.m.
Use the Monte Carlo method from the section Monte Carlo integration for complex-shaped domains to compute the area of a triangle with vertices at \( (-1,0) \), \( (1,0) \), and \( (3,0) \).
Filename: MC_triangle.m.