Chapter13.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:9ffa27748b4aa27853ec1fa96ce6d18b12d164a5d0db63a032fd06c7dbbaacb1"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Numerical Solution of Ordindary Differential Equations\n",
      "\n",
      "Ordinary differential equations describe many physical situations: spring-mass systems, resistor-capacitor-inductance circuits, bending of beams, chemical reactions, pendulums, etc. Their prominence in applied mathematics is due to the fact that most scientific laws are readily expressed in terms of rates of change. For example,\n",
      "\n",
      "$$ \\frac{du}{dt} = -0.27(u-60)^{5/4} $$\n",
      "\n",
      "is an equation describing the rate of change of temperature $u$ of a body losing heat by natural convection with constant-temperature surroundings. This is termed a first-order differential equation because the highest-order derivative is the first.\n",
      "\n",
      "If the equation contains derivatives of $n^{'th'}$ order, it is said to be an $n^{'th'}$ order differential equation. The solution to a differential equation is the function that satisfies the differential equation and that also satisfies certain initial conditions on the function. In solving differential equations analytically, one usually finds a general solution containing arbitrary coefficients and uses the initial conditions on the function to evaluate the coefficients. Typically, a first course in differential equations convers solution techniques for linear equations with constant equations. It can become very difficult to find solutions to nonlinear differential equations with analytical techniques, this is when we typically resolve to using numerical techniques.\n",
      "\n",
      "Numerical methods have no such limitation on the form of the differential equations, however we find the solution as a tabulation of the values of the function at various values of the independent variable, but not as the explicit function. We can always use an interpolation function for that if necessary. \n",
      "\n",
      "We will begin our discussion considering only first-order equations, and then to show how these same methods can be applied to systems of simultaneous first-order equations and to higher-order differential equations. We will use for our typical first-orer equation the form:\n",
      "\n",
      "$$\n",
      "\\frac{dy}{dx} = f(x,y)\\\\\n",
      "y(x_0) = y_0\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Taylor-Series Method\n",
      "\n",
      "This method is not strictly a numerical method, but it is sometimes used in conjunction with the numerical schemes, is of general applicability, and serves as an introductgion to the other techniques we will study. Consider the problem:\n",
      "\n",
      "$$\n",
      "\\frac{dy}{dx} = x + y, y(0) = 1 (*)\n",
      "$$\n",
      "\n",
      "We develop the relation between $y$ and $x$ by finding the coefficients of the Taylor series in which we expand $y$ about the point $x = x_0$\n",
      "\n",
      "$$\n",
      "y(x) = y(x_0) + y'(x_0)(x-x_0) + \\frac{y''(x_0)}{2!}(x-x_0)^2 + \\frac{y'''(x_0)}{3!}(x-x_0)^3 + ...\n",
      "$$\n",
      "\n",
      "Letting $x-x_0 = h$, we can write\n",
      "\n",
      "$$\n",
      "y(x) = y(x_0) + y'(x_0)h + \\frac{y''(x_0)}{2}h^2 + \\frac{y'''(x_0)}{6} h^3 + ...\n",
      "$$\n",
      "\n",
      "Since $y(x_0)$ is an intial condition, the first term is known from the initial condition $y(0)=1$. We get the coefficient of the second term by substituting $x=0, y=1$ into the equation for the first derivative, $(*)$\n",
      "\n",
      "$$\n",
      "y'(x_0) = y'(0) = 0 + 1 = 1\n",
      "$$\n",
      "\n",
      "We get equations for the second- and higher-order derivatives by successively differentiation of the first equation for the first derivative. Each of these is evaluated corresponding to $x=0$ to get the various coefficients:\n",
      "\n",
      "$$\n",
      "\\begin{matrix} y''(x) = 1+y' & y''(0) = 1+1=2,\\\\\n",
      "y'''(x) = y'' & y'''(0) = 2,\\\\\n",
      "y^{iv}(x) = y''' & y^{iv}(0) = 2, \\\\\n",
      "\\vdots\n",
      "y^{(n)}(x) = y^{(n-1)} & y^{(n)}(0) = 2\\end{matrix}\n",
      "$$\n",
      "\n",
      "(it's not always this easy to determine the coeffients)\n",
      "\n",
      "Now we can write the series approximation for $y$, letting $x=h$ be the value at which we wish to determine $y$:\n",
      "\n",
      "$$\n",
      "y(h) = 1 + h + h^2 + \\frac{1}{3}h^3 + \\frac{1}{12}h^4 + error\n",
      "$$\n",
      "\n",
      "The error term of a truncated Taylor series is simple to express. We take the next term and evaluate the derivative at the point $x = \\xi, 0<\\xi<h$, instead of at the point $x=x_0$. The error term of the Taylor series after $h^4$ term is \n",
      "\n",
      "$$\n",
      "Error = \\frac{y^{v}(\\xi)}{5!}h^5, 0<\\xi<h\n",
      "$$\n",
      "\n",
      "Knowing when to truncate the taylor series is more art than science. We normally truncate the Taylor series when the contribution of the last term is negligible to the number of decimal places to which we are working. A table of answers to the problem above is listed below:\n",
      "\n",
      "|x|y|y, analytical|\n",
      "|:-:|:-:|:-:|\n",
      "|0.0|1.000|1.0000|\n",
      "|0.1|1.1103|1.1103|\n",
      "|0.2|1.2428|1.2428|\n",
      "|0.3|1.3997|1.3997|\n",
      "|0.4|1.5835|1.5836|\n",
      "|0.5|1.7969|1.7974|\n",
      "\n",
      "Below is a plot showing the exact (analytical) solution in blue and the Taylor approximation to the solution in red. Manipulate the slider bar to see the affect adding more terms has ont he approximated solution."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Euler Methods\n",
      "\n",
      "We will now use what we've learned about Taylor's method to help derive Euler's method. Although Euler's method is seldom used in practice, the simplicity of its derivation can be used to illustrate the techniques involved in the construction of some of the more advanced techniques, without the cumbersome algebra that accompanies some of the other constructions.\n",
      "\n",
      "By inspection of the error term in Taylor's method we can see that this term will be small if the step size $h$ is small. In fact, if it is small enough, only a few terms of the series are needed for good accuracy. The Euler method may be thought of as following this idea to the extreme for first-order differential equations. Suppose we choose $h$ small enough that we may truncate after the first-derivative term. Then\n",
      "\n",
      "$$\n",
      "y(x_0 + h) = y(x_0) + hy'(x_0) + \\frac{y''(\\xi)}{2}h^2, x_0 < \\xi < x_0 + h\n",
      "$$\n",
      "\n",
      "In using this equation, the value of $y(x_0)$ is given by the initial condition and $y'(x_0)$ is evaluated from $f(x_0, y_0)$, given by the differential equation, $\\frac{dy}{dx} = f(x,y)$. It will of course be necessary to use this method iteratively, advancing the solution to $x = x_0 + 2h$ after $y(x_0 + h)$ has been found, then to $x=x_0 + 3h$, and so on. Utilizing a subscript notation for the successive $y-$values we have:\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_{n} + hy'_{n} + O(h^2)\n",
      "$$\n",
      "\n",
      "The trouble with this most simple method is its lack of accuracy, requiring an extremely small step size."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Psuedocode for Euler's Method\n",
      "\n",
      "To approximate the solution of the initial-value problem\n",
      "\n",
      "$$\n",
      "y'=f(t,y), a\\leq t \\leq b, y(a) = \\alpha\n",
      "$$\n",
      "\n",
      "  1. Set $h = \\frac{b-a}{N}$ where $N$ defines the number of points and therefore the step size on the interval $[a,b]$\n",
      "  1. Set $t=a$\n",
      "  1. Set $w=\\alpha$\n",
      "  1. For $i = 1, 2, ..., N$ do Steps A-B\n",
      "    1. Set $w = w + h f(t,w)$\n",
      "    1. Set $t= a + ih$\n",
      "  1. Output $w$ as approximation to $y$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Modified Euler Method\n",
      "\n",
      "In the simple Euler Method, we use the slope at the beginning of the interval, $y'_n$, to determine the increment to the function. This is always wrong. If the slope of the function were constant, the solution is an obvious linear relation. We need to use an average slope over the interval if we hope to estimate the change in $y$ with precision. Using the average of the slopes at the beginning and end of the interval:\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_{n} + h\\frac{y'_n + y'_{n+1}}{2}\n",
      "$$\n",
      "\n",
      "This will give an improved estimate for $y$ at $x_{n+1}$. However, we are unable to use this equation immediately, because $y'_{n+1}$ is unknown. Therefore we use the simple Euler method to estimate or \"predict\" a value of $y_{n+1}$ and then compute $y'_{n+1}$ based on this value. $y'_{n+1}$ is then used to go back and \"correct\" the original estimate of $y_{n+1}$. Since hte value of $y'_{n+1}$ was computed using the predicted value, of less than perfect accuracy, one is tempted to recorrect the $y_{n+1}$ as many times as will make a significant difference. This method is called the *modified Euler method*, or sometimes the *Euler predictor-corrector method*. The method is utilized below using *Python* to perform the algebra to solve the following differential equation.\n",
      "\n",
      "$$\n",
      "\\frac{dy}{dx} = x + y, y(0) = 1, h=0.02\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Compare with the analytic value of 1.1103. Modified Euler's method is slightly more accurate. We can find the error of the modified Euler method by comparing with the Taylor series.\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_n + y'_n h + \\frac{1}{2}y''_n h^2 + \\frac{y'''(\\xi)}{6} h^3, x_n < \\xi < x_n + h\n",
      "$$\n",
      "\n",
      "Replace the second derivative by the forward-difference approximation for $y'', {y'_{n+1} - y'_n}{h}$, which has an error of $O(h)$, and write the error term as $O(h^3)$:\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_n + h\\left (y'_n + \\frac{1}{2}\\left[\\frac{y'_{n+1} - y'_n}{h} + O(h)\\right]h\\right) + O(h^3) \\\\\n",
      "y_{n+1} = y_n + h(y'_n + \\frac{1}{2} y'_{n+1} - \\frac{1}{2}y'_n) + O(h^3) \\\\\n",
      "y_{n+1} = y_n + h(\\frac{y'_n + y'_{n+1}}{2}) + O(h^3)\n",
      "$$\n",
      "\n",
      "This shows that the error of one step of the modified Euler method is $O(h^3)$. This the \"local error\". There is an accumulation of errors from step to step, so that the error over the whole range of application, or \"global error\", is $O(h^2)$\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Runge-Kutta methods\n",
      "\n",
      "Probably one of the most widely used numerical techniques for solving differential equations are the Runge-Kutta methods. The development of these techniques is algebraically complicated, therefore we will not attempt to show every detail of the derivation of the popular $4^{th}$ order Runge-Kutta method; instead, to get an idea of how the Runge-Kutta methods are developed, we will study a second-order method. \n",
      "\n",
      "We start by writing the increment to $y$ as a weighted average of two estimates of $\\Delta y, k_1,$ and $k_2$. For the equation $\\frac{dx}{dy} = f(x,y)$.\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_n + ak_1 + bk_2  (*) \\\\\n",
      "k_1 = hf(x_n, y_n) \\\\\n",
      "k_2 = hf(x_n + \\alpha h, y_n + \\eta k_1)\n",
      "$$\n",
      "\n",
      "We can think of $k_1$ and $k_2$ as estimates of the change in $y$ when $x$ advances by $h$ because they are the product of the change in $x$ and a value for the slope of the curve, $\\frac{dy}{dx}$. The Runge-Kutta methods always use as the first estimate of $\\Delta y$ the simple Euler estimate; the other estimate is taken with $x$ and $y$ stepped up by the fraction $\\alpha$ and $\\beta$ of $h$ and of the earlier estimate of $\\Delta y, k_1$. We need to devise a scheme to choose the four parameters, $a, b, \\alpha, \\eta$. We do this by making the equations above agree as well as possible with the Taylor-series expansion, in which the $y$-derivatives are written in terms of $f$, from $\\frac{dy}{dx} = f(x,y)$,\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_n + h f(x_n, y_n) + (\\frac{h^2}{2})f'(x_n, y_n) + ...\n",
      "$$\n",
      "\n",
      "An equivalent form, since $\\frac{df}{dx} = f_x + f_y \\frac{dy}{dx} = f_x + f_y f$, is\n",
      "\n",
      "$$\n",
      "y_{n+1} y_n + h f_n + h^2(\\frac{1}{2} f_x + \\frac{1}{2} f_y f)_n  (**)\n",
      "$$\n",
      "\n",
      "If we rewrite $(*)$ with the $k_i$'s substituted in, we have:\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_n + \\alpha h f(x_n, y_n) + b h f [x_n + \\alpha h, y_n + \\beta h f(x_n, y_n)]  (***)\n",
      "$$\n",
      "\n",
      "we desire this equation to look as much like $(**)$ as possible. To make the last term comparable we expand $f(x,y)$ in a Taylor series in terms of $x_n, y_n$, remembering that $f$ is a function of two variables, retaining only first derivative terms.\n",
      "\n",
      "$$\n",
      "f[x_n + \\alpha h, y_n + \\beta h f(x_n, y_n)] = (f + f_x \\alpha h + f_y \\beta f)_n\n",
      "$$\n",
      "\n",
      "Substituting this into $(***)$, we have\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_n + \\alpha h f(x_n, y_n) + b h(f+f_x \\alpha h + f_y \\beta f)_n\n",
      "$$ \n",
      "\n",
      "or\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_n + (a+b)hf_n + h^2(\\alpha b f_x + \\beta b f_y f)_n\n",
      "$$\n",
      "\n",
      "This equation will be identical to $(**)$ if \n",
      "\n",
      "$$\n",
      "a + b = 1, \\alpha b = \\frac{1}{2}, \\beta b = \\frac{1}{2}\n",
      "$$\n",
      "\n",
      "Here we have 4 unknowns and three equations. If we choose one of them arbitrarily (with some restrictions). If we chose $\\alpha = \\frac{1}{2}$, the other variables are $b = \\frac{1}{2}, \\alpha 1, \\beta = 1$. Here we recover the modified Euler method previously discussed. The modified Euler method is a special case of a second-order Runge-Kutta.\n",
      "\n",
      "Fourth-order Runge-Kutta methods are the most widely used and are derived in a similar fashion. Of course since we have to compare terms through $h^4$ this is a much more complex derivation resulting in 11 equations and 13 unknowns. If we choose two of the unkowns arbitrarily we can resolve the rest of the unknowns. The most commonly used values lead to the algorithm:\n",
      "\n",
      "$$\n",
      "y_{n+1} = y_n + \\frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \\\\\n",
      "k_1 = hf(x_n,y_n) \\\\\n",
      "k_2 = hf(x_n + \\frac{1}{2}h, y_n + \\frac{1}{2}k_1) \\\\\n",
      "k_3 = hf(x_n + \\frac{1}{2}h, y_n + \\frac{1}{2}k_2) \\\\\n",
      "k_4 = hf(x_n + h, y_n + k_3)\n",
      "$$\n",
      "\n",
      "The local error term for the fourth-order Runge-Kutta is $O(h^5)$; the global error is $O(h^4)$. It is computationally more efficient than the modified Euler method because, while four evaluations of the function are required per step rather than two, the steps can be manyfold larger for the same accuracy."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Pseudocode for $4^{th}$ order Runge-Kutta method\n",
      "\n",
      "To approximate the solution of the initial-value problem\n",
      "\n",
      "$$\n",
      "y' = f(x,y), a \\leq t \\leq b, y(a) = \\alpha\n",
      "$$\n",
      "\n",
      "at $(N+1)$ equally spaced numbers in the interval $[a, b]$\n",
      "\n",
      "  1. Set $h = \\frac{b-a}{N}$\n",
      "  1. Set $t = a$\n",
      "  1. Set $w = \\alpha$\n",
      "  1. For $i = 1, 2, ..., N$ do Steps A-F\n",
      "    1. Set $K_1 = hf(t,w)$\n",
      "    1. Set $K_2 = h f (t + \\frac{h}{2}, w + \\frac{K_1}{2})$\n",
      "    1. Set $K_3 = h f (t + \\frac{h}{2}, w + \\frac{K_2}{2})$\n",
      "    1. Set $K_4 = h f (t+ h, w + K_3)$\n",
      "    1. Set $w = w + {K_1 + 2K_2 + 2K_3 + K_4}{6}$\n",
      "    1. Set $t = a + i h$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Runge-Kutta-Fehlberg method\n",
      "\n",
      "We saw that going from a second-order Runge-Kutta (in which Euler's method was a special case) to a fourth-order Runge-Kutta was a tremendous increase in efficiency (computational efficiency = fewest steps with the least work at each step). We might wish to use an even higher order formula (fifth, sixth, etc.). These higher-order formulas have been developed and can be used to advantage in determining a suitable step size, $h$. \n",
      "\n",
      "The Runge-Kutta-Fehlberg method uses a combination of fourth and fifth-order Runge-Kutta methods and compares estimates of $y_{n+1}$ that result from each of the two methods to get an error estimate. This error estimate is then used to correct the step size, $h$ at the next step. Both the fourth and fifth order Runge-Kutta use the same $k$'s, so only six functional evaluations are needed at every step and we get the added benefit of being able to increase $h$ depending of the error value. This algorithm is summarized below:\n",
      "\n",
      "$$\n",
      "\\begin{matrix}\n",
      "k_1  = hf(x_n, & y_n) & & & & & & \\\\\n",
      "k_2 = h f(x_n & + \\frac{h}{4},&  y_n & + \\frac{k_1}{4}) & & & &\\\\\n",
      "k_3 = h f(x_n & + \\frac{3h}{8}, & y_n & + \\frac{3k_1}{32} & + \\frac{9k_2}{32}) & & &\\\\\n",
      "k_4 = h f(x_n & + \\frac{12h}{13}, & y_n & + \\frac{1932k_1}{2197} & - \\frac{7200k_2}{2197} & + \\frac{7296k_3}{2197} & &\\\\\n",
      "k_5 = h f(x_n & + h, & y_n & + \\frac{439k_1}{216} & - 8k_2 & + \\frac{3680k_3}{513} & - \\frac{854k_4}{4104} &\\\\\n",
      "k_6 = h f(x_n & + \\frac{h}{2}, & y_n & - \\frac{8k_1}{27} & + 2k_2 &  - \\frac{3544k_3}{2565} & + \\frac{1859k_4}{4104} & - \\frac{11k_5}{40})\n",
      "\\end{matrix}\n",
      "$$\n",
      "\n",
      "$$\n",
      "\\hat{y}_{n+1} = y_n + \\left(\\frac{25k_1}{216} + \\frac{1408 k_3}{2565} + \\frac{2197k_4}{4104} - \\frac{k_5}{5}\\right) \\\\\n",
      "y_{n+1} = y_n + \\left( \\frac{16k_1}{135} + \\frac{6656k_3}{12825} + \\frac{28561k_4}{56430} - \\frac{9k_5}{50} + \\frac{2k_6}{55} \\right)\\\\\n",
      "E = \\frac{k_1}{360} - \\frac{128k_3}{4275} - \\frac{2197k_4}{75240} + \\frac{k_5}{50} + \\frac{2k_6}{55}\n",
      "$$\n",
      "\n",
      "where, $\\hat{y}_{n+1}$ is the fourth-order Runge-Kutta estimate, $y_{n+1}$ is the fifth-order Runge-Kutta estimate, and $E$ is the error between them."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Pseudocode for Runge-Kutta-Fehlburg method\n",
      "\n",
      "To approximate the solution of the initial-value problem\n",
      "\n",
      "$$\n",
      "y' = f(t,y), a \\leq t \\leq b, y(a) = \\alpha\n",
      "$$\n",
      "\n",
      "with local truncation error within a given tolerance:\n",
      "\n",
      "  1. Set $t = a$\n",
      "  1. Set $w = \\alpha$\n",
      "  1. Set $h = HMAX$\n",
      "  1. Set $FLAG = 1$\n",
      "  1. While ($FLAG = 1$) do Steps A-O\n",
      "    1. Set $K_1 = h f(t, w)$\n",
      "    1. Set $K_2 = h f(t + \\frac{h}{4}, w + \\frac{k_1}{4})$\n",
      "    1. Set $K_3 = h f(t + \\frac{3h}{8}, w + \\frac{3K_1}{32} + \\frac{9K_2}{32})$\n",
      "    1.  Set $K_4 = h f(t + \\frac{12h}{13}, w + \\frac{1932K_1}{2197} - \\frac{7200K_2}{2197} + \\frac{7296K_3}{2197})$\n",
      "    1. Set $K_5 = h f(t + h, w + \\frac{439K_1}{216} - 8K_2 + \\frac{3680K_3}{513} - \\frac{854K_4}{4104})$\n",
      "    1. Set $K_6 = h f\\left(t + \\frac{h}{2}, w - \\frac{8K_1}{27} + 2K_2 - \\frac{3544K_3}{2565} + \\frac{1859K_4}{4104} - \\frac{11K_5}{40}\\right)$\n",
      "    1. Set $R = \\frac{1}{h}\\left[ \\frac{K_1}{360} - \\frac{128K_3}{4275} - \\frac{2197K_4}{75240} + \\frac{K_5}{50} + \\frac{2K_6}{55}\\right]$\n",
      "    1. If $R \\leq TOL$ do Steps a-b\n",
      "      1. Set $t = t+h$\n",
      "      1. Set $w = w + \\left( \\frac{25K_1}{216} + \\frac{1408K_3}{2565} + \\frac{2197K_4}{4104} - \\frac{K_5}{5}\\right)$\n",
      " (Here is where we change the size of $h$, below is a common method)\n",
      "    1. Set $\\delta = 0.84(\\frac{TOL}{R})^{\\frac{1}{4}}$\n",
      "    1. If $\\delta \\leq 0.1$ do Step a\n",
      "      1. Set $h = 0.1 h$\n",
      "    1. If $\\delta \\geq 4$ do Step a\n",
      "      1. Set $h = 4h$\n",
      "    1. Else do Step a\n",
      "      1. Set $h = \\delta h$\n",
      "    1. If $h > HMAX$ do Step a\n",
      "      1. Set $h = HMAX$\n",
      "    1. If $t \\leq b$ do Step a\n",
      "      1. Set $FLAG = 0$\n",
      "    1. if $h < HMIN$ do Steps a-b\n",
      "      1. Set $FLAG = 0$\n",
      "      1. Print \"Minimum $h$ exceeded\""
     ]
    }
   ],
   "metadata": {}
  }
 ]
}