Runge–Kutta Methods and Order

Let $T_n$ and $U_n$ be the real evaluation functions of the Chebyshev polynomials defined by Chebyshev polynomials of the first and second kinds by their three-term recurrences. For every natural $n$,

$\mathrm{\mid }{T}_n(x)\mathrm{\mid }\le 1,\mathrm{\mid }{U}_n(x)\mathrm{\mid }\le n+1(-1\le x\le 1),$

and for $n\ge 1$, $\mathrm{\mid }{T}_n(x)\mathrm{\mid }>1$ whenever $\mathrm{\mid }x\mathrm{\mid }>1$. Also,

$T_n(1)=1,T_n(-1)=(-1{)}^n,U_n(1)=n+1,U_n(-1)=(-1{)}^n(n+1)\mathrm{.}$

For every natural $n\ge 1$, $T_n$ has degree $n$ and leading coefficient $2^{n-1}$, while $U_n$ has degree $n$ and leading coefficient $2^n$. This degree and leading-coefficient assertion is exactly the content of Degrees and leading coefficients of the Chebyshev polynomials; that item does not assert an extremum bound.

The evaluation functions are twice differentiable at $1$, with

$T_n^{\mathrm{\prime }}(1)=n^2,T_n^{\mathrm{\prime }\mathrm{\prime }}(1)=\frac{n^2(n^2-1)}{3}\mathrm{.}$

Here each natural number in a real formula denotes its canonical real image (The canonical natural $\iota (n)=n\cdot 1_F$ of a field), and powers use Integer powers $a^m$.

Let $d\ge 1$, let $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ be an open time-state domain, and let $f:D\to {\mathbb{R}}^d$ be the right-hand side of a first-order initial-value problem in the sense of First-order systems, initial value problems, and solutions on intervals. An explicit $s$-stage Runge–Kutta tableau consists of a positive integer $s$, real coefficients $a_{ij}$, $b_i$, and $c_i$ with $a_{ij}=0$ whenever $j\ge i$, and $c_i={\sum }_{j<i}{a}_{ij}(1\le i\le s)\mathrm{.}$ The finite row sum is empty, hence $c_1=0$.

For a step size $h>0$ and an admitted base pair $(t,y)\in D$, the stage graph consists of the finite tuples satisfying $Y_i=y+h{\sum }_{j<i}{a}_{ij}{K}_j,t_i=t+c_{i}h,K_i=f(t_i,Y_i)(i=1,\dots ,s)\mathrm{.}$ These equations characterize each stage from the preceding slopes because $a_{ij}=0$ for $j\ge i$. For a complete stage tuple, the candidate endpoint is $y^{+}=y+h{\sum }_{i=1}^s{b}_i{K}_i\mathrm{.}$

To include unavailable evaluations, let $(E_i{)}_{i=1}^s$ be evaluator response functions on $D$: at each point $E_i$ has value either $f(\tau ,u)$ or the symbol evaluation-failure for stage $i$. Stage-dependent responses allow availability to differ between repeated inputs. A successful prefix of length $r$, $0\le r\le s$, is a stage tuple through index $r$ whose stage pairs lie in $D$ and whose responses under $E_i$ equal the corresponding slopes $K_i$. The empty prefix is allowed. For any prefix of length $r<s$, the next pair $(t_{r+1},Y_{r+1})$ is characterized by the same displayed equations using its already specified slopes. There is a unique maximal successful prefix, since every admissible successful extension has the uniquely specified slope $f(t_i,Y_i)$.

The status-labelled stage map is the following mathematical map. For a maximal prefix of length $r<s$, its value is domain-failure with index $r+1$ when the next pair is outside $D$, and evaluation-failure with index $r+1$ when that pair lies in $D$ and its response under $E_{r+1}$ is evaluation-failure. For a complete prefix its value is domain-failure at output when $(t+h,y^{+})\mathrm{\notin }D$, and success with $(t+h,y^{+})$ otherwise. Thus no slope at an inadmissible pair is part of the stage graph. In the exact mathematical model $E_i=f$ on $D$ for every $i$; evaluation-failure represents only an unavailable numerical evaluation. On its success domain the state component $y^{+}$ defines a one-step map in the sense of One-step ODE methods, local defect, uniform order, and step stability.

