Show code cell content

Hide code cell content

import mmf_setup;mmf_setup.nbinit()
import logging;logging.getLogger('matplotlib').setLevel(logging.CRITICAL)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt

$\newcommand{\vect}[1]{\mathbf{#1}} \newcommand{\uvect}[1]{\hat{#1}} \newcommand{\abs}[1]{\lvert#1\rvert} \newcommand{\norm}[1]{\lVert#1\rVert} \newcommand{\I}{\mathrm{i}} \newcommand{\ket}[1]{\left|#1\right\rangle} \newcommand{\bra}[1]{\left\langle#1\right|} \newcommand{\braket}[1]{\langle#1\rangle} \newcommand{\Braket}[1]{\left\langle#1\right\rangle} \newcommand{\op}[1]{\mathbf{#1}} \newcommand{\mat}[1]{\mathbf{#1}} \newcommand{\d}{\mathrm{d}} \newcommand{\D}[1]{\mathcal{D}[#1]\;} \newcommand{\pdiff}[3][]{\frac{\partial^{#1} #2}{\partial {#3}^{#1}}} \newcommand{\diff}[3][]{\frac{\d^{#1} #2}{\d {#3}^{#1}}} \newcommand{\ddiff}[3][]{\frac{\delta^{#1} #2}{\delta {#3}^{#1}}} \newcommand{\floor}[1]{\left\lfloor#1\right\rfloor} \newcommand{\ceil}[1]{\left\lceil#1\right\rceil} % The following allows KaTeX to render these for significantly improved % performance on CoCalc. \newcommand{\DeclareMathOperator}[2]{\newcommand{#1}{#2}} \DeclareMathOperator{\Tr}{Tr} \DeclareMathOperator{\erf}{erf} \DeclareMathOperator{\erfi}{erfi} \DeclareMathOperator{\sech}{sech} \DeclareMathOperator{\sn}{sn} \DeclareMathOperator{\cn}{cn} \DeclareMathOperator{\dn}{dn} \DeclareMathOperator{\sgn}{sgn} \DeclareMathOperator{\order}{O} \DeclareMathOperator{\diag}{diag} % \newcommand{\mylabel}[1]{\label{#1}\tag{#1}} \newcommand{\degree}{\circ}$

This cell adds /home/docs/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:

• Choose "Trust Notebook" from the "File" menu.
• Re-execute this cell.
• Reload the notebook.

Runge-Kutta Truncation Error

Consider the Runge-Kutta RK4 formulation given in [Gezerlis, 2023] for evolving \(y'(t) = f(t, y)\):

\[\begin{align*} k_0 &= h f(t, y)\\ k_1 &= h f\Bigl(t+\tfrac{1}{2}h, y+\tfrac{1}{2}k_0\Bigr)\\ k_2 &= h f\Bigl(t+\tfrac{1}{2}h, y+\tfrac{1}{2}k_1\Bigr) \tag{8.64}\\ k_3 &= h f\Bigl(t+h, y+k_2\Bigr)\\ y(t+h) &\approx y + \frac{k_0 + 2k_1 + 2k_2 + k_3}{6}. \end{align*}\]

Tracing through, one can expand the final step \(y_{j+1} \approx y(t+h)\) to obtaining an expression in terms of \(y'(t) = f(t, y)\) and various derivatives. This expression can be simplified by noting:

\[\begin{align*} y'(t) &= f\\ y''(t) &= f_{,t} + f_{,y}y',\\ y'''(t) &= f_{,tt} + 2f_{,ty}y' + f_{,yy}(y')^2 + f_{,y}y'',\\ y''''(t) &= f_{,ttt} + 3f_{,tty}y' + 3f_{,tyy}(y')^2 + f_{,yyy}(y')^3 + 2f_{,ty}y'' + 3f_{,yy}y'y'' + f_{,ty} y'' + f_{,y}y''' \end{align*}\]

(Incomplete. Probably Mathematica or Maple can do a better job.)

In principle, the following code should be able to give the correct error formula, but SymPy has problems simplifying the various expressions.

from sympy import *
t, h = var('t, h', real=True)
f = Function('f', real=True)
y = Function('y', real=True)
dy = f(t, y(t))
d2y = dy.diff(t)
d3y = d2y.diff(t)
d4y = d3y.diff(t)
k0 = h*f(t, y(t))
k1 = h*f(t+h/2,  y(t)+k0/2)
k2 = h*f(t+h/2,  y(t)+k1/2)
k3 = h*f(t+h,  y(t)+k2)
y1 = y(t) + (k0 + 2*k1 + 2*k2 + k3)/6
err = ((y(t+h) - y1).series(h, 0, 5).simplify()
                    .replace(y(t).diff(t,t,t,t), d4y).simplify()
                    .replace(y(t).diff(t,t,t), d3y).simplify()
                    .replace(y(t).diff(t,t), d2y).simplify()
                    .replace(diff(y(t), t), dy).simplify())
display(err)
err.coeff(h, 3).doit().simplify()

\[\displaystyle \frac{h^{3} f{\left(t,y{\left(t \right)} \right)} \left. \frac{\partial^{2}}{\partial t\partial \xi} f{\left(t,\xi \right)} \right|_{\substack{ \xi=y{\left(t \right)} }}}{6} + \frac{h^{3} f{\left(t,y{\left(t \right)} \right)} \left. \frac{\partial^{2}}{\partial y{\left(t \right)}\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }}}{12} - \frac{h^{3} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{2}}{\partial t\partial _xi*f(t, y(t))/2 + y(t)} f{\left(t,_xi*f(t, y(t))/2 + y(t) \right)}}{24} - \frac{h^{3} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{2}}{\partial t\partial _xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t)} f{\left(t,_xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t) \right)}}{12} - \frac{h^{3} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{2}}{\partial y{\left(t \right)}\partial _t + xi/2} f{\left(_t + xi/2,y{\left(t \right)} \right)}}{24} - \frac{h^{3} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{2}}{\partial y{\left(t \right)}\partial _t + xi} f{\left(_t + xi,y{\left(t \right)} \right)}}{12} + \frac{h^{4} \left. \frac{\partial}{\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }} \left. \frac{\partial^{2}}{\partial t\partial \xi} f{\left(t,\xi \right)} \right|_{\substack{ \xi=y{\left(t \right)} }}}{12} - \frac{h^{4} \left. \frac{\partial}{\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }} \left. \frac{\partial^{2}}{\partial y{\left(t \right)}\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }}}{72} - \frac{h^{4} \frac{\partial^{2}}{\partial t\partial _xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t)} f{\left(t,_xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t) \right)} \left. \frac{\partial}{\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }}}{24} - \frac{h^{4} \frac{\partial^{2}}{\partial y{\left(t \right)}\partial _t + xi} f{\left(_t + xi,y{\left(t \right)} \right)} \left. \frac{\partial}{\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }}}{36} + \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \left. \frac{\partial^{3}}{\partial t^{2}\partial \xi_{2}} f{\left(t,\xi_{2} \right)} \right|_{\substack{ \xi_{2}=y{\left(t \right)} }}}{72} + \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \left. \frac{\partial^{3}}{\partial y{\left(t \right)}\partial \xi_{1}^{2}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }}}{12} - \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{3}}{\partial t^{2}\partial _xi*f(t, y(t))/2 + y(t)} f{\left(t,_xi*f(t, y(t))/2 + y(t) \right)}}{144} - \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{3}}{\partial t^{2}\partial _xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t)} f{\left(t,_xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t) \right)}}{48} - \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{3}}{\partial y{\left(t \right)}\partial _t + xi/2^{2}} f{\left(_t + xi/2,y{\left(t \right)} \right)}}{144} - \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{3}}{\partial y{\left(t \right)}\partial _t + xi^{2}} f{\left(_t + xi,y{\left(t \right)} \right)}}{36} + \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{d}{d y{\left(t \right)}} f{\left(t,y{\left(t \right)} \right)} \left. \frac{\partial^{2}}{\partial t\partial \xi} f{\left(t,\xi \right)} \right|_{\substack{ \xi=y{\left(t \right)} }}}{8} + \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{d}{d y{\left(t \right)}} f{\left(t,y{\left(t \right)} \right)} \left. \frac{\partial^{2}}{\partial y{\left(t \right)}\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }}}{36} - \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{d}{d y{\left(t \right)}} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{2}}{\partial t\partial _xi*f(t, y(t))/2 + y(t)} f{\left(t,_xi*f(t, y(t))/2 + y(t) \right)}}{48} - \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{d}{d y{\left(t \right)}} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{2}}{\partial t\partial _xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t)} f{\left(t,_xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t) \right)}}{12} - \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{d}{d y{\left(t \right)}} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{2}}{\partial y{\left(t \right)}\partial _t + xi/2} f{\left(_t + xi/2,y{\left(t \right)} \right)}}{48} - \frac{h^{4} f{\left(t,y{\left(t \right)} \right)} \frac{d}{d y{\left(t \right)}} f{\left(t,y{\left(t \right)} \right)} \frac{\partial^{2}}{\partial y{\left(t \right)}\partial _t + xi} f{\left(_t + xi,y{\left(t \right)} \right)}}{36} + \frac{h^{4} f^{2}{\left(t,y{\left(t \right)} \right)} \left. \frac{\partial^{3}}{\partial t\partial \xi^{2}} f{\left(t,\xi \right)} \right|_{\substack{ \xi=y{\left(t \right)} }}}{12} - \frac{h^{4} f^{2}{\left(t,y{\left(t \right)} \right)} \left. \frac{\partial^{3}}{\partial y{\left(t \right)}^{2}\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }}}{72} - \frac{h^{4} f^{2}{\left(t,y{\left(t \right)} \right)} \frac{\partial^{3}}{\partial t\partial _xi*f(t, y(t))/2 + y(t)^{2}} f{\left(t,_xi*f(t, y(t))/2 + y(t) \right)}}{144} - \frac{h^{4} f^{2}{\left(t,y{\left(t \right)} \right)} \frac{\partial^{3}}{\partial t\partial _xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t)^{2}} f{\left(t,_xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t) \right)}}{72} - \frac{h^{4} f^{2}{\left(t,y{\left(t \right)} \right)} \frac{\partial^{3}}{\partial y{\left(t \right)}^{2}\partial _t + xi/2} f{\left(_t + xi/2,y{\left(t \right)} \right)}}{144} - \frac{h^{4} f^{2}{\left(t,y{\left(t \right)} \right)} \frac{\partial^{3}}{\partial y{\left(t \right)}^{2}\partial _t + xi} f{\left(_t + xi,y{\left(t \right)} \right)}}{36} + O\left(h^{5}\right)\]

\[\displaystyle \frac{\left(- 2 \frac{\partial^{2}}{\partial y{\left(t \right)}\partial _t + xi} f{\left(_t + xi,y{\left(t \right)} \right)} - \frac{\partial^{2}}{\partial y{\left(t \right)}\partial _t + xi/2} f{\left(_t + xi/2,y{\left(t \right)} \right)} - 2 \frac{\partial^{2}}{\partial t\partial _xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t)} f{\left(t,_xi*f(t + xi/2, xi*f(t, y(t))/2 + y(t))/2 + y(t) \right)} - \frac{\partial^{2}}{\partial t\partial _xi*f(t, y(t))/2 + y(t)} f{\left(t,_xi*f(t, y(t))/2 + y(t) \right)} + 2 \left. \frac{\partial^{2}}{\partial y{\left(t \right)}\partial \xi_{1}} f{\left(\xi_{1},y{\left(t \right)} \right)} \right|_{\substack{ \xi_{1}=t }} + 4 \left. \frac{\partial^{2}}{\partial t\partial \xi} f{\left(t,\xi \right)} \right|_{\substack{ \xi=y{\left(t \right)} }}\right) f{\left(t,y{\left(t \right)} \right)}}{24}\]