{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Unstable iterations" ] }, { "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", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "%config InlineBackend.figure_format = 'svg'\n", "from numpy import zeros, float32, float64, abs\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first try with single precision. By default, Numpy uses double precision, so to use single precision, we have to explicitly specify this in Numpy functions." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "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" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "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]))\n", "\n", "plt.semilogy(range(n),x,'o-',range(n),y,'*--')\n", "plt.legend(('Iteration','Exact'))\n", "plt.xlabel('i');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The iterations start to blow up around iteration 9. Now we try with double precision." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "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" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "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]))\n", " \n", "plt.semilogy(range(n),abs(x),'o-',range(n),y,'*--')\n", "plt.legend(('Iteration','Exact'))\n", "plt.xlabel('i');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The iterations are stable for more steps but they eventually start to blow up." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 1 }