{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "example problem\n",
      "==="
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\\frac{dy}{dx} + 0.6y = 10exp(-{\\frac{(x-a)^2}{2(0.075)^2}})$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f(x,y):\n",
      "    return 10*np.exp(-(x-2)**2/(2*(0.075**2))) - 0.6*y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "RK4, step size = 0.1, range 0-4, 40 steps\n",
      "==="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "h = 0.1\n",
      "x0 = 0\n",
      "y0 = 0.5\n",
      "xlist = []\n",
      "xlist.append(x0)\n",
      "ylist = []\n",
      "ylist.append(y0)\n",
      "xi = x0\n",
      "yi = y0\n",
      "for i in xrange(40):\n",
      "    k1 = f(xi,yi)\n",
      "    k2 = f(xi+h/2,yi+k1*h/2)\n",
      "    k3 = f(xi+h/2,yi+k2*h/2)\n",
      "    k4 = f(xi+h,yi+k3*h)\n",
      "    xi += h\n",
      "    yi += (k1+2*k2+2*k3+k4)*h/6\n",
      "    xlist.append(xi)\n",
      "    ylist.append(yi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "RK4, step size = 0.01, range 0-4, 400 steps\n",
      "==="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "h = 0.01\n",
      "x0 = 0\n",
      "y0 = 0.5\n",
      "xlist = []\n",
      "xlist.append(x0)\n",
      "ylist = []\n",
      "ylist.append(y0)\n",
      "xi = x0\n",
      "yi = y0\n",
      "for i in xrange(400):\n",
      "    k1 = f(xi,yi)\n",
      "    k2 = f(xi+h/2,yi+k1*h/2)\n",
      "    k3 = f(xi+h/2,yi+k2*h/2)\n",
      "    k4 = f(xi+h,yi+k3*h)\n",
      "    xi += h\n",
      "    yi += (k1+2*k2+2*k3+k4)*h/6\n",
      "    xlist.append(xi)\n",
      "    ylist.append(yi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "RK4, step size = 0.001, range 0-4, 4000 steps\n",
      "==="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "h = 0.001\n",
      "x0 = 0\n",
      "y0 = 0.5\n",
      "rk4xlist = []\n",
      "rk4xlist.append(x0)\n",
      "rk4ylist = []\n",
      "rk4ylist.append(y0)\n",
      "xi = x0\n",
      "yi = y0\n",
      "for i in xrange(4000):\n",
      "    k1 = f(xi,yi)\n",
      "    k2 = f(xi+h/2,yi+k1*h/2)\n",
      "    k3 = f(xi+h/2,yi+k2*h/2)\n",
      "    k4 = f(xi+h,yi+k3*h)\n",
      "    xi += h\n",
      "    yi += (k1+2*k2+2*k3+k4)*h/6\n",
      "    rk4xlist.append(xi)\n",
      "    rk4ylist.append(yi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "RK5, step size = 0.001, range 0-4, 4000 steps\n",
      "==="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "h = 0.001\n",
      "x0 = 0\n",
      "y0 = 0.5\n",
      "rk5xlist = []\n",
      "rk5xlist.append(x0)\n",
      "rk5ylist = []\n",
      "rk5ylist.append(y0)\n",
      "xi = x0\n",
      "yi = y0\n",
      "for i in xrange(4000):\n",
      "    k1 = f(xi,yi)\n",
      "    k2 = f(xi+h/4,yi+k1*h/4)\n",
      "    k3 = f(xi+h/4,yi+k1*h/8+k2*h/8)\n",
      "    k4 = f(xi+h/2,yi-k2*h/2+k3*h)\n",
      "    k5 = f(xi+3*h/4,yi+3*k1*h/16+9*k4*h/16)\n",
      "    k6 = f(xi+h,yi+(-3*k1*h+2*k2*h+12*k3*h-12*k4*h+8*k5*h)/7)\n",
      "    xi += h\n",
      "    yi += (7*k1+32*k3+12*k4+32*k5+7*k6)*h/90\n",
      "    rk5xlist.append(xi)\n",
      "    rk5ylist.append(yi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "comparing the difference between rk4 and rk5 estimates\n",
      "==="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xtemp = np.asarray(rk5xlist)-np.asarray(rk4xlist)\n",
      "ytemp = np.asarray(rk5ylist)-np.asarray(rk4ylist)\n",
      "plt.plot(np.asarray(rk5xlist),ytemp)\n",
      "# difference in the estimation of the function value using rk4 and rk5 methods with same step size, 0.001"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}