Math 753/853 HW6 Quadrature (ungraded, don't turn in)

Problem 1

(a) Write a trapezoid(y,h) function that approximate definite integrals using the Trapezoid Rule

\begin{equation*} \int_a^b f(x) dx = h \left(\frac{1}{2} y_0 + \sum_{i=1}^{N-1} y_i + \frac{1}{2} y_N \right) + O(h^2) \end{equation*}

In this formula, the $y_i$ are values $y_i = f(x_i)$ on $N+1$ evenly spaced gridpoints between $a$ and $b$, where $x_i = a + i h$, and $h = (b-a)/N$.

Note that the natural mathematical notation for the Trapezoid rule uses 0-based indexing, but Julia uses 1-based indexing for its vectors. You will have to negotiate this difference when writing your code (i.e. adjust the formula to use 1-based indices).

In [ ]:

(b) Test your trapezoid(y,h) function on the integral $\int_0^1 x \, dx = 1/2$. To do this, choose a smallish value of $N$ (perhaps 10), construct a vector $y$ of length $N+1$, fill it with evelny spaced values of $y=f(x)$ between 0 and 1, and then run trapezoid(y,h) using the appropriate value of $h$. If you don't get 0.5 as an answer, something's wrong!

In [ ]:

Problem 2

Write a wrapper function trapezoid(f, a, b, N) that takes a function $f(x)$, and interval $a,b$ and a discretization number $N$, uses those arguments to construct a vector $y$ and a gridspacing $h$, and then calls trapezoid(y,h). Test it on $f(x) = x, a=0, b=1, N=10$.

In [ ]:

Problem 3

Use the trapezoid(f,a,b,N) function to approximate the integral

\begin{equation*} \int_0^1 x\, e^{2x} \, dx \end{equation*}

using $N=16$ and $N=32$. The exact value of the integral is $(e^2-1)/4$. What is the error in both cases? Does the error scale as expected, $O(h^2)$?

In [ ]:

Problem 4

(a) Write a simpson(y,h) function that approximates definite integrals using Simpson's rule

\begin{equation*} \int_a^b f(x) \, dx = \frac{h}{3} \left(y_0 + 4 \sum_{odd \, i=1}^{N-1} y_i + 2 \sum_{even \, i=2}^{N-2} y_i + y_N\right) + O(h^4) \end{equation*}

As in the Trapezoid Rule, the $y_i$ are values $y_i = f(x_i)$ on $N+1$ evenly spaced gridpoints between $a$ and $b$, where $x_i = a + i h$, and $h = (b-a)/N$.

Remember that for Simpson's Rule $N$ must be odd! And be careful, as you implement this formula written with 0-based indices using 1-based Julia indices, even and odd switch places!

In [ ]:

(b) Test your simpson(y,h) function on the definite integral $\int_0^1 x^2 \, dx = 1/3$ using a smallish value of $N$.

In [ ]:

Problem 5

As in problem 2, write a wrapper function simpson(f, a, b, N) that translates its arguments into appropriate values y,h and then calls simpson(y,h). Make sure it gets the right answer to $\int_0^1 x^2 \, dx = 1/3$.

In [ ]:

Problem 6

Use the simpson(f,a,b,N) function to approximate the integral

\begin{equation*} \int_0^1 x\, e^{2x} \, dx \end{equation*}

using $N=16$ and $N=32$. The exact value of the integral is $(e^2-1)/4$. What is the error in both cases? How do those errors compare to those of the Trapezoid Rule? Does the error scale as expected, $O(h^4)$?

In [2]:

Out[2]:
simpson (generic function with 2 methods)

Problem 7

Make a log-log plot of error versus $N$ for the integral of problem 6, for $N=2^n+1$ for $n=2$ through $10$. Show the error for Trapezoid in red and the error for Simpson in blue. Can you relate the slope of the log-log error lines to the expected errors $O(h^2)$ and $O(h^4)$?

In [ ]: