{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 13 - Numerical solution of ordinary differential equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_01 Pg No. 414"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " y(0.25) =  5.3125\n",
      " y(0.5) =  7.40685424949\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from math import sqrt\n",
    "#Taylor method\n",
    "\n",
    "def f(x,y):\n",
    "    F = x**2 + y**2\n",
    "    return F\n",
    "def d2y(x,y):\n",
    "    D2Y = 2*x + 2*y*f(x,y)\n",
    "    return D2Y\n",
    "def d3y(x,y):\n",
    "    D3Y = 2 + 2*y*d2y(x,y) + 2*f(x,y)**2\n",
    "    return D3Y\n",
    "def y(x):\n",
    "    Y = 1 + f(0,1)*x + d2y(0,1)*x**2/2 + d3y(0,1)*sqrt(x)\n",
    "    return Y\n",
    "print ' y(0.25) = ',y(0.25)\n",
    "print ' y(0.5) = ',y(0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_02 Pg No. 415"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Iteration-1\n",
      "\n",
      " dy(0) = 0.000000\n",
      " d2y(0) = 0.000000\n",
      " d3y(0) = 2.000000\n",
      " d4y(0) = 0.000000\n",
      " \n",
      "y(0) = 0.002667\n",
      "\n",
      "\n",
      " Iteration-2\n",
      "\n",
      " dy(0) = 0.040007\n",
      " d2y(0) = 0.400213\n",
      " d3y(0) = 2.005336\n",
      " d4y(0) = 0.106763\n",
      " \n",
      "y(0) = 0.021353\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from scipy.misc import factorial\n",
    "#Recursive Taylor Method\n",
    "\n",
    "def f(x,y):\n",
    "    F = x**2 + y**2\n",
    "    return F\n",
    "def d2y(x,y):\n",
    "    D2Y = 2*x + 2*y*f(x,y)\n",
    "    return D2Y\n",
    "def d3y(x,y):\n",
    "    D3Y = 2 + 2*y*d2y(x,y) + 2*f(x,y)**2\n",
    "    return D3Y\n",
    "def d4y(x,y):\n",
    "    D4Y = 6*f(x,y)*d2y(x,y) + 2*y*d3y(x,y)\n",
    "    return D4Y\n",
    "h = 0.2 #\n",
    "def y(x,y):\n",
    "    Y = y + f(x,y)*h + d2y(x,y)*h**2/2 + d3y(x,y)*h**3/6 + d4y(x,y)*h**4/factorial(4)\n",
    "    return Y\n",
    "x0 = 0#\n",
    "y0 = 0 #\n",
    "y_=[]\n",
    "for i in range(1,3):\n",
    "    y_.append(y(x0,y0))\n",
    "    print ' Iteration-%d\\n\\n dy(0) = %f\\n d2y(0) = %f\\n d3y(0) = %f\\n d4y(0) = %f\\n '%(i,f(x0,y0),d2y(x0,y0),d3y(x0,y0),d4y(x0,y0)) \n",
    "    x0 = x0 + i*h\n",
    "    y0 = y_[i-1]\n",
    "    print 'y(0) = %f\\n\\n'%(y_[i-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_03 Pg No. 417"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "for dy(x) = x**2 + y**2 the results are :\n",
      "y(0.1) =  0.000333334920635\n",
      "y(0.2) =  0.00266686984127\n",
      "y(1) =  0.349206349206\n",
      "for dy(x) = x*e**y the results are \n",
      "y(0.1) =  0.0050125208594\n",
      "y(0.2) =  0.0202013400268\n",
      "y(1) =  0.6487212707\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from math import exp \n",
    "#Picard's Method\n",
    "\n",
    "#y'(x) = x**2 + y**2,y(0) = 0\n",
    "#y(1) = y0 + integral(x**2 + y0**2,x0,x)\n",
    "#y(1) = x**3/3\n",
    "#y(2) = 0 + integral(xY2 + y1**2,x0,x)\n",
    "#     = integral(x**2 + x**6/9,0,x) = x**3/3 + x**7/63\n",
    "#therefore y(x) = x**3 /3 + x**7/63\n",
    "def y(x):\n",
    "    Y = x**3/3 + x**7/63 \n",
    "    return Y\n",
    "print 'for dy(x) = x**2 + y**2 the results are :'\n",
    "print 'y(0.1) = ',y(0.1)\n",
    "print 'y(0.2) = ',y(0.2)\n",
    "print 'y(1) = ',y(1)\n",
    "\n",
    "#y'(x) = x*e**y, y(0) = 0\n",
    "#y0 = 0 , x0 = 0\n",
    "#Y(1) = 0 + integral(x*e**0,0,x) = x**2/2\n",
    "#y(2) = 0 + integral( x*e**( x**2/2 ) ,0,x) = e**(x**2/2)-1\n",
    "#therefore y(x) = e**(x**2/2) - 1\n",
    "def y(x):\n",
    "    Y = exp(x**2/2) - 1 \n",
    "    return Y\n",
    "print 'for dy(x) = x*e**y the results are '                 \n",
    "print 'y(0.1) = ',y(0.1)\n",
    "print 'y(0.2) = ',y(0.2)\n",
    "print 'y(1) = ',y(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_04 Pg No. 417"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "for h = 0.500000\n",
      "\n",
      "y(1.500000) = 4.000000\n",
      "\n",
      "y(2.000000) = 7.875000\n",
      "\n",
      "\n",
      "for h = 0.250000\n",
      "\n",
      "y(1.250000) = 3.000000\n",
      "\n",
      "y(1.500000) = 4.421875\n",
      "\n",
      "y(1.750000) = 6.359375\n",
      "\n",
      "y(2.000000) = 8.906250\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "\n",
    "#Euler's Method\n",
    "\n",
    "def dy(x):\n",
    "    DY = 3*x**2 + 1\n",
    "    return DY\n",
    "x0 = 1\n",
    "y = [0,2] #\n",
    "#h = 0.5\n",
    "h = 0.5\n",
    "print 'for h = %f\\n'%(h)\n",
    "for i in range(2,4):\n",
    "    y.append(y[i-1] + h*dy(x0+(i-2)*h))\n",
    "    print 'y(%f) = %f\\n'%(x0+(i-1)*h,y[i])\n",
    "\n",
    "#h = 0.25\n",
    "h = 0.25\n",
    "print '\\nfor h = %f\\n'%(h)\n",
    "y = [0,2] #\n",
    "for i in range(2,6):\n",
    "    y.append(y[i-1] + h*dy(x0+(i-2)*h))\n",
    "    print 'y(%f) = %f\\n'%(x0+(i-1)*h,y[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_05 Pg No. 422"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " Step 1 \n",
      " x(1) = 1.500000\n",
      " y(1.500000) = 4.000000\n",
      "\n",
      "\n",
      " Et(1) = 0.875000\n",
      "\n",
      "\n",
      " Step 2 \n",
      " x(2) = 2.000000\n",
      " y(2.000000) = 7.875000\n",
      "\n",
      "\n",
      " Et(2) = 1.250000\n",
      "\n",
      "x      Est y   true y    Et      Globalerr\n",
      "1.5 \t4.0 \t4.875 \t0.875 \t0.875\n",
      "2.0 \t7.875 \t10.0 \t1.25 \t2.125\n"
     ]
    }
   ],
   "source": [
    "from numpy import nditer\n",
    "from __future__ import division\n",
    "#Error estimation in Euler's Method\n",
    "\n",
    "def dy(x):\n",
    "    DY = 3*x**2 + 1\n",
    "    return DY\n",
    "def d2y(x):\n",
    "    D2Y = 6*x\n",
    "    return D2Y\n",
    "def d3y(x):\n",
    "    D3Y = 6\n",
    "    return D3Y\n",
    "def exacty(x):\n",
    "    Exacty = x**3 + x\n",
    "    return Exacty\n",
    "x0 = 1\n",
    "y = 2\n",
    "h = 0.5\n",
    "X=[];ESTY=[];TRUEY=[];ET=[];GERR=[];\n",
    "for i in range(2,4):\n",
    "    x = x0 + (i-1)*h\n",
    "    y = y + h*dy(x0+(i-2)*h)\n",
    "    print '\\n Step %d \\n x(%d) = %f\\n y(%f) = %f\\n'%(i-1,i-1,x,x,y)\n",
    "    Et = d2y(x0+(i-2)*h)*h**2/2 +  d3y(x0+(i-2)*h)*h**3/6\n",
    "    print '\\n Et(%d) = %f\\n'%(i-1,Et)\n",
    "    truey = exacty(x0+(i-1)*h)\n",
    "    gerr = truey - y\n",
    "    X.append(x)\n",
    "    ESTY.append(y)\n",
    "    TRUEY.append(truey)\n",
    "    ET.append(Et)\n",
    "    GERR.append(gerr)\n",
    "\n",
    "#table = [x y(2:3) truey Et gerr]\n",
    "print 'x      Est y   true y    Et      Globalerr'    \n",
    "for a,b,c,d,e in nditer([X,ESTY,TRUEY,ET,GERR]):\n",
    "    print a,'\\t',b,'\\t',c,'\\t',d,'\\t',e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_06 Pg No. 427"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EULERS METHOD\n",
      "y(1.250000) = 3.000000 \n",
      " \n",
      "y(1.500000) = 4.200000 \n",
      " \n",
      "y(1.750000) = 5.600000 \n",
      " \n",
      "y(2.000000) = 7.200000 \n",
      " \n",
      "HEUNS METHOD\n",
      "\n",
      "Iteration 1 \n",
      " m1 = 4.000000\n",
      " ye(1.250000) = 3.000000 \n",
      " m2 = 4.800000 \n",
      " y(1.250000) = 3.100000 \n",
      "\n",
      "\n",
      "Iteration 2 \n",
      " m1 = 4.960000\n",
      " ye(1.500000) = 4.340000 \n",
      " m2 = 5.786667 \n",
      " y(1.500000) = 4.443333 \n",
      "\n",
      "\n",
      "Iteration 3 \n",
      " m1 = 5.924444\n",
      " ye(1.750000) = 5.924444 \n",
      " m2 = 6.770794 \n",
      " y(1.750000) = 6.030238 \n",
      "\n",
      "\n",
      "Iteration 4 \n",
      " m1 = 6.891701\n",
      " ye(2.000000) = 7.753163 \n",
      " m2 = 7.753163 \n",
      " y(2.000000) = 7.860846 \n",
      "\n",
      "x     Eulers   Heuns   \t\tAnalytical\n",
      "0.0 \t0.0 \t0.000000 \t0.0\n",
      "1.0 \t2.0 \t2.000000 \t2.0\n",
      "1.25 \t3.0 \t3.100000 \t2.2360679775\n",
      "1.5 \t4.2 \t4.443333 \t2.44948974278\n",
      "1.75 \t5.6 \t6.030238 \t2.64575131106\n",
      "2.0 \t7.2 \t7.860846 \t2.82842712475\n"
     ]
    }
   ],
   "source": [
    "from numpy import nditer,sqrt\n",
    "from __future__ import division\n",
    "\n",
    "#Heun's Method\n",
    "\n",
    "def f(x,y):\n",
    "    F = 2*y/x \n",
    "    return F\n",
    "def exacty(x):\n",
    "    Exacty = 2*sqrt(x)\n",
    "    return Exacty\n",
    "x = [0,1] #\n",
    "y = [0,2] #\n",
    "h = 0.25 #\n",
    "#Euler's Method\n",
    "print 'EULERS METHOD'\n",
    "for i in range(2,6):\n",
    "    x.append(x[i-1] + h )\n",
    "    y.append(y[i-1] + h*f(x[i-1],y[i-1]))\n",
    "    print 'y(%f) = %f \\n '%(x[i],y[i])\n",
    "\n",
    "eulery = y\n",
    "#Heun's Method\n",
    "print 'HEUNS METHOD'\n",
    "ye=[0,0]\n",
    "y = [0,2] #\n",
    "for i in range(2,6):\n",
    "    m1 = f(x[i-1],y[i-1]) #\n",
    "    ye.append(y[i-1] + h*f(x[(i-1)],y[(i-1)]))\n",
    "    m2 = f(x[(i)],ye[(i)]) \n",
    "    y.append(y[(i-1)] + h*(m1 + m2)/2)\n",
    "    print '\\nIteration %d \\n m1 = %f\\n ye(%f) = %f \\n m2 = %f \\n y(%f) = %f \\n'%(i-1,m1,x[(i)],ye[(i)],m2,x[(i)],y[(i)]) \n",
    "    \n",
    "    \n",
    "truey = exacty(x) \n",
    "print 'x     Eulers   Heuns   \\t\\tAnalytical'    \n",
    "for a,b,c,d in nditer([x,eulery,y,truey]):\n",
    "    print a,'\\t',b,'\\t%.6f'%c,'\\t',d\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_07 Pg NO. 433"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y(1.250000) = 3.111111 \n",
      " \n",
      "y(1.500000) = 4.468687 \n",
      " \n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "\n",
    "#Polygon Method\n",
    "\n",
    "def f(x,y):\n",
    "    F = 2*y/x\n",
    "    return F\n",
    "x=[1]\n",
    "y=[2]\n",
    "h = 0.25 #\n",
    "for i in range(1,3):\n",
    "    x.append(x[(i-1)] + h )\n",
    "    y.append(y[(i-1)] + h*f(  x[(i-1)]+ h/2  ,  y[(i-1)] + h*f( x[(i-1)]  ,  y[(i-1)] )/2  ))\n",
    "    print 'y(%f) = %f \\n '%(x[(i)],y[(i)])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_08 Pg No. 439"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Iteration - 1\n",
      " m1 = 0.000000\n",
      " m2 = 0.010000 \n",
      " m3 = 0.010001 \n",
      " m4 = 0.040004 \n",
      " y(0.200000) = 0.002667 \n",
      "\n",
      "\n",
      "Iteration - 2\n",
      " m1 = 0.040007\n",
      " m2 = 0.090044 \n",
      " m3 = 0.090136 \n",
      " m4 = 0.160428 \n",
      " y(0.400000) = 0.021360 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "\n",
    "#Classical Runge Kutta Method\n",
    "\n",
    "def f(x,y):\n",
    "    F = x**2 + y**2\n",
    "    return F\n",
    "h = 0.2\n",
    "x=[0]\n",
    "y=[0]\n",
    "\n",
    "for i in range(0,2):\n",
    "    m1 = f(  x[(i)]  ,  y[(i)]  ) #\n",
    "    m2 = f(  x[i] + h/2  ,  y[(i)] + m1*h/2  ) #\n",
    "    m3 = f(  x[(i)] + h/2  ,  y[(i)] + m2*h/2  ) #\n",
    "    m4 = f(  x[(i)] + h  ,  y[(i)] + m3*h  ) #\n",
    "    x.append(x[(i)] + h )\n",
    "    y.append(y[(i)] +  (m1 + 2*m2 + 2*m3 + m4)*h/6 )\n",
    "    \n",
    "    print '\\nIteration - %d\\n m1 = %f\\n m2 = %f \\n m3 = %f \\n m4 = %f \\n y(%f) = %f \\n'%(i+1,m1,m2,m3,m4,x[(i+1)],y[(i+1)])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_09 Pg No. 444"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "for four decimal places\n",
      "h =  0.0143668085653\n",
      "for six decimal places\n",
      "h =  0.00454318377739\n",
      "Note-We can use h = 0.01 for four decimal places and h = 0.004 for six decimal places\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "\n",
    "#Optimum step size\n",
    "\n",
    "x = 0.8 #\n",
    "h1 = 0.05 #\n",
    "y1 = 5.8410870 #\n",
    "h2 = 0.025 #\n",
    "y2 = 5.8479637 #\n",
    "\n",
    "#d = 4\n",
    "h = ((h1**4 - h2**4)*10**(-4)/(2*(y2 - y1)))**(1/4)\n",
    "print 'for four decimal places'\n",
    "print 'h = ',h \n",
    "\n",
    "#d = 6\n",
    "h = ((h1**4 - h2**4)*10**(-6)/(2*(y2 - y1)))**(1/4)\n",
    "print 'for six decimal places'\n",
    "print 'h = ' ,h\n",
    "print 'Note-We can use h = 0.01 for four decimal places and h = 0.004 for six decimal places'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_10 Pg NO. 446"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yp4 =  8.00914285714\n",
      "fp4 =  8.00914285714 yc4 =  8.00266666667\n",
      "f4 =  8.00266666667\n",
      "yc4 =  8.00212698413\n",
      "Note- the exact solution is y(2) = 8\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "\n",
    "#Milne-Simpson Predictor-Corrector method\n",
    "\n",
    "def f(x,y):\n",
    "    F = 2*y/x\n",
    "    return F\n",
    "x0 = 1 #\n",
    "y0 = 2 #\n",
    "h = 0.25 #\n",
    "#Assuming y1 ,y2 and y3(required for milne-simpson formula) are estimated using Fourth- order Runge kutta method\n",
    "x1 = x0 + h \n",
    "y1 = 3.13 #\n",
    "x2 = x1 + h \n",
    "y2 = 4.5 #\n",
    "x3 = x2 + h\n",
    "y3 = 6.13 #\n",
    "#Milne Predictor formula\n",
    "yp4 = y0 + 4*h*(2*f(x1,y1) - f(x2,y2) + 2*f(x3,y3))/3\n",
    "x4 = x3 + h \n",
    "fp4 = f(x4,yp4) #\n",
    "print 'yp4 = ',yp4\n",
    "print 'fp4 = ',fp4, \n",
    "#Simpson Corrector formula\n",
    "yc4 = y2 + h*( f(x2,y2) + 4*f(x3,y3) + fp4)/3 \n",
    "f4 = f(x4,yc4)\n",
    "print 'yc4 = ',yc4\n",
    "print 'f4 = ',f4 \n",
    "\n",
    "yc4 = y2 + h*( f(x2,y2) + 4*f(x3,y3) + f4)/3 \n",
    "print 'yc4 = ',yc4 \n",
    "print 'Note- the exact solution is y(2) = 8'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_11 Pg NO. 446"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adams Predictor formula\n",
      "yp4 =  8.01135714286\n",
      "fp4 =  8.01135714286\n",
      "\n",
      "Adams Corrector formula\n",
      "yc4 =  8.00727901786\n",
      "f4 =  8.00727901786\n",
      "\n",
      "refined-yc4 =  8.00689669364\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "\n",
    "#Adams-Bashforth-Moulton Method\n",
    "\n",
    "def f(x,y):\n",
    "    F = 2*y/x\n",
    "    return F\n",
    "x0 = 1 #\n",
    "y0 = 2 #\n",
    "h = 0.25 #\n",
    "x1 = x0 + h \n",
    "y1 = 3.13 #\n",
    "x2 = x1 + h \n",
    "y2 = 4.5 #\n",
    "x3 = x2 + h\n",
    "y3 = 6.13 #\n",
    "#Adams Predictor formula\n",
    "yp4 = y3 + h*(55*f(x3,y3) - 59*f(x2,y2) + 37*f(x1,y1) - 9*f(x0,y0))/24\n",
    "x4 = x3 + h \n",
    "fp4 = f(x4,yp4) \n",
    "print 'Adams Predictor formula'\n",
    "print 'yp4 = ',yp4 \n",
    "print 'fp4 = ',fp4 \n",
    "#Adams Corrector formula\n",
    "yc4 = y3 + h*( f(x1,y1) - 5*f(x2,y2) + 19*f(x3,y3) + 9*fp4 )/24 \n",
    "f4 = f(x4,yc4)\n",
    "print '\\nAdams Corrector formula'\n",
    "print 'yc4 = ',yc4 \n",
    "print 'f4 = ',f4 \n",
    "\n",
    "yc4 = y3 + h*( f(x1,y1) - 5*f(x2,y2) + 19*f(x3,y3) + 9*f4 )/24 \n",
    "print '\\nrefined-yc4 = ',yc4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_12 Pg No. 453"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " y4p = 0.557351\n",
      " f4p = -0.310640 \n",
      " y4c = 0.555396 \n",
      " Modifier Etc = 0.000067 \n",
      " Modified y4c = 0.555464 \n",
      "\n",
      "\n",
      " y5p = 0.500837\n",
      " f5p = -0.250837 \n",
      " y5c = 0.499959 \n",
      " Modifier Etc = 0.000030 \n",
      " Modified y5c = 0.499989 \n",
      "\n",
      "error =  1.11332736336e-05\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "\n",
    "#Milne-Simpson Method using modifier\n",
    "\n",
    "def f(y):\n",
    "    F = -y**2\n",
    "    return F\n",
    "x = [ 1 , 1.2 , 1.4 , 1.6 ] \n",
    "y = [ 1 , 0.8333333 , 0.7142857 , 0.625 ] \n",
    "h = 0.2 #\n",
    "\n",
    "for i in range(0,2):\n",
    "    yp = y[(i)] + 4*h*( 2*f( y[(i+1)] ) - f( y[(i+2)] ) + 2*f( y[(i+3)] ) )/3\n",
    "    fp = f(yp) #\n",
    "    yc = y[( i+2)] + h*(f( y[(i+2)] ) + 4*f( y[(i+3)] ) + fp )/3 \n",
    "    Etc = -(yc - yp)/29\n",
    "    y.append(yc + Etc)\n",
    "    print '\\n y%dp = %f\\n f%dp = %f \\n y%dc = %f \\n Modifier Etc = %f \\n Modified y%dc = %f \\n'%(i+4,yp,i+4,fp,i+4,yc,Etc,i+4,y[(i+4)])\n",
    "exactanswer = 0.5 #\n",
    "err = exactanswer - y[5] #\n",
    "print 'error = ',err"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_13 Pg No. 455"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1(1) = 0.000000\n",
      " m1(2) = 1.000000\n",
      " m2(1) = 0.200000\n",
      " m2(2) = 1.100000\n",
      " m(1) = 0.100000\n",
      " m(2) = 1.050000\n",
      " y1(0.1) = 1.010000\n",
      " y2(0.1) = -0.895000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "\n",
    "#System of differetial Equations\n",
    "\n",
    "def f1(x,y1,y2):\n",
    "    F1 = x + y1 + y2 \n",
    "    return F1\n",
    "def f2(x,y1,y2):\n",
    "    F2 = 1 + y1 + y2 \n",
    "    return F2\n",
    "\n",
    "x0 = 0 #\n",
    "y10 = 1 #\n",
    "y20 = -1 #\n",
    "h = 0.1 #\n",
    "m1= [f1( x0 , y10 , y20 )]\n",
    "m1.append(f2( x0 , y10 , y20 ))\n",
    "m2=[f1( x0+h , y10 + h*m1[0] , y20 + h*m1[1] )]\n",
    "m2.append(f2( x0+h , y10 + h*m1[0] , y20 + h*m1[1] ))\n",
    "m= [(m1[0] + m2[0])/2 ]\n",
    "m.append((m1[1] + m2[1])/2)\n",
    "\n",
    "y1_0_1 = y10 + h*m[0]\n",
    "y2_0_1 = y20 + h*m[1]\n",
    "\n",
    "print 'm1(1) = %f\\n m1(2) = %f\\n m2(1) = %f\\n m2(2) = %f\\n m(1) = %f\\n m(2) = %f\\n y1(0.1) = %f\\n y2(0.1) = %f\\n'%(m1[0],m1[1],m2[0],m2[1],m[0],m[1],y1_0_1,y2_0_1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example No. 13_14 Pg No. 457"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1(1) = 1.000000\n",
      " m1(2) = -2.000000\n",
      " m2(1) = 0.600000\n",
      " m2(2) = 0.600000\n",
      " m(1) = 0.800000\n",
      " m(2) = -0.700000\n",
      " y1(0.1) = 0.160000\n",
      " y2(0.1) = 0.860000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "#Higher Order Differential Equations\n",
    "\n",
    "x0 = 0\n",
    "y10 = 0\n",
    "y20 = 1\n",
    "h = 0.2\n",
    "m1 = [y20] #\n",
    "m1.append(6*x0 + 3*y10 - 2*y20)\n",
    "m2= [y20 + h*m1[1]]\n",
    "m2.append(6*(x0+h) + 3*(y10 + h*m1[0]) - 2*(y20 + h*m1[1]) )\n",
    "m = [(m1[0] + m2[0])/2 ]\n",
    "m.append((m1[1] + m2[1])/2)\n",
    "\n",
    "y1_0_2 = y10 + h*m[0]\n",
    "y2_0_2 = y20 + h*m[1]\n",
    "\n",
    "print 'm1(1) = %f\\n m1(2) = %f\\n m2(1) = %f\\n m2(2) = %f\\n m(1) = %f\\n m(2) = %f\\n y1(0.1) = %f\\n y2(0.1) = %f\\n'%(m1[0],m1[1],m2[0],m2[1],m[0],m[1],y1_0_2,y2_0_2) "
   ]
  }
 ],
 "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
}