{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Testing numerical methods" ] }, { "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": [ "Back in October, Greg Wilson [blogged about testing scientific software](http://software-carpentry.org/blog/2014/10/why-we-dont-teach-testing.html). He specifically points to [a lesson introducing the Euler method](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/01_03_PhugoidFullModel.ipynb) as part of the [Practical Numerical Methods with Python MOOC](http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about), and asks why a certain number (the measured convergence rate) is *close enough* to the expected value.\n", "\n", "The big question of testing in scientific software, how this should be done, and how it is being done, has led to the [*Close Enough for Scientific Work*](https://github.com/swcarpentry/close-enough-for-scientific-work) project, as launched and explained by Greg Wilson in [this blog post](http://software-carpentry.org/blog/2014/11/close-enough-for-scientific-work.html). If you're reading this, you should be following that project! As I've helped out on the MOOC, which is led by [Lorena Barba](http://lorenabarba.com/), I've been (over)thinking what sort of answer we could, or should, give. Well aware that we're heading for the [academic version of the balloon joke](https://www2.bc.edu/~radinr/Management_Humor/jokes.htm) (an answer that is correct, of little use, and took a long time to arrive), let's begin..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Background" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\newcommand{\\dt}{\\Delta t}\n", "\\newcommand{\\udt}[1]{u^{({#1})}(T)}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have a system of differential equations to solve: the solution will be $u(t)$. The numerical solution we get will depend on one parameter, the timestep $\\dt$. The resulting numerical approximation, at the point $t=T$, will be $\\udt{\\dt}$: the exact solution is $u(T)$. \n", "\n", "The method used is Euler's method, also [introduced in an earlier lesson in the MOOC](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/01_02_Phugoid_Oscillation.ipynb), and [explained in detail by Lorena in this screencast](https://www.youtube.com/watch?v=6i6qhqDCViA). The key result derived there is that Euler's method is *first order*, so that\n", "\n", "$$\n", "\\begin{equation}\n", " u(T) - \\udt{\\dt} = c_1 \\dt + c_2 \\dt^2 + {\\cal O}(\\dt^3) \\simeq c_1 \\dt.\n", "\\end{equation}\n", "$$\n", "\n", "For sufficiently small $\\dt$ the error is proportional to $\\dt$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we can't (usually) know the exact solution $u(T)$ ahead of time, a standard check is *grid convergence*: compute three solutions with different timesteps (eg, $\\dt, 2\\dt, 4\\dt$) and compare them. The comparison to use is\n", "\n", "$$\n", "\\begin{equation}\n", " s_m = \\log_2 \\left( \\frac{\\udt{4\\dt} - \\udt{2\\dt}}{\\udt{2\\dt} - \\udt{\\dt}} \\right) \\simeq 1.\n", "\\end{equation}\n", "$$\n", "\n", "The symbol $s_m$ is used as it's the *(measured) slope* of the best fit line through the errors, as shown on the [original lesson](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/01_03_PhugoidFullModel.ipynb):\n", "\n", "![Image](images/slope1.png)\n", "####Figure 1. Convergence test plot taken from the original lesson." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The [original lesson](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/01_03_PhugoidFullModel.ipynb) does precisely this comparison (using the phugoid model, with specific parameter values, and $\\dt = 0.001$), finding a measured slope $s_m = 1.014$. Is this *close enough* to 1?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## TL;DR: Answer" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answer**s**\n", "\n", "1. [*I don't care about the algorithm: I just want the answer to be right*](http://nbviewer.ipython.org/github/IanHawke/close-enough-balloons/blob/master/01-Close-Enough-Errorbars.ipynb): $0.585 \\lesssim s_m \\lesssim 1.585$ is close enough for Euler's method.\n", "2. [*I don't care about the answer: I just want the algorithm behaviour to be right*](http://nbviewer.ipython.org/github/IanHawke/close-enough-balloons/blob/master/02-Close-Enough-Slopes.ipynb): $s_m=1.014$ with $\\dt=0.001$ is close enough if $1 \\le s_m \\le 1.0093$ with $\\dt=0.0005$.\n", "3. [*I don't care about the behaviour: I just want to know I've implemented Euler's method*](http://nbviewer.ipython.org/github/IanHawke/close-enough-balloons/blob/master/03-Close-Enough-Just-Euler.ipynb): measure the local truncation error instead, and check *both* the convergence rate *and* the leading order constant. $s_m$ by itself is useless." ] } ], "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 }