Contents
RungeKutta methods¶
A RungeKutta method is a onestep method that computes the next time step solution as follows:
\begin{align*} y_i = & u^{n} + \Delta t \sum_{j=1}^{s} + a_{ij} f(y_j)) & (1\le j \le s) \\ u^{n+1} = & u^{n} + \Delta t \sum_{j=1}^{s} b_j f(y_j). \end{align*}
The simplest way to load a RungeKutta method is using the loadRKM function:
>> from nodepy import runge_kutta_method as rk
>> import numpy as np
>> rk44=rk.loadRKM('RK44')
Classical RK4
0.000 
0.500  0.500
0.500  0.000 0.500
1.000  0.000 0.000 1.000
_______________________________________
 0.167 0.333 0.333 0.167
Many wellknown methods are available through the loadRKM() function. Additionally, several classes of methods are available through the following functions:
See the documentation of these functions for more details.
More generally, any RungeKutta method may be instantiated by providing its Butcher coefficients, \(A\) and \(b\):
>> A=np.array([[0,0],[0.5,0]])
>> b=np.array([0,1.])
>> rk22=rk.RungeKuttaMethod(A,b)
Note that, because NumPy arrays are indexed from zero, the Butcher coefficient \(a_{21}\), for instance, corresponds to my_rk.a[1,0]. The abscissas \(c\) are automatically set to the row sums of \(A\) (this implies that every stage has has stage order at least equal to 1). Alternatively, a method may be specified in ShuOsher form, by coefficient arrays \(\alpha,\beta\):
>> rk22=rk.RungeKuttaMethod(alpha=alpha,beta=beta)
A separate subclass is provided for explicit RungeKutta methods: ExplicitRungeKuttaMethod. If a method is explicit, it is important to instantiate it as an ExplicitRungeKuttaMethod and not simply a RungeKuttaMethod, since the latter class has significantly less functionality in NodePy. Most significantly, time stepping is currently implemented for explicit methods, but not for implicit methods.
Examples:
>>> from nodepy.runge_kutta_method import *
Load a method:
>>> ssp104=loadRKM('SSP104')
Check its order of accuracy:
>>> ssp104.order() 4
Find its radius of absolute monotonicity:
>>> ssp104.absolute_monotonicity_radius() 5.999999999949068
Load a dictionary with many methods:
>>> RK=loadRKM() >>> sorted(RK.keys()) ['BE', 'BS5', 'BuRK65', 'CMR6', 'DP5', 'FE', 'Fehlberg45', 'GL2', 'GL3', 'Heun33', 'Lambert65', 'LobattoIIIA2', 'LobattoIIIA3', 'LobattoIIIC2', 'LobattoIIIC3', 'LobattoIIIC4', 'MTE22', 'Merson43', 'Mid22', 'NSSP32', 'NSSP33', 'PD8', 'RK44', 'RadauIIA2', 'RadauIIA3', 'SDIRK23', 'SDIRK34', 'SDIRK54', 'SSP104', 'SSP22', 'SSP22star', 'SSP33', 'SSP53', 'SSP54', 'SSP63', 'SSP75', 'SSP85', 'SSP95'] >>> print(RK['Mid22']) Midpoint RungeKutta 0  1/2  1/2 _______________  0 1
Many methods are naturally implemented in some ShuOsher form different from the Butcher form:
>>> ssp42 = SSPRK2(4) >>> ssp42.print_shu_osher() SSPRK(4,2)   1/3  1  1/3 2/3  1  1/3 1  1  1/3 _______________________________________________  1/4 3/4  1/4
 References:
Accuracy¶
The principal measure of accuracy of a RungeKutta method is its order of accuracy. By comparing the RungeKutta solution with the Taylor series for the exact solution, it can be shown that the local truncation error for small enough step size \(h\) is approximately
Error \(\approx Ch^p\),
where \(C\) is a constant independent of \(h\). Thus the expected asymptotic rate of convergence for small step sizes is \(p\). This error corresponds to the lowestorder terms that do not match those of the exact solution Taylor series. Typically, a higher order accurate method will provide greater accuracy than a lower order method even for practical step sizes.
In order to compare two methods with the same order of accuracy more detailed information about the accuracy of a method may be obtained by considering the relative size of the constant \(C\). This can be measured in various ways and is referred to as the principal error norm.
For example:
>> rk22.order()
2
>> rk44.order()
4
>> rk44.principal_error_norm()
0.014504582343198208
>> ssp104=rk.loadRKM('SSP104')
>> ssp104.principal_error_norm()
0.002211223747053554
Since the SSP(10,4) method has smaller principal error norm, we expect that it will provide better accuracy than the classical 4stage RungeKutta method for a given step size. Of course, the SSP method has 10 stages, so it requires more work per step. In order to determine which method is more efficient, we need to compare the relative accuracy for a fixed amount of work:
>> rk.relative_accuracy_efficiency(rk44,ssp104)
1.7161905294239843
This indicates that, for a desired level of error, the SSP(10,4) method will require about 72% more work.

RungeKuttaMethod.
principal_error_norm
(tol=1e13, mode='float')[source] The 2norm of the vector of leading order error coefficients.

RungeKuttaMethod.
error_metrics
()[source] Returns several measures of the accuracy of the RungeKutta method. In order, they are:
 \(A^{q+1}\): 2norm of the vector of leading order error coefficients
 \(A^{q+1}_{max}\): Maxnorm of the vector of leading order error coefficients
 \(A^{q+2}\) : 2norm of the vector of next order error coefficients
 \(A^{q+2}_{max}\): Maxnorm of the vector of next order error coefficients
 \(D\): The largest (in magnitude) coefficient in the Butcher array
Examples:
>>> from nodepy import rk >>> rk4 = rk.loadRKM('RK44') >>> rk4.error_metrics() (sqrt(1745)/2880, 1/120, sqrt(8531)/5760, 1/144, 1)
Reference: [kennedy2000]

RungeKuttaMethod.
stage_order
(tol=1e14)[source] The stage order of a RungeKutta method is the minimum, over all stages, of the order of accuracy of that stage. It can be shown to be equal to the largest integer k such that the simplifying assumptions \(B(\\xi)\) and \(C(\\xi)\) are satisfied for \(1 \\le \\xi \\le k\).
Examples:
>>> from nodepy import rk >>> rk4 = rk.loadRKM('RK44') >>> rk4.stage_order() 1 >>> gl2 = rk.loadRKM('GL2') >>> gl2.stage_order() 2
 References:
 Dekker and Verwer
 [butcher2003]
Classical (linear) stability¶

RungeKuttaMethod.
stability_function
(stage=None, mode='exact', formula='lts', use_butcher=False)[source] The stability function of a RungeKutta method is \(\\phi(z)=p(z)/q(z)\), where
$$p(z)=\det(I  z A + z e b^T)$$
$$q(z)=\det(I  z A)$$
The function can also be computed via the formula
$$`\phi(z) = 1 + b^T (IzA)^{1} e$$
where \(e\) is a column vector with all entries equal to one.
This function constructs the numerator and denominator of the stability function of a RungeKutta method.
For methods with rational coefficients, mode=’exact’ computes the stability function using rational arithmetic. Alternatively, you can set mode=’float’ to force computation using floating point, in case the exact computation is too slow.
For explicit methods, the denominator is simply \(1\) and there are three options for computing the numerator (this is the ‘formula’ option). These only affect the speed, and only matter if the computation is symbolic. They are:
 ‘lts’: SymPy’s lower_triangular_solve
 ‘det’: ratio of determinants
 ‘pow’: power series
For implicit methods, only the ‘det’ (determinant) formula is supported. If mode=’float’ is selected, the formula automatically switches to ‘det’.
The user can also select whether to compute the function based on Butcher or ShuOsher coefficients by setting \(use_butcher\).
 Output:
 p – Numpy poly representing the numerator
 q – Numpy poly representing the denominator
Examples:
>>> from nodepy import rk >>> rk4 = rk.loadRKM('RK44') >>> p,q = rk4.stability_function() >>> print(p) 4 3 2 0.04167 x + 0.1667 x + 0.5 x + 1 x + 1 >>> dc = rk.DC(3) >>> dc.stability_function(mode='exact') (poly1d([1/3888, 1/648, 1/24, 1/6, 1/2, 1, 1], dtype=object), poly1d([1], dtype=object)) >>> dc.stability_function(mode='float') (poly1d([ 2.57201646e04, 1.54320988e03, 4.16666667e02, 1.66666667e01, 5.00000000e01, 1.00000000e+00, 1.00000000e+00]), poly1d([ 1.])) >>> ssp3 = rk.SSPIRK3(4) >>> ssp3.stability_function() (poly1d([67/300 + 13*sqrt(15)/225, sqrt(15)/25 + 1/6, sqrt(15)/5 + 9/10, 1 + 2*sqrt(15)/5, 1], dtype=object), poly1d([2*sqrt(15)/25 + 31/100, 7/5 + 9*sqrt(15)/25, 3*sqrt(15)/5 + 12/5, 2 + 2*sqrt(15)/5, 1], dtype=object)) >>> ssp3.stability_function(mode='float') (poly1d([ 4.39037781e04, 1.17473328e02, 1.25403331e01, 5.49193338e01, 1.00000000e+00]), poly1d([ 1.61332303e04, 5.72599537e03, 7.62099923e02, 4.50806662e01, 1.00000000e+00])) >>> ssp2 = rk.SSPIRK2(1) >>> ssp2.stability_function() (poly1d([1/2, 1], dtype=object), poly1d([1/2, 1], dtype=object))

RungeKuttaMethod.
plot_stability_region
(N=200, color='r', filled=True, bounds=None, plotroots=False, alpha=1.0, scalefac=1.0, to_file=False, longtitle=True, fignum=None)[source] The region of absolute stability of a RungeKutta method, is the set
\(\{ z \in C : \phi (z)\le 1 \}\)
where \(\phi(z)\) is the stability function of the method.
 Input: (all optional)
 N – Number of gridpoints to use in each direction
 bounds – limits of plotting region
 color – color to use for this plot
 filled – if true, stability region is filled in (solid); otherwise it is outlined
 Example::
>>> from nodepy import rk >>> rk4 = rk.loadRKM('RK44') >>> rk4.plot_stability_region() <matplotlib.figure.Figure object at 0x...>
Nonlinear stability¶

RungeKuttaMethod.
absolute_monotonicity_radius
(acc=1e10, rmax=200, tol=3e16)[source] Returns the radius of absolute monotonicity (also referred to as the radius of contractivity or the strong stability preserving coefficient of a RungeKutta method.

RungeKuttaMethod.
circle_contractivity_radius
(acc=1e13, rmax=1000)[source] Returns the radius of circle contractivity of a RungeKutta method.
Example:
>>> from nodepy import rk >>> rk4 = rk.loadRKM('RK44') >>> rk4.circle_contractivity_radius() 1.000...
Reducibility of RungeKutta methods¶
Two kinds of reducibility (DJreducibility and HSreducibility) have been identified in the literature. NodePy contains functions for detecting both and transforming a reducible method to an equivalent irreducible method. Of course, reducibility is dealt with relative to some numerical tolerance, since the method coefficients are floating point numbers.

RungeKuttaMethod.
_dj_reducible_stages
(tol=1e13)[source] Determine whether the method is DJreducible.
A method is DJreducible if it contains any stage that does not influence the output.
Returns a list of unnecessary stages. If the method is DJirreducible, returns an empty list.
This routine may not work correctly for RK pairs, as it doesn’t check bhat.

RungeKuttaMethod.
dj_reduce
(tol=1e13)[source] Remove all DJreducible stages.
A method is DJreducible if it contains any stage that does not influence the output.
Examples:
Construct a reducible method: >>> from nodepy import rk >>> A=np.array([[0,0],[1,0]]) >>> b=np.array([1,0]) >>> rkm = rk.ExplicitRungeKuttaMethod(A,b) Check that it is reducible: >>> rkm._dj_reducible_stages() [1] Reduce it: >>> print(rkm.dj_reduce()) RungeKutta Method <BLANKLINE> 0  ______  1

RungeKuttaMethod.
_hs_reducible_stages
(tol=1e13)[source] Determine whether the method is HSreducible. A RungeKutta method is HSreducible if two rows of A are equal.
If the method is HSreducible, returns True and a pair of equal stages. If not, returns False and the minimum pairwise difference (in the maximum norm) between rows of A.
Examples:
Construct a reducible method: >>> from nodepy import rk >>> A=np.array([[1,0],[1,0]]) >>> b=np.array([0.5,0.5]) >>> rkm = rk.ExplicitRungeKuttaMethod(A,b) Check that it is reducible: >>> rkm._hs_reducible_stages() (True, [0, 1])
Composing RungeKutta methods¶
Butcher has developed an elegant theory of the group structure of RungeKutta methods. The RungeKutta methods form a group under the operation of composition. The multiplication operator has been overloaded so that multiplying two RungeKutta methods gives the method corresponding to their composition, with equal timesteps.
It is also possible to compose methods with nonequal timesteps using the compose() function.
Embedded RungeKutta Pairs¶

class
nodepy.runge_kutta_method.
ExplicitRungeKuttaPair
(A=None, b=None, bhat=None, alpha=None, beta=None, alphahat=None, betahat=None, name='RungeKutta Pair', shortname='RKM', description='', order=(None, None))[source] Class for embedded RungeKutta pairs. These consist of two methods with identical coefficients \(a_{ij}\) but different coefficients \(b_j\) such that the methods have different orders of accuracy. Typically the higher order accurate method is used to advance the solution, while the lower order method is used to obtain an error estimate.
An embedded RungeKutta Pair takes the form:
\begin{align*} y_i = & u^{n} + \Delta t \sum_{j=1}^{s} + a_{ij} f(y_j)) & (1\le j \le s) \\ u^{n+1} = & u^{n} + \Delta t \sum_{j=1}^{s} b_j f(y_j) \\ \hat{u}^{n+1} = & u^{n} + \Delta t \sum_{j=1}^{s} \hat{b}_j f(y_j). \end{align*}
That is, both methods use the same intermediate stages \(y_i\), but different weights. Typically the weights \(\\hat{b}_j\) are chosen so that \(\\hat{u}^{n+1}\) is accurate of order one less than the order of \(u^{n+1}\). Then their difference can be used as an error estimate.
The class also admits ShuOsher representations:
\begin{align*} y_i = & v_i u^{n} + \sum_{j=1}^s \alpha_{ij} y_j + \Delta t \sum_{j=1}^{s} + \beta_{ij} f(y_j)) & (1\le j \le s+1) \\ u^{n+1} = & y_{s+1} \hat{u}^{n+1} = & \hat{v}_{s+1} u^{n} + \sum_{j=1}^s \hat{\alpha}_{s+1,j} + \Delta t \sum_{j=1}^{s} \hat{\beta}_{s+1,j} f(y_j). \end{align*}
In NodePy, if rkp is a RungeKutta pair, the principal (usually higherorder) method is the one used if accuracy or stability properties are queried. Properties of the embedded (usually lowerorder) method can be accessed via rkp.embedded_method.
When solving an IVP with an embedded pair, one can specify a desired error tolerance. The step size will be adjusted automatically to achieve approximately this tolerance.
In addition to the ordinary RungeKutta initialization, here the embedded coefficients \(\hat{b}_j\) are set as well.
LowStorage RungeKutta methods¶
Typically, implementation of a RungeKutta method requires \(s \times N\) memory locations, where \(s\) is the number of stages of the method and \(N\) is the number of unknowns. Certain classes of RungeKutta methods can be implemented using substantially less memory, by taking advantage of special relations among the coefficients. Three main classes have been developed in the literature:
 2N (Williamson) methods
 2R (van der Houwen/Wray) methods
 2S methods
Each of these classes requires only \(2\times N\) memory locations. Additional methods have been developed that use more than two memory locations per unknown but still provide a substantial savings over traditional methods. These are referred to as, e.g., 3S, 3R, 4R, and so forth. For a review of lowstorage methods, see [ketcheson2010] .
In NodePy, lowstorage methods are a subclass of explicit RungeKutta methods (and/or explicit RungeKutta pairs). In addition to the usual properties, they possess arrays of lowstorage coefficients. They override the generic RK implementation of timestepping and use special memoryefficient implementations instead. It should be noted that, although these lowstorage algorithms are implemented, due to Python language restrictions an extra intermediate copy of the solution array will be created. Thus the implementation in NodePy is not really minimumstorage.
At the moment, the following classes are implemented:
2S : Methods using two registers (under Ketcheson’s assumption)
 2S* : Methods using two registers, one of which retains the previous
step solution
 3S* : Methods using three registers, one of which retains the previous
step solution
2S embedded pairs
3S* embedded pairs
2R embedded pairs
3R embedded pairs
Examples:
>>> from nodepy import lsrk, ivp
>>> myrk = lsrk.load_2R('DDAS47')
>>> print(myrk)
DDAS4()7[2R]
2R Method of Tselios \& Simos (2007)
0.000 
0.336  0.336
0.286  0.094 0.192
0.745  0.094 0.150 0.501
0.639  0.094 0.150 0.285 0.110
0.724  0.094 0.150 0.285 0.122 0.317
0.911  0.094 0.150 0.285 0.122 0.061 0.444
________________________________________________________
 0.094 0.150 0.285 0.122 0.061 0.346 0.187
>>> rk58 = lsrk.load_2R('RK58[3R]C')
>>> rk58.name
'RK5(4)8[3R+]C'
>>> rk58.order()
5
>>> problem = ivp.load_ivp('vdp')
>>> t,u = rk58(problem)
>>> u[1]
array([1.40278844, 1.23080499])
>>> import nodepy
>>> rk2S = lsrk.load_LSRK("{}/method_coefficients/582S_acc.txt".format(nodepy.__path__[0]),has_emb=True)
>>> rk2S.order()
5
>>> rk2S.embedded_method.order()
4
>>> rk3S = lsrk.load_LSRK(nodepy.__path__[0]+'/method_coefficients/583Sstar_acc.txt',lstype='3S*')
>>> rk3S.principal_error_norm()
0.00035742076...
2S/3S methods¶

class
nodepy.low_storage_rk.
TwoSRungeKuttaMethod
(betavec, gamma, delta, lstype, name='Lowstorage RungeKutta Method', description='', shortname='LSRK2S', order=None)[source] Class for lowstorage RungeKutta methods that use Ketcheson’s assumption (2S, 2S*, and 3S* methods).
This class cannot be used for embedded pairs. Use the class TwoSRungeKuttaPair instead.
The lowstorage coefficient arrays \(eta,\gamma,\delta\) follow the notation of [ketcheson2010] .
The argument lstype must be one of the following values:
 2S
 2S*
 3S*
Initializes the lowstorage method by storing the lowstorage coefficients and computing the Butcher coefficients.
2S/3S embedded pairs¶

class
nodepy.low_storage_rk.
TwoSRungeKuttaPair
(betavec, gamma, delta, lstype, bhat=None, name='Lowstorage RungeKutta Pair', description='', shortname='LSRK2S', order=None)[source] Class for lowstorage embedded RungeKutta pairs that use Ketcheson’s assumption (2S, 2S*, and 3S* methods).
This class is only for embedded pairs. Use the class TwoSRungeKuttaMethod for single 2S/3S methods.
The lowstorage coefficient arrays \(eta,\gamma,\delta\) follow the notation of [ketcheson2010] .
The argument lstype must be one of the following values:
 2S
 2S*
 2S_pair
 3S*
 3S*_pair
The 2S/2S*/3S* classes do not need an extra register for the error estimate, while the 2S_pair/3S_pair methods do.
Initializes the lowstorage pair by storing the lowstorage coefficients and computing the Butcher coefficients.
2R/3R methods¶

class
nodepy.low_storage_rk.
TwoRRungeKuttaMethod
(a, b, bhat=None, regs=2, name='2R RungeKutta Method', description='', shortname='LSRK2R', order=None)[source] Class for 2R/3R/4R lowstorage RungeKutta pairs.
These were developed by van der Houwen, Wray, and Kennedy et. al. Only 2R and 3R methods have been implemented so far.
 References:
Initializes the 2R method by storing the lowstorage coefficients and computing the Butcher array.
The coefficients should be specified as follows:
 For all methods, the weights \(b\) are used to fill in the appropriate entries in \(A\).
 For 2R methods, a is a vector of length \(s1\) containing the first subdiagonal of \(A\)
 For 3R methods, a is a \(2\times s1\) array whose first row contains the first subdiagonal of \(A\) and whose second row contains the second subdiagonal of \(A\).