This algorithm takes an autonomous vector field $f$ on an open state domain $V\subseteq {\mathbb{R}}^d$, an admitted initial state $y$, a step $h>0$, an integer stage count $s\ge 1$, and an evaluator that returns either $f(u)$ or an evaluation-failure status. In the exact mathematical model the evaluator returns $f(u)$ at every $u\in V$; the failure status records only an unavailable numerical evaluation. Define stage abscissae $c_j=j^2\mathrm{/}{s}^2$ for $0\le j\le s$. Set $Y_0=y$. Before each query check that its stage state belongs to $V$; if it does not, return domain-failure with that stage index. If the evaluator reports failure, return evaluation-failure with that stage index. Query at $Y_0$ and, on success, set $Y_1=Y_0+(h\mathrm{/}{s}^2)f(Y_0)$. For $j=2,\dots ,s$, query at $Y_{j-1}$ and, on success, set $Y_j=2Y_{j-1}-Y_{j-2}+\frac{2h}{s^2}f(Y_{j-1})\mathrm{.}$ If the final state $Y_s$ is outside $V$, return output-domain-failure; otherwise return success with state $Y_s$.

For a nonautonomous field $f:D\to {\mathbb{R}}^d$ on an open time-state domain $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$, input an admitted pair $(t,y)\in D$ and an evaluator that returns either $f(\tau ,u)$ or evaluation-failure. In the exact mathematical model it returns $f$ at every point of $D$. Use the same state recurrence, querying at the declared stage times: query at $(t+c_{0}h,Y_0)=(t,y)$ before setting $Y_1=Y_0+(h\mathrm{/}{s}^2)f(t,Y_0)$; for $j=2,\dots ,s$, query at $(t+c_{j-1}h,Y_{j-1})$ before setting $Y_j=2Y_{j-1}-Y_{j-2}+\frac{2h}{s^2}f(t+c_{j-1}h,Y_{j-1})\mathrm{.}$ Before each query, check that its time-state input lies in $D$; return domain-failure with that stage index if it does not, or evaluation-failure with that index if the evaluator fails. After the stages, return output-domain-failure if $(t+h,Y_s)\mathrm{\notin }D$; otherwise return success with $(t+h,Y_s)$. The scalar linear stability polynomial is established separately.

The stability analysis uses $T_s$ and $U_s$ from Chebyshev polynomials of the first and second kinds by their three-term recurrences; this item defines the executable recurrence and does not itself assert an order or stability theorem.

Input a finite-dimensional real vector space $V$, a linear endomorphism $A:V\to V$, an initial vector $y\in V$, a step $h>0$, and an integer stage count $s\ge 2$. Each application of $A$ is the exact value of this total linear map, so this mathematical procedure has no application-failure branch. Define $w_1=\frac{3}{s^2-1},b_s=\frac{s^2-1}{3s^2},a_s=1-b_s\mathrm{.}$ Set $V_0=y$ and $V_1=V_0+w_{1}hA(V_0)$. For $j=2,\dots ,s$, compute $V_j=2V_{j-1}+2w_{1}hA(V_{j-1})-V_{j-2}\mathrm{.}$ Return $a_s{V}_0+b_s{V}_s$. Each stage uses one application of $A$; no implicit system is solved. The output is the one-step update for the autonomous linear equation $y^{\mathrm{\prime }}=Ay$. Here $T_s$ denotes the Chebyshev polynomial defined by Chebyshev polynomials of the first and second kinds by their three-term recurrences; the equivalent polynomial in $hA$ and its certified real stability interval are proved separately.

