Butcher parameters of Runge-Kutta schemes
This functionality does not run in MATLAB.
numeric::butcher(EULER1 | RKF43 | xRKF43 | RK4 | RKF34 | xRKF34 | RKF54a | xRKF54a | RKF54b | xRKF54b | DOPRI54 | xDOPRI54 | CK54 | xCK54 | RKF45a | xRKF45a | RKF45b | xRKF45b | DOPRI45 | xDOPRI45 | CK45 | xCK45 | DOPRI65 | xDOPRI65 | DOPRI56 | xDOPRI56 | BUTCHER6 | RKF87 | xRKF87 | RKF78 | xRKF78 | DOPRI87 | xDOPRI87 | DOPRI78 | xDOPRI78 | GAUSS(s) | GAUSS = s, <digits>)
numeric::butcher(method) returns the Butcher parameters of the Runge-Kutta scheme named method.
An s-stage Runge-Kutta method for the numerical integration of a dynamical system with step size h is a map
The "intermediate stages" k1, …, ks are defined as the solutions of the algebraic equations
If the s×s "Butcher matrix" aij is strictly lower triangular, the method is called "explicit". In this case, the intermediate stages are computed recursively as:
Various numerical schemes arise from different choices of the Butcher parameters: the s×s-matrix aij, the weights b = [b1, …, bs] and the abscissae c = [c1, c2, …, cs].
Embedded pairs of Runge-Kutta methods consist of two methods that share the matrix aij and the abscissae ci, but use different weights bi.
The returned list [s, c, a, b1, b2, order1, order2] contains the Butcher data of the method: s is the number of stages, c is the list of abscissae, a is the Butcher matrix, b1 and b2 are lists of weights. The integers order1 and order2 are the orders of the scheme when using the weights b1 or b2, respectively, in conjunction with the matrix a and the abscissae c.
The methods EULER1 (order 1), RK4 (order 4) and BUTCHER6 (order 6) are single methods with b1 = b2 and order1 = order2. All other methods are embedded pairs of Runge-Kutta-Fehlberg (RKFxx), Dormand-Prince (DOPRIxx) or Cash-Karp (CKxx) type. The names indicate the orders of the subprocesses, e.g., CK45 is the Cash-Karp pair of orders 4 and 5. CK54 is the same pair with reversed ordering of the subprocesses. The second subprocess is used to produce a time step of the Runge-Kutta method, the first subprocess is used for estimating the error of the time step.
The methods GAUSS(s) or, equivalently, GAUSS = s are the implicit Gauss methods with s stages of order 2 s.
The data of all explicit methods are returned as exact rational numbers. The data of the Gauss methods are returned as floating-point numbers.
When computing the data for GAUSS(s), the function is sensitive to the environment variable DIGITS, which determines the numerical working precision.
The Butcher data of the classical 4 stage, 4th order Runge-Kutta scheme are:
Note that the weights b1 and b2 coincide: this classical method does not provide an embedded pair.
The Butcher data of the (implicit) 3 stage Gauss method:
DIGITS := 5: numeric::butcher(GAUSS(3)); delete DIGITS:
The Butcher data of the embedded Runge-Kutta-Fehlberg pair RKF34 of orders 3 and 4 are:
[s, c, a, b1, b2, order1, order2] := numeric::butcher(RKF34):
The number of stages s of the 4th order subprocess is 5, the abscissae c and the matrix a are given by:
s, c, a
Using these parameters with the weights
yields a numerical scheme of order 3 or 4, respectively:
delete s, c, a, b1, b2, order1, order2:
We plot the stability regions of the two sub-methods of DOPRI78. The stability function of a Runge-Kutta scheme with Butcher parameters (c, a, b) is given by
where e is the column vector (1, 1, …, 1)T. For an explicit s-stage scheme (the matrix a is strictly lower triangular), this stability function reduces to the polynomial
We compute the coefficients of the stability polynomials associated with the Butcher matrix a and the weights b1 and b2 of the sub-methods of DOPRI78:
[s, c, a, b1, b2, order1, order2] := numeric::butcher(DOPRI78): e := matrix([1 $ s]): a := float(matrix(a)): b1 := linalg::transpose(float(matrix(b1))): b2 := linalg::transpose(float(matrix(b2))): for i from 1 to s do c1[i] := (b1*a^(i-1)*e)[1, 1]; c2[i] := (b2*a^(i-1)*e)[1, 1]; end_for:
We define the stability polynomials:
z := x + I*y: p1 := 1 + _plus(c1[i]*z^i $ i = 1..s): p2 := 1 + _plus(c2[i]*z^i $ i = 1..s):
The boundary of the stability region is the curve defined by |p(z)| = 1. We plot these implicit curves associated with the stability polynomials p1(z) and p2(z) defined above:
plot(plot::Implicit2d(abs(p1) = 1, x = -6..1, y = 0..6, Color = RGB::Red, Legend = "Order 7"), plot::Implicit2d(abs(p2) = 1, x = -6..1, y = 0..6, Color = RGB::Blue, Legend = "Order 8"), Scaling = Constrained):
delete s, c, a, b1, b2, order1, order2, e, c1, c2, z, p1, p2:
The number of stages of the Gauss method: a positive integer
The number of significant digits with which the Butcher data of the methods GAUSS(s) are computed. The default value for digits is the current value of the environment variable DIGITS. This argument is only relevant for the Gauss methods.
The Butcher parameters provided in this original paper consist of rational approximations of solutions of the order equations of Runge-Kutta systems. The parameters provided by numeric::butcher are exact rational solutions of the order equations. The approximations given by Prince and Dormand coincide with the MuPAD® exact values through 16 decimal digits.
J.C. Butcher: The Numerical Analysis of Ordinary Differential Equations, Wiley, Chichester (1987).
E. Hairer, S.P. Norsett and G. Wanner: Solving Ordinary Differential Equations I, Springer, Berlin (1993).
The methods DOPRI87 and DOPRI78 correspond to the method RK8(7)13M published in:
P.J. Prince and J.R.Dormand: High order embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics 7(1), 1981.