{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trapezoidal rule when f'' is not finite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next function implements Trapezoid method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Performs Trapezoid quadrature\n",
    "def integrate(a,b,n,f):\n",
    "    h = (b-a)/n\n",
    "    x = np.linspace(a,b,n+1)\n",
    "    y = f(x)\n",
    "    res1 = 0.5*y[0] + np.sum(y[1:n]) + 0.5*y[n]\n",
    "    return h*res1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will integrate\n",
    "$$\n",
    "f(x) = x \\sqrt{x}, \\qquad x \\in [0,1]\n",
    "$$\n",
    "The exact integral is $2/5$ but $f''$ is not finite. Now we perform convergence test. The Peano kernel error formula gives the error estimate\n",
    "$$\n",
    "|E_n(f)| \\le \\frac{3 h^2}{16}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     2    -2.67766953e-02              0     4.68750000e-02\n",
      "     4    -7.01811086e-03         3.8154     1.17187500e-02\n",
      "     8    -1.81246480e-03         3.8721     2.92968750e-03\n",
      "    16    -4.63401302e-04         3.9112     7.32421875e-04\n",
      "    32    -1.17671210e-04         3.9381     1.83105469e-04\n",
      "    64    -2.97398625e-05         3.9567     4.57763672e-05\n",
      "   128    -7.49190899e-06         3.9696     1.14440918e-05\n",
      "   256    -1.88304417e-06         3.9786     2.86102295e-06\n",
      "   512    -4.72540682e-07         3.9849     7.15255737e-07\n",
      "  1024    -1.18449772e-07         3.9894     1.78813934e-07\n"
     ]
    }
   ],
   "source": [
    "f = lambda x: x * np.sqrt(x)\n",
    "a, b = 0.0, 1.0\n",
    "Ie = 2.0/5.0 # Exact integral\n",
    "n, N = 2, 10\n",
    "\n",
    "e = np.zeros(N)\n",
    "for i in range(N):\n",
    "    h = (b-a)/n\n",
    "    I = integrate(a,b,n,f)\n",
    "    e[i] = Ie - I\n",
    "    if i > 0:\n",
    "        print('%6d %18.8e %14.5g %18.8e'% (n,e[i],e[i-1]/e[i],3*h**2/16))\n",
    "    else:\n",
    "        print('%6d %18.8e %14.5g %18.8e'%(n,e[i],0,3*h**2/16))\n",
    "    n = 2*n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We still get $O(h^2)$ convergence and the error is smaller than the error estimate."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}