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

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

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

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

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

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

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

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