Homework 4

Student's Name: PLEASE INSERT YOUR NAME HERE (DOUBLE CLICK THIS BOX TO EDIT.)

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. Comparing approximations for the derivatives.

Consider the approximation $$f'(x_0)\approx \frac{f(x_0+h) - f(x_0)}{h}$$ for the case of $f(x)=\ln x$ at the point $x_0=1$. Compute the approximation with $h=10^{-k}$ for $k \in \{1,\ldots, 10\}$. Store the results in a list approx1. (So, approx1[3] should be the computation with $h=10^{-4}$.)

Do the same for the approximation $$f'(x_0)\approx \frac{f(x_0+h) - f(x_0-h)}{2h},$$ storing the result in approx2.

The absolute error is the difference of an approximation from the actual value. Store the list of errors in the first case in error1 and the list of errors in the second case in error2.

Discuss based on your computations what the best approximation is, and what the best choice of $h$ seems to be.

(This problem was based on Problems 1 and 4 from TAK § 3.2)

In [ ]:
# Tests:
assert type(approx1)==list, "Error: approx1 should be a list."
assert len(approx1)==10, "Error: approx1 should have 10 elements."
assert type(approx2)==list, "Error: approx2 should be a list."
assert len(approx2)==10, "Error: approx2 should have 10 elements."
assert type(error1)==list, "Error: error1 should be a list."
assert len(error1)==10, "Error: error1 should have 10 elements."
assert type(error2)==list, "Error: error2 should be a list."
assert len(error2)==10, "Error: error2 should have 10 elements."
assert abs(approx1[0] - m.log(1.1)/0.1) < 10**-8, "approx1[0] is incorrect."

2 Comparing formulas for the second derivative

Repeat exercise 1, this time comparing approximations for the second derivative. For approx1 use $$f''(x_0) \approx \frac{f(x_0+h) - 2 f(x_0) + f(x_0-h)}{h^2}.$$ for approx2 use $$f''(x_0) \approx \frac{-f(x_0+2h) + 16 f(x_0+h) - 30 f(x_0) + 16 f(x_0-h) - f(x_0-2h)}{12 h^2}.$$

Use the same function $f$ same point $x_0$ and same values of $h$ as before. Also compute the corresponding lists of errors error1 and error2.

Explain which approximation is better and what the best choice of $h$ is.

In [ ]:
# Tests:
assert type(approx1)==list, "Error: approx1 should be a list."
assert len(approx1)==10, "Error: approx1 should have 10 elements."
assert type(approx2)==list, "Error: approx2 should be a list."
assert len(approx2)==10, "Error: approx2 should have 10 elements."
assert type(error1)==list, "Error: error1 should be a list."
assert len(error1)==10, "Error: error1 should have 10 elements."
assert type(error2)==list, "Error: error2 should be a list."
assert len(error2)==10, "Error: error2 should have 10 elements."
assert abs(approx1[0] - -1.0050335853501344) < 10**-8, "approx1[0] is incorrect."

3. Tracking a drone

The time_position list below contains the time (in seconds) and position (in meters) of a drone. Estimate its signed velocity and acceleration for times $t \in \{2,4,6,8,10\}$. Place the pairs $(t,v)$ where $v$ is the estimated velocity in order in a list velocity. Similarly, place the pairs $(t,a)$ where $a$ is the estimated acceleration in a list acceleration.

(This is Problem 8 from TAK § 3.2.)

In [ ]:
time_position = [(0,0.7),(2,1.9),(4,3.7),(6,5.3),(8,6.3), \
                 (10,7.4),(12,8.0)]
for t,pos in time_position:
    print("At {} seconds, the drone is at position {} meters." \
         .format(t,pos))
In [ ]:
# insert code solving the problem here
In [ ]:
# For printing out your answers:
for t,v in velocity:
    print("At {} seconds, the drone's velocity is approximately {:.3f} m/s." \
         .format(t,v))
for t,a in acceleration:
    print("At {} seconds, the drone's acceleration is approximately {:.3f} m/s^2." \
         .format(t,a))
In [ ]:
# Tests
assert type(velocity)==list, "Error: velocity should be a list."
assert len(velocity)==5, "Error: velocity should have 5 elements."
assert len(velocity[0])==2, "Error: each member of velocity should be a pair."
assert type(acceleration)==list, "Error: acceleration should be a list."
assert len(acceleration)==5, "Error: acceleration should have 5 elements."
assert len(acceleration[0])==2, "Error: each member of acceleration should be a pair."

4. Farm-raised seafood

More and more seafood is being farm-raised these days. A model used for the rate of change for a fish population, $P(t)$ in farming ponds is given by $$P'(t)=b\left(1-\frac{P(t)}{P_M}\right)P(t) - h P(t),$$ where $b$ is the birth rate, $P_M$ is the maximum number of fish the pond can support, and $h$ is the rate the fish are harvested.

(a) Use forward differences to discretize this problem and write a script that takes in an initial population $P(0)$ as well as all model parameters and outputs (via a plot for example) the population over time.

(b) Demonstrate your model for a carrying capacity $P_M$ of $20,000$ fish, a birth rate of $b = 6\%$ and a harvesting rate of $h = 4\%$ if the pond starts with $5,000$ fish.

(c) Perform some numerical experiments to demonstrate how sensitive this resulting populations are to all of the model parameters in part (b).

Carry all this out in cells below, and describe your results in markdown cells.

(This is Problem 9 from TAK § 3.2.)

5. Finding a quadrature rule

Find the quadrature rule for the interval $[−2, 2]$ using the nodes $-1$, $0$, $1$. (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 [ ]:
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."

6. 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
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."

7. 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 problem 6.

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