Input a first-order right-hand side $f:D\to {\mathbb{R}}^d$ on an open time-state domain $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$, an admitted pair $(t,y)\in D$, a step $h>0$, an integer stage count $s\ge 2$, and a right-hand-side evaluator that returns either $f(\tau ,u)$ or an evaluation-failure status at each queried $(\tau ,u)\in D$. In the exact mathematical model the evaluator returns $f$ at every point of $D$; the failure status records only an unavailable numerical evaluation. Set $w_1=\frac{3}{s^2-1},b_0=b_1=b_2=\frac{1}{4},b_j=\frac{j^2-1}{3j^2}(3\le j\le s),a_j=1-b_j(0\le j\le s)\mathrm{.}$ For the stages define $c_0=0$, $c_1=b_1{w}_1$, and $c_j=b_j{w}_1{j}^2$ for $2\le j\le s$. Put $Y_0=y$. Before each query check that its time-state input belongs to the declared domain; if not, return domain-failure with that stage index. If the right-hand-side evaluator reports failure, return evaluation-failure with that stage index. Query $k_0=f(t,Y_0)$ and on success set $Y_1=Y_0+b_1{w}_{1}h{k}_0$. For each $j=2,\dots ,s$, define ${\mu }_j=\frac{2b_j}{b_{j-1}},{\nu }_j=-\frac{b_j}{b_{j-2}},m_j=\frac{2b_j{w}_1}{b_{j-1}},g_j=-a_{j-1}{m}_j,$ query $k_{j-1}=f(t+c_{j-1}h,Y_{j-1})$, and on success set $Y_j=(1-{\mu }_j-{\nu }_j)y+{\mu }_j{Y}_{j-1}+{\nu }_j{Y}_{j-2}+m_{j}h{k}_{j-1}+g_{j}h{k}_0\mathrm{.}$ If $(t+h,Y_s)$ is outside the declared domain, return output-domain-failure; otherwise return success with $(t+h,Y_s)$. All stages use only the initial state and already computed values, so the procedure contains no implicit solve.

For an autonomous field, omit the time arguments. Order, stability, and stage-amplification properties are established separately.

The coefficient $b_2$ agrees with the displayed formula at $j=2$; it is listed among the base values to make the recurrence unambiguous. The Chebyshev identities used to derive these coefficients are supplied by Chebyshev polynomials of the first and second kinds by their three-term recurrences.

The following are the explicit Runge–Kutta tableaus in the convention of Explicit Runge–Kutta tableaus, stage maps, and failure status.

Forward Euler has $s=1$, $A=(0)$, $b=(1)$, and $c=(0)$.

Explicit midpoint has $s=2$, $A=\left(\begin{array}{cc}0& 0\\ 1\mathrm{/}2& 0\end{array}\right)$, $b=(0,1)$, and $c=(0,1\mathrm{/}2)$.

The classical four-stage Runge–Kutta method has $s=4$, $A=\left(\begin{array}{cccc}0& 0& 0& 0\\ 1\mathrm{/}2& 0& 0& 0\\ 0& 1\mathrm{/}2& 0& 0\\ 0& 0& 1& 0\end{array}\right)$, $b=(1\mathrm{/}6,1\mathrm{/}3,1\mathrm{/}3,1\mathrm{/}6)$, and $c=(0,1\mathrm{/}2,1\mathrm{/}2,1)$.

Each displayed $c_i$ is the row sum required in the definition of an explicit tableau. These tableaus specify the stage equations; their order and scalar stability functions are proved separately.

Let $(A,b,c)$ be an $s$-stage explicit tableau in the sense of Explicit Runge–Kutta tableaus, stage maps, and failure status, with $c_i={\sum }_{j<i}{a}_{ij}$. For a vector $v\in {\mathbb{R}}^s$, write $v_i^{[r]}=v_i^r$. Define $(Ac{)}_i={\sum }_{j<i}{a}_{ij}{c}_j$, $(A{c}^{[2]}{)}_i={\sum }_{j<i}{a}_{ij}{c}_j^2$, and $(A^{2}c{)}_i={\sum }_{j<i}{a}_{ij}(Ac{)}_j$.

