{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Just Euler's method" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "css_file = 'https://raw.githubusercontent.com/ngcm/training-public/master/ipython_notebook_styles/ngcmstyle.css'\n", "HTML(url=css_file)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\newcommand{\\dt}{\\Delta t}\n", "\\newcommand{\\udt}[1]{u^{({#1})}(T)}\n", "\\newcommand{\\Edt}[1]{E^{({#1})}}\n", "\\newcommand{\\uone}[1]{u_{1}^{({#1})}}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the previous cases we've focused on the *behaviour* of the algorithm: whether it will give the correct answer in the limit, or whether it converges as expected. This is really what you want to do: you're trying to do science, to get an answer, and so implementing the precise algorithm should be secondary. If you are trying to implement a precise algorithm, it should be because of its (expected) behaviour, and so you should be testing for that!\n", "\n", "However, let's put that aside and see if we can work out how to test whether we've implemented exactly the algorithm we want: Euler's method. Checking convergence alone is not enough: the [Backwards Euler method](http://en.wikipedia.org/wiki/Backward_Euler_method) has identical convergence behaviour, as do whole families of other methods. We need a check that characterizes the method uniquely." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The *local truncation error* $\\Edt{\\dt}$ would be exactly such a check. This is the error produced by a single step from exact data, eg\n", "\n", "$$\n", "\\begin{equation}\n", " \\Edt{\\dt} = u_1 - u(\\dt).\n", "\\end{equation}\n", "$$\n", "\n", "For Euler's method we have\n", "\n", "$$\n", "\\begin{equation}\n", " u_{n+1} = u_n + \\dt f(t_n, u_n)\n", "\\end{equation}\n", "$$\n", "\n", "and so\n", "\n", "$$\n", "\\begin{equation}\n", " \\Edt{\\dt} = \\left| u_0 + \\dt f(0, u_0) - u(\\dt) \\right| = \\left| \\frac{\\dt^2}{2} \\left. u''\\right|_{t=0} \\right| + {\\cal O}(\\dt^3).\n", "\\end{equation}\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is all well and good, but we don't know the exact solution (in principle) at any point other than $t=0$, so cannot compute $u(\\dt)$, so cannot compute $\\Edt{\\dt}$. We only know $\\uone{\\dt}$ for whichever values of $\\dt$ we wish to compute." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use repeated Richardson extrapolation to get the solution $u(\\dt)$ to sufficient accuracy, however. On the *assumption* that the algorithm is first order (we can use the previous techniques to check this), we can use Richardson extrapolation to repeatedly remove the highest order error terms. We can thus find the local truncation errors." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from math import sin, cos, log, ceil\n", "import numpy\n", "from matplotlib import pyplot\n", "%matplotlib inline\n", "from matplotlib import rcParams\n", "rcParams['font.family'] = 'serif'\n", "rcParams['font.size'] = 16" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# model parameters:\n", "g = 9.8 # gravity in m s^{-2}\n", "v_t = 30.0 # trim velocity in m s^{-1} \n", "C_D = 1/40. # drag coefficient --- or D/L if C_L=1\n", "C_L = 1.0 # for convenience, use C_L = 1\n", "\n", "### set initial conditions ###\n", "v0 = v_t # start at the trim velocity (or add a delta)\n", "theta0 = 0.0 # initial angle of trajectory\n", "x0 = 0.0 # horizotal position is arbitrary\n", "y0 = 1000.0 # initial altitude" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def f(u):\n", " \"\"\"Returns the right-hand side of the phugoid system of equations.\n", " \n", " Parameters\n", " ----------\n", " u : array of float\n", " array containing the solution at time n.\n", " \n", " Returns\n", " -------\n", " dudt : array of float\n", " array containing the RHS given u.\n", " \"\"\"\n", " \n", " v = u[0]\n", " theta = u[1]\n", " x = u[2]\n", " y = u[3]\n", " return numpy.array([-g*sin(theta) - C_D/C_L*g/v_t**2*v**2,\n", " -g*cos(theta)/v + g/v_t**2*v,\n", " v*cos(theta),\n", " v*sin(theta)])" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def euler_step(u, f, dt):\n", " \"\"\"Returns the solution at the next time-step using Euler's method.\n", " \n", " Parameters\n", " ----------\n", " u : array of float\n", " solution at the previous time-step.\n", " f : function\n", " function to compute the right hand-side of the system of equation.\n", " dt : float\n", " time-increment.\n", " \n", " Returns\n", " -------\n", " u_n_plus_1 : array of float\n", " approximate solution at the next time step.\n", " \"\"\"\n", " \n", " return u + dt * f(u)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "T_values = numpy.array([0.001*2**(i) for i in range(10)])\n", "lte_values = numpy.zeros_like(T_values)\n", "for j, T in enumerate(T_values):\n", " dt_values = numpy.array([T*2**(i-8) for i in range(8)])\n", " v_values = numpy.zeros_like(dt_values)\n", " for i, dt in enumerate(dt_values):\n", " N = int(T/dt)+1\n", " t = numpy.linspace(0.0, T, N)\n", " u = numpy.empty((N, 4))\n", " u[0] = numpy.array([v0, theta0, x0, y0])\n", " for n in range(N-1):\n", " u[n+1] = euler_step(u[n], f, dt)\n", " v_values[i] = u[-1,0]\n", " v_next = v_values\n", " for s in range(1, len(v_values-1)):\n", " v_next = (2**s*v_next[1:]-v_next[0:-1])/(2**s-1)\n", " lte_values[j] = abs(v_values[0]-v_next)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 1.99954897e-09, 8.05573563e-09, 3.26825287e-08,\n", " 1.34407252e-07, 5.67045063e-07, 2.50349015e-06,\n", " 1.18961299e-05, 6.26360623e-05, 3.70836832e-04,\n", " 2.44291651e-03])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lte_values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have four values for the local truncation error. We can thus compute the convergence rate of the local truncation error itself (which should be two), and check that it is close enough to the expected value using the previous techniques:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Measured convergence rate (base dt 0.001) is 2.02375 (error is 0.02375).\n", "Measured convergence rate (base dt 0.002) is 2.04637 (error is 0.04637).\n", "Convergence error has reduced by factor 0.5121.\n" ] } ], "source": [ "s_m = numpy.zeros(2)\n", "for i in range(2):\n", " s_m[i] = log(abs((lte_values[2+i]-lte_values[1+i])/\n", " (lte_values[1+i]-lte_values[0+i]))) / log(2.0)\n", " print(\"Measured convergence rate (base dt {}) is {:.6g} (error is {:.4g}).\".format(\n", " T_values[i], s_m[i], abs(s_m[i]-2)))\n", "print(\"Convergence error has reduced by factor {:.4g}.\".format(\n", " abs(s_m[0]-2)/abs(s_m[1]-2)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the error has gone down considerably, and certainly $0.51 < 2/3$, so the convergence rate of the local truncation error is close enough to 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, that alone isn't enough to determine that this really is Euler's method: as noted above, the convergence rate of the local truncation error isn't the key point: the key point is that we can predict its *actual value* as\n", "\n", "$$\n", "\\begin{equation}\n", " \\Edt{\\dt} = \\frac{\\dt^2}{2} \\left| \\left. u''\\right|_{t=0} \\right| + {\\cal O}(\\dt^3) = \\frac{\\dt^2}{2} \\left| \\left( \\left. \\frac{\\partial f}{\\partial t} \\right|_{t=0} + f(0, u_0) \\left. \\frac{\\partial f}{\\partial u} \\right|_{t=0, u=u_0} \\right) \\right|.\n", "\\end{equation}\n", "$$\n", "\n", "For the specific problem considered here we have\n", "\n", "$$\n", "\\begin{equation}\n", " u = \\begin{pmatrix} v \\\\ \\theta \\\\ x \\\\ y \\end{pmatrix}, \\quad f = \\begin{pmatrix} -g\\sin \\theta - \\frac{C_D}{C_L} \\frac{g}{v_t^2} v^2 \\\\ -\\frac{g}{v}\\cos \\theta + \\frac{g}{v_t^2} v \\\\ v \\cos \\theta \\\\ v \\sin \\theta \\end{pmatrix}.\n", "\\end{equation}\n", "$$\n", "\n", "We note that $f$ does not explicitly depend on $t$ (so $\\partial f / \\partial t \\equiv 0$), and that the values of the parameters $g, C_D, C_L$ and $v_t$ are given above, along with the initial data $u_0 = (v_0, \\theta_0, x_0, y_0)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, let's find what the local truncation error should be." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import sympy\n", "sympy.init_printing()\n", "v, theta, x, y, g, CD, CL, vt, dt = sympy.symbols('v, theta, x, y, g, C_D, C_L, v_t, {\\Delta}t')\n", "u = sympy.Matrix([v, theta, x, y])\n", "f = sympy.Matrix([-g*sympy.sin(theta)-CD/CL*g/vt**2*v**2, \n", " -g/v*sympy.cos(theta)+g/vt**2*v, \n", " v*sympy.cos(theta), \n", " v*sympy.sin(theta)])\n", "dfdu = f.jacobian(u)\n", "lte=dt**2/2*dfdu*f" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": [ "iVBORw0KGgoAAAANSUhEUgAAAQMAAABlCAMAAAB+1XYKAAAAP1BMVEX///8AAAAAAAAAAAAAAAAA\n", "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRS\n", "TlMAMquZdlQQQO0wRIlmzd0i77t8bBwggJIAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAhwSURBVHgB\n", "7V1rg6MoECS+dk+Njzv//2+9fgCFpmNcMyazE/lgChoQygaVkhnnPic0+ZhXVncvk4TMsv2stKZ3\n", "Lrsmfeq0685dprygUCa2Hwo74uAy0SGEljueTZx6CWk/6re96U1XO9dPdJiF1uSgypu2SYhBHKjP\n", "83HUPEi0kKuKosnlxH3TNKP6G6CFnJWINNOcJEoXL1Mx66qPlNNyQrA5uFI7K3YcHxCPqOLpo1FO\n", "Y6KzkMu5loyq7BmVE18fQAvtMSdl+HzO5UPSW/Smu2HG5KAduI4mTh6IA2XC5sA5kWghV4pP9lRb\n", "MbEPTB0dAC20x5yUofqdq8t2agTxoSR/l9DIFfER/TE5GOUWAadBHGjgnriMpxIkWsg1I+esKH89\n", "0IioJuYW0EJ7zEkZPh97X+IIub+g7S0FzuRgEg5quWpcHeJAV3E04QCJFnL1NJLPNGGG0rHA1fph\n", "cQ/tMaNMTS6XOEKnLlHTT51MdNIMi4NqErLqMKcgDiSFneuICSRaiDJm09RGCirfGkoHtNAec1JG\n", "uuAdoc7GacxoHrhcy7LkKzILFge9jqNLGE6IA2klNWdBooU4YzNNVz1x33ThOQ3QQjTJh5wWMs1J\n", "Is0GfOLoCH46GOSRiC1psDkQEhMOQrxXD4mWjgcNEi1EGcqs7iaeCiQ0ETlAC+0xxzIymjEjhOnA\n", "t2H+Y3EAl9a8iAOJJRdukGghuhjcnmIaggsOETlAC+0x+zKluAE7gjTRYQDOuy8xiwM/B5aLOZHj\n", "Oul5S6H1J4kwA7lBnjMuVLoXNHK9gBbaY07K0BQUuqozQnXzbBjs/GtyMMrdrI2PGIgDkYszBdwD\n", "JFqo8nfmrHWT1CgcAFooybnZnJRxwQ2CI9TchuCH3O1ZMDnQRx2MIcSB3EW8oCAOkGgh7wdu7N0g\n", "3Hb80gJooT3mpAzcgGcEOhs/7pXLW2LkweTAdeSt1UCF9IEGcaC+o2f/PJMnpZgd5gTVck9o6dbU\n", "8jNCKbdcQAslOTebkzLlNY+h4xkhz+ilJfZ5CWwOqobeKIU3mfkd4hFd9c1bOIiJVka6K2e5f2eq\n", "qWmj3iAALeSsRKSZZiT6VQFtIQ+mfhQql533cZuDO5l/aPLJwZ37wg+93ve6dfrB6QfsG6cfnBx8\n", "sB/MtJbPHAtzreUzOZhrLZ/AQVjJxPPBXGtZ5wCSiZZHHCjRWugFrSm8NuPRpaOXBQ702hAhkJn4\n", "nJle9KKmI61+qLWscwDJRDlAPKJUa6mu9D4kr9wRtf7FhRY1AIHMxOfMsqqumo42+qHWssqBLgds\n", "1lrcyEsKso4SUU5v7xR4XQcQyEx8zpxoOnLqx1rLKgeQTKQ2U0xJtJY6KrpAuqwv4wMQyAFaaI85\n", "0XSk1XRdHmgtqxxAMlEOEAdKtJZRJDrOCiQFa6xfAAI5QAv9qXmu6WzQWtY4wDKxUoA4kFpEa3HT\n", "9dIUIjADiV2W0DQnIJADtNAfm2eajqz4rWstaxxAMtH2Iw6kFtFaaN2N70JDyStwHok5akyJ4gbt\n", "zU58UGbdDE1nk9ayzoGQGBWVOxIK9VNW3CpdAs6HCog5CAvLM/gg8Slzqun4RfYwI2CdmFvmw4yD\n", "6trFcL0kQqJmxggAEotqLSqq092uThDZCy/5zuCDxGfMqaYTFtlXtRbm4Nf0Szu5PCZCiZgQB+Ke\n", "ibuQN0hnW/rGAojMHQY5IFBiRyLQn5uh6WzUWv4lDu5+jwTJRNlBHCjRWjK5L7AfAPFQiKIPIFBi\n", "RyLQn5vDMCJNZ6PWMhsL2lMcIZloGuJAidZyERkpJyaAGHs3SeGDxKfM6ges6WzUWlY5gFCyTWvJ\n", "6DmnovsCnTwi0lQiB4BAiR2JQDvM0HQ2ai3rHEA82aS1uCbPVZpJUC+kiCMBAjlAC+0xR01no9ay\n", "zgGGxU9GJwfnmir79+kHJwenH+hMf44FHQtVcfdTHSXqex1n+shXNK3nrzf/qv0Lc33kKzj4+8bC\n", "XB/5CA5uBJK5PvJ3cAA1ZtleqDNAV5qbLrn/ypQK2AIJthUs69wTP34sRDUmNq/Uj3WhzgDRMjgF\n", "XWaX/LZAcrsXJVa+AxzOwVKnqbI80w9gsRMGyF2bjD/7DMEWSIy9KKHAnt/DOVjqNNTIQjmAOgOU\n", "rnpwd2jpISyHckxXJq29KJx5bzicA11Jw54YaqjnAOoM0IIDUyAx96Ls7T+XO5qDxQq0NNVzoM3m\n", "nTAJytq2yOOXxbICtRBI7L0ovo5dP0dzsFRjuJEpB6LOSMsV8UabuAvT3Ixi70XZ1Xlf6HgO5FpG\n", "nYZPm3Kga3ScCuRG+Qo6LomGGcEUSLjok+HLOXig03BzEw68OkOJQIRVwN4kkDzZfy7OHPz6/c8X\n", "1HSnipkao3nAQVBniBddfNZN6X77bNQl1BHWN6PcOf2W5P9+H/zOlKoxvkGRg7gTJu6JGeTul8s8\n", "GdyALxQztL4ZZUtn7+T58rGwPE+ixgRT4AA7YSJSb1CpLrrBps0oofIdv4dzcKPTxPkAO2ESxCS0\n", "ohRtFEh29HlZ5HgOljpNng3TlcVA7IQBos3rTaa7QDcKJMsO7Ygfz8GORr24yMnB8c/KL76ku053\n", "+sHpB+w4px+cHJx+oFPoORbOsXCOhW8zFu6LMNrEw4/fYD64FWEO7/X8BO/nYCnCzNv3itj7OTBE\n", "mFd0PDnH+zkwRJikfa+Ab+fAEmFe0fHkHG/nwBJhkva9An4DDmQddSbCvKLjyTnezsE5FuhqGCJM\n", "co1eAd/uB8lfGntFf61zvJ8DQ4SxGnpg2vs5gAhzYDdXq/4GHECEWW3pccZvwMFxndtY88nBuY7E\n", "rnL6wcnB6Qc6aepYkD/cknz3obYPOMb/w0F/S4dD+CPQH9D12EX5PxxF4f4HVD6e7pxnK1EAAAAA\n", "SUVORK5CYII=\n" ], "text/latex": [ "$$\\left[\\begin{matrix}0.00200083333333333 {\\Delta}t^{2}\\\\- 0.00266777777777778 {\\Delta}t^{2}\\\\- 0.1225 {\\Delta}t^{2}\\\\0\\end{matrix}\\right]$$" ], "text/plain": [ "⎡ 2 ⎤\n", "⎢0.00200083333333333⋅{\\Delta}t ⎥\n", "⎢ ⎥\n", "⎢ 2⎥\n", "⎢-0.00266777777777778⋅{\\Delta}t ⎥\n", "⎢ ⎥\n", "⎢ 2 ⎥\n", "⎢ -0.1225⋅{\\Delta}t ⎥\n", "⎢ ⎥\n", "⎣ 0 ⎦" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lte_0=lte.subs([(g,9.8),(vt,30.0),(CD,1.0/40.0),(CL,1.0),(v,30.0),(theta,0.0),(x,0.0),(y,1000.0)])\n", "lte_0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So let us check the local truncation error values, which are computed for `v`:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 0.00199955, 0.00201393, 0.00204266, 0.00210011, 0.00221502,\n", " 0.00244481, 0.00290433, 0.003823 , 0.00565852, 0.00931899])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lte_exact = float(lte_0[0]/dt**2)\n", "lte_values/T_values**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are indeed converging towards $0.002 \\dt^2$ as they should. To check this quantitatively, we use that our model is\n", "\n", "$$\n", "\\begin{equation}\n", " \\Edt{\\dt} = \\alpha \\dt^2 + {\\cal O}(\\dt^3),\n", "\\end{equation}\n", "$$\n", "\n", "with the exact value $\\alpha_e \\simeq 0.002$. So we can use our usual Richardson extrapolation methods applied to $\\Edt{\\dt}/\\dt^2$, to get a measured value for $\\alpha$ with an error interval:\n", "\n", "$$\n", "\\begin{equation}\n", " \\alpha_m = \\frac{8\\Edt{\\dt} - \\Edt{2\\dt} \\pm \\left| \\Edt{\\dt} - \\Edt{2\\dt} \\right|}{4\\dt^2}.\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Base dt=0.001: the measured alpha is in [0.00047112, 0.0034992]\n", "Does this contain the exact value? True\n", "Base dt=0.002: the measured alpha is in [0.00044604, 0.0035244]\n", "Does this contain the exact value? True\n", "Base dt=0.004: the measured alpha is in [0.00039575, 0.0035747]\n", "Does this contain the exact value? True\n", "Base dt=0.008: the measured alpha is in [0.00029522, 0.0036752]\n", "Does this contain the exact value? True\n", "Base dt=0.016: the measured alpha is in [9.4165e-05, 0.0038763]\n", "Does this contain the exact value? True\n", "Base dt=0.032: the measured alpha is in [-0.00030782, 0.0042784]\n", "Does this contain the exact value? True\n", "Base dt=0.064: the measured alpha is in [-0.0011113, 0.0050826]\n", "Does this contain the exact value? True\n", "Base dt=0.128: the measured alpha is in [-0.0027153, 0.0066903]\n", "Does this contain the exact value? True\n", "Base dt=0.256: the measured alpha is in [-0.0059063, 0.0099024]\n", "Does this contain the exact value? True\n" ] } ], "source": [ "for i in range(len(lte_values)-1):\n", " Edt = lte_values[i]\n", " E2dt = lte_values[i+1]\n", " dt = T_values[i]\n", " err1 = abs(Edt - E2dt)\n", " a_lo = (8.0*Edt - E2dt - err1)/(4.0*dt**2)\n", " a_hi = (8.0*Edt - E2dt + err1)/(4.0*dt**2)\n", " print(\"Base dt={:.4g}: the measured alpha is in [{:.5g}, {:.5g}]\".format(\n", " dt, a_lo, a_hi))\n", " print(\"Does this contain the exact value? {}\".format(\n", " a_lo <= lte_exact <= a_hi))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, to the limits that we can measure the local truncation error, we have implemented Euler's method." ] } ], "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.4.2" } }, "nbformat": 4, "nbformat_minor": 0 }