{
"cells": [
{
"cell_type": "code",
"execution_count": 316,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%precision 16\n",
"import numpy\n",
"from scipy.special import erf, fresnel\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Question 1 - Basic Quadrature Rules\n",
"\n",
"The Error Function $\\mathrm{erf}(x)$ is defined by the definite integral\n",
"\n",
"$$\n",
" \\mathrm{erf}(x) = \\frac{2}{\\pi}\\int_0^x e^{-t^2}dt\n",
"$$ \n",
"\n",
"We will use this function to test and understand the different errors associated with a range of both single interval and extended quadrature rules. \n",
"\n",
"Use the function definitions below to write a set of functions for approximating \n",
"$$\n",
" I[f] = \\int_a^b f(x) dx\n",
"$$ \n",
"which we will use for testing against $\\mathrm{erf}(x)$\n",
"\n",
"**(a)** Single step quadrature rules\n",
"\n",
"Write short python functions for calculating $I[f]$ for the following newton-Cotes quadrature rules on a single interval\n",
"\n",
"* Left-Rectangle rule\n",
"* Mid-point rule\n",
"* Trapezoidal rule\n",
"* Simpson's Rule"
]
},
{
"cell_type": "code",
"execution_count": 317,
"metadata": {},
"outputs": [],
"source": [
"def I_left(f, a, b):\n",
" \"\"\" Left rectangle rule\"\"\"\n",
" return (b-a)*f(a)\n",
" \n",
"def I_midpoint(f, a, b):\n",
" \"\"\" Newton-Cotes Midpoint rule\"\"\"\n",
" return (b-a)*f((a+b)/2)\n",
" \n",
"def I_trapezoid(f, a, b):\n",
" \"\"\" Newton-Cotes Midpoint rule\"\"\"\n",
" return (b-a)/2*(f(a)+f(b))\n",
" \n",
"def I_simpson(f, a, b):\n",
" \"\"\" Newton-Cotes Simpsons rule\"\"\"\n",
" return (b-a)/6*(f(a)+4*f((a+b)/2)+f(b))\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 318,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
   "x = 1.0, erf(x) = 0.8427007929497148\n",
   "\n",
   "success!\n"
]
}
],
"source": [
"# test against erf(1.)\n",
"\n",
"x = 1.0\n",
"I_true = erf(x)\n",
"print('x = {}, erf(x) = {}\\n'.format(x,I_true))\n",
"f = lambda x : 2./numpy.sqrt(numpy.pi) * numpy.exp(-x*x)\n",
"\n",

15. The initial value problem for the damped pendulum equation is...