runge_kutta_method

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.999999...
  

• Load a dictionary with many methods:

  >>> RK=loadRKM()
  >>> sorted(RK.keys()) 
  ['BE', 'BS3', 'BS5', 'BuRK65', 'CK5', 'CMR6', 'DP5', 'FE', 'Fehlberg43', 'Fehlberg45', 'GL2', 'GL3', 'HH5', 'HH5S', 'Heun22', 'Heun33', 'Lambert65', 'LobattoIIIA2', 'LobattoIIIA3', 'LobattoIIIC2', 'LobattoIIIC3', 'LobattoIIIC4', 'MTE22', 'Merson43', 'Mid22', 'NSSP32', 'NSSP33', 'PD8', 'RK44', 'RadauIIA2', 'RadauIIA3', 'SDIRK23', 'SDIRK34', 'SDIRK54', 'SS3', 'SSP104', 'SSP22', 'SSP22star', 'SSP33', 'SSP43', 'SSP53', 'SSP54', 'SSP63', 'SSP75', 'SSP85', 'SSP95', 'Soderlind43', ..., 'TR-BDF2', 'Tsit5', 'Zonneveld43']
  
  >>> print(RK['Mid22'])
  Midpoint Runge-Kutta
  
   0   |
   1/2 | 1/2
  _____|__________
       | 0    1
  

• Many methods are naturally implemented in some Shu-Osher 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:

• [But08]

• [HNorW93]

class nodepy.runge_kutta_method.RungeKuttaMethod(A=None, b=None, alpha=None, beta=None, name='Runge-Kutta Method', shortname='RKM', description='', mode='exact', order=None)

General class for implicit and explicit Runge-Kutta Methods. The method is defined by its Butcher array (\(A,b,c\)). It is assumed everywhere that \(c_i=\sum_j A_{ij}\).

A Runge-Kutta Method is initialized by providing either:

1. Butcher arrays \(A\) and \(b\) with valid and consistent dimensions; or

2. Shu-Osher arrays \(\alpha\) and \(\beta\) with valid and consistent dimensions

but not both.

The Butcher arrays are used as the primary representation of the method. If Shu-Osher arrays are provided instead, the Butcher arrays are computed by shu_osher_to_butcher().

Initialize a Runge-Kutta method. For explicit methods, the class ExplicitRungeKuttaMethod should be used instead.

TODO: make A a property and update c when it is changed

Now that we store (alpha, beta) as auxiliary data, maybe it’s okay to specify both (A, b) and (alpha, beta).

property p

Order of the method. This can be imposed and cached, which is advantageous to avoid issues with roundoff error and slow computation of the order conditions.

latex()

A LaTeX representation of the Butcher arrays.

Example:

>>> from nodepy import rk
>>> merson = rk.loadRKM('Merson43')
>>> print(merson.latex())
\begin{align}
\begin{array}{c|ccccc}
 &  &  &  &  & \\
\frac{1}{3} & \frac{1}{3} &  &  &  & \\
\frac{1}{3} & \frac{1}{6} & \frac{1}{6} &  &  & \\
\frac{1}{2} & \frac{1}{8} &  & \frac{3}{8} &  & \\
1 & \frac{1}{2} &  & - \frac{3}{2} & 2 & \\
\hline
 & \frac{1}{6} &  &  & \frac{2}{3} & \frac{1}{6}\\
 & \frac{1}{10} &  & \frac{3}{10} & \frac{2}{5} & \frac{1}{5}
\end{array}
\end{align}

print_shu_osher()

Pretty-prints the Shu-Osher arrays in the form:

  |        |
c | alpha  | beta
______________________
  | amp1   | bmp1

where amp1, bmp1 represent the last rows of \(\alpha, \beta\).

dj_reduce(tol=1e-13)

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

error_coefficient(tree, mode='exact')

Returns the coefficient in the Runge-Kutta method’s error expansion multiplying a single elementary differential, corresponding to a given tree.

Examples:

Construct an RK method and some rooted trees:
>>> from nodepy import rk, rt
>>> rk4 = rk.loadRKM('RK44')
>>> tree4 = rt.list_trees(4)[0]
>>> tree5 = rt.list_trees(5)[0]

The method has order 4, so this gives zero:
>>> rk4.error_coefficient(tree4)
0

This is non-zero, as the method doesn't
satisfy fifth-order conditions:
>>> rk4.error_coefficient(tree5)
-1/720

error_coeffs(p)

Returns the coefficients in the Runge-Kutta method’s error expansion multiplying all elementary differentials of the given order.

error_metrics(q=None, tol=1e-14)

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]

principal_error_norm(tol=1e-13, mode='float')

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

order(tol=1e-14, mode='float', extremely_high_order=False)

The order of a Runge-Kutta method.

Examples:

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> rk4.order()
4
>>> rk4.order(mode='exact')
4
>>> from nodepy.loadmethod import load_rk_from_file
>>> feagin14 = load_rk_from_file("rk1412.txt")
>>> feagin14.order(tol=1.0e-12, extremely_high_order=True)
This method has order at least 14.  Higher order conditions are not implemented in this routine.
14

mode == ‘float’: (default)

Check that conditions hold approximately, to within tolerance tol. Appropriate when coefficients are floating-point, or for faster checking of high-order methods.

mode == ‘exact’:

Check that conditions hold exactly. Appropriate when coefficients are specified as rational or algebraic numbers, but may be very slow for high order methods.

order_condition_residuals(p)

Generates and evaluates code to test whether a method satisfies the order conditions of order \(p\) (only).

effective_order(tol=1e-14)

Returns the effective order of a Runge-Kutta method. This may be higher than the classical order.

