Add a Runge-Kutta scheme

Current \(Q\)-generators based on Runge-Kutta schemes are implemented in the qmat.qcoeff.butcher submodule. Those are based on Butcher tables from classical schemes available in the literature, and the selected approach is to define one class for one scheme.

Standard scheme

In order to add a new RK, search first for its section in the butcher.py file, depending on its type (explicit or implicit) and its order. Then add a new class at the bottom of this section following this template :

@registerRK
class NewRK(RK):
    """Some new RK method from ..."""
    A = ... # TODO
    b = ... # TODO
    c = ... # TODO

    @property
    def order(self): return ... # TODO

Here the registerRK decorators interfaces the classical register decorator for QGenerator classes, but also :

check if the dimensions of the A, b and c are consistent
register the generator in a specific category with all RK-type generators

💡 You can use either the built-in list or Numpy nd.array to add the class attributes A, b and c.

Tip : for large Butcher table, you can also use this approach (from the CashKarp class) :

A = np.zeros((6, 6))
A[1, 0] = 1.0 / 5.0
A[2, :2] = [3.0 / 40.0, 9.0 / 40.0]
A[3, :3] = [0.3, -0.9, 1.2]
A[4, :4] = [-11.0 / 54.0, 5.0 / 2.0, -70.0 / 27.0, 35.0 / 27.0]
A[5, :5] = [1631.0 / 55296.0, 175.0 / 512.0, 575.0 / 13824.0, 44275.0 / 110592.0, 253.0 / 4096.0]

Convergence testing

To test your scheme … you don’t have to do anything 🥳 : all RK schemes are automatically tested thanks to the registration mechanism, that checks (in particular) the convergence order of each scheme (global truncation error).

⚠️ Depending on the implemented RK scheme, convergence test may fail … in that case no worries 😉 you’ll just have to adapt your scheme to the test, as explained below :

All convergence tests are done on the following Dahlquist problem :

u0 = 1        # unitary initial solution
lam = 1j      # purely imaginary lambda
T = 2*np.pi   # one time period

They use three numbers of time-steps for the convergence analysis, depending on the order of the method (see here …).

But this automatic time-step size selection may not be adapted for methods with high error constant that require finer time-steps to actually see the theoretical order. In that case, simply add a CONV_TEST_NSTEPS class attribute storing a list with higher numbers of time-steps in increasing order, high enough so the convergence test passes.

📜 See SDIRK2_2 implementation for an example …

Embedded scheme

Butcher tables with additional coefficients for embedded method can also be added to new or existing RK schemes. For that, simply define a b2 class attribute :

@registerRK
class NewRK(RK):
    """Some new RK method from ..."""
    ## previous coefficients ... 
    b2 = ... # embedded coefficients

Per default, \(Q\)-generators define the order of the embedded method (using those additional coefficient) as one order less than the method’s order (that is, returned by the order property). If this is not the case, then you should override the weightEmbedded property from the base class :

@registerRK
class NewRK(RK):
    # ...
    @property
    def weightsEmbedded(self):
        return ...  # effective embedded order