cds-numerical-methods/Week 2/6 Composite Numerical Integration: Trapezoid and Simpson Rules.ipynb

CDS: Numerical Methods Assignments¶

• See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.

• Solutions must be submitted via the Jupyter Hub.

• Make sure you fill in any place that says YOUR CODE HERE or "YOUR ANSWER HERE".

Submission¶

1. Name all team members in the the cell below
2. make sure everything runs as expected
3. restart the kernel (in the menubar, select Kernel$\rightarrow$Restart)
4. run all cells (in the menubar, select Cell$\rightarrow$Run All)
5. Check all outputs (Out[*]) for errors and resolve them if necessary
6. submit your solutions in time (before the deadline)

team_members = ""

Composite Numerical Integration: Trapezoid and Simpson Rules¶

In the following we will implement the composite trapezoid and Simpson rules to calculate definite integrals. These rules are defined by

\begin{align} \int_a^b \, f(x)\, dx &\approx \frac{h}{2} \left[ f(a) + 2 \sum_{j=1}^{n-1} f(x_j) + f(b) \right] &\text{trapezoid} \\ &\approx \frac{h}{3} \left[ f(a) + 2 \sum_{j=1}^{n/2-1} f(x_{2j}) + 4 \sum_{j=1}^{n/2} f(x_{2j-1}) + f(b) \right] &\text{Simpson} \end{align}

with $a = x_0 < x_1 < \dots < x_{n-1} < x_n = b$ and $x_k = a + kh$. Here $k = 0, \dots, n$ and $h = (b-a) / n$ is the step size.

In [ ]:

# Import packages here ...

# YOUR CODE HERE
raise NotImplementedError()

Task 1¶

Implement both integration schemes as Python functions $\text{trapz(yk, dx)}$ and $\text{simps(yk, dx)}$. The argument $\text{yk}$ is an array of length $n+1$ representing $y_k = f(x_k)$ and $\text{dx}$ is the step size $h$. Compare your results with Scipy's functions $\text{scipy.integrate.trapz(yk, xk)}$ and $\text{scipy.integrate.simps(yk, xk)}$ for a $f(x_k)$ of your choice.

Try both even and odd $n$. What do you see? Why?

Hint: go to the Scipy documentation. How are even and odd $n$ handled there?

In [ ]:

def trapz(yk, dx):
    # YOUR CODE HERE
    raise NotImplementedError()
    
def simps(yk, dx):
    # YOUR CODE HERE
    raise NotImplementedError()

In [ ]:

# Compare your results here ...

# YOUR CODE HERE
raise NotImplementedError()

Task 2¶

Implement at least one test function for each of your integration functions.

In [ ]:

def test_trapz():
    # YOUR CODE HERE
    raise NotImplementedError()
    
def test_simps():
    # YOUR CODE HERE
    raise NotImplementedError()
    
test_trapz()
test_simps()

Task 3¶

Study the accuracy of these integration routines by calculating the following integrals for a variety of step sizes $h$:

• $\int_0^1 \, x\, dx$
• $\int_0^1 \, x^2\, dx$
• $\int_0^1 \, x^\frac{1}{2}\, dx$

The integration error is defined as the difference (not the absolute difference) between your numerical results and the exact results. Plot the integration error as a function of $h$ for both integration routines and all listed functions. Comment on the comparison between both integration routines. Does the sign of the error match your expectations? Why?

In [ ]:

# YOUR CODE HERE
raise NotImplementedError()

In [ ]: