{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Homework 5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"INSERT YOUR NAME HERE"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"__Directions:__ Add work to this notebook to solve the problems below.\n",
"\n",
"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.\n",
"\n",
"Before you submit your work, make sure it runs correctly without errors. Click “Kernel > Restart & Run All” to restart and run everything.\n",
"\n",
"Some problems are derived from problems in these books:\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": "markdown",
"metadata": {},
"source": [
"### Standard imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import math as m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Finding a quadrature rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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 \n",
"$$Q(f)=c_- f(-1) + c_0 f(0) + c_+ f(1)$$ \n",
"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": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"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": [
"## 2. Different estimates for the integral"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"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": [
"## 3. Gaussian Quadrature"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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 the previous problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"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.\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Composite Trapezoid Rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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**10**-8:\n",
" print(\"Your trapezoid rule returns an incorrect answer \" + \\\n",
" \"when f(x)=sin(x), a=0, b=pi, and N={}.\".format(N))\n",
" print(\"The correct answer should be {}.\".format(sol))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5. Example of use of the composite trapezoid rule"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use your `trapezoid_rule` function from the previous problem to estimate \n",
"$$\\ln 5 = \\int_1^5 \\frac{1}{x}~dx.$$\n",
"\n",
"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\\}.$\n",
"\n",
"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.)\n",
"\n",
"(This problem is close to TAK exercise 4 from § 3.4.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"assert len(trap_rule_ln_5)==11, 'The list trap_rule_ln_5 should have 11 elements.'\n",
"for elt in trap_rule_ln_5:\n",
" assert type(elt)==float, 'The elements of trap_rule_ln_5 should be floats'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6. Using data from sensors"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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. \n",
"\n",
"Use numerical integration to determine the distance she covered in meters. Store the result of your calculation in the variable `distance`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"t_v_table = [(0.0, 0.0),\n",
" (0.5, 7.51),\n",
" (1.0, 8.10),\n",
" (1.5, 8.93),\n",
" (2.0, 9.32),\n",
" (2.5, 9.76),\n",
" (3.0, 10.22),\n",
" (3.5, 10.56),\n",
" (4.0, 11.01),\n",
" (4.5, 11.22),\n",
" (5.0, 11.22)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Tests for your code:**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"assert type(distance) == float, 'distance should be a float'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 7. Error control"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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 \n",
"$$E_N = T_N - \\int_a^b f(x)~dx.$$\n",
"Recall that our error formula guarantees that \n",
"$$E_N =\\frac{-(b-a)^3 f''(y)}{12 N^2}$$\n",
"for some $y \\in (a,b)$. Therefore \n",
"$$|E_N| < \\frac{(b-a)^3 M}{12 N^2},$$\n",
"if $M$ satisfies $|f''(x)|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)|**