The first eighth-order Runge--Kutta methods

Jim Verner

In October, 1966, I entered a Ph.D. program at the University of London Institute of Computer Science intending to study methods for the numerical solution of differential equations. During the early months, I attended classes in elementary numerical analysis, learned about solving linear equations, and encountered a view by Leslie Fox [L. Fox, 1964], that experience "indicates that linear algebra is involved, wholely or in part, in about 75 percent of all scientific problems". This linearity formed the basis of many approaches I have since used to derive new Runge–Kutta and related methods for solving initial value problems in ordinary differential equations.

I also learned about methods for other scientific problems, and later I began studies of recent papers on numerical ordinary differential equations under the tutelage of Graeme Cooper. Graeme encouraged me to study John Butcher's paper on the convergence of A-B methods [J.C. Butcher, 1966], for such problems and then to give a seminar based on these studies. Next year, he suggested I read some of John's papers on Runge–Kutta methods [J.C. Butcher, 1963], [J.C. Butcher, 1964] as well as some of his own [G.J. Cooper, 1966], [G.J. Cooper, 1967], based on quadrature methods.

During the following months, Graeme accepted a position at the University of Edinburgh, and invited me to continue my studies with him. I accepted his invitation, and continued working through the various known approaches to multistep and Runge–Kutta methods, and working on a paper of my own [J.H. Verner, 1969], during the end of 1967 and into 1968.

In the middle of 1968, Graeme announced that he was departing for a month's holiday, and asked if I could continue studying a design he had for some explicit 11-stage Runge–Kutta methods based on Lobatto quadrature rules of order eight. For such a method to have order eight for initial value problems, there were 55 unknown coefficients, and 11 weights that needed to be chosen to satisfy a set of 200 polynomial order conditions in these unknowns. In the early 1960's, John had explored these polynomial order conditions in general, and had shown that these equations could be conveniently sequenced using bijective mappings of the 200 rooted trees having up to eight nodes.

Graeme based his design on two lower order methods found by John for which the nodes and weights were those of Lobatto quadrature rules of orders 4 (namely $\{0, \frac12,1\}$), and 6 (namely $\{0, \frac{5-\sqrt{5}}{10}, \frac{5+\sqrt{5}}{10},1\}$) respectively. In addition, certain coefficients and weights were to be selected equal to 0, and the nodes were subdivided in subsets with coefficients constrained to satisfy increasing 'internal order' conditions. This design allowed reduction of the larger system to a smaller set of 46 polynomial equations of which any selection of coefficients satisfying these constraints satisfied 40 of these equations. There remained some flexibility of several arbitrary parameters, and as well, the 11 abscissae and an appropriate sequencing thereof had yet to be selected from the five Lobatto quadrature nodes.

Even though this reduced system of polynomial equations had fewer unknowns than equations, John had already shown that the 37 equations for a 7-stage Runge–Kutta method of order 6 could be solved using 28 parameters The first sixth order methods [J.C. Butcher, 1964], and the 85 equations for a 9-stage method of order 7 could be solved using 45 parameters. Indeed, these systems of order equations led a charmed existence: there were strong internal dependencies among them that admitted their solution with fewer variables than equations. In his derivation of methods, John had shown that there existed some 'simplifying conditions' that reduced the form and number of polynomial equations that had to be satisfied. Thus, Graeme's hope that this system of 200 equations in 66 variables had a solution might be realizable.

