from IPython.core.display import HTML
css_file = 'https://raw.githubusercontent.com/ngcm/training-public/master/ipython_notebook_styles/ngcmstyle.css'
HTML(url=css_file)
Back in October, Greg Wilson blogged about testing scientific software. He specifically points to a lesson introducing the Euler method as part of the Practical Numerical Methods with Python MOOC, and asks why a certain number (the measured convergence rate) is close enough to the expected value.
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 project, as launched and explained by Greg Wilson in this blog post. 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, 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 (an answer that is correct, of little use, and took a long time to arrive), let's begin...
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)$.
The method used is Euler's method, also introduced in an earlier lesson in the MOOC, and explained in detail by Lorena in this screencast. The key result derived there is that Euler's method is first order, so that
$$ \begin{equation} u(T) - \udt{\dt} = c_1 \dt + c_2 \dt^2 + {\cal O}(\dt^3) \simeq c_1 \dt. \end{equation} $$For sufficiently small $\dt$ the error is proportional to $\dt$.
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
$$ \begin{equation} s_m = \log_2 \left( \frac{\udt{4\dt} - \udt{2\dt}}{\udt{2\dt} - \udt{\dt}} \right) \simeq 1. \end{equation} $$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:
The original lesson 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?
Answers