{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Slopes"
]
},
{
"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{\\sm}[1]{s_m^{({#1})}}\n",
"$$\n",
"Let's suppose that instead of caring about the answer $u$ we care about the algorithm behaving in the right fashion, $\\udt{\\dt} = u(T) + c_1 \\dt^{s_e}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given three pieces of data we can measure the slope:\n",
"\n",
"$$\n",
"\\begin{align}\n",
" s_m & = \\log_2 \\left| \\frac{\\udt{4\\dt}-\\udt{2\\dt}}{\\udt{2\\dt}-\\udt{\\dt}} \\right| \\\\ \n",
" & = \\log_2 \\left| \\frac{u(T) + c_1 4^{s_e} \\dt^{s_e} - u(T) + c_1 2^{s_e} \\dt^{s_e} + {\\cal O}(\\dt^{s_e+1})}{u(T) + c_1 2^{s_e} \\dt^{s_e} - u(T) + c_1 \\dt^{s_e} + {\\cal O}(\\dt^{s_e+1)}} \\right| \\\\\n",
" & = \\log_2 \\left| 2^{s_e} + {\\cal O}(\\dt) \\right| \\\\\n",
" & = s_e + {\\cal O}(\\dt).\n",
"\\end{align}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So now we can apply the same Richardson extrapolation tricks to our measured slope. If the measured slope is denoted $\\sm{\\dt}$ then we can use Richardson extrapolation as before to compute\n",
"\n",
"$$\n",
"\\begin{align}\n",
" s_R &= 2\\sm{\\dt} - \\sm{2\\dt}, \\\\ E_R & = s_R - \\sm{\\dt}.\n",
"\\end{align}\n",
"$$\n",
"\n",
"We then state that our algorithm is *close enough* to correct if $s_e$ lies in the interval $s_R \\pm E_R$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"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": true
},
"outputs": [],
"source": [
"T = 100.0\n",
"dt_values = numpy.array([0.001*2**i for i in range(-1,3)])\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]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Measured convergence rate (base dt 0.0005) is 1.011621.\n",
"Measured convergence rate (base dt 0.001) is 1.023266.\n"
]
}
],
"source": [
"s_m = numpy.zeros(2)\n",
"for i in range(2):\n",
" s_m[i] = log(abs((v_values[2+i]-v_values[1+i])/(v_values[1+i]-v_values[0+i]))) / log(2.0)\n",
" print(\"Measured convergence rate (base dt {}) is {:.7g}.\".format(dt_values[i], s_m[i]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then use Richardson extrapolation applied to the measured convergence rates (at finite resolution) to get the \"*actual*\" convergence rate and its accuracy:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The convergence rate lies in the interval [0.9883297,1.011621]\n"
]
}
],
"source": [
"s_exact = 2*s_m[0]-s_m[1]\n",
"error_e = abs(s_exact - s_m[0])\n",
"print(\"The convergence rate lies in the interval [{:.7g},{:.7g}]\".format(\n",
" s_exact-error_e, s_exact+error_e))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that $\\sm{\\dt} \\simeq 1.023$ is *close enough* when $\\sm{\\dt/2} \\simeq 1.012$, as the interval is consistent with $s_e=1$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finding the limits"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Once again, we want to know what the limiting values are. In this case we're interested in fixing $s_e$ and $\\sm{\\dt}$, and finding the value of $\\sm{\\dt/2}$ that is close enough.\n",
"\n",
"So the model we'll use is that $s_m = s_e + c_1 \\dt$, as above, so we need\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" 2 \\sm{\\dt/2} - \\sm{\\dt} - | \\sm{\\dt/2} - \\sm{\\dt} | \\le s_e \\le 2 \\sm{\\dt/2} - \\sm{\\dt} + | \\sm{\\dt/2} - \\sm{\\dt} |.\n",
"\\end{equation}\n",
"$$\n",
"\n",
"In the case where $\\sm{\\dt/2} < \\sm{\\dt}$ this gives\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" 3 \\sm{\\dt/2} - 2\\sm{\\dt} \\le s_e \\le \\sm{\\dt/2} \\quad \\implies \\quad \\sm{\\dt/2} \\le \\tfrac{1}{3} \\left( 2 \\sm{\\dt} + s_e \\right).\n",
"\\end{equation}\n",
"$$\n",
"\n",
"In the case where $\\sm{\\dt/2} > \\sm{\\dt}$ this gives\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" \\sm{\\dt/2} \\le s_e \\le 3 \\sm{\\dt/2} - 2 \\sm{\\dt} \\quad \\implies \\quad \\sm{\\dt/2} \\ge \\tfrac{1}{3} \\left( 2 \\sm{\\dt} + s_e \\right).\n",
"\\end{equation}\n",
"$$\n",
"\n",
"We can write this in terms of the *error* in the slope as\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" | \\sm{\\dt/2} - s_e | \\le \\frac{2}{3} | \\sm{\\dt} - s_e |;\n",
"\\end{equation}\n",
"$$\n",
"\n",
"that is, the error must go down by at least a factor of $2/3$ when the step size is cut in half.\n",
"\n",
"So, if we fix $s_e=1$ (as for Euler's method), we can see the *upper* bound on $\\sm{\\dt/2}$ given $s_m = 1.023$ is:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Upper bound when s_e=1 and s_m (coarse resolution)=1.023: 1.0155107.\n"
]
}
],
"source": [
"s_e = 1.0\n",
"s_mdthalf_max = (2.0*s_m[1] + s_e)/3.0\n",
"print(\"Upper bound when s_e=1 and s_m (coarse resolution)=1.023: {:.7f}.\".format(s_mdthalf_max))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In general, for Euler's method, the bound varies linearly:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAgQAAAGWCAYAAAAHRJtDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
"AAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm4XWV59/HvnQBhjFZiQoIQwmWZZBREGaoBCq3FoUId\n",
"WgUVB8Rq61trQdQacHjBYkEtOLQg1rEgiOCAUEy0gmABFRB9wWZCAkmwwQRQA8n9/rHWSfbZ7DPu\n",
"ae29v5/rOtc+Z+21du6zPLJ/+1n3ep7ITCRJ0mCb0u0CJElS9xkIJEmSgUCSJBkIJEkSBgJJkoSB\n",
"QJIkUeFAEBGzI+LaiNjYxGt8MCI2RsRrW1mbJEn9ppKBICJOAG4EdgMmNVFCRDwD+LvyeCdbkCRp\n",
"FJUMBMA7gaOBm4GY5Gt8GLihieMlSRoYVQ0ER2bm0skeHBEHA4cDn2hZRZIk9bFKBoJsfj7ljwJn\n",
"AutbUI4kSX2vkoGgGRHx58C0zLys27VIktQrtuh2Aa0UEVsC5wCndLsWSZJ6Sb+NEJwG3JWZN3W7\n",
"EEmSeknfjBBExFOBM4DnN3p6hGO8HVGSNHAy80nvi30TCIDnAU8Al0ds+j23Lx/Pjoh3AFdk5gdr\n",
"D2p0UqosIhZk5oJu19HPPMft5znuDM9z+/XiOR7pw3AvBIIRP8VHxCxgVRauBXate/4FwELgfZn5\n",
"7+0tU5Kk3tULPQQjDfcfAawALhzHsT01CiBJUqdVMhBExEURsQQ4EciIWBIRi8u7CIasA9ZQhIL6\n",
"43cqj/8yxQjDeeVrHNqJ+ttsUbcLGACLul3AAFjU7QIGxKJuFzAAFnW7gFaJ5ucA6l0Rkb3WQyBJ\n",
"UjNGeu+r5AiBJEnqLAOBJEkyEEiSJAOBJEnCQCBJkjAQSJIkDASSJAkDgSRJwkAgSZIwEEiSJAwE\n",
"kiQJA4EkScJAIEmSMBBIkiQMBJIkCQOBJEnCQCBJkjAQSJIkDASSJAkDgSRJwkAgSZKALbpdgCRJ\n",
"ar+I2JPtOGak5w0EkiT1sYiYzg6cyF7M5zHg0cb7GQgkSepDETGFqRzKzpzEAUzj2Szji+w80v4G\n",
"AkmS+kxEzOZpnMQ8nsVhrGQGq8c6xkAgSVKfiIit2IbjmMfLeA7r2ZslxPiONRBIktQHImJPZnIK\n",
"ezOTQ1nBdjw+keMNBJIk9bBNTYN7M5/n8TBzWTaZ16lsIIiI2cBngeMy0/kSJEmq0bBpcAs2Tvb1\n",
"KhkIIuIE4DxgPZATOG474O3A8cA25ddG4JLMPL8NpUqS1HE1TYP7chgPjqdpcCyVDATAO4GjgQXA\n",
"HhM4bh7wAeClmfktgIj4M+DqiNg2Mz/U6kIlSeqUBk2Di8fbNDiWqgaCIzMzIyb8Wz4CfHIoDABk\n",
"5rci4k7gBMBAIEnqSc02DY6lkoEgM8d9maDuuKXA3zR46inA/2umJkmSuqFVTYNjqWQgaJWyp+Dd\n",
"FL0E7+5yOZIkjVtd0+DWzTYNjqVvA0FE/AB4DvBL4MTM/EmXS5IkaVwmM9Ngs/o2EGTmkRExFfgr\n",
"4IaIeFdm/ku365IkaSQRsRVbT26mwWb1bSAAyMwNwOcj4gjgvIi4PDNX1u4TEQtqflyUmYs6WKIk\n",
"SUAbmwZvZTcWsxsAq5k+0m59FQgiYktgYxkEat0BbAXsCwwLBJm5oDPVSZL0ZMOWJz6sDU2Dh7CU\n",
"Q1gKwOfYhUc4qNFuvRAIRrzjICJmAatq7kp4D/A74Jy6XXcrH3/d8uokSZqETU2DcziJA5ufabBZ\n",
"vTAlcMOrJ+VlgBXAhTWbE/jriNi7Zr8jgbcA37OxUJJUBRGxE0/jnRzIaZzAoxzKr7oZBqCiIwQR\n",
"cRHwQmAGkBGxhOLNfs/MHLqmsg5YQxEKhnyO4hbDL0Uxq9EWwAbgw8DHO1S+JEkNlU2DxzKPEzrd\n",
"NDiWmOQcQH0hIjIzK/I/hSSpn0XEHszkDe2aaXBcPscuLOGURu99lRwhkCSpX0TEdKZzAntzFIex\n",
"hl3bM9NgswwEkiS1Qdk0+Bx25uRWLE/cbgYCSZJarGwaHFqeuCMzDTbLQCBJUouUyxNXsmlwLAYC\n",
"SZJaoN3LE7ebgUCSpCYMW564wk2DYzEQSJI0CXXLE1e+aXAsBgJJkiaoXJ74NcxjPw7nQXasftPg\n",
"WAwEkiSNU9k0WLs88eJeaRoci4FAkqRx6PWmwbEYCCRJGkW/NA2OxUAgSVID/dY0OBYDgSRJdcqm\n",
"waGZBh/shZkGm2UgkCSp1M9Ng2MxEEiSRP83DY7FQCBJGmiD0jQ4FgOBJGkgDVrT4FgMBJKkgVOz\n",
"PPF+g9I0OBYDgSRpYDRYnnhgmgbHYiCQJA2EiNiDmZzCXuzEc7l/0JoGx2IgkCT1tYiYznROYG+O\n",
"4nk8zFyWdrumKjIQSJL6kk2DE2MgkCT1nZqZBp/FYay0aXBsBgJJUt+IiK3YethMg0tsGhwfA4Ek\n",
"qS8M+kyDzTIQSJJ6Ws1Mgy/gMB4e1JkGm2UgkCT1pAZNg8ttGpw8A4EkqefUzDS4r02DrVHZQBAR\n",
"s4HPAsdl5pRu1yNJ6r4GMw3aNNgilQwEEXECcB6wHsgJHDcbeAvw58AUit/vbuD9mXlXG0qVJHVI\n",
"OdPgG2wabI+qfvJ+J3A0cDNMKPu9H3gV8GeZuR9wILABuCUi9m15lZKktouI6fGUeB17cyYvYhpH\n",
"scww0HqVHCEAjszMjJjwOFAC52bm/QCZ+fuIOAP4C+DNwN+0tkxJUrvYNNhZlQwEmTnuywR13gZP\n",
"+mN5oHx86uQrkiR1kssTd14lA8FkZeaGBpv3KB8XdbAUSdIkuDxx9/RVIBjBm4G7gM93uxBJ0sic\n",
"abC7+joQRMQxwCuAP8rMhn9YEbGg5sdFmbmoA6VJkko1Mw3O5zDWONNgi93KbixmNwBWM32k3WLy\n",
"l+vbLyIuBU6ezDwEEXEA8A3g5Zl58wj7ZGY6GCVJXbCpaXCnTU2DK2wabLPPsQtLOKXRe19fjhBE\n",
"xP7A14BXjhQGJEndU7M88b42DVZDLwSCEYcwImIWsKr2roQyDFwFvCYzbyq3zaaYnOgt7S5WkjSy\n",
"smmwdnlimwYrohcCQcM/lYg4Avg+8GngreW2/YAbgK8Cu0fE7uXuM4A921+qJGkkNg1WWyUDQURc\n",
"BLyQ4o08I2IJxUjBnjXNgeuANcCKmkMXAE8DTi2/ai1qY8mSpBHYNNgbKt1U2G42FUpS+9g0WEGD\n",
"1lQoSequmqbBZ7k8cW8wEEiSWqZB06DLE/cIA4EkqSVsGuxtBgJJUlNsGuwPBgJJ0qQ0WJ54mU2D\n",
"vctAIEmaMJcn7j8GAknSuA1bnvgQHmcfZxrsFwYCSdK42DTY3wwEkqRR2TQ4GAwEkqSGbBocLAYC\n",
"SdKTONPg4DEQSJI2iYit2NqZBgeRgUCSBNg0OOgMBJI04GqaBl/AYTxs0+BgMhBI0oBq0DS43KbB\n",
"wWUgkKQBVDPT4L42DQoMBJI0UMqmwWKmQZsGVcNAIEkDwqZBjcZAIEl9zqZBjYeBQJL6lE2DmggD\n",
"gST1IZsGNVEGAknqI8OWJ7ZpUBNgIJCkPmHToJphIJCkHjdseeLn8TBzbRrUxBkIJKlHuTyxWslA\n",
"IEk9qGZ54n05jAdtGlSzDASS1EPKpsHa5YkX2zSoVjAQSFKPsGlQ7TSl2wWMJiJmR8S1ETHha2IR\n",
"MS0izouIDRHx/HbUJ0mdEBHTY3q8nr15Ny9iGkexzDCgVqvsCEFEnACcB6wHcoLH7gt8HtgIDqZJ\n",
"6k3ONKhOqmwgAN4JHA0sAPaY4LGnl8fvAny2tWVJUvvVNA0+y5kG1QlVDgRHZmZGTOoD/mszc2NE\n",
"vK7FNUlSWzVoGnSmQXVEZQNBZk7oMkHdsQ6pSeo5Ng2qmyobCCRpUAybafAw1rg8sbrBQCBJXeJM\n",
"g6qSgQ8EEbGg5sdFmbmoS6VIGiA1yxPv50yDaqtb2Y3F7AbAaqaPtNvAB4LMXNDtGiQNjgbLEzvT\n",
"oNrrEJZyCEsB+By78AgHNdpt4AOBJHVKROzBTN5g06CqqFcCwYh3HETELGBVM3clSFI7RcR0pnMC\n",
"e3OUyxOrqnolEDQcUIuII4DvA58G3jrR4yWpncqmweewMyfbNKiqq2wgiIiLgBcCM4CMiCUUIwV7\n",
"ZubQMNs6YA2wou7Yk4GzgO3LY74SEb8DTs7M/+rQryBpgNU0De7rTIPqBTHII+0RkZnp6IGkltnU\n",
"NLjTpqbBlY5RqjI+xy4s4ZRG732VHSGQpF5j06B6mYFAkppk06D6gYFAkibJmQbVTwwEkjQJ5fLE\n",
"r7FpUP3CQCBJE+DyxOpXBgJJGieXJ1Y/MxBI0hhcnliDwEAgSSOwaVCDxEAgSQ2UTYNDMw26PLH6\n",
"noFAkmq4PLEGlYFAkko2DWqQGQgkDTybBiUDgaQBZtOgtJmBQNJAsmlQGs5AIGmgNJhp0KZBCQOB\n",
"pAGyqWlwL2bxXO63aVDabEKBICKmAXsA2wOrgcWZ6fU2SZU2rGmwWJ54abdrkqpmXIEgIp4DnArs\n",
"ADwIrAOeAsyIiFXARzNzeduqlKRJsGlQGr8xA0FEnAw8DpyamRsaPL8NcFJE/DIzv9uGGiVpwmqa\n",
"Bp/l8sTS2EYNBOWb/Xcz81cj7ZOZvwU+ExHPbHVxkjRREbEVW7s8sTRRowaC8s1+xDBQt+8vW1KR\n",
"JE1S2TT4evZmljMNShMzZawdIuKAiPhpRPwmIr4SEU8vt786Ir7d/hIlaXQRMT2mx+vZizM5nm04\n",
"imWGAWlixtNUuAD4R+CXwBHAlyPipMz8YkRc0M7iJGk0Ng1KrTOeQPCNzPx6+f3PIuIy4L0RcV4b\n",
"65KkUUXETjUzDdo0KDVpPIEgI2I/itsO35uZD0fE6cBbgK3bWp0k1SmbBmuXJ7ZpUGqBMQNBZl4S\n",
"EX8K3AM8Um7bAFxYzkEgSR0REXswkze4PLHUemPddvhs4DiKyYj+PTOfqH0+My9vY22SBJQzDU7n\n",
"hHJ54oddnlhqvbFuO7wduD0iZgIvj4inAYuBqzPz950oUNLgatA0uNymQak9xjV1cWauAv4VICLm\n",
"AW8qJy26E7i+0QyGzYqI2cBngeMyc8zbIyX1l5qmwf1cnlhqvwmvdpiZS4B/AYiI/YG/jYipwA8y\n",
"84etKCoiTgDOA9YDOcFjt6S4TfIvgCeAtcA/ZOaNrahNUnuVyxPXNg26PLHUAU0tf5yZdwB3REQA\n",
"z25NSQC8EziaYg6EPSZ47CeA+cARmfnriHgDcF1EHJ6ZP21hjZJabNPyxDYNSh03oaH4iDiq5vtj\n",
"I2ILgCzc1sK6jszMpRM9KCL2BN4EnJOZvy5ruxhYAnyohfVJaqFNMw3uzZm8iGnONCh13niXP74G\n",
"uBfYGBG3ZeZa4EbgVcAXWl1UZk7oMkGNlwEBLKzbvhA4NSK2zczHmipOUss406BUHeO9ZPAy4LkU\n",
"1+a/GhHbAbdSvPm2PBA0YX9gA7C8bvsSit91H4q6JXVZzfLE+3I4D7KjTYNSN433LoMngBsj4p8z\n",
"8zsRsRVwMFTu/8AzgMcajDCsLR93rD8gIg5qe1WSam3FFhzLLuzNzmzHrqxkJdNZyfRuFyb1vY0j\n",
"v+9PqKkwM79TPq4HWnJHQdfN4uJN3/8BDzCLB7tYjdTfHuLpPMxvuJ/FrOG7/IbHuLvbRUl97vfs\n",
"zRPsA0DyxEi7Tfoug4iYA7wrM//PZF+jDR4CtouIqBslGPrk8esnHXEaV3eiMGmgrWUatzCHVfyE\n",
"1Xy2vH1ZUhdExCsbbZ/whD8RsXNEfAK4Gzit2cJa7KcUv9MuddvnAY+Dn0WkjtpIcAezuYIZ3MYX\n",
"WM3ZhgGpmsYdCCLiGRFxIfAzimbCfelMD8GIdxxExKxyDoQhXyv3P6pu16OA67zDQOqglWzP15nH\n",
"Qu5mGWfmb/P6+vVQJFXHmIGgDAIXAXcBG4F9MvNtmfmrtldXljBCXUcAK4ALh7Zl5j3AZ4B3R8SO\n",
"5X6vpxgheE/7S5XEeqZyE7tyFVvxMy5gDR/PzKo1IEuqM54egiuB+yiCwIo21wNAGUBeSHHXQEbE\n",
"EopP/ntm5tBkJeuANRShoNbbgfdT3BXxOMUdBseVsypKaqfF7MjNTGcF1/IIV2fmo90uSdL4xFhz\n",
"AEXEFOClFNflr8rM5TXP3ZeZ9dfrG73GYcCtNW/mlRARyQLO6nYdUs9byzRuZjb3stymQanayr77\n",
"J42+jzlCkJkbga+V1+pfHBEvAa7JzFHXIy+7GJ8L3ATcQDGr4ecnU7ykitpIcBc7cStTWMUX+R0L\n",
"7ROQetO4bzssb+O7ugwGx5fBYJtRDpkNfAV4AUWz32IMBFL/WMX23MRMlnEba/iSfQJSb5vM8scJ\n",
"fKMMBqONEizPzB8BPwL+aZL1Saqa9UzlNuZwB4+wio+xgdubWH9EUkWMepdBRGwdES9q9Fy5wuHV\n",
"5X7RYKKDp0TEv0XE/HKqY0m9bjE7cjm7cCP/yQOcmU/kbYYBqT+MOkKQmb+LiHsi4lzgB8D3ypUO\n",
"AYiImcCxwIHAuXWHzwO+CbwI+FBErMzME1pavaTO2Nw0eB+rOT8zF3e7JEmtNeZdBrDpToOXlF9z\n",
"gC0pVhW8F7gyM29ocMwJmXllzc9bZ+bvWlV4K3iXgTSGoabB25jCSi6zaVDqfZO+ywA23WlwVfk1\n",
"XjtFxBnA5Zn5P1ULA5LGsLJsGlxu06A0CCa9uNE47Eix2NCHI+KZwI8ys2prH0iqt56p3MrO3Mkj\n",
"rOYCnuDH9glI/a+dgeAGYNvM/CRAROzaxn9LUis406A0sJoOBBFxPPCXwFuBacCRwI2ZeVPtfrUz\n",
"HEqqmOEzDf6zMw1Kg6cVIwR7AP8FbAvcQrEC4sqIeH9m3tqC15fULsObBp1pUBpgLblkkJmfjoiX\n",
"AjMpbkH8DXAOYCCQqsqZBiXVaEUgmBMRBwOvA27JzDUAEeElAqmKamcatGlQUqkVgeBc4FPAbsCp\n",
"ABHxQsBmJKlqFrMjP2QHHuA7Ng1KqtV0IMjMh4C/GPo5Iv4A+CquXyBVx1qmcQuzucemQUmNtfy2\n",
"w8xcExFPz8zHWv3akibI5YkljVNb5iEwDEgV4EyDkiagnRMTSeqGzTMNrrNpUNJ4GQikfuJMg5Im\n",
"yUAg9QOXJ5bUJAOB1MtsGpTUIgYCqVfZNCiphQwEUq+xaVBSGxgIpF7iTIOS2sRAIPWCYqbBOdzD\n",
"MmcalNQOBgKpyoYvT/wFmwYltYuBQKoqmwYldZCBQKoamwYldUElA0FEzATOBw4uN90JvCMz7x/H\n",
"sXOBc4DnAk8AvwE+kJlXt6lcqXVsGpTUJVO6XUC9iNgKuJ4irOxTfj0KLIyI7cY4diZwE7AVsEdm\n",
"7gGcB1wZEce3tXCpGWuZxnXsxrdZw72clevyy4YBSZ1UuUAAvBbYDzg9Mzdm5kbgdGB34LQxjn0b\n",
"MBs4Y6jxKjP/A7gZ+Kf2lSxN0kaCO5jNV5nB7XyJ1ZztHQSSuqGKgeBEYFlmLh3akJkrgbvL50Zz\n",
"CLA+M++t234nsFdE/GErC5WaspLt+TrzWMjdLOfM/G1e5x0Ekrqlij0E+wO/aLB9KXD0GMc+CkSD\n",
"7RvLxz2B+rAgdZZNg5IqqIqBYAawrsH2tcC2ETEtM38/wrG3AydGxP6ZeUfN9gPLx+ktrFOaOJsG\n",
"JVVUFQNBM5+U/gX4a+D8iHg58DDwBoqeBIDfNlmbNDmblyde7kyDkqqoioHgIWCHBtunA4+OMjpA\n",
"Zq6LiCOBsykaCR8DFgF/B3wGuO9JB13G/E3f785SDmHppCuX6g2fadDliSV1XETMh5r3upH2q9ql\n",
"y4j4NrBXZs6r234nsC4zD5/Ea/4D8B7g6Zm5vmZ7soCzmq1ZasiZBiVVUERkZj6p366KdxlcCcwt\n",
"JxgCICJmAXsBV9TuGBGzIiJqft4mIv64wWu+CPhCbRiQ2mY9U7mJXbmKrfgZF7CGjxsGJFVdFQPB\n",
"pRS3CZ4bEVMjYgrFzIOLgU8O7RQRRwArgAtrjp0FXBMRB5b7TImIvwd2At7XmfI10BazI19lV27i\n",
"eh7g3fl43u4dBJJ6QeV6CDLz8Yg4lmLq4rspmgzvBI7OzMdqdl0HrKEIBUPWANcAV0XEWorbDW8C\n",
"Ds/M/+1E/RpQNg1K6nGV6yHoJHsI1LThTYOX2TQoqepG6iGo3AiB1DNWlU2Dy2walNT7DATSRK1n\n",
"Krcxhzt4xJkGJfULA4E0Ec40KKlPGQik8bBpUFKfMxBIoxlqGryVKaxypkFJ/ctAII2kmGnw6Szn\n",
"dpsGJfU7A4FUb/PyxI+wio+zAScXktT3DARSLZsGJQ0oA4EENg1KGngGAg02mwYlCTAQaJC5PLEk\n",
"bWIg0OCpbRp0pkFJAgwEGjQ2DUpSQwYCDQabBiVpVAYC9TebBiVpXAwE6l82DUrSuBkI1H82Nw2u\n",
"s2lQksbHQKD+spgduZkdWGHToCRNhIFA/WFz0+B9rOb8zFzc7ZIkqZcYCNTbbBqUpJYwEKh3rSqb\n",
"BpfZNChJzTIQqPesZyq3MYc7nGlQklrFQKDeUjQNbs8KrrNpUJJax0Cg3rCWadzCHO5huU2DktR6\n",
"BgJV21DT4G1MYSVfsGlQktrDQKDqsmlQkjrGQKDqsWlQkjrOQKBqcXliSeoKA4GqYXPT4DKXJ5ak\n",
"zqtkIIiImcD5wMHlpjuBd2Tm/eM4djbwAeAIYD3F7/hl4NzMfLw9FWvShs80aNOgJHVJ5QJBRGwF\n",
"XA/8Atin3HwJsDAiDhptCDkipgDfAqYCh2fmmog4ELgJmAG8o63Fa2KGL0/8xcx8qNslSdKgmtLt\n",
"Ahp4LbAfcHpmbszMjcDpwO7AaWMcuxdwAHBxZq4ByMyfANcBr2pfyZqQ9UzlJnblKrbkZ1zAGj5u\n",
"GJCk7qrcCAFwIrAsM5cObcjMlRFxd/nceaMcOzTUvGXd9i0BLxdUgU2DklRJVQwE+1NcLqi3FDh6\n",
"tAMz856I+BJwakRcnpnLIuLo8ri/bnmlGr/NyxMvt2lQkqqnioFgBrCuwfa1wLYRMS0zfz/K8a+l\n",
"aEi8NyJWA9sCb8vMS1pfqsY0fKZBlyeWpIqqYiCY9AQ0EbE1RUPiVGBuZj4QEQcB10TE/pn5t60q\n",
"UuMwvGnQmQYlqcKqGAgeAnZosH068OgYowOnUNxu+EeZ+QBAZv44Ij4CXBARX8rMW4YdcRnzN32/\n",
"O0s5hKVNVa+iafBWduZO1jnToCR1V0TMh5r3uhFUMRDcQXG3QL15FPMRjGa/8vHeuu1DP+8PDA8E\n",
"r2DRxMrTqGwalKRKycxFsPm9LiLe32i/Kt52eCUwNyLmDm2IiFkUIeGK2h0jYlZERM2mleXjXIab\n",
"W/e8Wm0t07iO3fg2D3MvZ+e6/LJhQJJ6R1RtJDcitgRuBX4OvJqip+Bi4HDgoMx8rNzvCOD7wKcz\n",
"863ltt2AnwI3Aydm5iMRsSvwXYpbEvfPzPU1/1aygLM69Kv1p+EzDV5m06AkVVtEZGZG/fbKjRCU\n",
"0wsfC2wA7i6/tgeOHgoDpXXAGmBFzbFLgUOB/wX+OyJ+ClwLfBM4sjYMqAVWsj1fZx4LuZvlnJm/\n",
"zesNA5LUmyo3QtBJjhBM0uamwUdYzaU2DUpS7xhphKCKTYWqssXsyM1MZwXX2jQoSf3DQKDx2bw8\n",
"8XJWc35mLu52SZKk1jEQaHTDZxp0eWJJ6lMGAo1sJdvzQ2ayzJkGJanfGQj0ZLVNg6v4GBu43aZB\n",
"SepvBgINZ9OgJA0kA4EKRdPg7LJp0OWJJWnAGAgGncsTS5IwEAy2VeXyxDYNStLAMxAMovVM5Tbm\n",
"cAePujyxJAkMBIPH5YklSQ0YCAbFWqZxM7O516ZBSdKTGQj63fDliW0alCQ1ZCDoZyvLpsHlNg1K\n",
"kkZnIOhHm2caXGfToCRpPAwE/caZBiVJk2Ag6Bc2DUqSmmAg6HXONChJagEDQS+zaVCS1CIGgl5U\n",
"uzyxTYOSpBYwEPQamwYlSW1gIOgVxfLEc7iHZTYNSpJazUBQdcObBr9g06AkqR0MBFVm06AkqUMM\n",
"BFXkTIOSpA4zEFSNyxNLkrrAQFAVNg1KkrrIQNBtw5cntmlQktQVBoJusmlQklQRlQwEETETOB84\n",
"uNx0J/COzLx/jONeB5wDPFD31JbAPsAxmbmwtdVOgjMNSpIqpnKBICK2Aq4HfkHxJg5wCbAwIg4a\n",
"o8kugYsy8+y61/xL4CPAotZXPEHFTIM7sMKmQUlSdVQuEACvBfYDXpqZGwEi4nTgfuA04LxRjv0+\n",
"sH2D7W8ELunqp3CbBiVJFVbFQHAisCwzlw5tyMyVEXF3+dyIgaDRm2xE7A48H3hdyysdD2calCT1\n",
"gCoGgv0pLhfUWwocPYnXeyPwncy8r5miJsWmQUlSj6hiIJgBrGuwfS2wbURMy8zfj+eFImIqcDLw\n",
"1hbWNzabBiVJPaaKgaCVb5zHl4/faOFrjs7liSVJPaiKgeAhYIcG26cDj453dKA01Ey4ccQ9LmP+\n",
"pu93ZymHsHQCr7/ZWqZxM7O5l+U2DUqSqiIi5kPNe91I+1VtJDsivg3slZnz6rbfCazLzMPH+Tqz\n",
"gSXAH47UPxARyQLOaqrg4TMNXmbToCSpyiIiMzPqt1dxhOBK4NMRMTczlwFExCxgL+CM2h3L7atG\n",
"uD7/euCGtjYT2jQoSeoTU7pdQAOXUsxMeG5ETI2IKRSzDy4GPjm0U0QcAawALqx/gYgI4BTgM22p\n",
"cD1TuYlduYot+RkXsIaPGwYkSb2sciMEmfl4RBxLMXXx3RRNhncCR2fmYzW7rgPWUISCekcBWwPX\n",
"tLzAzcsTX8cjfN2mQUlSP6hcD0EnTaiHYHPT4H2s5rOZubjN5UmS1HK91ENQLcObBr9o06AkqR8Z\n",
"CEazqmwaXGbToCSpvxkIGlnPVG5jDnfwCKv4GBu43ZkGJUn9zEBQb3PToMsTS5IGhoFgSLE88Wzu\n",
"caZBSdLgMRAMX57YpkFJ0kAyEFzNPJsGJUmDznkItuBglyeWJA2KkeYhGPhA0OikSJLUr0Z676vi\n",
"WgaSJKnDDASSJMlAIEmSDASSJAkDgSRJwkAgSZIwEEiSJAwEkiQJA4EkScJAIEmSMBBIkiQMBJIk\n",
"CQOBJEm4abafAAAJ6UlEQVTCQCBJkjAQSJIkDASSJAkDgSRJwkAgSZIwEEiSJAwEkiSJigaCiJgZ\n",
"EV+MiF+UX5dHxM4TOP6AiPh6RNweET8vX+PcdtYsSVIvq1wgiIitgOuBLYB9yq9HgYURsd04jj8c\n",
"+E/ggsx8dmbuDXwceHn7qpYkqbdFZna7hmEi4k3Ap4HdM3NpuW0WcD9wRmaeN8qxAdwNfDUz31ez\n",
"fQvgmMz8Tt3+mZnR+t9CkqRqGum9r3IjBMCJwLKhMACQmSsp3uhPHOPYI4E9gW/UbszMJ+rDQK+K\n",
"iPndrqHfeY7bz3PcGZ7n9uunc1zFQLA/sKTB9qXAfmMce3j5uG1EXBERd0XEnRHxoYjYupVFdtH8\n",
"bhcwAOZ3u4ABML/bBQyI+d0uYADM73YBrbJFtwtoYAawrsH2tRRv9NMy8/cjHLtL+fjvwCsy84cR\n",
"sS/wbeAQ4E9aXq0kSX2giiMEzTQ1DI0CfD4zfwiQmXcB5wLHRsTzmy1OkqR+VMURgoeAHRpsnw48\n",
"OsroAGweWfhJ3fahnw8Bvl/7RERUq6tyHCLi/d2uod95jtvPc9wZnuf265dzXMVAcAewV4Pt84A7\n",
"xzj25+Vj/cjHhkbbvcNAkqRCFS8ZXAnMjYi5QxvK2w73Aq6o3TEiZpW3Gg75FsWb//51r7lv+fjf\n",
"rS9XkqTeV8V5CLYEbqX4tP9qip6CiynuIDgoMx8r9zuCYvj/05n51prjPwqcBLwgM38eEXOA7wH/\n",
"k5l/2tFfRpKkHlG5SwaZ+XhEHAucTzH3QFJcKjh6KAyU1gFrgBV1L/Euij6EqyPiCYrf8QqgL67x\n",
"SJLUDpUbIZDUuyJiNvBZ4LjMrOIlyb7geW6/VpzjiPggcCbw+sz8XCvrawf/kCrAxZzar5lzHBGz\n",
"I+LfynP704j4WUS8t7y8pVJEnADcCOzGBG8fjogtI+ID5Tm+MyJuLC8Lqs5kz3P5d3xW+Td8Z3mu\n",
"ryjnalGNZv6Wa17jGcDflcf3xCdvA0GXuZhT+zVzjiNiCkWz6qHA4Zl5AEVvy5nAP7Wz7h70TuBo\n",
"4GZgonfwfILib/bIzNwPuAS4LiIOaG2JfWGy5/n9wKuAPyvP8YEUTdi3GAqepJm/5SEfBm5o4viO\n",
"MxB032sppmQ+PTM3ZuZG4HRgd+C00Q4s77C4GPhUZi6seeozYx07YCZ9jinubjkAuDgz1wBk5k+A\n",
"6yj+46rNjqxdg2S8ImJP4E3AOZn5a4DMvJhiCvMPtbTC/jCp80zxKfXczLwfoJzT5QxgG+DNrSuv\n",
"L0z2HAMQEQdTNMJ/omUVdYCBoPtczKn9mjnHT5SP9ZcHtgQeb1WB/SAn35D0MopPUQvrti8EjouI\n",
"bZsqrM80cZ7fRnFNvNYD5eNTJ19R/2niHA/5KMUo4voWlNMxBoLuczGn9pv0Oc7Me4AvAacOzY0R\n",
"EUdTDCd650pr7E8xdL28bvsSNl/mUZMyc0ODN7o9ysdFHS6nb0XEnwPTMvOybtcyUZW77XAAuZhT\n",
"+zVzjqG45HA+cG9ErAa2Bd6WmZe0vtSBNAN4rMGb1drycccO1zNI3gzcBXy+24X0g7LR+BzglG7X\n",
"MhmOEHSfizm136TPcTnSsogiYM3NzJ0pRgfOioiPtaY8qfMi4hjgFRQfJrz81RqnAXdl5k3dLmQy\n",
"DATd1+7FnNTcOT4FOAJ4V2Y+AJCZPwY+Arw9Ip7b6mIH0EPAdnXTkEPxvw/ArztcT98r7964FHhx\n",
"Zv6iy+X0hYh4KkWT5hmNnu5wOZNiIOi+OygWbqrX8sWcBlgz53iox+Deuu1DP9evm6GJ+ynF3+ou\n",
"ddvnUTRu3t3xivpYROwPfA14ZWbe3O16+sjzKJqQL4+IH0fEj4F/LZ87u9z23u6VNzbfMLrPxZza\n",
"r5lzvLJ8nMtwc+ue13AjXqZpcI6/Vu5/VN2uRwHX1U1ZruEmcp6HwsBVwGuGhrXLCYs+1d4ye9q4\n",
"znFmXpuZu2bmQUNfwBvLXd9XbvtgJwqeLANB911K8Sn13IiYWk6Ecw6wGPjk0E7lrG0rgAuHtmXm\n",
"rygmIXpjROxd7jcH+AeK/5B+r1O/RMVdyiTPcXnsOuCDEbF9ud+uwN8D9wDXdqD+XtRwiHSEv+N7\n",
"KObOeHdE7Fju93qKEYL3tL/Unjbu8xwR+1FMlPMdYPeIeE1EvAZ4JcXty2ps3Od4lGN74pKBdxl0\n",
"mYs5tV8z5zgzl0bEocAC4L8jYj3FHATfBD6QmT11n3E7RcRFwAsp7hrIiFhCca73rGlaG+nv+O0U\n",
"f7M3RsTjFHcYHJeZd3Sk+B7SxHleADwNOLX8qrWojSX3nCb/lomInYAfUjR+J3BeRCyguEzzo/b/\n",
"BpPj4kaSJMlLBpIkyUAgSZIwEEiSJAwEkiQJA4EkScJAIEmSMBBIkiQMBJIkCQOBpA4pl5Ke7LHb\n",
"tLIWSU9mIJDUdhGxB/DHddumRMSPIuIZ43iJHSPi5PZUJwkMBJLarFxM6g2Z+Y26p14K3AecOcJx\n",
"W0TENbBpIa/HI+I5bS1WGmAGAknt9hLgB7UbyiVjd6FY1OivRhgleC5wb83Pl/HkRXkktYiBQFK7\n",
"vQSoX4r7xcBVmbkC+Dzw7tonI+JPgH8ENkbEUQCZuQHYISKmtb9kafAYCCS125zMXDv0Qzk6sGtm\n",
"Li83/V/g1RGx89A+mfkdijXkz87MhTWv9SCwRwdqlgaOgUBSu21Z9/PxwDVDP5SjBF8AzhjaFhFb\n",
"AdvVBonSOmD7NtUpDTQDgaR22zD0TTk6MC8zl9Xtcw5wUkTMKX8+GLg1IraNiGNr9nsKsLqt1UoD\n",
"ykAgqd3ujYiZ5fcvBt4fEffVfgE/BLYG3lXut5riksEJQO0lgxnA4g7VLQ2UyMxu1yCph0XEDOB0\n",
"YEW56Z7M/GbN88cAO2XmF5v8d7YBzsnMv23mdSQ15giBpGZ9CvhKZp4PPAb8Ze2TmXkDcGB5uaAZ\n",
"JwEfbfI1JI3AEQJJTYmIHwNLgUuBW4DH6psBy1GEYzLzPyb5bzyT4s6E7zZXraSRGAgkNSUiDgU+\n",
"QjGR0CrgkMy08U/qMV4ykDRpEfEHmfmjzJwPzKFoBjyyu1VJmgwDgaRJiYhtgV9FxAvLTb8BVgI3\n",
"dq8qSZPlJQNJkxYRC4BlwHbArsBlmXlrRBwOHAP8mmJ2wQOBJRShYVfggcn2E0hqjy26XYCk3pWZ\n",
"C0Z4alvgIeApmXlRRCRwfGa+sew5eBVgIJAqxEsGklouM/8TmA9cXm46gmJ6YoAXADd1oSxJozAQ\n",
"SGqXZ2bm/5TfH8bmEHA8cENEHNadsiQ1YiCQ1HIRsSvwo/L7rYGVmbm+fHox8CLg1i6VJ6kBmwol\n",
"SZIjBJIkyUAgSZIwEEiSJAwEkiQJA4EkScJAIEmSMBBIkiQMBJIkCfj/LVJRjcnQte0AAAAASUVO\n",
"RK5CYII=\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"s_m_values = numpy.linspace(0.5, 1.5)\n",
"bound = (2.0*s_m_values + 1.0)/3.0\n",
"pyplot.figure(figsize=(8,6))\n",
"pyplot.fill_between(s_m_values, numpy.ones_like(s_m_values), bound,\n",
" facecolor='green', alpha=0.5)\n",
"#pyplot.ylim(1.0, 2.5)\n",
"pyplot.xlim(0.5, 1.5)\n",
"pyplot.xlabel(r\"$s_m^{(\\Delta t)}$\")\n",
"pyplot.ylabel(r\"$s_m^{(\\Delta t / 2)}$\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we want the algorithm to behave in the way expected of Euler's method, then measuring a single convergence rate is useful but not sufficient. Instead we must measure the convergence rate twice, once with base timestep $\\dt$ and once with base timestep $\\dt/2$. The difference between the measured slope $\\sm{\\dt}$ and the expected slope $s_e=1$ must go down by at least $2/3$ when the timestep is reduced to $\\dt/2$."
]
}
],
"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
}