We turn now to describe the derivation of the 11-stage Runge–Kutta methods of order eight. In this Lobatto quadrature design, the abscissae were to be confined to choices from the five values $$ \left\{\begin{array}{ccccc} 0, & \frac{7-\sqrt{21}}{14}, & \frac12,& \frac{7+\sqrt{21}}{14}, & 1 \end{array} \right\}; $$ the vector of weights were confined to corresponding weights $\{\frac{1}{20}, \frac{49}{180}, \frac{16}{45}, \frac{49}{180}, \frac{1}{20}\}$ of this quadrature rule or else zero. Early studies suggested that the choices $c_1=0, c_{11}=1$, and the remaining nine internal nodes of the method be selected from the three internal nodes of the Lobatto quadrature. However, the sequencing of an appropriate choice was not known. (In a different design, Alan Curtis [A.R. Curtis, 1970], found that two abscissae corresponding to weights equal to zero and two other non-zero weights could be chosen almost arbitrarily to yield a parametric family of solutions.) In addition, the patterns from the methods of orders 4 and 6 suggested a choice of 15 coefficients being zero: \begin{align*} a_{i,2}&=0,& i&=5,6,\dots,11, \\ a_{i,j}&=0,& i&=8,9,10,11,\quad j=3,4. \end{align*} This would leave $55-15=40$ elements of the lower triangular matrix $A$ yet to be selected. As well, each of the weights $\{b_2,b_3,\dots,b_7\}$ were also set equal to zero so that any particular choice of all the abscissae would uniquely determine all weights of the method.

The plan was then to restrict the coefficients of $A$ remaining so that the internal stages had increasing stage order up to order 4, and then to use the remaining freedom among the coefficients and of the selection of the internal abscissae to satisfy the remaining order equations. These 30 choices to satisfy the 'internal stage orders' and the constraints on the weights reduced the system to one of 46 polynomial equations in the 40 unknowns. These 40 variables could be selected in various ways to satisfy at least 40 of the 46 equations in the reduced system; from these solutions we wanted to isolate some solutions that satisfied more equations.

With this model, Graeme established that were these assumptions used together with yet unknown choices to satisfy the 46 equations, then all 200 of the required order conditions would be satisfied, and the method obtained would have order 8.

After describing this scenario to me in detail, Graeme departed, encouraging me to study this reduced system, and expressing a hope that I would be able to determine additional contraints on the coefficients that would further reduce this system of 46 algebraic equations to a smaller number that could be solved for the 40 unknown coefficients of the method.

Here are the assumptions that were to be imposed on the coefficients (allowing each stage to have a specified 'internal order').
The set of eleven nodes is partitioned sequentially into five subsets as $\bigcup_{i=0}^4 C_i$, where \begin{align*} C_0 &= \{c_0=1\},\\ \quad C_i &= \{c_j: 2+i(i-1)/2 \le j \le 1+i(i+1)/2\}. \end{align*} For each stage $j$ in subset $C_i$, $i$ constraints were imposed as follows: for $q_{j,k}=\sum_{l=1}^{j-1}a_{jl}c_l^{k-1}-c_j^k/k$, we require \begin{align*} C_1:& &q_{j,k}&=0, & j&=2, & k&=1,\\ C_2:& &q_{j,k}&=0, & j&=3,4,& k&=1,2,\\ C_3:& &q_{j,k}&=0, & j&=5,6,7,& k&=1,2,3,\\ C_4:& &q_{j,k}&=0, & j&=8,9,10,&\quad k&=1,2,3,4. \end{align*} While we could also have required that $q_{11,k}=0, \ k=1,..,4,$ it turns out that requiring $$ b_j(1-c_j)=\sum_{i=j+1}^{11}b_ia_{ij},\quad j=1,5,..,10, $$ provides 7 linear constraints on the non-zero coefficients $a_{11,j},\ j=1,5,..,10,$ a choice that leads to $\{q_{11,k}=0,\ k=1,..,4\},$ and also simplifies other order conditions.

It will be observed that for requiring subsets of distinct abscissae, these 33 constraints can be satisfied by solving subsets of linear equations after particular choices are made for the abscissae. In addition to these constraints, there remain 13 additional algebraic equations to solve in order to achieve order 8, yet to satisfy these, there remain only 7 arbitrary choices of the 40 non-zero coefficients. The plan was to find some dependence among these additional algebraic equations, and show that all 13 could be satisfied by solving some equations that sequentially become linear. Each of these found equations might be interleaved with those already solved to form fully determinable subsets of linear equations.

