{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 12 - Numerical integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_01 Pg No. 373"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "case(a)\n",
      "It =  5\n",
      "Ett =  1.0\n",
      "Iexact =  4.75\n",
      "True error =  0.25 Here Error bound is an overestimate of true error\n",
      "case(b)\n",
      "It =  1.59375\n",
      "Ett =  0.09375\n",
      "Iexact =  1.515625\n",
      "True error =  0.078125 Here Error bound is an overestimate of true error\n"
     ]
    }
   ],
   "source": [
    "from scipy.misc import derivative\n",
    "from mpmath import quad\n",
    "#Trapezoidal Rule\n",
    "\n",
    "def f(x):\n",
    "    F = x**3 + 1\n",
    "    return F\n",
    "\n",
    "#case(a)\n",
    "a = 1#\n",
    "b = 2 #\n",
    "h = b - a #\n",
    "It = (b-a)*(f(a)+f(b))/2\n",
    "d2f = derivative(f,2,n=2)\n",
    "Ett = h**3*d2f/12.0\n",
    "Iexact = quad(f,[1,2])\n",
    "Trueerror = It - Iexact\n",
    "print 'case(a)'\n",
    "print 'It = ',It\n",
    "print 'Ett = ',Ett\n",
    "print 'Iexact = ',Iexact\n",
    "print 'True error = ',Trueerror,\n",
    "print 'Here Error bound is an overestimate of true error'\n",
    "\n",
    "#case(b)\n",
    "a = 1#\n",
    "b = 1.5 #\n",
    "h = b - a #\n",
    "It = (b-a)*(f(a)+f(b))/2\n",
    "Ett = h**3*derivative(f,1.5,n=2)/12\n",
    "Iexact = quad(f,[1,1.5])\n",
    "Trueerror = It - Iexact\n",
    "print 'case(b)'\n",
    "print 'It = ',It\n",
    "print 'Ett = ',Ett\n",
    "print 'Iexact = ',Iexact\n",
    "print 'True error = ',Trueerror,\n",
    "print 'Here Error bound is an overestimate of true error'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_02 Pg No. 376"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intergral for case(a),Ia =  2.54308063482\n",
      "intergral for case(b),Ib =  0.0\n",
      "exact integral,Iexact =  2.3504\n",
      "n = 4 case is better than n = 2 case\n"
     ]
    }
   ],
   "source": [
    "from math import exp\n",
    "#Tapezoidal rule\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    F = exp(x)\n",
    "    return F\n",
    "a = -1 #\n",
    "b = 1 #\n",
    "\n",
    "#case(a)\n",
    "n = 2\n",
    "h = (b-a)/n \n",
    "I = 0\n",
    "for i in range(1,n+1):\n",
    "    I = I + f(a+(i-1)*h)+f(a+i*h)\n",
    "\n",
    "I = h*I/2 #\n",
    "print 'intergral for case(a),Ia = ',I\n",
    "\n",
    "#case(b)\n",
    "n = 4\n",
    "h = (b-a)/n \n",
    "I = 0\n",
    "for i in range(1,n+1):\n",
    "    I = I + f(a+(i-1)*h)+f(a+i*h)#\n",
    "\n",
    "I = h*I/2 #\n",
    "Iexact = 2.35040\n",
    "print 'intergral for case(b),Ib = ',I\n",
    "print 'exact integral,Iexact = ',Iexact\n",
    "print 'n = 4 case is better than n = 2 case'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_03 Pg No. 381"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h =  1\n",
      "Integral for case(a) , Is1 =  2.36205375654\n",
      "h =  0.785398163397\n",
      "Integral for case(b),Is1 =  1.14238405466\n"
     ]
    }
   ],
   "source": [
    "from math import pi,sin,exp,sqrt\n",
    "#Simpon's 1/3 rule\n",
    "\n",
    "#case(a)\n",
    "def f(x):\n",
    "    F = exp(x)\n",
    "    return F\n",
    "a = -1#\n",
    "b = 1#\n",
    "h = (b-a)/2 \n",
    "x1 = a+h\n",
    "Is1 = h*( f(a) + f(b) + 4*f(x1) )/3 \n",
    "print 'h = ',h\n",
    "print 'Integral for case(a) , Is1 = ',Is1\n",
    "\n",
    "#case(b)\n",
    "def f(x):\n",
    "    F = sqrt(sin(x))\n",
    "    return F\n",
    "a = 0\n",
    "b = pi/2\n",
    "h = (b-a)/2 \n",
    "x1 = a+h\n",
    "Is1 = h*( f(a) + f(b) + 4*f(x1) )/3\n",
    "print 'h = ',h\n",
    "print 'Integral for case(b),Is1 = ',Is1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_04 Pg No.382"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h =  0.392699081699\n",
      "Integral value for n = 4 is 1.17822754435\n",
      "h =  0.261799387799\n",
      "Integral value for n = 6 is 1.18728120346\n"
     ]
    }
   ],
   "source": [
    "from math import pi,sin,exp,sqrt\n",
    "#Simpon's 1/3 rule\n",
    "\n",
    "def f(x):\n",
    "    F = sqrt( sin(x) )\n",
    "    return F\n",
    "x0 = 0 #\n",
    "xa = pi/2 #\n",
    "\n",
    "#case(a) n = 4\n",
    "n = 4 #\n",
    "h = ( xa-x0 )/n\n",
    "I = 0 \n",
    "for i in range(1,n/2+1):\n",
    "    I = I + f(x0 + (2*i-2)*h) + 4*f(x0 + (2*i-1)*h) + f(x0 + 2*i*h) #\n",
    "I = h*I/3\n",
    "print 'h = ',h\n",
    "print 'Integral value for n = 4 is',I \n",
    "\n",
    "#case(b) n = 6\n",
    "n = 6\n",
    "h = ( xa-x0 )/n\n",
    "I = 0 \n",
    "for i in range(1,n/2+1):\n",
    "    I = I + f(x0 + (2*i-2)*h) + 4*f(x0 + (2*i-1)*h) + f(x0 + 2*i*h) #\n",
    "\n",
    "I = h*I/3\n",
    "print 'h = ',h\n",
    "print 'Integral value for n = 6 is',I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_05 Pg No. 386"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral of x**3 +1 from 1 to 2 :  0\n",
      "Integral of sqrt(sin(x)) from 0 to pi/2: 1.16104132669\n"
     ]
    }
   ],
   "source": [
    "from math import pi,sin,exp,sqrt\n",
    "\n",
    "#Simpson's 3/8 rule\n",
    "\n",
    "#case(a)\n",
    "def f(x):\n",
    "    F = x**3 + 1\n",
    "    return F\n",
    "a = 1 #\n",
    "b = 2 #\n",
    "h = (b-a)/3 \n",
    "x1 = a + h \n",
    "x2 = a + 2*h\n",
    "Is2 = 3*h*( f(a) + 3*f(x1) + 3*f(x2) + f(b) )/8 #\n",
    "print 'Integral of x**3 +1 from 1 to 2 : ',Is2\n",
    "#case(b)\n",
    "def f(x):\n",
    "    F = sqrt( sin(x) )\n",
    "    return F\n",
    "a = 0 #\n",
    "b = pi/2 #\n",
    "h = (b-a)/3 \n",
    "x1 = a + h \n",
    "x2 = a + 2*h\n",
    "Is2 = 3*h*( f(a) + 3*f(x1) + 3*f(x2) + f(b) )/8 #\n",
    "print 'Integral of sqrt(sin(x)) from 0 to pi/2:',Is2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_06 Pg No. 387"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ib =  1.18061711033\n"
     ]
    }
   ],
   "source": [
    "from math import sin,sqrt,pi\n",
    "#Booles's Five-Point formula\n",
    "\n",
    "def f(x):\n",
    "    F = sqrt( sin(x) )\n",
    "    return F\n",
    "x0 = 0#\n",
    "xb = pi/2 #\n",
    "n = 4 #\n",
    "h = (xb - x0)/n\n",
    "Ib = 2*h*(7*f(x0) + 32*f(x0+h) + 12*f(x0 + 2*h) + 32*f(x0+3*h) + 7*f(x0+4*h))/45#\n",
    "print 'Ib = ',Ib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_07 Pg No. 391"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R(0,0) =  0.75\n",
      "\n",
      "R(1,0) = 0.333333 \n",
      "\n",
      "\n",
      "R(2,0) = 0.342857 \n",
      "\n",
      "\n",
      "R(1,1) = 0.194444 \n",
      "\n",
      "\n",
      "R(2,1) = 0.346032 \n",
      "\n",
      "\n",
      "R(2,2) = 0.356138 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from numpy import zeros\n",
    "def f(x):\n",
    "    F = 1.0/x\n",
    "    return F\n",
    "#since we can't have (0,0) element in matrix we start with (1,1)\n",
    "a = 1 ;\n",
    "b = 2 ;\n",
    "h = b-a ;\n",
    "R=zeros([3,3])\n",
    "R[0,0] = h*(f(a)+f(b))/2 \n",
    "print 'R(0,0) = ',R[0,0]\n",
    "\n",
    "h=[h]\n",
    "for i in range(2,4):\n",
    "    h.append((b-a)/2**(i-1))\n",
    "    s = 0\n",
    "    for k in range(1,2**(i-2)+1):\n",
    "        s = s + f(a + (2*k - 1)*h[i-1]);\n",
    "    \n",
    "    R[i-1,0] = R[i-1,0]/2 + h[i-1]*s;\n",
    "    print '\\nR(%i,0) = %f \\n'%(i-1,R[i-1,0])\n",
    "\n",
    "for j in range(2,4):\n",
    "    for i in range(j,4):\n",
    "        R[i-1,j-1] = (4**(j-1)*R[i-1,j-2] - R[i-2,j-2])/(4**(j-1)-1);\n",
    "        print '\\nR(%i,%i) = %f \\n'%(i-1,j-1,R[i-1,j-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_08 Pg No. 397"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x1 =  -0.57735026919\n",
      "x2 =  0.57735026919\n",
      "I =  2.34269608791\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt,exp \n",
    "#Two Point Gauss -Legedre formula\n",
    "\n",
    "def f(x):\n",
    "    F = exp(x)\n",
    "    return F\n",
    "x1 = -1/sqrt(3)\n",
    "x2 = 1/sqrt(3)\n",
    "I = f(x1) + f(x2)\n",
    "print 'x1 = ',x1\n",
    "print 'x2 = ',x2\n",
    "print 'I = ',I, "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 12_09 Pg No. 398"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " g(z) = exp(-(2.000000*z + 0.000000)/2) \n",
      " C = 2.000000 \n",
      " Ig = 4.685392 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt,exp\n",
    "#Gaussian two point formula\n",
    "\n",
    "a = -2 #\n",
    "b = 2 #\n",
    "def f(x):\n",
    "    F = exp(-x/2)\n",
    "    return F\n",
    "A = (b-a)/2 \n",
    "B = (a+b)/2\n",
    "C = (b-a)/2\n",
    "def g(z):\n",
    "    G = exp(-1*(A*z+B)/2)\n",
    "    return G\n",
    "w1 = 1 #\n",
    "w2 = 1 #\n",
    "z1 = -1/sqrt(3)\n",
    "z2 = 1/sqrt(3)\n",
    "Ig = C*( w1*g(z1) + w2*g(z2) )\n",
    "print ' g(z) = exp(-(%f*z + %f)/2) \\n C = %f \\n Ig = %f \\n'%(A,B,C,Ig)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}