{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "baf256be",
   "metadata": {},
   "source": [
    "Quiz 9 MCS 471, on Friday 28 October 2022, due 10:50am"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2f9c208",
   "metadata": {},
   "source": [
    "Consider ${\\displaystyle \\int_{0}^{\\pi/2} \\cos(x) dx}$.\n",
    "\n",
    "1. Apply the composite trapezoidal rule to approximate this integral, \n",
    "   using 1, 2, 4, and 8 subintervals of $[0,\\pi/2]$.\n",
    "\n",
    "2. Extrapolate to improve the approximations.\n",
    "\n",
    "3. How many decimal places are correct in your extrapolations?  Justify your answer."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f8c27d2",
   "metadata": {},
   "source": [
    "# answer to question 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c236fd3e",
   "metadata": {},
   "source": [
    "We copy the function for the composite Trapezoidal rule from the lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ff7ee88d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "adaptrap"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    function adaptrap(f::Function,a::Float64,b::Float64,n::Int64)\n",
    "\n",
    "Returns a vector of n approximations for the definite integral of\n",
    "the function f over the interval [a,b], the i-th entry in the vector\n",
    "uses 2^i function evaluations.\n",
    "\n",
    "Example:\n",
    "\n",
    "    t = adaptrap(cos,0,pi/2,10)\n",
    "\"\"\"\n",
    "function adaptrap(f::Function,a::Float64,b::Float64,n::Int64)\n",
    "    t = zeros(n)\n",
    "    h = (b-a)        # size of subinterval\n",
    "    m = 1            # number of subintervals\n",
    "    t[1] = (f(a) + f(b))*h/2\n",
    "    for i = 2:n\n",
    "        h = h/2\n",
    "        for j=0:m-1\n",
    "            t[i] = t[i] + f(a+h+j*2*h)\n",
    "        end;\n",
    "        t[i] = t[i-1]/2 + h*t[i]\n",
    "        m = 2*m\n",
    "    end\n",
    "    return t\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ad9dcc1e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Vector{Float64}:\n",
       " 0.7853981633974483\n",
       " 0.9480594489685199\n",
       " 0.9871158009727754\n",
       " 0.9967851718861697"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T = adaptrap(cos,0.0,pi/2,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd1a3daf",
   "metadata": {},
   "source": [
    "# answer to question 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd0ebd8d",
   "metadata": {},
   "source": [
    "Extrapolation to improve the accuracy of approximations to integrals is Romberg integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0be339f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "romberg"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    function romberg(t::Array{Float64,1})\n",
    "\n",
    "Applies extrapolation to the approximations in t,\n",
    "returned by the composite Trapezoidal rule.\n",
    "\n",
    "Example:\n",
    "\n",
    "    t = adaptrap(cos,0,pi/2,6)\n",
    "    et = romberg(t);\n",
    "\"\"\"\n",
    "function romberg(t::Array{Float64,1})\n",
    "    n = length(t)\n",
    "    et = zeros(n,n)\n",
    "    et[:,1] = t\n",
    "    for i = 2:n\n",
    "        for j = 2:i\n",
    "            r = 4^(j-1)\n",
    "            et[i,j] = (r*et[i,j-1] - et[i-1,j-1])/(r - 1);\n",
    "       end\n",
    "    end\n",
    "    return et\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a0374106",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Matrix{Float64}:\n",
       " 0.785398  0.0      0.0       0.0\n",
       " 0.948059  1.00228  0.0       0.0\n",
       " 0.987116  1.00013  0.999992  0.0\n",
       " 0.996785  1.00001  1.0       1.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = romberg(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad7fe63a",
   "metadata": {},
   "source": [
    "# answer to question 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35c2c981",
   "metadata": {},
   "source": [
    "The exact answer equals 1.  We compare the last element in the table to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "84c591c2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0000000081440208"
     ]
    },
    {
     "data": {
      "text/plain": [
       "8.144020791078788e-9"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(R[end,end])\n",
    "R[end,end] - 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c585dc5",
   "metadata": {},
   "source": [
    "We have about 9 decimal places of accuracy."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.0",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}