{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 11 - Numerical differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 11_01 Pg No. 348"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h\ty\terr\n",
      "0.2 \t2.2 \t0.2\n",
      "0.1 \t2.1 \t0.1\n",
      "0.05 \t2.05 \t0.05\n",
      "0.01 \t2.01 \t0.01\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from scipy.misc import derivative\n",
    "\n",
    "#First order forward difference\n",
    "\n",
    "def f(x):\n",
    "    F=x**2\n",
    "    return F\n",
    "                 \n",
    "\n",
    "def df(x,h):\n",
    "    DF=(f(x+h)-f(x))/h\n",
    "    return DF\n",
    "\n",
    "dfactual = derivative(f,1)\n",
    "h = [0.2, 0.1, 0.05, 0.01 ]\n",
    "print 'h\\ty\\terr'\n",
    "for hh in h:\n",
    "    y = df(1,hh)\n",
    "    err = y - dfactual\n",
    "    print hh,'\\t',y,'\\t',err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 11_02 Pg No. 350"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h\ty\terr\n",
      "0.2 \t2.0 \t-4.4408920985e-16\n",
      "0.1 \t2.0 \t4.4408920985e-16\n",
      "0.05 \t2.0 \t4.4408920985e-16\n",
      "0.01 \t2.0 \t1.7763568394e-15\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from scipy.misc import derivative\n",
    "\n",
    "#Three-Point Formula\n",
    "\n",
    "def f(x):\n",
    "    F = x**2\n",
    "    return F\n",
    "def df(x,h):\n",
    "    DF = (f(x+h)-f(x-h))/(2*h)\n",
    "    return DF\n",
    "dfactual = derivative(f,1)\n",
    "h = [0.2, 0.1, 0.05, 0.01 ]\n",
    "print 'h\\ty\\terr'\n",
    "for hh in h:\n",
    "    y = df(1,hh)\n",
    "    err = y - dfactual\n",
    "    print hh,'\\t',y,'\\t',err\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 11_03 Pg No. 353"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h\ty\t\terr\n",
      "0.01 \t0.898257285429 \t-0.00218981692372\n",
      "0.015 \t0.897151155627 \t-0.00329594672573\n",
      "0.02 \t0.896037563392 \t-0.00440953896077\n",
      "0.025 \t0.894916522709 \t-0.00553057964367\n",
      "0.03 \t0.893788047675 \t-0.00665905467759\n",
      "0.035 \t0.892652152498 \t-0.00779494985432\n",
      "0.04 \t0.891508851498 \t-0.00893825085448\n",
      "\n",
      "M2 = 2.061e-08 & hopt = 4.706e-01\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from scipy.misc import derivative\n",
    "from numpy import arange,sin,cos,sqrt\n",
    "x = 0.45#\n",
    "def f(x):\n",
    "    F = sin(x)\n",
    "    return F\n",
    "def df(x,h):\n",
    "    DF = (f(x+h) - f(x))/h\n",
    "    return DF\n",
    "dfactual =  cos(x)#\n",
    "h = arange(0.01,0.005+0.04,0.005)\n",
    "n = len(h)#\n",
    "print 'h\\ty\\t\\terr'\n",
    "for hh in h:\n",
    "    y = df(x,hh)\n",
    "    err = y - dfactual \n",
    "    print hh,'\\t',y,'\\t',err\n",
    "\n",
    "#using 16 significant digits so the bound for roundoff error is 0.5*10**(-16)\n",
    "e = 0.5*10**(-16)\n",
    "M2 = max(sin(x+h))#\n",
    "hopt = 2*sqrt(e/M2)#\n",
    "print '\\nM2 = %.3e & hopt = %.3e'%(hopt,M2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 11_04 Pg No. 354"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " d2fexact = -0.732 \n",
      " err = -6.097e-06 \n",
      " y = -0.732 \n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from numpy import arange,sin,cos,sqrt\n",
    "\n",
    "#Approximate second derivative\n",
    "\n",
    "x = 0.75#\n",
    "h = 0.01#\n",
    "def f(x):\n",
    "    F = cos(x)\n",
    "    return F\n",
    "def d2f(x,h):\n",
    "    D2F = ( f(x+h) - 2*f(x) + f(x-h) )/h**2\n",
    "    return D2F\n",
    "y = d2f(0.75,0.01)#\n",
    "d2fexact = -cos(0.75)\n",
    "err = d2fexact - y #\n",
    "print ' d2fexact = %.3f \\n err = %.3e \\n y = %.3f '%(d2fexact,err,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 11_05 Pg No. 358"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "v(5) =  4.25\n",
      "v(7) =  5.5\n",
      "v(9) =  6.75\n"
     ]
    }
   ],
   "source": [
    "#Differentiation of tabulated data\n",
    "\n",
    "T = range(5,10)\n",
    "s = [10,  14.5,  19.5,  25.5 , 32 ]\n",
    "h = T[1]-T[0]\n",
    "n = len(T)\n",
    "def v(t):\n",
    "    if t in T and T.index(t) == 0:\n",
    "        V = (-3*s[T.index(t) ] + 4*s[T.index(t+h)]  - s[T.index(t+2*h) ] ) /(2*h) #Three point forward difference formula\n",
    "    elif t in T and T.index(t) == n-1:\n",
    "        V = ( 3*s[T.index(t)] - 4*s[ T.index(t-h)]  + s[T.index(t-2*h) ]  )/(2*h) #Backward difference formula\n",
    "    else:\n",
    "        V = ( s[T.index(t+h)]  - s[T.index(t-h)]  )/(2*h) #Central difference formula\n",
    "    return V\n",
    "\n",
    "v_5 = v(5)\n",
    "v_7 = v(7)\n",
    "v_9 = v(9)\n",
    "\n",
    "print 'v(5) = ',v_5\n",
    "print 'v(7) = ',v_7\n",
    "print 'v(9) = ',v_9\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 11_06 Pg No. 359"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a(7) =  0.5\n"
     ]
    }
   ],
   "source": [
    "T = range(5,10)\n",
    "s = [10,  14.5 , 19.5 , 25.5,  32 ]\n",
    "h = T[1]-T[0]\n",
    "def a(t):\n",
    "    A = ( s[T.index(t+h) ] - 2*s[ T.index(t) ] + s[T.index(t-h)  ] )/h**2\n",
    "    return A\n",
    "a_7 = a(6)\n",
    "\n",
    "print 'a(7) = ',a_7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 11_7 Pg No. 359"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solving the above equations we get \n",
      "\n",
      " y1 = y(0.25) = -0.117562 \n",
      " y2 = y(0.5) = -0.168475 \n",
      " y3 = y(0.75) = -0.139088 \n",
      " \n"
     ]
    }
   ],
   "source": [
    "from numpy import mat,divide,linalg\n",
    "h = 0.25 #\n",
    "#y''(x) = e**(x**2)\n",
    "#y(0) = 0 , y(1) = 0\n",
    "# y''(x) = y(x+h) - 2*y(x) + y(x-h)/h**2  = e**(x**2)\n",
    "#( y(x + 0.25) - 2*y(x) + y(x-0.25) )/0.0625 = e**(x**2)\n",
    "#y(x+0.25) - 2*y(x) + y(x - 0.25) = 0.0624*e**(x**2)\n",
    "#y(0.5) - 2*y(0.25) + y(0) = 0.0665\n",
    "#y(0.75) - 2*y(0.5) + y(0.25) = 0.0803\n",
    "#y(1) - 2*y(0.75) + y(0.5) = 0.1097\n",
    "#given y(0) = y(1) = 0\n",
    "#\n",
    "#0 + y2 - 2y1 = 0.06665\n",
    "#y3 - 2*y2 + y1 = 0.0803\n",
    "#-2*y3 + y2 + 0 = 0.1097\n",
    "#Therefore\n",
    "A = mat([[0, 1, -2],[1, -2, 1],[-2, 1, 0 ]])\n",
    "B = mat([[ 0.06665],[0.0803] ,[0.1097 ]])\n",
    "C = linalg.solve(A,B)\n",
    "print 'solving the above equations we get \\n\\n y1 = y(0.25) = %f \\n y2 = y(0.5) = %f \\n y3 = y(0.75) = %f \\n '%(C[2],C[1],C[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 11_08 Pg No. 362"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "richardsons technique -\n",
      "\n",
      "df(0.5) =  1.64850499031\n",
      "D(rh) = D(0.25) =  1.71828182846\n",
      "D(0.5) =  1.66594919985\n",
      "Exact df(0.5) =  1.93757920531\n",
      "The result by richardsons technique is much better than other results\n",
      "for r = 2\n",
      "df(x) =  1.64518270284\n",
      "D(rh) =  1.93757920531\n"
     ]
    }
   ],
   "source": [
    "from numpy import arange,exp\n",
    "from scipy.misc import derivative\n",
    "\n",
    "x = arange(-0.5,0.25+1.5,0.25)\n",
    "h = 0.5 ;\n",
    "r = 1.0/2 ;\n",
    "\n",
    "def f(x):\n",
    "    F = exp(x)\n",
    "    return F\n",
    "def D(x,h):\n",
    "    D = (f(x + h) - f(x-h) )/(2*h) \n",
    "    return D\n",
    "def df(x,h,r):\n",
    "    Df = (D(x,r*h) - r**2*D(x,h))/(1-r**2)\n",
    "    return Df\n",
    "\n",
    "df_05 = df(0.5,0.5,1.0/2);\n",
    "print 'richardsons technique -\\n\\ndf(0.5) = ',df_05\n",
    "print 'D(rh) = D(0.25) = ',D(0.5,0.5)\n",
    "print 'D(0.5) = ',D(0.5,0.25)\n",
    "\n",
    "dfexact = derivative(f,0.5)\n",
    "print 'Exact df(0.5) = ',dfexact\n",
    "print 'The result by richardsons technique is much better than other results'\n",
    "\n",
    "#r = 2\n",
    "print 'for r = 2'\n",
    "print 'df(x) = ',df(0.5,0.5,2)\n",
    "print 'D(rh) = ',D(0.5,2*0.5)"
   ]
  }
 ],
 "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
}