Tutorial-RK_Butcher_Table_Validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Validating Runge Kutta Butcher tables using Truncated Taylor Series\n",
    "## Authors: Zach Etienne, Brandon Clark, & Gabriel M Steward\n",
    "\n",
    "\n",
    "## This tutorial notebook is designed to validate the Butcher tables contained within the Butcher dictionary constructed in the [RK Butcher Table Dictionary](Tutorial-RK_Butcher_Table_Dictionary.ipynb) NRPy+ module.  \n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [MoLtimestepping/RK_Butcher_Table_Validation.py](../edit/MoLtimestepping/RK_Butcher_Table_Validation.py) stores the `Validate` function for validating convergence orders for Runge Kutta methods.\n",
    "* [MoLtimestepping/RK_Butcher_Table_Dictionary.py](../edit/MoLtimestepping/RK_Butcher_Table_Dictionary.py) [\\[**tutorial**\\]](Tutorial-RK_Butcher_Table_Dictionary.ipynb) accesses the Butcher table dictionary `Butcher_dict` for known explicit Runge Kutta methods.\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "Starting with the ODE (ordinary differential equation) initial value problem:\n",
    "$$\n",
    "y'(t) = f(y,t)\\ \\ \\ y\\left(t=0\\right)=y_0,\n",
    "$$\n",
    "for various choices of $f(y,t)$, this module validates the Runge Kutta (RK) methods coded in [RK_Butcher_Table_Dictionary.py](../edit/MoLtimestepping/RK_Butcher_Table_Dictionary.py) [**tutorial notebook**](Tutorial-RK_Butcher_Table_Dictionary.ipynb) as follows.\n",
    "\n",
    "Given $y_0$ and a smooth $f(y,t)$, all explicit RK methods provide an estimate for $y_1 = y\\left(\\Delta t\\right)$, with an error term that is proportional to $\\left(\\Delta t\\right)^m$, where $m$ is an integer typically greater than zero. This error term corresponds to the *local* truncation error. For RK4, for example, while the *total accumulated truncation error* (i.e., the accumulated error at a fixed final time $t_f$) is proportional to $\\left(\\Delta t\\right)^4$, the *local* truncation error (i.e., the error after one arbitrarily chosen timestep $\\Delta t$) is proportional to $\\left(\\Delta t\\right)^5$.\n",
    "\n",
    "If the exact solution $y(t)$ is known as a closed-form expression, then $y\\left(\\Delta t\\right)$ can be *separately* written as a Taylor expansion about $y(t=0)$:\n",
    "\n",
    "$$\n",
    "y\\left(\\Delta t\\right) = \\sum_{n=0}^\\infty \\frac{y^{(n)}(t=0)}{n!} \\left(\\Delta t\\right)^n,\n",
    "$$\n",
    "where $y^{(n)}(t=0)$ is the $n$th derivative of $y(t)$ evaluated at $t=0$.\n",
    "\n",
    "The above expression will be known exactly. Furthermore, if one chooses a numerical value for $y_0$ *and leaves $\\Delta t$ unspecified*, any explicit RK method will provide an estimate for $y\\left(\\Delta t\\right)$ of the form\n",
    "\n",
    "$$\n",
    "y\\left(\\Delta t\\right) = \\sum_{n=0}^\\infty a_n \\left(\\Delta t\\right)^n,\n",
    "$$\n",
    "where $a_n$ *must* match the Taylor expansion of the *exact* solution at least up to and including terms proportional to $\\left(\\Delta t\\right)^m$, where $m$ is the order of the local truncation error. If this is *not* the case, then the Butcher table was almost certainly *not* typed correctly.\n",
    "\n",
    "Therefore, comparing the numerical result with unspecified $\\Delta t$ against the exact Taylor series provides a convenient (though not perfectly robust) means to verify that the Butcher table for a given RK method was typed correctly. Multiple typos in the Butcher tables were found using this approach.\n",
    "\n",
    "**Example from Z. Etienne's MATH 521 (Numerical Analysis) lecture notes:**\n",
    "\n",
    "Consider the ODE\n",
    "$$\n",
    "y' = y - 2 t e^{-2t},\\quad y(0)=y(t_0)=0.\n",
    "$$\n",
    "\n",
    "* Solve this ODE exactly, then Taylor expand the solution about $t=0$ to\n",
    "approximate the solution at $y(t=\\Delta t)$ to the fifth order in $\\Delta\n",
    "t$.\n",
    "* Next, solve this ODE using Heun's method (second order in total accumulated truncation error, third order in local truncation error) {\\it by hand} with a step size of\n",
    "$\\Delta t$ to find $y(\\Delta t)$. Confirm that the solution obtained\n",
    "when using Heun's method has an error term that is at worst\n",
    "$\\mathcal{O}\\left((\\Delta t)^3\\right)$. If the dominant error is\n",
    "proportional to a higher power of $\\Delta t$, explain the discrepancy.\n",
    "\n",
    "* Finally, solve this ODE using the Ralston method {\\it by hand}\n",
    "  with a step size of $\\Delta t$ to find $y(\\Delta t)$. Is the\n",
    "  coefficient on the dominant error term closer to the exact solution\n",
    "  than Heun's method?\n",
    "\n",
    "We can solve this equation via the method of integrating factors,\n",
    "which states that ODEs of the form:\n",
    "$$\n",
    "y'(t) + p(t) y(t) = g(t)\n",
    "$$\n",
    "are solved via \n",
    "$$\n",
    "y(t) = \\frac{1}{\\mu(t)} \\left[ \\int \\mu(s) g(s) ds + c \\right],\n",
    "$$\n",
    "where the integrating factor $\\mu(t)$ is given by\n",
    "$$\n",
    "\\mu(t) = \\exp\\left(\\int p(t) dt\\right).\n",
    "$$\n",
    "\n",
    "Here, $p(t)=-1$ and $g(t) = - 2 t e^{-2t}$. Then\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu(t) = \\exp\\left(-\\int dt\\right) = e^{-t+c} = k e^{-t}\n",
    "\\end{equation}\n",
    "and\n",
    "\n",
    "\\begin{align}\n",
    "y(t) &= e^t/k  \\left[ \\int k e^{-s} (- 2 s e^{-2s}) ds + c \\right] = -2 e^t \\left[ \\int s e^{-3s} ds + c' \\right] \\\\\n",
    "&= -2 e^t \\left[ e^{-3 t} \\left(-\\frac{t}{3}-\\frac{1}{9}\\right) + c' \\right] = -2 e^{-2t} \\left(-\\frac{t}{3}-\\frac{1}{9}\\right) -2 c' e^t \\\\\n",
    "&= e^{-2t} \\left(2\\frac{t}{3}+\\frac{2}{9}\\right) + c'' e^t. \\\\\n",
    "\\end{align}\n",
    "\n",
    "If $y(0)=0$ then we can compute the integration constant $c''$, and\n",
    "$y(t)$ becomes\n",
    "$$\n",
    "y(t) = \\frac{2}{9} e^{-2 t} \\left(3 t + 1 - e^{3 t}\\right).\n",
    "$$\n",
    "\n",
    "The Taylor Series expansion of the exact solution about $t=0$\n",
    "evaluated at $y(\\Delta t)$ yields\n",
    "$$\n",
    "y(\\Delta t) = -(\\Delta t)^2+(\\Delta t)^3-\\frac{3 (\\Delta t)^4}{4}+\\frac{23 (\\Delta\n",
    "  t)^5}{60}-\\frac{19 (\\Delta t)^6}{120}+O\\left((\\Delta t)^7\\right).\n",
    "$$\n",
    "\n",
    "Next, we evaluate $y(\\Delta t)$ using Heun's method. We know $y(0)=y_0=0$ and\n",
    "$f(y,t)=y - 2 t e^{-2t}$, so\n",
    "\\begin{align}\n",
    "k_1 &= \\Delta t f(y(0),0) \\\\\n",
    "    &= \\Delta t \\times 0 \\\\\n",
    "    &= 0 \\\\\n",
    "k_2 &= \\Delta t f(y(0)+k_1,0+\\Delta t) \\\\\n",
    "   &= \\Delta t f(y(0)+0,0+\\Delta t) \\\\\n",
    "   &= \\Delta t (-2 \\Delta t e^{-2\\Delta t}) \\\\\n",
    "   &= -2 (\\Delta t)^2 e^{-2\\Delta t} \\\\\n",
    "y(\\Delta t) &= y_0 + \\frac{1}{2} (k_1 + k_2) + \\mathcal{O}\\left((\\Delta t)^3\\right) \\\\\n",
    "&= 0 - (\\Delta t)^2 e^{-2\\Delta t} \\\\\n",
    "&= - (\\Delta t)^2 ( 1 - 2 \\Delta t + 2 (\\Delta t)^2 + ...) \\\\\n",
    "&= - (\\Delta t)^2 + 2 (\\Delta t)^3 + \\mathcal{O}\\left((\\Delta t)^4\\right).\n",
    "\\end{align}\n",
    "\n",
    "Thus the coefficient on the $(\\Delta t)^3$ term is wrong, but\n",
    "this is completely consistent with the fact that our stepping\n",
    "scheme is only third-order accurate in $\\Delta t$.\n",
    "\n",
    "In the below approach, the RK result is subtracted from the exact Taylor series result, as a check to determine whether the RK Butcher table was coded correctly; if it was not, then the odds are good that the RK results will not match the expected local truncation error order. Multiple $f(y,t)$ are coded below to improve the robustness of this test.\n",
    "\n",
    "As NRPy+'s butcher tables also contain the information to run Adams-Bashforth (AB) methods, those too shall also be validated. Fortunately, the exact same validation method that we use for the RK methods also functions for the AB methods, all that needs to change is the specific implementation code. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): Initialize needed Python/NRPy+ modules\n",
    "1. [Step 2](#table_validate) Validate Convergence Order of Butcher Tables\n",
    "    1. [Step 2.a](#rhs): Defining the right-hand side of the ODE\n",
    "    1. [Step 2.b](#validfunc): Defining a Validation Function\n",
    "    1. [Step 2.c](#rkvalid): Validating RK Methods against ODEs\n",
    "    1. [Step 2.d](#arkvalid): Validating Inherently Adaptive RK Methods against ODEs\n",
    "    1. [Step 2.e](#abvalid): Validating AB Methods against ODEs\n",
    "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize needed Python/NRPy+ modules [Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "Let's start by importing all the needed modules from Python/NRPy+:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:16:35.516909Z",
     "iopub.status.busy": "2021-03-07T17:16:35.516113Z",
     "iopub.status.idle": "2021-03-07T17:16:36.042355Z",
     "shell.execute_reply": "2021-03-07T17:16:36.042843Z"
    }
   },
   "outputs": [],
   "source": [
    "import sympy as sp              # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import numpy as np              # NumPy: A numerical methods module for Python\n",
    "from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='table_validate'></a>\n",
    "\n",
    "# Step 2:  Validate Convergence Order of Butcher Tables [Back to [top](#toc)\\]\n",
    "$$\\label{table_validate}$$\n",
    "\n",
    "\n",
    "Each Butcher table/Runge Kutta method is tested by solving an ODE. Comparing the Taylor series expansions of the exact solution and the numerical solution as discussed in the **Introduction** above will confirm whether the method converges to the appropriate order. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rhs'></a>\n",
    "\n",
    "## Step 2.a:  Defining the right-hand side of the ODE [Back to [top](#toc)\\]\n",
    "$$\\label{rhs}$$\n",
    "\n",
    "Consider the form of ODE $y'=f(y,t)$. The following begins to construct a dictionary `rhs_dict` of right-hand side functions for us to validate explicit Runge Kutta methods. The most up-to-date catalog of functions stored in `rhs_dict` can be found in the [RK_Butcher_Table_Validation.py](../edit/MoLtimestepping/RK_Butcher_Table_Validation.py) module. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:16:36.049940Z",
     "iopub.status.busy": "2021-03-07T17:16:36.049261Z",
     "iopub.status.idle": "2021-03-07T17:16:36.051301Z",
     "shell.execute_reply": "2021-03-07T17:16:36.051767Z"
    }
   },
   "outputs": [],
   "source": [
    "def fypt(y,t): # Yields expected convergence order for all cases\n",
    "    #            except DP6 which converge to higher order (7, respectively)\n",
    "    return y+t\n",
    "\n",
    "def fy(y,t): # Yields expected convergence order for all cases\n",
    "    return y\n",
    "\n",
    "def feypt(y,t): # Yields expected convergence order for all cases\n",
    "    return sp.exp(1.0*(y+t))\n",
    "\n",
    "def ftpoly6(y,t): # Yields expected convergence order for all cases, L6 has 0 error\n",
    "    return 2*t**6-389*t**5+15*t**4-22*t**3+81*t**2-t+42\n",
    "rhs_dict = {'ypt':fypt, 'y':fy, 'eypt':feypt, 'tpoly6':ftpoly6}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='validfunc'></a>\n",
    "\n",
    "## Step 2.b: Defining a Validation Function  [Back to [top](#toc)\\]\n",
    "$$\\label{validfunc}$$\n",
    "\n",
    "To validate each Butcher table we compare the exact solutions to ODEs with the numerical solutions using the Runge Kutta scheme built into each Butcher table. The following is a function that\n",
    "\n",
    "1. solves the ODE exactly,\n",
    "2. solves the ODE numerically for a given Butcher table, and\n",
    "3. compares the two solutions and checks for the order of convergence by returning their difference.\n",
    "\n",
    "The `Validate()` function inputs a specified `Butcher_key`, the starting guess solution and time `y_n`, `t_n`, and the right-hand side of the ODE corresponding to a specified initial value problem, `rhs_key`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:16:36.063763Z",
     "iopub.status.busy": "2021-03-07T17:16:36.062979Z",
     "iopub.status.idle": "2021-03-07T17:16:36.066476Z",
     "shell.execute_reply": "2021-03-07T17:16:36.065801Z"
    }
   },
   "outputs": [],
   "source": [
    "from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
    "def Validate(Butcher_key, yn, tn, rhs_key):\n",
    "    # 1. First we solve the ODE exactly\n",
    "    y = sp.Function('y')\n",
    "    sol = sp.dsolve(sp.Eq(y(t).diff(t), rhs_dict[rhs_key](y(t), t)), y(t)).rhs\n",
    "    constants = sp.solve([sol.subs(t,tn)-yn])\n",
    "    exact = sol.subs(constants)\n",
    "\n",
    "    # 2. Now we solve the ODE numerically using specified Butcher table\n",
    "\n",
    "    # Access the requested Butcher table\n",
    "    Butcher = Butcher_dict[Butcher_key][0]\n",
    "    # Determine number of predictor-corrector steps\n",
    "    L = len(Butcher)-1\n",
    "    # Set a temporary array for update values\n",
    "    k = np.zeros(L, dtype=object)\n",
    "    # Initialize intermediate variable\n",
    "    yhat = 0\n",
    "    # Initialize the updated solution\n",
    "    ynp1 = 0\n",
    "    for i in range(L):\n",
    "        #Initialize and approximate update for solution\n",
    "        yhat = yn\n",
    "        for j in range(i):\n",
    "            # Update yhat for solution using a_ij Butcher table coefficients\n",
    "            yhat += Butcher[i][j+1]*k[j]\n",
    "            if Butcher_key == \"DP8\" or Butcher_key == \"L6\":\n",
    "                yhat = 1.0*sp.N(yhat,20) # Otherwise the adding of fractions kills performance.\n",
    "        # Determine the next corrector variable k_i using c_i Butcher table coefficients\n",
    "        k[i] = dt*rhs_dict[rhs_key](yhat, tn + Butcher[i][0]*dt)\n",
    "        # Update the solution at the next iteration ynp1 using Butcher table coefficients\n",
    "        ynp1 += Butcher[L][i+1]*k[i]\n",
    "    # Finish determining the solution for the next iteration\n",
    "    ynp1 += yn\n",
    "\n",
    "    # Determine the order of the RK method\n",
    "    order = Butcher_dict[Butcher_key][1]+2\n",
    "    # Produces Taylor series of exact solution at t=tn about t = 0 with the specified order\n",
    "    exact_series = sp.series(exact.subs(t, dt),dt, 0, order)\n",
    "    num_series = sp.series(ynp1, dt, 0, order)\n",
    "    diff = exact_series-num_series\n",
    "    return diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rkvalid'></a>\n",
    "\n",
    "## Step 2.c: Validating RK Methods against ODEs  [Back to [top](#toc)\\]\n",
    "$$\\label{rkvalid}$$\n",
    "\n",
    "The following makes use of the `Validate()` function above to demonstrate that each method within the Butcher table dictionary converges to the expected order for the given right-hand side expression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:16:36.111715Z",
     "iopub.status.busy": "2021-03-07T17:16:36.080988Z",
     "iopub.status.idle": "2021-03-07T17:16:58.312416Z",
     "shell.execute_reply": "2021-03-07T17:16:58.311759Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RK method: \"Euler\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^2) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  2    ⎛  3⎞\n",
      "dt  + O⎝dt ⎠\n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"RK2 Heun\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^3) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  3         \n",
      "dt     ⎛  4⎞\n",
      "─── + O⎝dt ⎠\n",
      " 3          \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"RK2 MP\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^3) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  3         \n",
      "dt     ⎛  4⎞\n",
      "─── + O⎝dt ⎠\n",
      " 3          \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"RK2 Ralston\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^3) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  3         \n",
      "dt     ⎛  4⎞\n",
      "─── + O⎝dt ⎠\n",
      " 3          \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"RK3\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^4) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  4         \n",
      "dt     ⎛  5⎞\n",
      "─── + O⎝dt ⎠\n",
      " 12         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"RK3 Heun\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^4) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  4         \n",
      "dt     ⎛  5⎞\n",
      "─── + O⎝dt ⎠\n",
      " 12         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"RK3 Ralston\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^4) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  4         \n",
      "dt     ⎛  5⎞\n",
      "─── + O⎝dt ⎠\n",
      " 12         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"SSPRK3\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^4) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  4         \n",
      "dt     ⎛  5⎞\n",
      "─── + O⎝dt ⎠\n",
      " 12         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"RK4\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^5) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  5         \n",
      "dt     ⎛  6⎞\n",
      "─── + O⎝dt ⎠\n",
      " 60         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"DP5\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^6) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "    6          \n",
      "  dt      ⎛  7⎞\n",
      "- ──── + O⎝dt ⎠\n",
      "  1800         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"DP5alt\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^6) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "     6         \n",
      "13⋅dt     ⎛  7⎞\n",
      "────── + O⎝dt ⎠\n",
      "231000         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"CK5\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^6) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  6          \n",
      "dt      ⎛  7⎞\n",
      "──── + O⎝dt ⎠\n",
      "3600         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"DP6\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^7) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      " ⎛  8⎞\n",
      "O⎝dt ⎠\n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"L6\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^7) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "                              5                               6               \n",
      "- 1.0587911840678754238e-22⋅dt  - 2.6469779601696885596e-23⋅dt  + 0.0013227513\n",
      "\n",
      "               7    ⎛  8⎞\n",
      "227513227513⋅dt  + O⎝dt ⎠\n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"DP8\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^9) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "                                                           2                  \n",
      "3.6854403535034607753e-18⋅dt + 5.8394451383711465375e-18⋅dt  + 3.7764963953332\n",
      "\n",
      "             3                              4                               5 \n",
      "980617e-18⋅dt  + 9.542884942003761195e-19⋅dt  + 1.2718729098615353529e-19⋅dt  \n",
      "\n",
      "                              6                               7               \n",
      "+ 3.9082629581905451582e-20⋅dt  + 4.8075737201581968464e-21⋅dt  + 5.1688448526\n",
      "\n",
      "                8                              9    ⎛  10⎞\n",
      "907316834e-22⋅dt  + 7.2078645877627939543e-9⋅dt  + O⎝dt  ⎠\n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "t, dt = sp.symbols('t dt')\n",
    "# Set initial conditions\n",
    "t0 = 0\n",
    "y0 = 1\n",
    "# Set RHS of ODE\n",
    "function = 'ypt'# This can be changed, just be careful that the initial conditions are satisfied\n",
    "for key,value in Butcher_dict.items():\n",
    "    if key not in {\"AHE\", \"ABS\", \"ARKF\", \"ACK\", \"ADP5\", \"ADP8\", \"AB\"}:\n",
    "        print(\"RK method: \\\"\"+str(key)+\"\\\".\")\n",
    "        y = sp.Function('y')\n",
    "        print(\" When solving y'(t) = \"+str(rhs_dict[function](y(t),t))+\", y(\"+str(t0)+\")=\"+str(y0)+\",\")\n",
    "        local_truncation_order = list(value)[1]+1\n",
    "        print(\" the first nonzero term should have local truncation error proportional to O(dt^\"+str(local_truncation_order)+\") or a higher power of dt.\")\n",
    "        print(\"Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\")\n",
    "        sp.pretty_print(Validate(key, y0, t0, function))\n",
    "    #     print(\"\\n\")\n",
    "        print(\" (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\\n\")\n",
    "        if key == \"DP8\":\n",
    "            break # Keep this code from trying to validate Adaptive and Adams-Bashforth methods\n",
    "            # They need to be read differently."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='arkvalid'></a>\n",
    "\n",
    "## Step 2.d: Validating Inherently Adaptive RK Methods against ODEs  [Back to [top](#toc)\\]\n",
    "$$\\label{arkvalid}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code validates inherently adaptive RK methods. Each method has two validation checks as each method has, individually, two different methods inside of it. One should have the first error term's order be equal to that of the method's order, and the other check should have the first error term's order be the method's order plus one. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ValidateARK1(Butcher_key, yn, tn, rhs_key):\n",
    "    # 1. First we solve the ODE exactly\n",
    "    y = sp.Function('y')\n",
    "    sol = sp.dsolve(sp.Eq(y(t).diff(t), rhs_dict[rhs_key](y(t), t)), y(t)).rhs\n",
    "    constants = sp.solve([sol.subs(t,tn)-yn])\n",
    "    exact = sol.subs(constants)\n",
    "\n",
    "    # 2. Now we solve the ODE numerically using specified Butcher table\n",
    "\n",
    "    # Access the requested Butcher table\n",
    "    Butcher = Butcher_dict[Butcher_key][0]\n",
    "    # Determine number of predictor-corrector steps\n",
    "    L = len(Butcher)-1\n",
    "    # Set a temporary array for update values\n",
    "    k = np.zeros(L, dtype=object)\n",
    "    # Initialize intermediate variable\n",
    "    yhat = 0\n",
    "    # Initialize the updated solution\n",
    "    ynp1 = 0\n",
    "    for i in range(L-1):\n",
    "        #Initialize and approximate update for solution\n",
    "        yhat = yn\n",
    "        for j in range(i):\n",
    "            # Update yhat for solution using a_ij Butcher table coefficients\n",
    "            yhat += Butcher[i][j+1]*k[j]\n",
    "            if Butcher_key == \"DP8\" or Butcher_key == \"L6\":\n",
    "                yhat = 1.0*sp.N(yhat,20) # Otherwise the adding of fractions kills performance.\n",
    "        # Determine the next corrector variable k_i using c_i Butcher table coefficients\n",
    "        k[i] = dt*rhs_dict[rhs_key](yhat, tn + Butcher[i][0]*dt)\n",
    "        # Update the solution at the next iteration ynp1 using Butcher table coefficients\n",
    "        ynp1 += Butcher[L][i+1]*k[i]\n",
    "    # Finish determining the solution for the next iteration\n",
    "    ynp1 += yn\n",
    "\n",
    "    # Determine the order of the RK method\n",
    "    order = Butcher_dict[Butcher_key][1]+2\n",
    "    # Produces Taylor series of exact solution at t=tn about t = 0 with the specified order\n",
    "    exact_series = sp.series(exact.subs(t, dt),dt, 0, order)\n",
    "    num_series = sp.series(ynp1, dt, 0, order)\n",
    "    diff = exact_series-num_series\n",
    "    return diff\n",
    "\n",
    "def ValidateARK2(Butcher_key, yn, tn, rhs_key):\n",
    "    # 1. First we solve the ODE exactly\n",
    "    y = sp.Function('y')\n",
    "    sol = sp.dsolve(sp.Eq(y(t).diff(t), rhs_dict[rhs_key](y(t), t)), y(t)).rhs\n",
    "    constants = sp.solve([sol.subs(t,tn)-yn])\n",
    "    exact = sol.subs(constants)\n",
    "\n",
    "    # 2. Now we solve the ODE numerically using specified Butcher table\n",
    "\n",
    "    # Access the requested Butcher table\n",
    "    Butcher = Butcher_dict[Butcher_key][0]\n",
    "    # Determine number of predictor-corrector steps\n",
    "    L = len(Butcher)-1\n",
    "    # Set a temporary array for update values\n",
    "    k = np.zeros(L, dtype=object)\n",
    "    # Initialize intermediate variable\n",
    "    yhat = 0\n",
    "    # Initialize the updated solution\n",
    "    ynp1 = 0\n",
    "    for i in range(L-1):\n",
    "        #Initialize and approximate update for solution\n",
    "        yhat = yn\n",
    "        for j in range(i):\n",
    "            # Update yhat for solution using a_ij Butcher table coefficients\n",
    "            yhat += Butcher[i][j+1]*k[j]\n",
    "            if Butcher_key == \"DP8\" or Butcher_key == \"L6\":\n",
    "                yhat = 1.0*sp.N(yhat,20) # Otherwise the adding of fractions kills performance.\n",
    "        # Determine the next corrector variable k_i using c_i Butcher table coefficients\n",
    "        k[i] = dt*rhs_dict[rhs_key](yhat, tn + Butcher[i][0]*dt)\n",
    "        # Update the solution at the next iteration ynp1 using Butcher table coefficients\n",
    "        ynp1 += Butcher[L-1][i+1]*k[i]\n",
    "    # Finish determining the solution for the next iteration\n",
    "    ynp1 += yn\n",
    "\n",
    "    # Determine the order of the RK method\n",
    "    order = Butcher_dict[Butcher_key][1]+2\n",
    "    # Produces Taylor series of exact solution at t=tn about t = 0 with the specified order\n",
    "    exact_series = sp.series(exact.subs(t, dt),dt, 0, order)\n",
    "    num_series = sp.series(ynp1, dt, 0, order)\n",
    "    diff = exact_series-num_series\n",
    "    return diff\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RK method: \"AHE\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first calculation's first nonzero term should have local truncation error proportional to O(dt^2) or a higher power of dt.\n",
      " the second calculation's first nonzero term should have local truncation error proportional to O(dt^3) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "        3         \n",
      "  2   dt     ⎛  4⎞\n",
      "dt  + ─── + O⎝dt ⎠\n",
      "       3          \n",
      "  3         \n",
      "dt     ⎛  4⎞\n",
      "─── + O⎝dt ⎠\n",
      " 3          \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"ABS\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first calculation's first nonzero term should have local truncation error proportional to O(dt^3) or a higher power of dt.\n",
      " the second calculation's first nonzero term should have local truncation error proportional to O(dt^4) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "    3     4         \n",
      "  dt    dt     ⎛  5⎞\n",
      "- ─── + ─── + O⎝dt ⎠\n",
      "   24    24         \n",
      "  4         \n",
      "dt     ⎛  5⎞\n",
      "─── + O⎝dt ⎠\n",
      " 12         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"ARKF\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first calculation's first nonzero term should have local truncation error proportional to O(dt^5) or a higher power of dt.\n",
      " the second calculation's first nonzero term should have local truncation error proportional to O(dt^6) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "    5     6         \n",
      "  dt    dt     ⎛  7⎞\n",
      "- ─── + ─── + O⎝dt ⎠\n",
      "  390   360         \n",
      "     6         \n",
      "17⋅dt     ⎛  7⎞\n",
      "────── + O⎝dt ⎠\n",
      " 9360          \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"ACK\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first calculation's first nonzero term should have local truncation error proportional to O(dt^5) or a higher power of dt.\n",
      " the second calculation's first nonzero term should have local truncation error proportional to O(dt^6) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "        5          6         \n",
      "  277⋅dt    4541⋅dt     ⎛  7⎞\n",
      "- ─────── + ──────── + O⎝dt ⎠\n",
      "   614400   7372800          \n",
      "  6          \n",
      "dt      ⎛  7⎞\n",
      "──── + O⎝dt ⎠\n",
      "3600         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"ADP5\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first calculation's first nonzero term should have local truncation error proportional to O(dt^5) or a higher power of dt.\n",
      " the second calculation's first nonzero term should have local truncation error proportional to O(dt^6) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "       5        6         \n",
      "  97⋅dt    17⋅dt     ⎛  7⎞\n",
      "- ────── + ────── + O⎝dt ⎠\n",
      "  60000    180000         \n",
      "    6          \n",
      "  dt      ⎛  7⎞\n",
      "- ──── + O⎝dt ⎠\n",
      "  1800         \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "RK method: \"ADP8\".\n",
      " When solving y'(t) = t + y(t), y(0)=1,\n",
      " the first calculation's first nonzero term should have local truncation error proportional to O(dt^8) or a higher power of dt.\n",
      " the second calculation's first nonzero term should have local truncation error proportional to O(dt^9) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "                                                  2                          3\n",
      "-7.71097267488672e-19⋅dt + 2.91023006428162e-18⋅dt  + 5.29778138968287e-19⋅dt \n",
      "\n",
      "                         4                          5                         \n",
      " - 1.5349356222021e-19⋅dt  + 1.65528504999405e-19⋅dt  + 4.19247608610368e-20⋅d\n",
      "\n",
      " 6                          7                         8                       \n",
      "t  + 6.27844437078994e-21⋅dt  - 4.85333183539141e-7⋅dt  + 3.49344710134931e-7⋅\n",
      "\n",
      "  9    ⎛  10⎞\n",
      "dt  + O⎝dt  ⎠\n",
      "                                                 2                          3 \n",
      "3.68531467298237e-18⋅dt + 5.84980541251108e-18⋅dt  + 3.78072846284061e-18⋅dt  \n",
      "\n",
      "                         4                          5                         \n",
      "+ 9.53606673846474e-19⋅dt  + 1.26972675407126e-19⋅dt  + 3.90778437343018e-20⋅d\n",
      "\n",
      " 6                          7                          8                      \n",
      "t  + 4.80748915688373e-21⋅dt  + 5.17189453562972e-22⋅dt  + 7.20786458776279e-9\n",
      "\n",
      "   9    ⎛  10⎞\n",
      "⋅dt  + O⎝dt  ⎠\n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "t, dt = sp.symbols('t dt')\n",
    "# Set initial conditions\n",
    "t0 = 0\n",
    "y0 = 1\n",
    "# Set RHS of ODE\n",
    "toggle = 0 # This is a bookkeeping device for knowing when we reached the Inherently Adaptive methods\n",
    "for key,value in Butcher_dict.items():\n",
    "    if key in {\"AHE\", \"ABS\", \"ARKF\", \"ACK\", \"ADP5\", \"ADP8\"}:\n",
    "        if (key == \"AHE\"):\n",
    "            toggle = 1\n",
    "        if (toggle == 1): #only do anything once we actually hit adaptive methods.\n",
    "            print(\"RK method: \\\"\"+str(key)+\"\\\".\")\n",
    "            y = sp.Function('y')\n",
    "            print(\" When solving y'(t) = \"+str(rhs_dict[function](y(t),t))+\", y(\"+str(t0)+\")=\"+str(y0)+\",\")\n",
    "            local_truncation_order = list(value)[1]+1\n",
    "            print(\" the first calculation's first nonzero term should have local truncation error proportional to O(dt^\"+str(local_truncation_order-1)+\") or a higher power of dt.\")\n",
    "            print(\" the second calculation's first nonzero term should have local truncation error proportional to O(dt^\"+str(local_truncation_order)+\") or a higher power of dt.\")\n",
    "            print(\"Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\")\n",
    "            if key == \"ADP8\": #ADP8 takes up too much of the screen in analytic, we need to print as decimals.\n",
    "                sp.pretty_print(ValidateARK1(key, y0, t0, function).evalf())\n",
    "                sp.pretty_print(ValidateARK2(key, y0, t0, function).evalf())\n",
    "            else:\n",
    "                sp.pretty_print(ValidateARK1(key, y0, t0, function))\n",
    "                sp.pretty_print(ValidateARK2(key, y0, t0, function))\n",
    "        #     print(\"\\n\")\n",
    "            print(\" (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='abvalid'></a>\n",
    "\n",
    "## Step 2.e: Validating AB Methods against ODEs  [Back to [top](#toc)\\]\n",
    "$$\\label{abvalid}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code validates the AB methods, all of them from order 1 to 19. Their first remaining error term should be precisely one order higher than the method order.\n",
    "\n",
    "Since AB methods have to validate using multiple past values, the validation function is hard-coded in, since we need exact solutions at those past points for this to work. Fortunately $y'=y$ is suitable for this purpose, as it sufficiently tests all the orders up to 19. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is taken almost directly from the butcher validation from NRPy+\n",
    "\n",
    "def ValidateAB(ABorder, yn, tn, rhs_key): # custom function for validating AB methods.\n",
    "\n",
    "    order = ABorder\n",
    "\n",
    "    # 1. First we solve the ODE exactly\n",
    "    y = sp.Function('y')\n",
    "    sol = sp.dsolve(sp.Eq(y(t).diff(t), fy(y(t), t)), y(t)).rhs\n",
    "    constants = sp.solve([sol.subs(t,tn)-yn])\n",
    "    exact = sol.subs(constants)\n",
    "    exact_series = sp.series(exact.subs(t, dt),dt, 0, order+2)\n",
    "\n",
    "    # 2. Now we solve the ODE numerically using specified Butcher table\n",
    "\n",
    "    # Access the requested Butcher table\n",
    "    Butcher = Butcher_dict.get(\"AB\")[0]\n",
    "    # Determine number of predictor-corrector steps\n",
    "    L = len(Butcher)-1 #set to -2 for adaptive methods.\n",
    "    # Set a temporary array for update values\n",
    "    k = np.zeros(L, dtype=object)\n",
    "    # Initialize intermediate variable\n",
    "    yhat = 1\n",
    "    # Initialize the updated solution\n",
    "    ynp1 = 0\n",
    "    for i in range(ABorder):\n",
    "        # Initialize and approximate update for solution\n",
    "        # Update yhat for solution using \"past values\",\n",
    "        # Which means evaluating the exact answer at x=(-i*dt).\n",
    "        # If the user wishes to validate another ODE, the respective section\n",
    "        # In the below line will have to be changed.\n",
    "        yhat += dt*Butcher[ABorder-1][i]*sp.exp((-i)*dt)\n",
    "    # Finish determining the solution for the next iteration\n",
    "    num_series = sp.series(yhat, dt, 0, order+2)\n",
    "\n",
    "    # Produces Taylor series of exact solution at t=tn about t = 0 with the specified order\n",
    "    diff = exact_series-num_series\n",
    "    return diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AB method order: \"1\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^2) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "  2         \n",
      "dt     ⎛  3⎞\n",
      "─── + O⎝dt ⎠\n",
      " 2          \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"2\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^3) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "    3         \n",
      "5⋅dt     ⎛  4⎞\n",
      "───── + O⎝dt ⎠\n",
      "  12          \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"3\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^4) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "    4         \n",
      "3⋅dt     ⎛  5⎞\n",
      "───── + O⎝dt ⎠\n",
      "  8           \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"4\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^5) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "      5         \n",
      "251⋅dt     ⎛  6⎞\n",
      "─────── + O⎝dt ⎠\n",
      "  720           \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"5\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^6) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "     6         \n",
      "95⋅dt     ⎛  7⎞\n",
      "────── + O⎝dt ⎠\n",
      " 288           \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"6\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^7) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "        7         \n",
      "19087⋅dt     ⎛  8⎞\n",
      "───────── + O⎝dt ⎠\n",
      "  60480           \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"7\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^8) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "       8         \n",
      "5257⋅dt     ⎛  9⎞\n",
      "──────── + O⎝dt ⎠\n",
      " 17280           \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"8\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^9) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "          9          \n",
      "1070017⋅dt     ⎛  10⎞\n",
      "─────────── + O⎝dt  ⎠\n",
      "  3628800            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"9\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^10) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "        10          \n",
      "25713⋅dt      ⎛  11⎞\n",
      "────────── + O⎝dt  ⎠\n",
      "  89600             \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"10\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^11) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "           11          \n",
      "26842253⋅dt      ⎛  12⎞\n",
      "───────────── + O⎝dt  ⎠\n",
      "   95800320            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"11\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^12) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "          12          \n",
      "4777223⋅dt      ⎛  13⎞\n",
      "──────────── + O⎝dt  ⎠\n",
      "  17418240            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"12\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^13) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "               13          \n",
      "703604254357⋅dt      ⎛  14⎞\n",
      "───────────────── + O⎝dt  ⎠\n",
      "  2615348736000            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"13\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^14) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "               14          \n",
      "106364763817⋅dt      ⎛  15⎞\n",
      "───────────────── + O⎝dt  ⎠\n",
      "   402361344000            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"14\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^15) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "                15          \n",
      "1166309819657⋅dt      ⎛  16⎞\n",
      "────────────────── + O⎝dt  ⎠\n",
      "  4483454976000             \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"15\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^16) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "           16          \n",
      "25221445⋅dt      ⎛  17⎞\n",
      "───────────── + O⎝dt  ⎠\n",
      "   98402304            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"16\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^17) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "                   17          \n",
      "8092989203533249⋅dt      ⎛  18⎞\n",
      "───────────────────── + O⎝dt  ⎠\n",
      "  32011868528640000            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"17\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^18) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "                 18          \n",
      "85455477715379⋅dt      ⎛  19⎞\n",
      "─────────────────── + O⎝dt  ⎠\n",
      "  342372925440000            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"18\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^19) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n",
      "                       19          \n",
      "12600467236042756559⋅dt      ⎛  20⎞\n",
      "───────────────────────── + O⎝dt  ⎠\n",
      "   51090942171709440000            \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n",
      "AB method order: \"19\".\n",
      " When solving y'(t) = y, y(0)=1,\n",
      " the first nonzero term should have local truncation error proportional to O(dt^20) or a higher power of dt.\n",
      "Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                      20          \n",
      "1311546499957236437⋅dt      ⎛  21⎞\n",
      "──────────────────────── + O⎝dt  ⎠\n",
      "  5377993912811520000             \n",
      " (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "t, dt = sp.symbols('t dt')\n",
    "# Set initial conditions\n",
    "t0 = 0\n",
    "y0 = 1\n",
    "\n",
    "i = 1\n",
    "while i < 20:\n",
    "    print(\"AB method order: \\\"\"+str(i)+\"\\\".\")\n",
    "    y = sp.Function('y')\n",
    "    print(\" When solving y'(t) = y, y(0)=1,\") # make this adaptable to all possible inputs.\n",
    "    print(\" the first nonzero term should have local truncation error proportional to O(dt^\"+str(i+1)+\") or a higher power of dt.\")\n",
    "    print(\"Subtracting the numerical result from the exact Taylor expansion, we find a local truncation error of:\")\n",
    "    sp.pretty_print(ValidateAB(i, y0, t0, function))\n",
    "#     print(\"\\n\")\n",
    "    print(\" (Coefficients of order 1e-15 or less may generally be ignored, as these are at roundoff error.)\\n\")\n",
    "    i = i+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 3: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-RK_Butcher_Table_Validation.pdf](Tutorial-RK_Butcher_Table_Validation.pdf). (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:16:58.316815Z",
     "iopub.status.busy": "2021-03-07T17:16:58.316182Z",
     "iopub.status.idle": "2021-03-07T17:17:02.223642Z",
     "shell.execute_reply": "2021-03-07T17:17:02.222688Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-RK_Butcher_Table_Validation.tex, and compiled LaTeX file\n",
      "    to PDF file Tutorial-RK_Butcher_Table_Validation.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-RK_Butcher_Table_Validation\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}