It is evident that some subsets of the 13 additional equations could be satisfied independently. I ventured to extract subsets that looked appropriate, and wrote code to solve these on the assumption that some results would lead me further. Simultaneously, I tried to reformulate these remaining equations algebraically, again on a quest for further dependencies. Just before Graeme was to return, I had found that 5 of the equations could be eliminated leaving only 8 (rather than 13) additional equations, or a total of 41 constraints on the 40 variables. By this point, I had subsets of linear equations that could be solved sequentially to yield values of all 40 variable coefficients for any distribution of the three internal Lobatto nodes among the nine internal nodes of the method. Accordingly, I wrote an Algol program to solve for the 40 nonzero coefficients, and selected a few sequences of the nodes to determine values of the remaining constraint. Unfortunately, no sequence of the nodes that I chose led to the hoped-for result. (Because each test required overnight computing, there was insufficient time to test all possible sequencings of the nodes before Graeme was to return.)

I explained to Graeme in detail what I had discovered. He was pleased with the fact that there had been some progress, and suggested that I return the next day after he had reviewed my results in detail. I had told him all of the possible distributions of the nodes, and of the few selections I had made and tested although none of the particular sequences of nodes I had selected for testing satisfied all 200 order conditions.

The next day, he observed to me that of the trial choices I had made, each pair of adjacent nodes was distinct, and suggested I resequence the nodes so that up to three adjacent pairs were equal. He suggested that, should some such choice lead to the anticipated result, he would invite me to the Club to share a pint. It is very likely you can guess what happened next! Indeed - that pint showed that me that the completion of my program was in sight. Graeme's suggestion did lead to a method of order 8. Thereafter, my recollection is that we had found twelve methods, and current computations show this probably happened. Indeed, for each method obtained, another method is obtained by exchanging $\sqrt{21}$ with its negative value, and a complete search of other permutations of the internal nodes leads to a total of exactly 72 distinct methods with this design. These results first appear in [James H. Verner, 1969b], and later within a more general derivation of high order explicit Runge–Kutta methods [G.J. Cooper and J.H. Verner, 1972].

A complete method may be expressed in a Butcher tableau consisting of a strictly lower triangular $11 \times 11$ matrix together with a vector of 11 weights. The abscissae are then obtained as the row sums of the matrix. Here, we record equations that exactly specify values of the coefficients by rows of the strict lower triangle of the $11 \times 11$ matrix $a$ to be found in a Butcher tableau. For each $i, i=2..11,$ there are $i-1$ coefficients determined.

Solution to Row 2 (Internal order is 1) $$a_{2,1}=c_2.$$

Solution to Row 3 (Internal order is 2) \begin{align*} a_{3,2}&=\frac{c_3^2}{2c_2},\\ a_{3,1}&=c_3-a_{3,2}. \end{align*} To prescribe a constant ratio of $q_{i,3}$ to $a_{i,2}c_2^2,\ i=3,4,$ define $$A_{31}=\frac{a_{3,2}c_2^2 -\frac13c_3^3}{a_{3,2}c_2^2}. $$

Equations for Row 4 (Internal order is 2) \begin{alignat*}2 q_{4,k}&=0, &\enspace k=1,2,\\ q_{4,3}&=A_{31}a_{4,2}c_2^2.&\\ \end{alignat*}

Equations for Row 5 (Internal order is 3) \begin{align*} a_{5,2}&=0, &\\ q_{5,k}&=0, &k=1,2,3.\\ \end{align*} To prescribe a constant ratio of $q_{i,4}$ to $\Sigma_{j=3}^{i-1} a_{i,j}a_{j,2}c_2^2,\ i=5,6,7,\ $ define $$A_{51}=\displaystyle{\frac{q_{5,4}}{(a_{5,3}a_{3,2}+a_{5,4}a_{4,2})c_2^2}}.$$

Equations for Row 6 (Internal order is 3) \begin{alignat*}2 a_{6,2}&=0, &\\ q_{6,k}&=0, &\enspace k=1,2,3,\\ q_{6,4}&=A_{51}(a_{6,3}a_{3,2} + a_{6,4}a_{4,2})c_2^2. \end{alignat*}

