Homework 5


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

Be careful to name the objects you create as described in the problem. Some problems may include some tests you can run to test your code.

Before you submit your work, make sure it runs correctly without errors. Click “Kernel > Restart & Run All” to restart and run everything.

Some problems are derived from problems in these books:

  • 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.

Standard imports

In [ ]:
import math as m

1. Finding a quadrature rule

Find the quadrature rule for the interval $[−2, 2]$ using the nodes $-1$, $0$, $1$ with the greatest degree of precision. (That is, we are looking for an integral approximation of the form $$Q(f)=c_- f(-1) + c_0 f(0) + c_+ f(1)$$ whose degree of precision is as larges as possible.)

Write a function quad_rule1(f) which takes a function $f:[-2,2] \to {\mathbb R}$ and returns the integral approximation $Q(f)$.

What this quadrature rule's degree of precision? Store the result in the variable quad_rule1_deg.

(This is a modification of Question 2 from TAK § 3.3.)

In [ ]:

Tests for your code:

In [ ]:
assert str(type(quad_rule1)).split("'")[1]=="function", \
    "Error: quad_rule1 is not a function"
assert type(quad_rule1_deg)==int, "Error: quad_rule1_deg is not an integer."
assert abs(quad_rule1(lambda x:1)-4)<10**-8, \
    "Error, quad_rule1 is not exact for constant functions."
assert abs(quad_rule1(lambda x:x))<10**-8, \
    "Error, quad_rule1 is not exact for linear functions."
assert abs(quad_rule1(lambda x:x**2)-16/3)<10**-8, \
    "Error, quad_rule1 is not exact for quadratic functions."

2. Different estimates for the integral

Write functions midpoint_rule(f, a, b), trapezoid_rule(f, a, b), which simpsons_rule(f, a, b) take as input a function $f$, and floats $a$ and $b$. The functions should return the corresponding integral approximation of $f$ over the interval $[a,b]$.

Apply the rules to approximate $\int_0^1 \frac{1}{1+x^2}~dx$ and compare to the actual value of the integral: $\frac{\pi}{4}$. Store the errors (approximation minus actual value) in the variables midpoint_error, trapezoid_error, and simpsons_error. Compare the results.

(Similar to Question 4 from TAK § 3.3.)

In [ ]:

Tests for your code:

In [ ]:
assert str(type(midpoint_rule)).split("'")[1]=="function", \
    "Error: midpoint_rule is not a function"
assert abs(midpoint_rule(lambda x:x**2,-1,1)) < 10**-8, \
    "Error: midpoint_rule(lambda x:x**2,-1,1) should be zero."
assert str(type(trapezoid_rule)).split("'")[1]=="function", \
    "Error: trapezoid_rule is not a function"
assert abs(trapezoid_rule(lambda x:x**2,-1,1)-2) < 10**-8, \
    "Error: trapezoid_rule(lambda x:x**2,-1,1) should be two."
assert str(type(simpsons_rule)).split("'")[1]=="function", \
    "Error: simpsons_rule is not a function"
assert abs(simpsons_rule(lambda x:x**2,-1,1)-2/3) < 10**-8, \
    "Error: simpsons_rule(lambda x:x**2,-1,1) should be two thirds."

3. Gaussian Quadrature

Find the Gaussian quadrature rule using three nodes in the interval $[0,1]$. Write a function gaussian_rule_3(f) which takes as input a function $f$ and returns the value of the Gaussian quadrature rule applied to $f$.

Check that your rule has degree of precision 5.

Apply this rule to approximate $\int_0^1 \frac{1}{1+x^2}~dx$. Store the error in gaussian_error. Compare the result to the approximations found in the previous problem.

In [ ]:

Tests for your code:

In [ ]:
# Tests
assert str(type(gaussian_rule_3)).split("'")[1]=="function", \
    "Error: gaussian_rule_3 is not a function"
assert abs(gaussian_rule_3(lambda x: x)-0.5) < 10**-8, \
    "Error: gaussian_rule_3 is not exact for the function x."

4. Composite 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 [ ]:

Tests for your code:

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))

5. Example of use of the composite trapezoid rule

Use your trapezoid_rule function from the previous problem 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 absolute errors and observe that the absolute errors decrease by a factor of approximately $1/4$ each time $k$ increases by one. (The absolute error is the absolute value of the difference between the actual value and your approximation.)

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

In [ ]:

Tests for your code:

In [11]:
assert len(trap_rule_ln_5)==11, 'The list trap_rule_ln_5 should have 11 elements.'
for elt in trap_rule_ln_5:
    assert type(elt)==float, 'The elements of trap_rule_ln_5 should be floats'

6. 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 [2]:
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 [ ]:

Tests for your code:

In [ ]:
assert type(distance) == float, 'distance should be a float'

7. 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 Trapezoid 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 [ ]:

Tests for your code:

In [ ]:
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."

8. Estimating pi over 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 the previous 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 earlier in this assignment 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.

In [ ]: