Contents

Runge-Kutta methods

A Runge-Kutta method is a one-step 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 Runge-Kutta 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 well-known methods are available through the loadRKM() function. Additionally, several classes of methods are available through the following functions:

  • Optimal strong stability preserving methods: SSPRK2(s), SSPRK3(s), SSPIRK2(), etc.

  • Integral deferred correction methods: DC(s)

  • Extrapolation methods: extrap(s)

  • Runge-Kutta Chebyshev methods: RKC1(s), RKC2(s)

See the documentation of these functions for more details.

More generally, any Runge-Kutta 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 Shu-Osher form, by coefficient arrays \(\alpha, \beta\):

>> rk22=rk.RungeKuttaMethod(alpha=alpha,beta=beta)

A separate subclass is provided for explicit Runge-Kutta 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 *
References:

Accuracy

The principal measure of accuracy of a Runge-Kutta method is its order of accuracy. By comparing the Runge-Kutta 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 lowest-order 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 4-stage Runge-Kutta 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=1e-13, mode='float')[source]

The 2-norm of the vector of leading order error coefficients.

RungeKuttaMethod.error_metrics(q=None, tol=1e-14)[source]

Returns several measures of the accuracy of the Runge-Kutta method. In order, they are:

  • \(A^{q+1}\): 2-norm of the vector of leading order error coefficients

  • \(A^{q+1}_{max}\): Max-norm of the vector of leading order error coefficients

  • \(A^{q+2}\) : 2-norm of the vector of next order error coefficients

  • \(A^{q+2}_{max}\): Max-norm 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()
main method has order 4
(sqrt(1745)/2880, 1/120, sqrt(8531)/5760, 1/144, 1)

Reference: [KCL00]

RungeKuttaMethod.stage_order(tol=1e-14)[source]

The stage order of a Runge-Kutta 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:

Classical (linear) stability

RungeKuttaMethod.stability_function(stage=None, mode='exact', formula='lts', use_butcher=False)[source]

The stability function of a Runge-Kutta 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 (I-zA)^{-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 Runge-Kutta 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 Shu-Osher 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.57201646e-04, 1.54320988e-03, 4.16666667e-02, 1.66666667e-01,
       5.00000000e-01, 1.00000000e+00, 1.00000000e+00]), poly1d([1.]))
>>> ssp3 = rk.SSPIRK3(4)
>>> ssp3.stability_function()
(poly1d([-67/300 + 13*sqrt(15)/225, 1/6 - sqrt(15)/25, 9/10 - sqrt(15)/5,
       -1 + 2*sqrt(15)/5, 1], dtype=object), poly1d([31/100 - 2*sqrt(15)/25, -7/5 + 9*sqrt(15)/25, 12/5 - 3*sqrt(15)/5,
       -2 + 2*sqrt(15)/5, 1], dtype=object))

>>> ssp3.stability_function(mode='float')  
(poly1d([4.39037781e-04, 1.17473328e-02, 1.25403331e-01, 5.49193338e-01,
       1.00000000e+00]), poly1d([ 1.61332303e-04, -5.72599537e-03, 7.62099923e-02, -4.50806662e-01,
       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 Runge-Kutta 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() 
<Figure size...

Nonlinear stability

RungeKuttaMethod.absolute_monotonicity_radius(acc=1e-10, rmax=200, tol=3e-16)[source]

Returns the radius of absolute monotonicity (also referred to as the radius of contractivity or the strong stability preserving coefficient of a Runge-Kutta method.

RungeKuttaMethod.circle_contractivity_radius(acc=1e-13, rmax=1000)[source]

Returns the radius of circle contractivity of a Runge-Kutta method.

Example:

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> rk4.circle_contractivity_radius() 
1.000...

Reducibility of Runge-Kutta methods

Two kinds of reducibility (DJ-reducibility and HS-reducibility) 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=1e-13)[source]

Determine whether the method is DJ-reducible.

A method is DJ-reducible if it contains any stage that does not influence the output.

Returns a list of unnecessary stages. If the method is DJ-irreducible, returns an empty list.

This routine may not work correctly for RK pairs, as it doesn’t check bhat.

RungeKuttaMethod.dj_reduce(tol=1e-13)[source]

Remove all DJ-reducible stages.

A method is DJ-reducible 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())
Runge-Kutta Method
<BLANKLINE>
 0 |
___|___
   | 1
RungeKuttaMethod._hs_reducible_stages(tol=1e-13)[source]

Determine whether the method is HS-reducible. A Runge-Kutta method is HS-reducible if two rows of A are equal.

If the method is HS-reducible, 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 Runge-Kutta methods

Butcher has developed an elegant theory of the group structure of Runge-Kutta methods. The Runge-Kutta methods form a group under the operation of composition. The multiplication operator has been overloaded so that multiplying two Runge-Kutta methods gives the method corresponding to their composition, with equal timesteps.

It is also possible to compose methods with non-equal timesteps using the compose() function.

Embedded Runge-Kutta Pairs

class nodepy.runge_kutta_method.ExplicitRungeKuttaPair(A=None, b=None, bhat=None, alpha=None, beta=None, alphahat=None, betahat=None, name='Runge-Kutta Pair', shortname='RKM', description='', order=(None, None))[source]

Class for embedded Runge-Kutta 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 Runge-Kutta 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 Shu-Osher 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 Runge-Kutta pair, the principal (usually higher-order) method is the one used if accuracy or stability properties are queried. Properties of the embedded (usually lower-order) 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 Runge-Kutta initialization, here the embedded coefficients \(\hat{b}_j\) are set as well.

Low-Storage Runge-Kutta methods

Typically, implementation of a Runge-Kutta 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 Runge-Kutta 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 low-storage methods, see [Ket10] .

In NodePy, low-storage methods are a subclass of explicit Runge-Kutta methods (and/or explicit Runge-Kutta pairs). In addition to the usual properties, they possess arrays of low-storage coefficients. They override the generic RK implementation of time-stepping and use special memory-efficient implementations instead. It should be noted that, although these low-storage 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 minimum-storage.

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

  • 2N methods pairs

  • 2R methods and embedded pairs

  • 3R methods and embedded pairs

Examples:

>>> from nodepy import lsrk, ivp
>>> myrk = lsrk.load_low_storage('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_low_storage('RK58[3R]C').__num__()
>>> 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/58-2S_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/58-3Sstar_acc.txt',lstype='3S*')
>>> rk3S.principal_error_norm() 
0.00035742076...

2S/3S methods

class nodepy.low_storage_rk.TwoSRungeKuttaMethod(betavec, gamma, delta, lstype, name='Low-storage Runge-Kutta Method', description='', shortname='LSRK2S', order=None)[source]

Class for low-storage Runge-Kutta 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 low-storage coefficient arrays \(\beta,\gamma,\delta\) follow the notation of [Ket10].

The argument lstype must be one of the following values:

  • 2S

  • 2S*

  • 3S*

Initializes the low-storage method by storing the low-storage coefficients and computing the Butcher coefficients.

2S/3S embedded pairs

class nodepy.low_storage_rk.TwoSRungeKuttaPair(betavec, gamma, delta, lstype, bhat=None, name='Low-storage Runge-Kutta Pair', description='', shortname='LSRK2S', order=None)[source]

Class for low-storage embedded Runge-Kutta 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 low-storage coefficient arrays \(\beta,\gamma,\delta\) follow the notation of [Ket10] .

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 low-storage pair by storing the low-storage coefficients and computing the Butcher coefficients.

2R/3R methods

class nodepy.low_storage_rk.TwoRRungeKuttaMethod(a, b, bhat=None, regs=2, name='2R Runge-Kutta Method', description='', shortname='LSRK2R', order=None)[source]

Class for 2R/3R/4R low-storage Runge-Kutta 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 low-storage 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 \(s-1\) containing the first subdiagonal of \(A\)

  • For 3R methods, a is a \(2\times s-1\) array whose first row contains the first subdiagonal of \(A\) and whose second row contains the second subdiagonal of \(A\).