Equations for Row 7 (Internal order is 3) \begin{alignat*}2 a_{7,2}&=0, &\\ q_{7,k}&=0, &\enspace k=1,2,3,\\ q_{7,4}&=A_{51}(a_{7,3}a_{3,2} + a_{7,4}a_{4,2})c_2^2,&\\ a_{7,3}\big(a_{5,4}(c_6-c_7)(a_{6,3}a_{3,2} + a_{6,4}a_{4,2})&\\ \qquad-a_{6,4}(c_5-c_7)(a_{5,3}a_{3,2} + a_{5,4}a_{4,2})\big)&\\ \enspace -a_{7,4}\big(a_{5,3}(c_6-c_7)(a_{6,3}a_{3,2} + a_{6,4}a_{4,2})&\\ \qquad-a_{6,3}(c_5-c_7)(a_{5,3}a_{3,2} + a_{5,4}a_{4,2})\big)&=0. \end{alignat*}

Equations for Row 8 (Order is 4) \begin{align*} a_{8,j}&=0, & j&=2,3,4,\\ q_{8,k}&=0, & k&=1,2,3,4. \end{align*}

Equations for Row 9 (Order is 4): assuming $c_8=c_7$: \begin{alignat*}2 a_{9,j}&=0, & j&=2,3,4,\\ q_{9,k}&=0, & k&=1,2,3,4,\\ b_9(c_9-1)(c_9-c_{10})a_{9,8}a_{8,7} &=b_8\big((3-4c_8+c_8^4)/12 &\\ &\quad \quad -\ (1+c_{10})(2-3c_8 +c_8^3)/6 + c_{10}(1-2c_8+c_8^2)/2\big). \end{alignat*}

Equations for Row 10 (Order is 4) \begin{alignat*}2 a_{10,j}&=0, & j&=2,3,4,\\ q_{10,k}&=0, & k&=1,2,3,4,\\ b_{9}(c_9-1)a_{9,8}a_{8,7}+ b_{10}(c_{10}-1)(a_{10,9}(a_{9,7}+a_{9,8}) +a_{10,8}a_{8,7})&=b_8(c_8-1)^3/6,&\\ b_{10}(c_{10}-1)a_{10,9}(c_8-c_9)(a_{9,7}+a_{9,8})&=b_{8}(c_8-1)^4/24. \end{alignat*}

Solution to Row 11 (Order is 4) \begin{alignat*}2 a_{11,j}&=0, &\quad j&=2,3,4,\\ a_{11,j}&=\frac{b_{j}(1-c_j)-\sum_{i=j+1}^{10}b_{i}a_{i,j}}{b_{11}}, &\quad j&=1,5,..,10. \end{alignat*}

If these equations are solved sequentially as given, each row $i$ contains a set of $i-1$ equations that are linear in $a_{i,j},\ j=1..i-1,$ and can be solved uniquely for appropriate choices of the nodes. Each set $C_i, \ i=1,2,3,4,$ must have distinct non-zero values. Nodes $c_1=0,\ c_{11}=1$, and the remaining nodes must be selected from the three 'internal nodes' of the Lobatto quadrature rule. Node $c_2$ can be any of the three internal Lobatto nodes, $C_2\subset C_3$, and $c_7=c_8$ cannot be in $C_2$. There are 72 sequences of the nodes that satisfy these choices with $c_8$ being any of the three internal nodes. If the coefficients of an 11-stage method are obtained from the equations above, the method obtained will have order 8, and the coefficients of the method will exactly satisfy the 200 order equations.

Simultaneously, yet unknown to this author, a four parameter family of (11-8) Runge–Kutta methods were derived by Alan Curtis [A.R. Curtis, 1970], a researcher at the National Physical Laboratory in the United Kingdom. On investigation, it may be found that these two sets of methods are distinct, but possibly may have a common intersection. By comparing these two sets of methods, we found that more unknown methods of this type exist. For example, it turns out that in the Cooper-Verner family, the node $c_2$ can be chosen almost arbitrarily to give a one-parameter family of methods.