Documentation |
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" k_{1}, …, k_{s} are defined as the solutions of the algebraic equations
.
If the s×s "Butcher matrix" a_{ij} 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 a_{ij}, the weights b = [b_{1}, …, b_{s}] and the abscissae c = [c_{1}, c_{2}, …, c_{s}].
Embedded pairs of Runge-Kutta methods consist of two methods that share the matrix a_{ij} and the abscissae c_{i}, but use different weights b_{i}.
The returned list [s, c, a, b_{1}, b_{2}, order_{1}, order_{2}] contains the Butcher data of the method: s is the number of stages, c is the list of abscissae, a is the Butcher matrix, b_{1} and b_{2} are lists of weights. The integers order_{1} and order_{2} are the orders of the scheme when using the weights b_{1} or b_{2}, 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 b_{1} = b_{2} and order_{1} = order_{2}. 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.
The Butcher data are called by the routines numeric::odesolve, numeric::odesolve2, and numeric::odesolveGeometric.
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:
numeric::butcher(RK4)
Note that the weights b_{1} and b_{2} 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
b1, b2
yields a numerical scheme of order 3 or 4, respectively:
order1, order2
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 b_{1} and b_{2} 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 p_{1}(z) and p_{2}(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:
s |
The number of stages of the Gauss method: a positive integer |
digits |
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.