Let $f:U\to {\mathbb{R}}^d$ be autonomous and of class $C^4$ on an open convex set $U\subseteq {\mathbb{R}}^d$. Assume an exact solution $y$ and all stages from exact starting values lie in a common compact convex $K\subset U$, and that the derivatives of $f$ through order four are bounded on $K$. For a component $f_i$, let $J=Df(y)$, $H=D^{2}f(y)$, and $G=D^{3}f(y)$ denote its coordinate derivative arrays: $(Jv{)}_i={\sum }_j{\mathrm{\partial }}_j{f}_i(y)v_j$, $H_i(u,v)={\sum }_{j,k}{\mathrm{\partial }}_j{\mathrm{\partial }}_k{f}_i(y)u_j{v}_k$, and $G_i(u,v,w)={\sum }_{j,k,\mathrm{\ell }}{\mathrm{\partial }}_j{\mathrm{\partial }}_k{\mathrm{\partial }}_{\mathrm{\ell }}{f}_i(y)u_j{v}_k{w}_{\mathrm{\ell }}$. These arrays act componentwise on vectors; $J^{2}v$ and $J^{3}v$ mean repeated matrix application, and $JH(u,v)$ means $J(H(u,v))$. The following coefficient equations are sufficient for the indicated uniform local orders on $K$:

• order at least one: ${\sum }_i{b}_i=1$;
• order at least two: also ${\sum }_i{b}_i{c}_i=1\mathrm{/}2$;
• order at least three: also ${\sum }_i{b}_i{c}_i^2=1\mathrm{/}3$ and ${\sum }_i{b}_i(Ac{)}_i=1\mathrm{/}6$;
• order at least four: also ${\sum }_i{b}_i{c}_i^3=1\mathrm{/}4$, ${\sum }_i{b}_i(A{c}^{[2]}{)}_i=1\mathrm{/}12$, ${\sum }_i{b}_i{c}_i(Ac{)}_i=1\mathrm{/}8$, and ${\sum }_i{b}_i(A^{2}c{)}_i=1\mathrm{/}24$.

Here local order $p$ means a uniform exact-start defect of norm $O(h^{p+1})$ as in One-step ODE methods, local defect, uniform order, and step stability. The conditions are sufficient coefficient matches; this item asserts no necessity claim.

For a nonautonomous field $f(t,y)$, adjoin the time equation $t^{\mathrm{\prime }}=1$ and apply the same autonomous argument to $(t,y{)}^{\mathrm{\prime }}=(1,f(t,y))$. Since $c_i={\sum }_{j<i}{a}_{ij}$, the time component of stage $i$ is exactly $t+c_{i}h$; thus this extension uses the actual tableau abscissae.

Input an initial pair $(t,y)$, a step $h>0$, an open time-state domain $D$, and a right-hand-side evaluator for $f:D\to {\mathbb{R}}^d$. Before each query, check that its time-state pair is in $D$; return domain-failure with that stage index if it is not. If the evaluator reports failure, return evaluation-failure with that stage index. Execute:

1. Query $f(t,y)$. On success set $Y_1=y+hf(t,y)$.
2. Query $f(t+h,Y_1)$. On success set $Y_2=\frac{3}{4}y+\frac{1}{4}(Y_1+hf(t+h,Y_1))$.
3. Query $f(t+h\mathrm{/}2,Y_2)$. On success set $y^{+}=\frac{1}{3}y+\frac{2}{3}(Y_2+hf(t+h\mathrm{/}2,Y_2))$.
4. If $(t+h,y^{+})\mathrm{\notin }D$, return output-domain-failure; otherwise return success with $(t+h,y^{+})$.

The procedure uses only explicit function evaluations; it has no implicit stage solve. If a numerical evaluator can fail even at a domain point, its failure is returned unchanged with the corresponding stage index. The exact mathematical method assumes all queried values exist.

Under the smoothness and compact-tube hypotheses of Sufficient explicit Runge–Kutta order conditions through order four, the forward Euler, explicit midpoint, and classical four-stage tableaus of The explicit Euler, midpoint, and classical fourth-order Runge–Kutta tableaus have uniform local orders at least one, two, and four, respectively.

For the scalar test equation $y^{\mathrm{\prime }}=\lambda y$, set $z=h\lambda$. The corresponding stability functions, defined by $y^{+}=R(z)y$, are $R_{\mathrm{E}}(z)=1+z,R_{\mathrm{M}}(z)=1+z+\frac{z^2}{2},R_{\mathrm{R}\mathrm{K}4}(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}\mathrm{.}$ The absolute stability region of a listed method is $\{z\in \mathbb{C}:\mathrm{\mid }R(z)\mathrm{\mid }\le 1\}$. On the negative real axis, the Euler and midpoint regions are exactly $[-2,0]$. No maximal real interval is asserted here for RK4.