Homework 5

Student's Name: PLEASE INSERT YOUR NAME HERE

Directions: Add work to this notebook to solve the problems below.

Check your work. I give some tests that should be passed if you get it correct.

Your notebook should run without printing errors and without user input. If this is not the case, points will be deducted from your grade.

Problem Sources:

  • LL: Programming for Computations - Python by Svein Linge and Hans Petter Langtangen, 2nd edition.
  • L: A Primer on Scientific Programming with Python by Hans Petter Langtangen, 2nd edition.
  • TAK: Applied Scientific Computing With Python by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh.
In [ ]:
# Standard imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv

1. Trapezoid rule

Write a function trapezoid_rule(f, a, b, N) which takes as input a function $f$ defined on $[a,b]$ (where $a$ and $b$ are floats with $a<b$) and an integer $N\geq 1$. The function should return the estimate for $\int_a^b f(x) ~dx$ obtained by applying the trapezoid rule to each interval obtained by dividing $[a,b]$ into $N$ equal sized intervals.

In [ ]:
In [ ]:
# Tests
for N,sol in [(1, 1.9236706937217898e-16),
             (2, 1.5707963267948966),
             (3, 1.8137993642342178),
             (4, 1.8961188979370398),
             (10, 1.9835235375094544),
             (100, 1.999835503887444)]:
    if abs(trapezoid_rule(m.sin, 0, m.pi, N)-sol)>10**-8:
        print("Your trapezoid rule returns an incorrect answer " + \
              "when f(x)=sin(x), a=0, b=pi, and N={}.".format(N))
        print("The correct answer should be {}.".format(sol))

2. Example of use of trapezoid rule

Use your trapezoid_rule function to estimate $$\ln 5 = \int_1^5 \frac{1}{x}~dx.$$

Create a list trap_rule_ln_5 whose entry trap_rule_ln_5[k] contains the result obtained by applying the trapezoid rule with $N=2^k$ equal-sized subintervals of $[1,5]$. The list should have length $11$ and contains values trap_rule_ln_5[k] for $k \in \{0,1,2, \ldots, 10\}.$

Compute errors and observe that the errors decrease by a factor of approximately $1/4$ each time $k$ increases by one.

(This problem is close to TAK exercise 4 from § 3.4.)

In [ ]:

3. Using the three rules

Use the midpoint, trapezoid, and Simpson's rule to estimate $$N_{0,1}(2) = \frac{1}{2} + \frac{1}{\sqrt{2 \pi}} \int_0^2 \exp(\tfrac{-x^2}{2})~dx$$ using $N=10$ subintervals. Store the result in the variables N_mid, N_trap, and N_simp.

You may write your own implementations of these methods or use versions in TAK or in Prof. Hooper's notebook.

(The actual integral evaluates to approximately $0.9772498680518$. Which is most accurate?)

In [ ]:

4. Using data from sensors

A sensor was used to track the speed of a bicyclist during the first five seconds of a race. Measurements times and velocities are stored as a list of pairs $(t,v)$ which are stored in the list t_v_table below. Here times is measured in seconds, and velocity in meters per second.

Use numerical integration to determine the distance she covered in meters. Store the result of your calculation in the variable distance.

In [ ]:
t_v_table = [(0.0, 0.0),
 (0.5, 7.51),
 (1.0, 8.10),
 (1.5, 8.93),
 (2.0, 9.32),
 (2.5, 9.76),
 (3.0, 10.22),
 (3.5, 10.56),
 (4.0, 11.01),
 (4.5, 11.22),
 (5.0, 11.22)]
In [ ]:

5. Error control

Let $f$ be a function with a continuous second derivative on $[a,b]$. Let $T_N$ be the result obtained by applying the Trapesoid rule to estimate the integral of $f$ over $[a,b]$ using $N$ equal-sized subintervals. The error in this integral estimate is then $$E_N = T_N - \int_a^b f(x)~dx.$$ Recall that our error formula guarantees that $$E_N =\frac{-(b-a)^3 f''(y)}{12 N^2}$$ for some $y \in (a,b)$. Therefore $$|E_N| < \frac{(b-a)^3 M}{12 N^2},$$ if $M$ satisfies $|f''(x)|<M$ for all $x \in (a,b)$.

Use this to write a function trap_needed_N(a, b, M, epsilon) which takes as input endpoints of an interval $[a,b]$, an $M>0$, and an $\epsilon>0$, all given as floats. The function should returns the smallest $N$ such that it can be guaranteed that the error satisfies $|E_N|<\epsilon$ for any function $f$ such that $|f''(x)|<M$ for every $x \in (a,b)$.

In [ ]:
In [ ]:
# TESTS:
assert trap_needed_N(-1, 1, 2, 0.1) == 4, \
    "The answer to trap_needed_N(-1, 1, 2, 0.1) should be 4."
assert trap_needed_N(0,10,5,0.01) == 205, \
    "The answer to trap_needed_N(0,10,5,0.01) should be 205."

6: Estimating $\frac{\pi}{12}$.

Consider the unit circle in the plane, $x^2+y^2=1$. A sector representing one twelveth of the unit circle is bounded by the $y$-axis and the line $y=x\sqrt{3}$. Therefore, we have $$\frac{\pi}{12} = \int_{0}^{\frac{1}{2}} \sqrt{1-x^2} - x\sqrt{3}~dx.$$

Use your trap_needed_N function from problem 5 to figure out how many subintervals are needed so that the trapezoid rule yields the correct value of $\frac{\pi}{12}$ up to an error of less than $\epsilon = 0.00001$. Store the result in the variable pi_12_N.

Use your trapezoid_rule function from problem 1 to compute $\frac{\pi}{12}$ to an error less than $\epsilon$. Store the result in the variable pi_12_approx. Check that your answer is as accurate as it should be using m.pi.

To check that code worked, we compute the error in our estimate using m.pi.

In [ ]: