{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework 4\n",
"\n",
"Student's Name: PLEASE INSERT YOUR NAME HERE (DOUBLE CLICK THIS BOX TO EDIT.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Directions:__ Add work to this notebook to solve the problems below. \n",
"\n",
"Check your work. I give some tests that should be passed if you get it correct.\n",
"\n",
"**Your notebook should run without printing errors and without user input. If this is not the case, points will be deducted from your grade.**\n",
"\n",
"Problem Sources:\n",
"* __LL__: *Programming for Computations - Python* by Svein Linge and Hans Petter Langtangen, 2nd edition.\n",
"* __L__: *A Primer on Scientific Programming with Python* by Hans Petter Langtangen, 2nd edition.\n",
"* __TAK__: *Applied Scientific Computing With Python* by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Standard imports\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import math as m\n",
"from mpmath import mp, iv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1. Comparing approximations for the derivatives.\n",
"\n",
"Consider the approximation\n",
"$$f'(x_0)\\approx \\frac{f(x_0+h) - f(x_0)}{h}$$\n",
"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\\}$.\n",
"Store the results in a list `approx1`. (So, `approx1[3]` should be the computation with $h=10^{-4}$.)\n",
"\n",
"Do the same for the approximation \n",
"$$f'(x_0)\\approx \\frac{f(x_0+h) - f(x_0-h)}{2h},$$\n",
"storing the result in `approx2`.\n",
"\n",
"The absolute error is the difference of an approximation from the actual value. Store the list of errors in the first case in `error1`\n",
"and the list of errors in the second case in `error2`.\n",
"\n",
"Discuss based on your computations what the best approximation is, and what the best choice of $h$ seems to be.\n",
"\n",
"(This problem was based on Problems 1 and 4 from TAK § 3.2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Tests:\n",
"assert type(approx1)==list, \"Error: approx1 should be a list.\"\n",
"assert len(approx1)==10, \"Error: approx1 should have 10 elements.\"\n",
"assert type(approx2)==list, \"Error: approx2 should be a list.\"\n",
"assert len(approx2)==10, \"Error: approx2 should have 10 elements.\"\n",
"assert type(error1)==list, \"Error: error1 should be a list.\"\n",
"assert len(error1)==10, \"Error: error1 should have 10 elements.\"\n",
"assert type(error2)==list, \"Error: error2 should be a list.\"\n",
"assert len(error2)==10, \"Error: error2 should have 10 elements.\"\n",
"assert abs(approx1[0] - m.log(1.1)/0.1) < 10**-8, \"approx1[0] is incorrect.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2 Comparing formulas for the second derivative\n",
"\n",
"Repeat exercise 1, this time comparing approximations for the second derivative. For `approx1` use\n",
"$$f''(x_0) \\approx \\frac{f(x_0+h) - 2 f(x_0) + f(x_0-h)}{h^2}.$$\n",
"for `approx2` use \n",
"$$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}.$$\n",
"\n",
"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`.\n",
"\n",
"Explain which approximation is better and what the best choice of $h$ is."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Tests:\n",
"assert type(approx1)==list, \"Error: approx1 should be a list.\"\n",
"assert len(approx1)==10, \"Error: approx1 should have 10 elements.\"\n",
"assert type(approx2)==list, \"Error: approx2 should be a list.\"\n",
"assert len(approx2)==10, \"Error: approx2 should have 10 elements.\"\n",
"assert type(error1)==list, \"Error: error1 should be a list.\"\n",
"assert len(error1)==10, \"Error: error1 should have 10 elements.\"\n",
"assert type(error2)==list, \"Error: error2 should be a list.\"\n",
"assert len(error2)==10, \"Error: error2 should have 10 elements.\"\n",
"assert abs(approx1[0] - -1.0050335853501344) < 10**-8, \"approx1[0] is incorrect.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. Tracking a drone"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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`.\n",
"\n",
"(This is Problem 8 from TAK § 3.2.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"time_position = [(0,0.7),(2,1.9),(4,3.7),(6,5.3),(8,6.3), \\\n",
" (10,7.4),(12,8.0)]\n",
"for t,pos in time_position:\n",
" print(\"At {} seconds, the drone is at position {} meters.\" \\\n",
" .format(t,pos))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# insert code solving the problem here\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# For printing out your answers:\n",
"for t,v in velocity:\n",
" print(\"At {} seconds, the drone's velocity is approximately {:.3f} m/s.\" \\\n",
" .format(t,v))\n",
"for t,a in acceleration:\n",
" print(\"At {} seconds, the drone's acceleration is approximately {:.3f} m/s^2.\" \\\n",
" .format(t,a))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Tests\n",
"assert type(velocity)==list, \"Error: velocity should be a list.\"\n",
"assert len(velocity)==5, \"Error: velocity should have 5 elements.\"\n",
"assert len(velocity[0])==2, \"Error: each member of velocity should be a pair.\"\n",
"assert type(acceleration)==list, \"Error: acceleration should be a list.\"\n",
"assert len(acceleration)==5, \"Error: acceleration should have 5 elements.\"\n",
"assert len(acceleration[0])==2, \"Error: each member of acceleration should be a pair.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4. Farm-raised seafood\n",
"\n",
"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\n",
"$$P'(t)=b\\left(1-\\frac{P(t)}{P_M}\\right)P(t) - h P(t),$$\n",
"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.\n",
"\n",
"**(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.\n",
"\n",
"**(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.\n",
"\n",
"**(c)** Perform some numerical experiments to demonstrate how sensitive this resulting populations are to all of the model parameters in part (b).\n",
"\n",
"**Carry all this out in cells below, and describe your results in markdown cells.**\n",
"\n",
"(This is Problem 9 from TAK § 3.2.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5. Finding a quadrature rule\n",
"\n",
"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.)\n",
"\n",
"Write a function `quad_rule1(f)` which takes a function $f:[-2,2] \\to {\\mathbb R}$ and returns the integral approximation $Q(f)$. \n",
"\n",
"What this quadrature rule's degree of precision? Store the result in the variable `quad_rule1_deg`.\n",
"\n",
"(This is a modification of Question 2 from TAK § 3.3.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"assert str(type(quad_rule1)).split(\"'\")[1]==\"function\", \\\n",
" \"Error: quad_rule1 is not a function\"\n",
"assert type(quad_rule1_deg)==int, \"Error: quad_rule1_deg is not an integer.\"\n",
"assert abs(quad_rule1(lambda x:1)-4)<10**-8, \\\n",
" \"Error, quad_rule1 is not exact for constant functions.\"\n",
"assert abs(quad_rule1(lambda x:x))<10**-8, \\\n",
" \"Error, quad_rule1 is not exact for linear functions.\"\n",
"assert abs(quad_rule1(lambda x:x**2)-16/3)<10**-8, \\\n",
" \"Error, quad_rule1 is not exact for quadratic functions.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 6. Different estimates for the integral\n",
"\n",
"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]$.\n",
"\n",
"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.\n",
"\n",
"(Similar to Question 4 from TAK § 3.3.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Tests\n",
"assert str(type(midpoint_rule)).split(\"'\")[1]==\"function\", \\\n",
" \"Error: midpoint_rule is not a function\"\n",
"assert abs(midpoint_rule(lambda x:x**2,-1,1)) < 10**-8, \\\n",
" \"Error: midpoint_rule(lambda x:x**2,-1,1) should be zero.\"\n",
"assert str(type(trapezoid_rule)).split(\"'\")[1]==\"function\", \\\n",
" \"Error: trapezoid_rule is not a function\"\n",
"assert abs(trapezoid_rule(lambda x:x**2,-1,1)-2) < 10**-8, \\\n",
" \"Error: trapezoid_rule(lambda x:x**2,-1,1) should be two.\"\n",
"assert str(type(simpsons_rule)).split(\"'\")[1]==\"function\", \\\n",
" \"Error: simpsons_rule is not a function\"\n",
"assert abs(simpsons_rule(lambda x:x**2,-1,1)-2/3) < 10**-8, \\\n",
" \"Error: simpsons_rule(lambda x:x**2,-1,1) should be two thirds.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 7. Gaussian Quadrature\n",
"\n",
"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$.\n",
"\n",
"Check that your rule has degree of precision 5.\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Tests\n",
"assert str(type(gaussian_rule_3)).split(\"'\")[1]==\"function\", \\\n",
" \"Error: gaussian_rule_3 is not a function\"\n",
"assert abs(gaussian_rule_3(lambda x: x)-0.5) < 10**-8, \\\n",
" \"Error: gaussian_rule_3 is not exact for the function x.\""
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.9"
}
},
"nbformat": 4,
"nbformat_minor": 2
}