1 files changed, 518 insertions, 0 deletions
diff --git a/Numerical_Methods_by_E._Balaguruswamy/chapter12.ipynb b/Numerical_Methods_by_E._Balaguruswamy/chapter12.ipynb
new file mode 100644
index 00000000..9bd54b00
--- /dev/null
+++ b/Numerical_Methods_by_E._Balaguruswamy/chapter12.ipynb
@@ -0,0 +1,518 @@
+{
+ "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
+}