Example:

>>> from nodepy import rk
>>> RK4 = rk.loadRKM('RK44')
>>> RK4.effective_order()
Effective order is at least 4.  Higher effective order conditions not yet implemented.
4

effective_order_condition_residuals(q)

Generates and evaluates code to test whether a method satisfies the effective order \(q\) conditions (only).

Similar to order_condition_residuals(), but at the moment works only for \(q \le 4\). (enough to find Explicit SSPRK)

stage_order(tol=1e-14)

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:

• [DV84]

• [But08]

stability_function_unexpanded()

Compute the stability function expression but don’t simplify it. This can be useful for performance reasons.

Example:

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> rk4.stability_function_unexpanded()
z*(z/2 + 1)/3 + z*(z*(z/2 + 1)/2 + 1)/3 + z*(z*(z*(z/2 + 1)/2 + 1) + 1)/6 + z/6 + 1
>>> rk4.stability_function_unexpanded().simplify()
z**4/24 + z**3/6 + z**2/2 + z + 1

stability_function(stage=None, mode='exact', formula='lts', use_butcher=False)

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))

E_polynomial()

Return the E-polynomial of the method.

plot_stability_function(bounds=[-20, 1])

Plot the value of the stability function along the negative real axis.

Example:

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

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)

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...

plot_order_star(N=200, bounds=[-5, 5, -5, 5], plotroots=False, color=('w', 'b'), filled=True, fignum=None)

The order star of a Runge-Kutta method is the set

$$ \{ z \in C : | \phi(z)/\exp(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, order star is filled in (solid); otherwise it is outlined

Example::

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> rk4.plot_order_star() 
<Figure size...

circle_contractivity_radius(acc=1e-13, rmax=1000)

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...

absolute_monotonicity_radius(acc=1e-10, rmax=200, tol=3e-16)

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.

linear_monotonicity_radius(acc=1e-10, tol=1e-15, tol2=1e-08)

Computes Horvath’s monotonicity radius of the stability function.

TODO: clean this up.

is_algebraically_stable()

Checks whether the Runge-Kutta method is algebraically stable, i.e. whether the matrix

\[B A + A^T B - b b^T\]

is positive semidefinite and all weights \(b_i \geq 0\).

Examples:

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> rk4.is_algebraically_stable() 
False

>>> from nodepy import rk
>>> lobatto = rk.loadRKM('LobattoIIIC4')
>>> lobatto.is_algebraically_stable() 
True

See [But08].

optimal_shu_osher_form()

Gives a Shu-Osher form in which the SSP coefficient is evident (i.e., in which \(\alpha_{ij}, \beta_{ij} \ge 0\) and \(\alpha_{ij} / \beta_{ij} = c\) for every \(\beta_{ij} \ne 0\)).

Input:

• A RungeKuttaMethod

Output:

• alpha, beta – Shu-Osher arrays

The ‘optimal’ Shu-Osher arrays are given by

$$\alpha= K(I+cA)^{-1}$$ $$\beta = c \alpha$$

where \(K=[ A; b^T]\).

Example:

>>> from nodepy import rk
>>> rk2 = rk.loadRKM('MTE22')
>>> rk2.optimal_shu_osher_form()
(array([[0, 0, 0],
       [1.00000000000000, 0, 0],
       [0.625000000060027, 0.374999999939973, 0]], dtype=object), array([[0, 0, 0],
       [0.666666666666667, 0, 0],
       [4.00177668780088e-11, 0.750000000000000, 0]], dtype=object))

See [Hig05].

canonical_shu_osher_form(r)

Returns \(d, P\) where \(P\) is the matrix \(P=r(I+rK)^{-1}K\) and \(d\) is the vector \(d=(I+rK)^{-1}e=(I-P)e\).

Note that this can be computed for any value of \(r\), including values for which \(d, P\) may have negative entries.

lp_perturb(r, tol=None)

Find a perturbation via linear programming.

Use linear programming to determine if there exists a perturbation of this method with radius of absolute monotonicity at least \(r\).

The linear program to be solved is

\[\begin{split}\begin{align} (I-2\alpha^{down}_r)\alpha_r + \alpha^{down}_r & = (\alpha^{up}_r ) \ge 0 \\ (I-2\alpha^{down}_r)v_r & = \gamma_r \ge 0. \end{align}\end{split}\]

This function requires cvxpy.

ssplit(r, P_signs=None, delta=None)

Sympy exact version of split.

If \(P_{signs}\) is passed, use that as the sign pattern of the \(P\) matrix. This is useful if \(r\) is symbolic (since then in general the signs of elemnts of \(P\) are unknown).

resplit(r, tol=1e-15, max_iter=5)

is_splittable(r, tol=1e-15)

optimal_perturbed_splitting(acc=1e-12, rmax=50.01, tol=1e-13, algorithm='split')

Return the optimal downwind splitting of the method along with the optimal downwind SSP coefficient.

The default algorithm (split with iteration) is not provably correct. The LP algorithm is. See the paper (Higueras & Ketcheson) for more details.

Example:

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> r, d, alpha, alphatilde = rk4.optimal_perturbed_splitting(algorithm='split')
>>> print(r) 
0.68501606...

propagation_matrix(L, dt)

Returns the solution propagation matrix for the linear autonomous system with RHS equal to the matrix \(L\), i.e. it returns the matrix \(G\) such that when the Runge-Kutta method is applied to the system \(u'(t)=Lu\) with stepsize \(dt\), the numerical solution is given by \(u^{n+1} = G u^n\).

Input:

• self – a Runge-Kutta method

• L – the RHS of the ODE system

• dt – the timestep

The formula for \(G\) is (if \(L\) is a scalar): \(G = 1 + b^T L (I-A L)^{-1} e\)

where \(A\) and \(b\) are the Butcher arrays and \(e\) is the vector of ones. If \(L\) is a matrix, all quantities above are replaced by their Kronecker product with the identity matrix of size \(m\), where \(m\) is the number of stages of the Runge-Kutta method.

is_explicit()

is_zero_stable()

is_FSAL()

True if method is “First Same As Last”.

nodepy.runge_kutta_method.sign_split(M)

Given a matrix \(M\), return two matrices. The first contains the positive entries of \(M\); the second contains the negative entries of \(M\), multiplied by \(-1\).

nodepy.runge_kutta_method.redistribute_gamma(gamma, alpha_up, alpha_down)

class nodepy.runge_kutta_method.ExplicitRungeKuttaMethod(A=None, b=None, alpha=None, beta=None, name='Runge-Kutta Method', shortname='RKM', description='', mode='exact', order=None)

Class for explicit Runge-Kutta methods. Mostly identical to RungeKuttaMethod, but also includes time-stepping and a few other functions.

Initialize a Runge-Kutta method. For explicit methods, the class ExplicitRungeKuttaMethod should be used instead.

TODO: make A a property and update c when it is changed

Now that we store (alpha, beta) as auxiliary data, maybe it’s okay to specify both (A, b) and (alpha, beta).

imaginary_stability_interval(mode='exact', eps=1e-14)

Length of imaginary axis half-interval contained in the method’s region of absolute stability.

Examples:

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> rk4.imaginary_stability_interval() 
2.8284271247461...

real_stability_interval(mode='exact', eps=1e-14)

Length of negative real axis interval contained in the method’s region of absolute stability.

Examples:

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> I = rk4.real_stability_interval()
>>> print("{:.10f}".format(I))
2.7852935634

linear_absolute_monotonicity_radius(acc=1e-10, rmax=50, tol=3e-16)

Returns the radius of absolute monotonicity of the stability function of a Runge-Kutta method.

TODO: implement this functionality for implicit methods.

is_explicit()

work_per_step()

Number of function evaluations required for one step.

num_seq_dep_stages()

Number of sequentially dependent stages.

Number of sequential function evaluations that must be made.

Examples:

Extrapolation methods are parallelizable:
>>> from nodepy import rk
>>> ex4 = rk.extrap(4)
>>> len(ex4)
7
>>> ex4.num_seq_dep_stages()
4

So are deferred correction methods:
>>> dc4 = rk.DC(4)
>>> len(dc4)
17
>>> dc4.num_seq_dep_stages()
8

Unless `\theta` is non-zero:
>>> rk.DC(4,theta=1).num_seq_dep_stages()
20

internal_stability_polynomials(stage=None, mode='exact', formula='lts', use_butcher=False)

The internal stability polynomials of a Runge-Kutta method depend on the implementation and must therefore be constructed base on the Shu-Osher form used for the implementation. By default the Shu-Osher coefficients are used. The Butcher coefficients are used if use_butcher=True or if Shu-Osher coefficients are not defined.

The formula for the polynomials is: Modified Shu-Osher form: \((alphastarmp1+z betastarmp1)(I-alphastar-z betastar)^{-1}\) Butcher array: \(z b^T(I-zA)^{-1}\).

Note that in the first stage no perturbation is introduced because for an explicit method the first stage is equal to the solution at the current time level. Therefore, the first internal polynomial is set to zero.

For symbolic computation, this routine has been significantly modified for efficiency relative to particular classes of methods. Two formulas are implemented, one based on SymPy’s Matrix.lower_triangular_solve() and the other using a power series for the inverse. Different choices of these two are more efficient for different classes of methods (this only matters for methods with very many stages).

Options

• use_butcher

Output:

• numpy array of internal stability polynomials

Examples:

>>> from nodepy import rk
>>> rk4 = rk.loadRKM('RK44')
>>> theta = rk4.internal_stability_polynomials()
>>> for p in theta:
...     print(p)
         3          2
0.08333 x + 0.1667 x + 0.3333 x
        2
0.1667 x + 0.3333 x

0.1667 x

internal_stability_plot(bounds=None, N=200, use_butcher=False, formula='lts', levels=[1, 100, 500, 1000, 1500, 10000])

Plot internal stability regions.

Plots the \(\epsilon\)-internal-stability region contours.

By default the Shu-Osher coefficients are used. The Butcher coefficients are used if use_butcher=True or if Shu-Osher coefficients are not defined.

Examples:

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

maximum_internal_amplification(N=200, use_butcher=False, formula='lts')

The maximum amount by which any stage error is amplified, assuming the step size is taken so that the method is absolutely stable:

\(\max_{z \in S,j} |\theta_j(z)|\)

where \(S = \{z \in C : |R(z)|\le 1.\)

Here \(R(z)\) is the stability function and \(\theta_j(z)\) are the internal stability functions.

By default the Shu-Osher coefficients are used. The Butcher coefficients are used if use_butcher=True or if Shu-Osher coefficients are not defined.

Examples:

>>> from nodepy import rk
>>> ssp2 = rk.SSPRK2(6)
>>> ssp2.maximum_internal_amplification() 
(1.097405..., 0.83333333...)
>>> ssp2.maximum_internal_amplification(use_butcher=True) 
(2.037051..., 0.0)

class nodepy.runge_kutta_method.ContinuousRungeKuttaMethod(A=None, b=None, alpha=None, beta=None, b_dense=None, name='Continuous Runge-Kutta Method', shortname='CRKM', description='', mode='exact', order=None)

Initialize a Runge-Kutta method. For explicit methods, the class ExplicitRungeKuttaMethod should be used instead.

TODO: make A a property and update c when it is changed

Now that we store (alpha, beta) as auxiliary data, maybe it’s okay to specify both (A, b) and (alpha, beta).

class nodepy.runge_kutta_method.ContinuousExplicitRungeKuttaMethod(A=None, b=None, alpha=None, beta=None, b_dense=None, name='Continuous Runge-Kutta Method', shortname='CRKM', description='', mode='exact', order=None)

Initialize a Runge-Kutta method. For explicit methods, the class ExplicitRungeKuttaMethod should be used instead.

TODO: make A a property and update c when it is changed

Now that we store (alpha, beta) as auxiliary data, maybe it’s okay to specify both (A, b) and (alpha, beta).

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))

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.

property main_method

Return the main method of the pair (usually the higher-order one).

property embedded_method

Always recompute the embedded method on the fly. This may be inefficient.

error_metrics(q=None, p=None, tol=1e-14)

Return full set of error metrics for an embedded RK pair. See [KCL00] p. 181

Example:

>>> from nodepy import rk
>>> bs5 = rk.loadRKM('BS5')
>>> bs5.error_metrics()
main method has order 5
embedded method has order 4
(43*sqrt(83011)/558835200, 43/3386880, sqrt(29695176594765489880490334265)/810521680634265600, 1451/15966720, sqrt(870269901055795)/277898765760, 10147/131580855, sqrt(51577359825120524319571156056057595)/219308015066060340, 26201089/40912704, sqrt(5250600078722255566247933273951710555)/2193080150660603400, 305343067/400035328, 482048/414219, 5987277*sqrt(72241974756542598745)/243675572295622600, 5987277/36366848)

is_FSAL()

True if method is “First Same As Last”.

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)

Plot the absolute stability region of an RK pair. By default, the region of the main method is filled in red and the region of the embedded method is outlined in black.

Example:

>>> from nodepy import rk
>>> bs5 = rk.loadRKM('BS5')
>>> bs5.plot_stability_region() 
<Figure size...

plot_I_controller_stability(beta1=1.0, N=200, color='r', filled=True, bounds=None, plotroots=False, alpha=1.0, scalefac=1.0, longtitle=True, to_file_region=False, to_file_controller=False, fignum_region=None, fignum_controller=None, plot_region=True)

Plot the absolute stability region and the function characterizing stepsize control stability for an I controller of an RK pair, cf. [HW96]. The I controller is of the form

\[\begin{equation} h_{n+1} = \left( \frac{\mathrm{TOL}}{\mathrm{err}_{n}} \right)^{\beta_1/k} h_{n}, \end{equation}\]

where \(h\) is the stepsize, \(TOL\) the tolerance, and \(err = O(h^k)\) the error estimate.

By default, the region of the main method is filled in red and the region of the embedded method is outlined in black.

Example:

>>> from nodepy import rk
>>> bs5 = rk.loadRKM('BS5')
>>> bs5.plot_I_controller_stability() 
(<Figure size...

plot_PI_controller_stability(beta1=0.6666666666666666, beta2=-0.3333333333333333, N=200, color='r', filled=True, bounds=None, plotroots=False, alpha=1.0, scalefac=1.0, longtitle=True, to_file_region=False, to_file_controller=False, fignum_region=None, fignum_controller=None, plot_region=True)

Plot the absolute stability region and the function characterizing stepsize control stability for a PI controller of an RK pair, cf. [HW96]. The PI controller is of the form

\[\begin{equation} h_{n+1} = \left( \frac{\mathrm{TOL}}{\mathrm{err}_{n}} \right)^{\beta_1/k} \left( \frac{\mathrm{TOL}}{\mathrm{err}_{n-1}} \right)^{\beta_2/k} h_{n}, \end{equation}\]

where \(h\) is the stepsize, \(TOL\) the tolerance, and \(err = O(h^k)\) the error estimate.

By default, the region of the main method is filled in red and the region of the embedded method is outlined in black.

Example:

>>> from nodepy import rk
>>> bs5 = rk.loadRKM('BS5')
>>> bs5.plot_PI_controller_stability() 
(<Figure size...

plot_PID_controller_stability(beta1=0.49, beta2=-0.34, beta3=0.1, N=200, color='r', filled=True, bounds=None, plotroots=False, alpha=1.0, scalefac=1.0, longtitle=True, to_file_region=False, to_file_controller=False, fignum_region=None, fignum_controller=None, plot_region=True)

Plot the absolute stability region and the function characterizing stepsize control stability for a PID controller of an RK pair. The PID controller is of the form

\[\begin{equation} h_{n+1} = \left( \frac{\mathrm{TOL}}{\mathrm{err}_{n}} \right)^{\beta_1/k} \left( \frac{\mathrm{TOL}}{\mathrm{err}_{n-1}} \right)^{\beta_2/k} \left( \frac{\mathrm{TOL}}{\mathrm{err}_{n-2}} \right)^{\beta_3/k} h_{n}, \end{equation}\]

where \(h\) is the stepsize, \(TOL\) the tolerance, and \(err = O(h^k)\) the error estimate.

By default, the region of the main method is filled in red and the region of the embedded method is outlined in black.

Example:

>>> from nodepy import rk
>>> bs5 = rk.loadRKM('BS5')
>>> bs5.plot_PID_controller_stability() 
(<Figure size...

nodepy.runge_kutta_method.elementary_weight(tree)

Constructs Butcher’s elementary weights for a Runge-Kutta method

Currently doesn’t work right; note that two of the 5th-order weights appear identical. The _str version below works correctly and produces NumPy code. But it would be nice to have this version working so that we could symbolically simplify the expressions.

In order to do things correctly, we need a symbolic system that includes support for either:

• Two different types of multiplication; or

• Full tensor expressions

The latter is now available in Sympy, and I’ve started a test implementation. The main issue now is that things like

$$ A \times A^2 $$

don’t get parentheses when they really mean

$$ (A \times A)^2. $$

It’s not really a bug since \(A\times(A^2)\) does show parentheses, but it will make it harder to parse into code.

Examples:

>>> from nodepy import rk, rt
>>> tree = rt.list_trees(2)[0]
>>> tree
'{T}'
>>> rk.elementary_weight(tree)
b*c

See [But08].

nodepy.runge_kutta_method.elementary_weight_str_jump(tree, parent_label_index=None, node_label_index=0, method_type='ERK')

Constructs elementary weights for a Runge-Kutta method as strings suitable for setting constraints in JuMP.

Here the indexing is set up for explicit methods.

Examples:

>>> from nodepy import rk, rt
>>> tree = rt.list_trees(5)[3]
>>> rk.elementary_weight_str_jump(tree)
'sum(b[i]*sum(A[i,j]*c[j]*sum(A[j,k]*c[k] for k in 1:j-1) for j in 1:i-1) for i in 1:s)'

nodepy.runge_kutta_method.elementary_weight_str_matlab(tree, root=True)

Constructs elementary weights for a Runge-Kutta method as strings suitable for MATLAB execution.

Examples:

>>> from nodepy import rk, rt
>>> tree = rt.list_trees(5)[3]
>>> rk.elementary_weight_str_matlab(tree)
"(b'*((A*((A*c).*c))))"
>>> rk.elementary_weight_str_matlab(rt.RootedTree('{T^10}'))
"(b'*(c.^10))"
>>> rk.elementary_weight_str_matlab(rt.RootedTree('{{T^11}T}'))
"(b'*((A*(c.^11))))"

nodepy.runge_kutta_method.elementary_weight_str(tree, style='python')

Constructs Butcher’s elementary weights for a Runge-Kutta method as strings suitable for numpy execution.

For MATLAB code or JuMP code, use the corresponding dedicated function instead of this one.

Examples:

>>> from nodepy import rk, rt
>>> tree = rt.list_trees(5)[0]
>>> rk.elementary_weight_str(tree)
'dot(b,dot(A,c**3))'

>>> rk.elementary_weight_str(rt.RootedTree('{T^10}'))
'dot(b,c**10)'
>>> rk.elementary_weight_str(rt.RootedTree('{{T^11}T}'))
'dot(b,dot(A,c**11))'

nodepy.runge_kutta_method.RKeta(tree)

nodepy.runge_kutta_method.RKeta_str(tree)

Computes eta(t) for Runge-Kutta methods

nodepy.runge_kutta_method.discrete_adjoint(meth)

Returns the discrete adjoint of a Runge-Kutta method

nodepy.runge_kutta_method.is_absolutely_monotonic_poly(r, tol, p)

Returns \(1\) if the polynomial \(p\) is absolutely monotonic at \(z=-r\).

nodepy.runge_kutta_method.shu_osher_change_alpha_ij(alpha, beta, i, j, val)

Input:

• alpha, beta: Shu-Osher arrays

• i,j: indices

• val – real number

Output: Shu-Osher arrays alpha, beta with

alpha[i,j]=alphaa[i,j] + val.

nodepy.runge_kutta_method.shu_osher_zero_alpha_ij(alpha, beta, i, j)

Input:

• Shu-Osher arrays alpha, beta

• Indices i,j

Output: Shu-Osher arrays alpha, beta with alpha[i,j]=0.

nodepy.runge_kutta_method.shu_osher_zero_beta_ij(alpha, beta, i, j)

Input:

• Shu-Osher arrays alpha, beta

• Indices i,j

Output:

• Shu-Osher arrays alpha, beta with beta[i,j]=0.

nodepy.runge_kutta_method.shu_osher_to_butcher(alpha, beta)

Accepts a Shu-Osher representation of a Runge-Kutta

method and returns the Butcher coefficients.

\begin{align*} A = & (I-\alpha_0)^{-1} \beta_0 \\ b = & \beta_1 + \alpha_1 \end{align*}

See [GKS09].

nodepy.runge_kutta_method.loadRKM(which='All')

Load a set of standard Runge-Kutta methods for testing. The following methods are included:

Explicit methods:

• ‘FE’: Forward (explicit) Euler

• ‘RK44’: Classical 4-stage 4th-order

• ‘SSP22’: Trapezoidal rule 2nd-order [SO88]

• ‘MTE22’: Minimal truncation error 2-stage 2nd-order

• ‘Mid22’: Explicit midpoint 2-stage 2nd-order

• ‘SSP33’: Optimal 3rd-order SSP method of Shu & Osher [SO88] with embedded method of [CFS18]

• ‘SSP43’: Optimal 3rd-order SSP method of Kraaijevanger [Kra91] with embedded method of [CFS18]

• ‘Heun33’: Third-order method of Heun [Heu00]

• ‘SSP22star’: Optimal 2nd-order downwind SSP

• ‘NSSP32’: [WS07]

• ‘NSSP33’: [WS07]

• ‘SSP104’: Optimal 10-stage, 4th-order SSP method [Ket08] with embedded method of [CFS18]

• ‘Merson43’ Merson 4(3) pair [HNorW93] pg. 167

• ‘DP5’: Dormand-Prince 5th-order [DP80]

• ‘PD8’: Prince-Dormand 8th-order and 7th-order pair [PD81]

• ‘CMR6’: Calvo et al.’s 6(5) pair [CMRandez90]

• ‘Zonneveld43’: 4(3) pair of Zonneveld [HNorW93] Table 4.2

• ‘Soderlind43’: 4(3) pair of Söderlind (also available as ‘Söderlind43’)

• ‘Fehlberg43’: 4(3) pair of Fehlberg [Feh69]

• ‘Fehlberg45’: 5(4) pair of Fehlberg [Feh69]

• ‘Lambert65’:

• ‘Tsit5’: 5(4) pair of Tsitouras [Tsi11]

• ‘HH5’: 5(4) pair of Highham and Hall with lower error constant [HH90]

• ‘HH5S’: 5(4) pair of Highham and Hall with increased stepsize stability [HH90]

• ‘BuRK65’: 6-stage, 5th-order method of Butcher

• ‘CK5’: Cash-Karp 5(4)6 [CK90]

• ‘BS3’: Bogacki-Shampine 3(2)4 pair [BS89]

• ‘BS5’: Bogacki-Shampine 5(4)8 pair [BS96]

• ‘SS3’: Sharp-Smart 3(2)4 pair [SS93]

• ‘SSP75’: Ruuth-Spiteri optimal downwind SSP [RS04]

• ‘SSP85’: Ruuth-Spiteri optimal downwind SSP [RS04]

• ‘SSP95’: Ruuth-Spiteri optimal downwind SSP [RS04]

• ‘SSP54’: Ruuth-Spiteri optimal downwind SSP [RS04]

• ‘SSP53’: Ruuth-Spiteri optimal downwind SSP [RS04]

• ‘SSP63’: Ruuth-Spiteri optimal downwind SSP [RS04]

Diagonally Implicit methods:

• ‘BE’: Backward Euler

• ‘SDIRK23’: [Nor74]

• ‘SDIRK34’: [Nor74]

• ‘SDIRK54’: [HW96]

• ‘TR-BDF2’: [BCF+85]

Fully Implicit methods:

• ‘GL2’: 2-stage Gauss-Legendre; see e.g. [But64]

• ‘GL3’: 3-stage Gauss-Legendre; see e.g. [But64]

• ‘LobattoIIIA2’: [Ehl69]

• ‘LobattoIIIA3’: [Ehl69]

• ‘LobattoIIIC2’: [Chi71]

• ‘LobattoIIIC3’: [Chi71]

• ‘LobattoIIIC4’: [Chi71]

• ‘RadauIIA2’: [Ehl69]

• ‘RadauIIA3’: [Ehl69]

nodepy.runge_kutta_method.RK22_family(gamma)

Construct a 2-stage second order Runge-Kutta method

Input: gamma – family parameter

Output: An ExplicitRungeKuttaMethod

Examples:

>>> from nodepy import rk
>>> print(rk.RK22_family(-1))
Runge-Kutta Method

 0    |
 -1/2 | -1/2
______|____________
      | 2     -1

nodepy.runge_kutta_method.RK44_family(w)

Construct a 4-stage fourth order Runge-Kutta method

Input: w – family parameter

Output: An ExplicitRungeKuttaMethod

Examples:

>>> from nodepy import rk
>>> print(rk.RK44_family(1))
Runge-Kutta Method

 0    |
 1/2  | 1/2
 1/2  | 1/3   1/6
 1    |       -2    3
______|________________________
      | 1/6   -1/3  1     1/6

nodepy.runge_kutta_method.SSPRK2(m)

Construct the optimal m-stage, second order SSP Explicit Runge-Kutta method (\(m \ge 2\)).

Input: m – number of stages

Output: A ExplicitRungeKuttaMethod

Examples:

Load the 4-stage method:
>>> SSP42=SSPRK2(4)
>>> print(SSP42)
SSPRK(4,2)
<BLANKLINE>
 0   |
 1/3 | 1/3
 2/3 | 1/3  1/3
 1   | 1/3  1/3  1/3
_____|____________________
     | 1/4  1/4  1/4  1/4

>>> SSP42.absolute_monotonicity_radius()
2.999999999974534

See [Ket08].

nodepy.runge_kutta_method.SSPRK3(m)

Construct the optimal m-stage third order SSP Runge-Kutta method (\(m = n^2\), \(n \ge 2\))

Input: m – number of stages

Output: A RungeKuttaMethod

Examples:

Load the 4-stage method:
>>> SSP43=SSPRK3(4)
>>> print(SSP43)
SSPRK43
<BLANKLINE>
 0   |
 1/2 | 1/2
 1   | 1/2  1/2
 1/2 | 1/6  1/6  1/6
_____|____________________
     | 1/6  1/6  1/6  1/2

>>> SSP43.absolute_monotonicity_radius()
1.9999999999527063

See [Ket08].

nodepy.runge_kutta_method.SSPRKm(m)

Construct the optimal m-stage, linearly mth order SSP Explicit Runge-Kutta method (\(m \ge 2\)).

Input: m – number of stages

Output: A ExplicitRungeKuttaMethod

Examples:

Load the 4-stage method:
>>> SSP44=SSPRKm(4)
>>> print(SSP44)
SSPRK44
<BLANKLINE>
 0    |
 1    | 1
 2    | 1     1
 3    | 1     1     1
______|________________________
      | 5/8   7/24  1/24  1/24


>>> SSP44.absolute_monotonicity_radius()
0.9999999999308784

See [GST01].

nodepy.runge_kutta_method.SSPIRK1(m)

Construct the m-stage, first order unconditionally SSP Implicit Runge-Kutta method with smallest coefficient of \(z^2\) (in the stability polynomial)

Input: m – number of stages

Output: A RungeKuttaMethod

Examples:

Load the 4-stage method:
>>> ISSP41=SSPIRK1(4)
>>> print(ISSP41)
SSPIRK41
<BLANKLINE>
 1/4 | 1/4
 1/2 | 1/4  1/4
 3/4 | 1/4  1/4  1/4
 1   | 1/4  1/4  1/4  1/4
_____|____________________
     | 1/4  1/4  1/4  1/4

nodepy.runge_kutta_method.SSPIRK2(m)

Construct the optimal m-stage, second order SSP Implicit Runge-Kutta method (\(m \ge 2\)).

Input: m – number of stages

Output: A RungeKuttaMethod.

Examples:

Load the 4-stage method:
>>> ISSP42=SSPIRK2(4)
>>> print(ISSP42)
SSPIRK42
<BLANKLINE>
 1/8 | 1/8
 3/8 | 1/4  1/8
 5/8 | 1/4  1/4  1/8
 7/8 | 1/4  1/4  1/4  1/8
_____|____________________
     | 1/4  1/4  1/4  1/4

>>> ISSP42.absolute_monotonicity_radius() # doctest: +ELLIPSIS
7.99...

See [KMG09].

nodepy.runge_kutta_method.SSPIRK3(m)

Construct the optimal m-stage, third order SSP Implicit Runge-Kutta method (\(m \ge 2\)).

Input: m – number of stages

Output: A RungeKuttaMethod

Examples:

Load the 4-stage method:
>>> ISSP43=SSPIRK3(4)
>>> print(ISSP43)  # doctest: +NORMALIZE_WHITESPACE
SSPIRK43
<BLANKLINE>
 1/2 - sqrt(15)/10 | 1/2 - sqrt(15)/10
 1/2 - sqrt(15)/30 | sqrt(15)/15        1/2 - sqrt(15)/10
 sqrt(15)/30 + 1/2 | sqrt(15)/15        sqrt(15)/15        1/2 - sqrt(15)/10
 sqrt(15)/10 + 1/2 | sqrt(15)/15        sqrt(15)/15        sqrt(15)/15        1/2 - sqrt(15)/10
___________________|____________________________________________________________________________
                   | 1/4                1/4                1/4                1/4

>>> x=ISSP43.absolute_monotonicity_radius()
>>> print("{:.5f}".format(x))
6.87298

See [KMG09].

nodepy.runge_kutta_method.RKC1(m, epsilon=0)

Construct the m-stage, first order explicit Runge-Kutta-Chebyshev methods of Verwer (\(m \ge 1\)).

‘epsilon’ is a damping parameter used to avoid tangency of the stability region boundary to the negative real axis.

Input: m – number of stages

Output: A ExplicitRungeKuttaMethod

Examples:

Load the 4-stage method:
>>> RKC41=RKC1(4)
>>> print(RKC41)
Runge-Kutta-Chebyshev (4,1)
<BLANKLINE>
 0    |
 1/16 | 1/16
 1/4  | 1/8   1/8
 9/16 | 3/16  1/4   1/8
______|________________________
      | 1/4   3/8   1/4   1/8

See [VS04].

nodepy.runge_kutta_method.RKC2(m, epsilon=0)

Construct the m-stage, second order Explicit Runge-Kutta-Chebyshev methods of Verwer (\(m \ge 2\)).

Inputs:

• m – number of stages

• epsilon – damping factor

Output: A ExplicitRungeKuttaMethod

Examples:

Load the 4-stage method:
>>> RKC42=RKC2(4)
>>> print(RKC42)
Runge-Kutta-Chebyshev (4,2)
<BLANKLINE>
 0      |
 1/5    | 1/5
 1/5    | 1/10    1/10
 8/15   | -8/45   32/135  64/135
________|________________________________
        | -51/64  3/8     1       27/64

See [VS04].

nodepy.runge_kutta_method.dcweights(x)

Takes a set of abscissae \(x\) and an index \(i\), and returns the quadrature weights for the interval \([x_i, x_{i+1}]\). Used in construction of deferred correction methods.

nodepy.runge_kutta_method.DC_pair(s, theta=0.0, grid='eq')

Spectral deferred correction embedded pairs. See also the help for DC().

Examples:

>>> from nodepy import rk
>>> DC2 = rk.DC_pair(2)
>>> print(DC2)
Picard 3(2)

 0     |
 1/2   | 1/2
 1     | 1/2    1/2
 1/2   | 5/24   1/3    -1/24
 1     | 1/6    2/3    1/6
_______|___________________________________
       | 1/6    0      0      2/3    1/6
       | 1/6    2/3    1/6

nodepy.runge_kutta_method.DC(s, theta=0, grid='eq', num_corr=None)

Spectral deferred correction methods. For now, based on explicit Euler and equispaced points. For theta=0, this is Picard iteration.

Input: s – number of grid points and number of correction iterations

Output: A ExplicitRungeKuttaMethod

Note that the number of stages is NOT equal to \(s\). The order is equal to \(s+1\).

Examples:

>>> from nodepy import rk
>>> dc3 = rk.DC(3)
>>> dc3.order()
4
>>> dc3.principal_error_norm() 
0.0069444...
>>> dc3_cheb = rk.DC(3,grid='cheb')
>>> dc3_cheb.order()
4
>>> dc3_cheb.principal_error_norm() 
0.0066478...

See [DGR00, GKS09].

nodepy.runge_kutta_method.extrap(k, base='euler', seq='harmonic', embedded=False, shuosher=False)

Construct extrapolation methods as Runge-Kutta methods.

Input:

• k – number of grid points & number of extrapolation iterations

• base – the base method to be used (‘euler’, ‘midpoint’, or ‘implicit euler’)

• seq – extrapolation sequence

Output: A ExplicitRungeKuttaMethod (or

RungeKuttaMethod if base==’implicit euler’)

Examples:

>>> from nodepy import rk
>>> ex3 = rk.extrap(3)
>>> print(ex3)
Ex-Euler 3

 0   |
 1/2 | 1/2
 1/3 | 1/3
 2/3 | 1/3       1/3
_____|____________________
     | 0    -2   3/2  3/2

>>> ex3.num_seq_dep_stages()
3
>>> ex3.principal_error_norm()
0.04606423319938055
>>> ex3.principal_error_norm(mode='exact')
sqrt(11)/72

>>> ex4 = rk.extrap(2,'midpoint')
>>> print(ex4)
Ex-Midpoint 2

 0    |
 1/2  | 1/2
 1/4  | 1/4
 1/2  |             1/2
 3/4  | 1/4               1/2
______|______________________________
      | 0     -1/3  2/3   0     2/3

>>> ex4.order()
4

See [HNorW93] chapter II.9.

nodepy.runge_kutta_method.extrap_pair(p, base='euler')

Returns an embedded RK pair. If the base method is Euler, the prinicpal method has order \(p\) and the embedded method has order \(p-1\). If the base method is midpoint, the orders are \(2p, 2(p-1)\).

Examples:

>>> from nodepy import rk
>>> ex32 = rk.extrap_pair(3,base='Euler')
>>> ex32.order()
3
>>> ex32.embedded_method.order()
2
>>> ex42 = rk.extrap_pair(2,base='midpoint')
>>> ex42.order()
4
>>> ex42.embedded_method.order()
2

nodepy.runge_kutta_method.runge_kutta_order_conditions(p, ind='all')

This is the current method of producing the code on-the-fly to test order conditions for RK methods. May be deprecated soon.

nodepy.runge_kutta_method.RKOCstr2code(ocstr)

Converts output of runge_kutta_order_conditions() to numpy-executable code.

nodepy.runge_kutta_method.compose(RK1, RK2, h1=1, h2=1)

The method obtained by applying RK2, followed by RK1, each with half the timestep.

Output:

The method
     c_2 | A_2  0
   1+c_1 | b_2 A_1
   _____________
         | b_2 b_1

but with everything divided by two.
The `b_2` matrix block consists of `m_1` (row) copies of `b_2`.

Examples:

What method is obtained by two successive FE steps?
>>> from nodepy import rk
>>> fe=rk.loadRKM('FE')
>>> print(fe*fe)
Runge-Kutta Method
<BLANKLINE>
 0     |
 0.500 | 0.500
_______|______________
       | 0.500  0.500

TODO: Generalize this for any number of inputs

nodepy.runge_kutta_method.python_to_fortran(code)

nodepy.runge_kutta_method.python_to_matlab(code)

Convert Python code string (order condition) to MATLAB code string Doesn’t really work yet. We need to do more parsing.

nodepy.runge_kutta_method.relative_accuracy_efficiency(rk1, rk2, mode='float', tol=1e-14)

Compute the accuracy efficiency of method rk1 relative to that of rk2, for two methods with the same order of accuracy.

The relative accuracy efficiency is

$$ \eta = \frac{s_2}{s_1} \left(\frac{A_2}{A_1}\right)^{1/p+1} $$

where \(s_1,s_2\) are the number of stages of the two methods and \(A_1,A_2\) are their principal error norms.

If the result is \(>1\), method 1 is more efficient.

Examples:

Compare Fehlberg's method with Dormand-Prince
>>> from nodepy import rk
>>> dp5 = rk.loadRKM('DP5')
>>> f45 = rk.loadRKM('Fehlberg45')
>>> rk.relative_accuracy_efficiency(dp5,f45) # doctest: +ELLIPSIS
1.22229116499...

nodepy.runge_kutta_method.accuracy_efficiency(rk1, parallel=False, mode='float', tol=1e-14, p=None)

Compute the accuracy efficiency of method rk1.

The accuracy efficiency is

$$ \eta = \frac{1}{s_1} \left(\frac{1}{A_1}\right)^{1/p+1} $$

where \(s_1\) are the number of stages of the the method and \(A_1\) is its principal error norms.

Examples:

Accuracy efficiency of Dormand-Prince
>>> from nodepy import rk
>>> dp5 = rk.loadRKM('DP5')
>>> rk.accuracy_efficiency(dp5) # doctest: +ELLIPSIS
0.5264921944121...

nodepy.runge_kutta_method.linearly_stable_step_size(rk, L, acc=1e-07, tol=1e-13, plot=1)

Determine the maximum linearly stable step size for Runge-Kutta method rk applied to the IVP \(u' = Lu\), by computing the eigenvalues of \(L\) and determining the values of the stability function of rk at the eigenvalues.

Note that this analysis is not generally appropriate if \(L\) is non-normal.

Examples:

>>> from nodepy import rk, semidisc

4th-order Runge-Kutta scheme:
>>> rk44=rk.loadRKM('RK44')

Centered differences on a grid with spacing 1/100:
>>> L1=semidisc.centered_diffusion_matrix(100)
>>> L2=semidisc.centered_advection_diffusion_matrix(1.,1./500,100)

>>> print("{:.5f}".format(rk.linearly_stable_step_size(rk44,L1,plot=0)))
0.00007
>>> print("{:.5f}".format(rk.linearly_stable_step_size(rk44,L2,plot=0)))
0.02423

>>> sd = semidisc.load_semidisc('spectral difference advection',order=1)
>>> print("{:.5f}".format(rk.linearly_stable_step_size(rk44,sd.L,plot=0)))
0.01393