{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 21 Newtin-cotes integration formula"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex : 21.1   Pg : 612"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Error Et= 1.468\n",
      "The percent relative error et= 89.467 %\n",
      "The approximate error estimate without using the true value= 2.56\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "def f(x):\n",
    "    y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)\n",
    "    return y\n",
    "tval=1.640533#\n",
    "a=0#\n",
    "b=0.8#\n",
    "fa=f(a)#\n",
    "fb=f(b)#\n",
    "l=(b-a)*((fa+fb)/2)#\n",
    "Et=tval-l##error\n",
    "et=Et*100/tval##percent relative error\n",
    "\n",
    "#by using approximate error estimate\n",
    "\n",
    "#the second derivative of f\n",
    "def g(x):\n",
    "    y=-400+4050*x-10800*x**2+8000*x**3\n",
    "    return y\n",
    "\n",
    "from sympy.mpmath import quad\n",
    "f2x=quad(g,[0,0.8])/(b-a)##average value of second derivative\n",
    "Ea=-(1/12)*(f2x)*(b-a)**3#\n",
    "print \"The Error Et=\",round(Et,3)\n",
    "print \"The percent relative error et=\",round(et,3),\"%\"\n",
    "print \"The approximate error estimate without using the true value=\",Ea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex : 21.2   Pg : 613"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Error Et= 0.572\n",
      "The percent relative error et= 34.85 %\n",
      "The approximate error estimate without using the true value= 0.64\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)\n",
    "    return y\n",
    "a=0#\n",
    "b=0.8#\n",
    "tval=1.640533#\n",
    "n=2#\n",
    "h=(b-a)/n#\n",
    "fa=f(a)#\n",
    "fb=f(b)#\n",
    "fh=f(h)#\n",
    "l=(b-a)*(fa+2*fh+fb)/(2*n)#\n",
    "Et=tval-l##error\n",
    "et=Et*100/tval##percent relative error\n",
    "\n",
    "#by using approximate error estimate\n",
    "\n",
    "#the second derivative of f\n",
    "def g(x):\n",
    "    y=-400+4050*x-10800*x**2+8000*x**3\n",
    "    return y\n",
    "f2x=quad(g,[0,0.8])/(b-a)##average value of second derivative\n",
    "Ea=-(1/12)*(f2x)*(b-a)**3/(n**2)#\n",
    "print \"The Error Et=\",round(Et,3)\n",
    "print \"The percent relative error et=\",round(et,3),\"%\"\n",
    "print \"The approximate error estimate without using the true value=\",Ea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex :21.3    Pg : 614"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No. of segments= 10\n",
      "Segment size= 1.0\n",
      "Estimated d= 288.749146143 m\n",
      "0.237014701487 et(%)\n",
      "---------------------------------------------------------\n",
      "No. of segments= 20\n",
      "Segment size= 0.5\n",
      "Estimated d= 289.263574224 m\n",
      "0.0592795228803 et(%)\n",
      "---------------------------------------------------------\n",
      "No. of segments= 50\n",
      "Segment size= 0.2\n",
      "Estimated d= 298.382319223 m\n",
      "et(%) -3.09125177877\n",
      "---------------------------------------------------------\n",
      "No. of segments= 100\n",
      "Segment size= 0.1\n",
      "Estimated d= 293.915596452 m\n",
      "et(%) -1.54799665905\n",
      "---------------------------------------------------------\n",
      "No. of segments= 100\n",
      "Segment size= 0.1\n",
      "Estimated d= 293.915596452 m\n",
      "et(%) -1.54799665905\n",
      "---------------------------------------------------------\n",
      "No. of segments= 200\n",
      "Segment size= 0.05\n",
      "Estimated d= 289.43343055 m\n",
      "et(%) 0.000594070904571\n",
      "---------------------------------------------------------\n",
      "No. of segments= 200\n",
      "Segment size= 0.05\n",
      "Estimated d= 289.43343055 m\n",
      "et(%) 0.000594070904571\n",
      "---------------------------------------------------------\n",
      "No. of segments= 500\n",
      "Segment size= 0.02\n",
      "Estimated d= 290.332334709 m\n",
      "et(%) -0.309977799375\n",
      "---------------------------------------------------------\n",
      "No. of segments= 1000\n",
      "Segment size= 0.01\n",
      "Estimated d= 289.883809248 m\n",
      "et(%) -0.155012011658\n",
      "---------------------------------------------------------\n",
      "No. of segments= 2000\n",
      "Segment size= 0.005\n",
      "Estimated d= 289.435129352 m\n",
      "et(%) 7.13401428866e-06\n",
      "---------------------------------------------------------\n",
      "No. of segments= 2000\n",
      "Segment size= 0.005\n",
      "Estimated d= 289.435129352 m\n",
      "et(%) 7.13401428866e-06\n",
      "---------------------------------------------------------\n",
      "No. of segments= 5000\n",
      "Segment size= 0.002\n",
      "Estimated d= 289.435143766 m\n",
      "et(%) 2.15393877364e-06\n",
      "---------------------------------------------------------\n",
      "No. of segments= 5000\n",
      "Segment size= 0.002\n",
      "Estimated d= 289.435143766 m\n",
      "et(%) 2.15393877364e-06\n",
      "---------------------------------------------------------\n",
      "No. of segments= 10000\n",
      "Segment size= 0.001\n",
      "Estimated d= 289.480018962 m\n",
      "et(%) -0.0155022506708\n",
      "---------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "g=9.8##m/s**2# acceleration due to gravity\n",
    "m=68.1##kg\n",
    "c=12.5##kg/sec# drag coefficient\n",
    "def f(t):\n",
    "    from numpy import exp\n",
    "    v=g*m*(1-exp(-c*t/m))/c\n",
    "    return v\n",
    "tval=289.43515##m\n",
    "a=0#\n",
    "b=10#\n",
    "fa=f(a)#\n",
    "fb=f(b)#\n",
    "from numpy import arange, exp\n",
    "for i in arange(10,21,10):\n",
    "    n=i#\n",
    "    h=(b-a)/n#\n",
    "    print \"No. of segments=\",i\n",
    "    print \"Segment size=\",h\n",
    "    j=a+h#\n",
    "    s=0#\n",
    "    while j<b:\n",
    "        s=s+f(j)#\n",
    "        j=j+h#\n",
    "    \n",
    "    l=(b-a)*(fa+2*s+fb)/(2*n)#\n",
    "    Et=tval-l##error\n",
    "    et=Et*100/tval##percent relative error\n",
    "    print \"Estimated d=\",l,\"m\"\n",
    "    print et,\"et(%)\"\n",
    "    print \"---------------------------------------------------------\"\n",
    "\n",
    "for i in arange(50,101,50):\n",
    "    n=i#\n",
    "    h=(b-a)/n#\n",
    "    print \"No. of segments=\",i\n",
    "    print \"Segment size=\",h\n",
    "    j=a+h#\n",
    "    s=0#\n",
    "    while j<b:\n",
    "        s=s+f(j)#\n",
    "        j=j+h#\n",
    "    \n",
    "    l=(b-a)*(fa+2*s+fb)/(2*n)#\n",
    "    Et=tval-l##error\n",
    "    et=Et*100/tval##percent relative error\n",
    "    print \"Estimated d=\",l,\"m\"\n",
    "    print \"et(%)\",et\n",
    "    print \"---------------------------------------------------------\"\n",
    "\n",
    "for i in arange(100,201,100):\n",
    "    n=i#\n",
    "    h=(b-a)/n#\n",
    "    print \"No. of segments=\",i\n",
    "    print \"Segment size=\",h\n",
    "    j=a+h#\n",
    "    s=0#\n",
    "    while j<b:\n",
    "        s=s+f(j)#\n",
    "        j=j+h#\n",
    "    \n",
    "    l=(b-a)*(fa+2*s+fb)/(2*n)#\n",
    "    Et=tval-l##error\n",
    "    et=Et*100/tval##percent relative error\n",
    "    print \"Estimated d=\",l,\"m\"\n",
    "    print \"et(%)\",et\n",
    "    print \"---------------------------------------------------------\"\n",
    "\n",
    "for i in arange(200,501,300):\n",
    "    n=i#\n",
    "    h=(b-a)/n#\n",
    "    print \"No. of segments=\",i\n",
    "    print \"Segment size=\",h\n",
    "    j=a+h#\n",
    "    s=0#\n",
    "    while j<b:\n",
    "        s=s+f(j)#\n",
    "        j=j+h#\n",
    "    \n",
    "    l=(b-a)*(fa+2*s+fb)/(2*n)#\n",
    "    Et=tval-l##error\n",
    "    et=Et*100/tval##percent relative error\n",
    "    print \"Estimated d=\",l,\"m\"\n",
    "    print \"et(%)\",et\n",
    "    print \"---------------------------------------------------------\"\n",
    "\n",
    "for i in arange(1000,2001,1000):\n",
    "    n=i#\n",
    "    h=(b-a)/n#\n",
    "    print \"No. of segments=\",i\n",
    "    print \"Segment size=\",h\n",
    "    j=a+h#\n",
    "    s=0#\n",
    "    while j<b:\n",
    "        s=s+f(j)#\n",
    "        j=j+h#\n",
    "    \n",
    "    l=(b-a)*(fa+2*s+fb)/(2*n)#\n",
    "    Et=tval-l##error\n",
    "    et=Et*100/tval##percent relative error\n",
    "    print \"Estimated d=\",l,\"m\"\n",
    "    print \"et(%)\",et\n",
    "    print \"---------------------------------------------------------\"\n",
    "\n",
    "for i in arange(2000,5001,3000):\n",
    "    n=i#\n",
    "    h=(b-a)/n#\n",
    "    print \"No. of segments=\",i\n",
    "    print \"Segment size=\",h\n",
    "    j=a+h#\n",
    "    s=0#\n",
    "    while j<b:\n",
    "        s=s+f(j)#\n",
    "        j=j+h#\n",
    "    \n",
    "    l=(b-a)*(fa+2*s+fb)/(2*n)#\n",
    "    Et=tval-l##error\n",
    "    et=Et*100/tval##percent relative error\n",
    "    print \"Estimated d=\",l,\"m\"\n",
    "    print \"et(%)\",et\n",
    "    print \"---------------------------------------------------------\"\n",
    "\n",
    "for i in arange(5000,10001,5000):\n",
    "    n=i#\n",
    "    h=(b-a)/n#\n",
    "    print \"No. of segments=\",i\n",
    "    print \"Segment size=\",h\n",
    "    j=a+h#\n",
    "    s=0#\n",
    "    while j<b:\n",
    "        s=s+f(j)#\n",
    "        j=j+h#\n",
    "    \n",
    "    l=(b-a)*(fa+2*s+fb)/(2*n)#\n",
    "    Et=tval-l##error\n",
    "    et=Et*100/tval##percent relative error\n",
    "    print \"Estimated d=\",l,\"m\"\n",
    "    print \"et(%)\",et\n",
    "    print \"---------------------------------------------------------\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex : 21.4   Pg : 618"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "l= 1.37\n",
      "The Error Et= 0.27\n",
      "The percent relative error et= 16.645 %\n",
      "The approximate error estimate without using the true value= 0.273\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)\n",
    "    return y\n",
    "a=0#\n",
    "b=0.8#\n",
    "tval=1.640533#\n",
    "n=2#\n",
    "h=(b-a)/n#\n",
    "fa=f(a)#\n",
    "fb=f(b)#\n",
    "fh=f(h)#\n",
    "l=(b-a)*(fa+4*fh+fb)/(3*n)#\n",
    "print\"l=\", round(l,2)\n",
    "Et=tval-l##error\n",
    "et=Et*100/tval##percent relative error\n",
    "\n",
    "#by using approximate error estimate\n",
    "\n",
    "#the fourth derivative of f\n",
    "def g(x):\n",
    "    y=-21600+48000*x\n",
    "    return y\n",
    "\n",
    "f4x=quad(g,[0,0.8])/(b-a)##average value of fourth derivative\n",
    "Ea=-(1/2880)*(f4x)*(b-a)**5#\n",
    "print \"The Error Et=\",round(Et,2)\n",
    "print \"The percent relative error et=\",round(et,3),\"%\"\n",
    "print \"The approximate error estimate without using the true value=\",round(Ea,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex : 23.5   Pg : 620"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "l= 1.62\n",
      "The Error Et= 0.02\n",
      "The percent relative error et= 1.04 %\n",
      "The approximate error estimate without using the true value= 0.017\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)\n",
    "    return y\n",
    "a=0#\n",
    "b=0.8#\n",
    "tval=1.640533#\n",
    "n=4#\n",
    "h=(b-a)/n#\n",
    "fa=f(a)#\n",
    "fb=f(b)#\n",
    "j=a+h#\n",
    "s=0#\n",
    "count=1#\n",
    "while j<b:\n",
    "    if (-1)**count==-1:\n",
    "        s=s+4*f(j)#\n",
    "    else:\n",
    "        s=s+2*f(j)#\n",
    "    \n",
    "    count=count+1#\n",
    "    j=j+h#\n",
    "\n",
    "l=(b-a)*(fa+s+fb)/(3*n)#\n",
    "print\"l=\", round(l,2)\n",
    "Et=tval-l##error\n",
    "et=Et*100/tval##percent relative error\n",
    "\n",
    "#by using approximate error estimate\n",
    "\n",
    "#the fou:rth derivative of f\n",
    "def g(x):\n",
    "    y=-21600+48000*x\n",
    "    return y\n",
    "f4x=quad(g,[0,0.8])/(b-a)##average value of fourth derivative\n",
    "Ea=-(1/(180*4**4))*(f4x)*(b-a)**5#\n",
    "print \"The Error Et=\",round(Et,2)\n",
    "print \"The percent relative error et=\",round(et,3),\"%\"\n",
    "print \"The approximate error estimate without using the true value=\",round(Ea,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex :23.6    Pg : 625"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Part A:\n",
      "l= 1.519\n",
      "The Error Et= 0.12\n",
      "The percent relative error et= 7.398 %\n",
      "The approximate error estimate without using the true value= 0.121\n",
      "---------------------------------------------------\n",
      "Part B:\n",
      "l= 1.645\n",
      "The Error Et= -0.005\n",
      "The percent relative error et= -0.277 %\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)\n",
    "    return y\n",
    "a=0#\n",
    "b=0.8#\n",
    "tval=1.640533#\n",
    "#part a\n",
    "n=3#\n",
    "h=(b-a)/n#\n",
    "fa=f(a)#\n",
    "fb=f(b)#\n",
    "j=a+h#\n",
    "s=0#\n",
    "count=1#\n",
    "while j<b:\n",
    "    s=s+3*f(j)#\n",
    "    count=count+1#\n",
    "    j=j+h#\n",
    "\n",
    "l=(b-a)*(fa+s+fb)/(8)#\n",
    "print \"Part A:\"\n",
    "print \"l=\",round(l,3)\n",
    "Et=tval-l##error\n",
    "et=Et*100/tval##percent relative error\n",
    "\n",
    "#by using approximate error estimate\n",
    "\n",
    "#the fourth derivative of f\n",
    "def g(x):\n",
    "    y=-21600+48000*x\n",
    "    return y\n",
    "f4x=quad(g,[0,0.8])/(b-a)##average value of fourth derivative\n",
    "Ea=-(1/6480)*(f4x)*(b-a)**5#\n",
    "print \"The Error Et=\",round(Et,2)\n",
    "print \"The percent relative error et=\",round(et,3),\"%\"\n",
    "print \"The approximate error estimate without using the true value=\",round(Ea,3)\n",
    "#part b\n",
    "n=5#\n",
    "h=(b-a)/n#\n",
    "l1=(a+2*h-a)*(fa+4*f(a+h)+f(a+2*h))/6#\n",
    "l2=(a+5*h-a-2*h)*(f(a+2*h)+3*(f(a+3*h)+f(a+4*h))+fb)/8#\n",
    "l=l1+l2#\n",
    "print \"---------------------------------------------------\"\n",
    "print \"Part B:\"\n",
    "print \"l=\", round(l,3)\n",
    "Et=tval-l##error\n",
    "et=Et*100/tval##percent relative error\n",
    "print \"The Error Et=\", round(Et,3)\n",
    "print \"The percent relative error et=\", round(et,3), \"%\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex : 23.7   Pg : 626"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "l= 1.59480096\n",
      "The Error Et= 0.04573204\n",
      "The percent relative error et= 2.78763304365 %\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)\n",
    "    return y\n",
    "tval=1.640533#\n",
    "x=[0, 0.12, 0.22, 0.32, 0.36, 0.4 ,0.44 ,0.54 ,0.64 ,0.7 ,0.8]\n",
    "func=[]\n",
    "for i in range(0,11):\n",
    "    func.append(f(x[i]))#\n",
    "\n",
    "l=0#\n",
    "for i in range(0,10):\n",
    "    l=l+(x[i+1]-x[i])*(func[i]+func[i+1])/2#\n",
    "\n",
    "print \"l=\",l\n",
    "Et=tval-l##error\n",
    "et=Et*100/tval##percent relative error\n",
    "print \"The Error Et=\",Et\n",
    "print \"The percent relative error et=\",et,\"%\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex : 23.8   Pg : 230"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "l= 1.60364091733\n",
      "The Error Et= 0.0368920826667\n",
      "The percent relative error et= 2.2487863802 %\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    y=(0.2+25*x-200*x**2+675*x**3-900*x**4+400*x**5)\n",
    "    return y\n",
    "tval=1.640533#\n",
    "x=[0, 0.12, 0.22, 0.32, 0.36, 0.4 ,0.44 ,0.54, 0.64, 0.7, 0.8]\n",
    "func =[]\n",
    "for i in range(0,11):\n",
    "    func.append(f(x[i]))\n",
    "\n",
    "l1=(x[1]-x[0])*((f(x[0])+f(x[1]))/2)#\n",
    "l2=(x[3]-x[1])*(f(x[3])+4*f(x[2])+f(x[1]))/6#\n",
    "l3=(x[6]-x[3])*(f(x[3])+3*(f(x[4])+f(x[5]))+f(x[6]))/8#\n",
    "l4=(x[8]-x[6])*(f(x[6])+4*f(x[7])+f(x[8]))/6\n",
    "l5=(x[9]-x[8])*((f(x[9])+f(x[8]))/2)#\n",
    "l6=(x[10]-x[9])*((f(x[10])+f(x[9]))/2)#\n",
    "l=l1+l2+l3+l4+l5+l6#\n",
    "print \"l=\",l\n",
    "Et=tval-l##error\n",
    "et=Et*100/tval##percent relative error\n",
    "print \"The Error Et=\",Et\n",
    "print \"The percent relative error et=\",et,\"%\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex : 23.9   Pg : 629"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The average termperature is= 53.0\n"
     ]
    }
   ],
   "source": [
    "def f(x,y):\n",
    "    t=2*x*y+2*x-x**2-2*y**2+72\n",
    "    return t\n",
    "Len=8##m,length\n",
    "wid=6##m,width\n",
    "a=0#\n",
    "b=Len#\n",
    "n=2#\n",
    "h=(b-a)/n#\n",
    "a1=0#\n",
    "b1=wid#\n",
    "h1=(b1-a1)/n#\n",
    "\n",
    "fa=f(a,0)#\n",
    "fb=f(b,0)#\n",
    "fh=f(h,0)#\n",
    "lx1=(b-a)*(fa+2*fh+fb)/(2*n)#\n",
    "\n",
    "fa=f(a,h1)#\n",
    "fb=f(b,h1)#\n",
    "fh=f(h,h1)#\n",
    "lx2=(b-a)*(fa+2*fh+fb)/(2*n)#\n",
    "\n",
    "fa=f(a,b1)#\n",
    "fb=f(b,b1)#\n",
    "fh=f(h,b1)#\n",
    "lx3=(b-a)*(fa+2*fh+fb)/(2*n)#\n",
    "\n",
    "l=(b1-a1)*(lx1+2*lx2+lx3)/(2*n)#\n",
    "\n",
    "avg_temp=l/(Len*wid)#\n",
    "print\"The average termperature is=\", avg_temp"
   ]
  }
 ],
 "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
}