{
 "metadata": {
  "name": "",
  "signature": "sha256:fc18dc0db333034491fb24c6d4a925206386fe68260cc2e8bfe29433e758df85"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Consider the iterations \n",
      "$$\n",
      "x_0 = 1, \\qquad x_1 = \\frac{1}{3}, \\qquad x_i = \\frac{13}{3}x_{i-1} - \\frac{4}{3} x_{i-2}, \\qquad i=2,3,\\ldots\n",
      "$$\n",
      "The exact solution is\n",
      "$$\n",
      "x_i = \\frac{1}{3^i}\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import zeros, float32, float64"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We first try with single precision."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 20\n",
      "x = zeros(n, dtype=float32)\n",
      "y = zeros(n, dtype=float32)\n",
      "x[0], y[0] = 1.0, 1.0\n",
      "x[1], y[1] = 1.0/3.0, 1.0/3.0\n",
      "\n",
      "for i in range(2,n):\n",
      "   x[i] = (13.0/3.0)*x[i-1] - (4.0/3.0)*x[i-2]\n",
      "   y[i] = y[i-1]/3.0\n",
      "   print \"%6d %20.10e %20.10e\" % (i, x[i], y[i])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "     2     1.1111115664e-01     1.1111111194e-01\n",
        "     3     3.7037219852e-02     3.7037037313e-02\n",
        "     4     1.2346410193e-02     1.2345679104e-02\n",
        "     5     4.1181510314e-03     4.1152262129e-03\n",
        "     6     1.3834409183e-03     1.3717421098e-03\n",
        "     7     5.0404260401e-04     4.5724736992e-04\n",
        "     8     3.3959673601e-04     1.5241578512e-04\n",
        "     9     7.9952902161e-04     5.0805261708e-05\n",
        "    10     3.0118301511e-03     1.6935087842e-05\n",
        "    11     1.1985225603e-02     5.6450294323e-06\n",
        "    12     4.7920204699e-02     1.8816764396e-06\n",
        "    13     1.9167391956e-01     6.2722546090e-07\n",
        "    14     7.6669335365e-01     2.0907515363e-07\n",
        "    15     3.0667726994e+00     6.9691715510e-08\n",
        "    16     1.2267090797e+01     2.3230571244e-08\n",
        "    17     4.9068363190e+01     7.7435240442e-09\n",
        "    18     1.9627345276e+02     2.5811746074e-09\n",
        "    19     7.8509381104e+02     8.6039153580e-10\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The iterations start to blow up around iteration 9. Now we try with double precision."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 30\n",
      "x = zeros(n, dtype=float64)\n",
      "y = zeros(n, dtype=float64)\n",
      "x[0], y[0] = 1.0, 1.0\n",
      "x[1], y[1] = 1.0/3.0, 1.0/3.0\n",
      "\n",
      "for i in range(2,n):\n",
      "   x[i] = (13.0/3.0)*x[i-1] - (4.0/3.0)*x[i-2]\n",
      "   y[i] = y[i-1]/3.0\n",
      "   print \"%6d %20.10e %20.10e\" % (i, x[i], y[i])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "     2     1.1111111111e-01     1.1111111111e-01\n",
        "     3     3.7037037037e-02     3.7037037037e-02\n",
        "     4     1.2345679012e-02     1.2345679012e-02\n",
        "     5     4.1152263374e-03     4.1152263374e-03\n",
        "     6     1.3717421124e-03     1.3717421125e-03\n",
        "     7     4.5724737062e-04     4.5724737083e-04\n",
        "     8     1.5241578946e-04     1.5241579028e-04\n",
        "     9     5.0805260180e-05     5.0805263425e-05\n",
        "    10     1.6935074827e-05     1.6935087808e-05\n",
        "    11     5.6449773443e-06     5.6450292695e-06\n",
        "    12     1.8814687225e-06     1.8816764232e-06\n",
        "    13     6.2639467164e-07     6.2722547439e-07\n",
        "    14     2.0575194713e-07     2.0907515813e-07\n",
        "    15     5.6398875392e-08     6.9691719376e-08\n",
        "    16    -2.9940802813e-08     2.3230573125e-08\n",
        "    17    -2.0494197938e-07     7.7435243751e-09\n",
        "    18    -8.4816084023e-07     2.5811747917e-09\n",
        "    19    -3.4021076685e-06     8.6039159724e-10\n",
        "    20    -1.3611585443e-05     2.8679719908e-10\n",
        "    21    -5.4447393362e-05     9.5599066360e-11\n",
        "    22    -2.1778992398e-04     3.1866355453e-11\n",
        "    23    -8.7115981276e-04     1.0622118484e-11\n",
        "    24    -3.4846392900e-03     3.5407061615e-12\n",
        "    25    -1.3938557173e-02     1.1802353872e-12\n",
        "    26    -5.5754228696e-02     3.9341179572e-13\n",
        "    27    -2.2301691478e-01     1.3113726524e-13\n",
        "    28    -8.9206765914e-01     4.3712421747e-14\n",
        "    29    -3.5682706366e+00     1.4570807249e-14\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The iterations are stable for more steps but they eventually start to blow up."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}