1 Introduction↩︎

Runge Kutta (RK) methods for the numerical solution of IVPs in systems of differential equations \[\label{A1} \dfrac{d y(t)}{dt}= f(t, y(t)), \quad y(0)= y_0 \in \mathbb{R}^d, \quad t \in [0,T]\tag{1}\] are among the most popular algorithms for the numerical solution of (1 ). In particular explicit \(s\)–stage RK methods advance the numerical solution of 1 from \((t_n=n h, y_n)\) to \((t_{n+1}= t_n+h, y_{n+1})\) with the formulas \[\label{A2} y_{n+1} = y_n + h \; \left[ b_1 \; f_{n,1} + b_2 \; f_{n,2}+ \ldots + b_n \; f_{n,s} \right],\tag{2}\] where \[\begin{align} f_{n,1} &=& f( t_n , y_n) \nonumber \\ f_{n,2} &=& f( t_n + c_2 h , y_n + h \; a_{21} f_{n,1}) \nonumber \\ \vdots & \vdots \label{A3} \\ f_{n,s} &=& f( t_n + c_s h, y_n + h a_{s,1} f_{n,1}+ \ldots +h a_{s,s-1} f_{n, s-1} ). \nonumber \end{align}\tag{3}\] It is common to collect the available parameters of (2 ),(3 ) , \(b_j, c_j, a_{jk}\) in the form of the Butcher array [1], [2]

\[\label{A4} \begin{array}{cccccc} c_1=0 & \vline & 0 & & & \\ c_2 & \vline & a_{21} & 0 & & \\ \vdots & \vline & \vdots & & & \\ \vdots & \vline & \vdots & & & \\ c_s & \vline & a_{s1} & \ldots & a_{s,s-1} & 0 \\ \hline & \vline & b_1 & b_2 & \ldots & b_s \end{array} \qquad \begin{array}{ccc} c & \vline & A \\ \hline 0 & \vline & b^T \end{array}\tag{4}\] with the condition \[c_i = \sum_{j=1}^{i-1} \; a_{ij} \qquad i=1,..,s\]

Usually the free parameters of (4 ) are selected in such a way that the method achieves the highest accuracy possible for all IVP (1 ) with \(f(t,y)\) sufficiently smooth. Denoting by \(y(t; t_n,y_n)\) the local solution of differential problem \(y'(t) = f(t, y(t))\) that passes by the point \((t_n, y_n)\), that is, \(y(t_n)=y_n\), the accuracy is measured by the local error, defined as the error after one step \(\Vert y_{n+1} - y(t_n+h; t_n,y_n) \Vert\). Then the method has order \(p\) [2] if \[\label{A5} \Vert y_{n+1} - y(t_n+h; t_n,y_n) \Vert= \mathcal{O} (h^{p+1}), \quad h \to 0\tag{5}\] for all IVP (1 ). This condition implies that the global error of (1 ) behaves as \(\Vert y(T) - y_N \Vert = \mathcal{O} (h^{p})\) with \(N h = T\). In this setting it is well known that for \(s \le 4\) it is possible to get order \(p=s\) whereas for \(s>4\) the maximum order \(p < s\).

Clearly the computational cost of an explicit RK method is associated to the number of function evaluations at each step i.e. the number of stages. With the aim to reduce this cost we examine here explicit RK methods with \(c_s=1\) in which the first evaluation of the \(n\)-th step, \(f_{n,1}=f(t_n,y_n)\), is substituted by the last evaluation of the \((n-1)\)-th step, \(f_{n-1,s}\) i.e. \(f_{n,1}=f( t_{n-1}+ c_s \; h, y_{n-1}+h\sum_{j=1}^{s-1} a_{s,j} f_{n-1,j})\). These methods will be called explicit RK methods reusing last stage, and we will denote them by RKRLS. Thus instead of (2 ),(3 ) we have the method

\[\label{A6} y_{n+1} = y_n + h \; \left[ b_1 \; f_{n,1} + b_2 \; f_{n,2}+ \ldots + b_n \; f_{n,s} \right],\tag{6}\] with \[\begin{align} f_{n,1} &=& f_{n-1,s} = f( t_{n-1}+ c_s \; h, y_{n-1} + h \sum_{j=1}^{s-1} a_{s,j} f_{n-1,j}) \nonumber \\ f_{n,2} &=& f( t_n + c_2 h , y_n + h \; a_{21} f_{n,1}) \nonumber \\ \vdots & \vdots \label{A7} \\ f_{n,s} &=& f( t_n + c_s h, y_n + h a_{s,1} f_{n,1}+ \ldots +h a_{s,s-1} f_{n, s-1} ). \nonumber \end{align}\tag{7}\]

In this setting (except for the first step) we only have \((s-1)\) function evaluations per step. However, even in the case that the explicit RK method (2 ),(3 ) has order \(p\), the associated RKRLS method (6 )(7 ) will not have in general the same convergence order.

The method (6 )(7 ) advances the numerical solution from \(t_n\) to \(t_{n+1}=t_n+h\) by considering at \(t_n\) two approximations to the solution: \(y_n \simeq y(t_n)\) and \(f_{n,1}\equiv f_{n-1,s}= f(t_{n-1} +h, y_{n-1}+ h \sum_{j=1}^{s-1} a_{s,j} f_{n-1,j}) \simeq y'(t_n)\), obtaining with them the next two approximations \(y_{n+1}\) and \(f_{n+1,1}\). The local error cannot be defined as in standard one step methods based only in the approximation \(y_n\), but a two-approximations, one step formulation can be used to apply the theory in [2] (Section III.8, p. 436). Here for technical reasons, instead of \(f_{n-1,s}\) we will use the quantity \(y_{n,1}\) defined by \[y_{n,1} \equiv y_{n-1} + h \sum_{j=1}^{s-1} a_{s,j} f_{n-1,j}.\] Then the method (6 )(7 ) can be written as a two-approximations, one step method in the form \[\label{asonestep} \begin{pmatrix} y_{n+1} \\ y_{n+1,1} \end{pmatrix} = \begin{pmatrix} y_{n} + h \sum_{i=1}^s b_i f_{n,i} \\[6pt] y_{n} + h \sum_{i=1}^{s-1} a_{si} f_{n,i} \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} y_{n} \\ y_{n,1}\end{pmatrix} + h \Phi(t_n, y_n, y_{n,1},h)\tag{8}\] where \(f_{n,1}\equiv f(t_n, y_{n,1})= f_{n-1,s}\) depends on the second carried vector \(y_{n,1}\) and is in fact the last stage of the previous step. Since this value has been computed in the previous step, it does not need to be computed in this step.

The vector \(y_{n+1,1}\) is an approximation to the exact value \(y(t_{n} + h)\). Thus, the order of the method can be studied by means of standard Butcher series as in [2] (Section III.8, p. 436). Unfortunately, \(y_{n+1,1}\) usually approximates \(y(t_{n} + h)\) with order \(q<p\) and this theory will ensure that the method has order \(\min(q+1,p)\). To get order \(p\), \(y_{n+1,1}\) should have order \(p-1\) at least and this can be a very restrictive condition.

Two-step Runge-Kutta methods were considered for the first time in [3], where they were named pseudo-Runge-Kutta methods. They have been considered later, for example in [4] and [5]. In both cases the new approximation \(y_{n+1}\) can be expressed as \[\label{manystepsmethod} y_{n+1}= y_n + h \Phi(t_n,y_n, v_{1n}, v_{2n}, \ldots)\tag{9}\] where \(v_{jn}\) are approximations to the solution \(y(t)\) of 1 at certain previous points \(t_{n-1} + a_j h\), of order \(p\). In [6] and [7] the computation of \(y_{n+1}\) uses the stages of the previous step so that it can be expressed as 9 but now \(v_j\) are approximations to \(y'(t_{n-1} + a_j h)\) of order \(p-1\). In all these cases the theory in [2] gives order \(p\) for the resulting methods.

In [8] general two-steps Runge-Kutta methods were studied. Since the approach in [2] could not be directly applied, the authors turned to an approach due to Albrecht [9], [10] to derive order conditions for the method. Nevertheless, in [11] it was shown that the results in [2] can also be applied to general two-step Runge-Kutta methods formulated as 9 by assuming that \(v_{in}\) approximates certain Butcher series \(B(\Phi_i, h)\) evaluated at the solution \(y(t)\) at \(t_n\). Here, we will study the accuracy order of RKRLS methods using this last approach. Starting with a standard RK method order \(p\) we derive additional conditions on their coefficients that ensure a global error of the RKRLS of order \(p\) even though the approximation \(y_{n+1,1}\) has order \(q < p-1\).

The paper is organized as follows: In Section 2 we present the theoretical results on the order of convergence of Runge-Kutta methods reusing the last stage. We define the local order of the methods and we prove that local order \(p\) implies global convergence of order \(p\). Then we obtain conditions on the coefficients of the RK to have RKRLS method of order \(p\), for \(p=2,3,4,5\). In Section 3 we obtain RKRLS methods with orders from 2 to 5, showing that there exist methods with order \(p\) and \(s=p\) stages (\(p-1\) effective stages) for \(p \le 4\) and order \(5\) with \(s=6\) stages. In Section 4 we verify the theoretical results by means of some numerical experiments. Finally, in Section 5 we present some conclusions and future work.

2 Order and convergence↩︎

2.1 Order theory↩︎

For a Runge-Kutta method \((A,b)\) with \(s\) stages, \(n\) consecutive steps of size \(h\) are equivalent to one step of size \(H=nh\) with the \(n\times s\) stages composed method \((A^{(n)}, b^{(n)})\) with coefficients \[A^{(n)}=\dfrac{1}{n} \begin{pmatrix} A & \\ eb^T & A \\ \cdots & \cdots & \ddots \\ e b^T & \cdots & e b^T & A \end{pmatrix}, \qquad b^{(n)}= \dfrac{1}{n} \begin{pmatrix} b \\ b \\ \vdots \\ b \end{pmatrix},\] where \(e= (1, \ldots, 1)^T \in \mathbb{R}^s\). Note that if the RK base method \((A, b)\) has order 1, since \(Ae=c\), the composed method satisfies \[c^{(n)} = A^{(n)} e^{(n)}=\dfrac{1}{n} \begin{pmatrix} c \\ e+c \\ \vdots \\ (n-1)e +c \end{pmatrix}.\] where \(e^{(n)}= (1, \ldots, 1)^T \in \mathbb{R}^{ns}\). Here, a superindex \((n)\) represents vectors and matrices of dimension \(sn\).

It can be easily proved that the RK method \((A,b)\) has local order \(p\) if and only if the composed method \((A^{(n)}, b^{(n)})\) has also the same order for any \(n\ge2\). Then, the local error of \((A,b)\) behaves as \(\mathcal{O}(h^{p+1})\), and the local error of the n-step composed method \((A^{(n)}, b^{(n)})\) behaves as \(\mathcal{O}(H^{p+1})\).

For a RKRLS method, the first stage is substituted by the last one of the previous step, that is, at the step from \(t_n\) to \(t_{n+1}=t_n+h\), \(f_{n,1}= f_{n-1,s} = h f(y_{n-1} + ha_{s,1}f_{n-1,1} + \ldots + h a_{s,s-1} f_{n-1,s-1})\). The first two steps of the method are equivalent to one step of the \(2s\) stage method with coefficients \[\widehat A^{(2)}=\dfrac{1}{2} \begin{pmatrix} 0 & & & & & | &\\ a_{21} & 0 & & & & |\\ \cdots & \cdots & \ddots &\ddots & & |\\ a_{s1} & a_{s2}& \cdots & a_{s,s-1} & 0 & | & 0 \\[3pt] \hline a_{s1} & a_{s2}& \cdots & a_{s,s-1} & 0 & | & 0 & 0\\ b_1 & b_2 & \cdots & b_{s-1} & b_s & | & a_{21} &0 \\ \cdots & \cdots & \ddots &\ddots & & | \\ b_1 & b_2 & \cdots & b_{s-1} & b_s & | & a_{s1} & \cdots & a_{s,s-1} & 0 \\ \end{pmatrix}, \; \widehat b^{(2)}= \dfrac{1}{2} \begin{pmatrix} b \\ b \\ \end{pmatrix}.\] that we will name 2-step composed method. Observe that the rows \(s\) and \(s+1\) of the matrix \(\widehat A^{(2)}\) are identical, which means that the method is reducible and the number of function evaluations required is in practice \(2s-1\), in agreement with the saving of one evaluation at each step.

In general, \(n\) consecutive steps of size \(h\) of the RKRLS method are equivalent to one step of size \(H=nh\) with the corresponding \(n\)-step composed method \((\widehat A^{(n)}, \widehat b^{(n)})\) given by \[\label{nreused} \widehat A^{(n)}=\dfrac{1}{n} \begin{pmatrix} A & \\ E & A \\ e b^T & E & A \\ \cdots & \cdots & \ddots &\ddots \\ e b^T & \cdots & e b^T & E & A \end{pmatrix}, \; \widehat b^{(n)}= \dfrac{1}{n} \begin{pmatrix} b \\ b \\ \vdots \\ b \end{pmatrix},\tag{10}\] where \(E=e b^T+M\), \(M=e_1 (a-b)^T\), \(a=( a_{s,1}, \ldots, a_{s,s-1}, 0)^T\) is the last row of the matrix \(A\) and \(e_1= (1, 0, \ldots, 0)^T\in \mathbb{R}^s\). Clearly, since \(a^T e=1\), \(E e=e\) and \(\widehat c^{(n)} = \widehat A^{(n)} e^{(n)}= c^{(n)}\).

The standard definition of local error of order \(p\) can not be directly applied to RKRLS methods because it does not allow us to obtain local order \(p\). Then, using the composed methods we define the local error of these methods in the following way

Let us see now that, with the above definition, if a RKRLS method has local order \(p\), then it is convergent of order \(p\), that is, the global error satisfies \(y_n - y(t_n)= \mathcal{O}(h^p)\).

Hereafter we will assume, without loss of generality, that the differential problem is autonomous, that is \(f=f(y)\).

Assume that the method has order \(p\) and let us consider the composed method with \(r\) steps (we will use the case \(r=p-1\) later). Starting from the initial point \((t_0, y_0)=(t_0, y(t_0))\) we get \(y_{r}\), approximation of order \(p\) to \(y(t_0+rh)\) and \(y_{r,1}=y_{r-1} + h a_{s1}f_{r-1,s} + h a_{s2} f_{r,2} + \ldots + h a_{s,s-1} f_{r,s-1}\) that approximates \(y(t_0+rh)\) in general with order \(1\). This means that \(y_r- y(t_0+rh) = \mathcal{O}(h^{p+1})\) and \(y_{r,1} - y(t_0+rh) = \mathcal{O}(h^2)\). Using that the approximation \(y_{r,1}\) is in fact the approximation provided by the method \((A^{(r)}, (b^T,\ldots,b^T, a^T)^T)\) we can obtain easily the coefficients of the Butcher series of \(y_{2,1}\) around the initial point \((t_0, y_0)\) given by \[y_{r,1} = y_0 + \dfrac{1}{r}(b^T e + a^T e) f(y_0) (rh)+ \dfrac{1}{2!} 2\dfrac{1}{r^2} (b^Tc +a^T (e+c)) f'(y_0) f(y_0) (rh)^2 + \cdots\] To study the effect of this approximation in the next step, it is more convenient to have the expansion of \(y_{r,1}\) around the current point \((t_r,y_r)\). This can be obtained by means of the adjoint of the method \((A^{(r)}, b^{(r)})\) (see [2], p. 219). This method is given by \[\label{adjointmethod} \begin{array}{l} \bar A^{(r)} = e^{(r)} {b^{(r)}}^T - A^{(r)}= \dfrac{1}{r} \begin{pmatrix} e b^T -A & e b^T & & \cdots & e b^T \\[5pt] - M & e b^T -A & & \cdots& e b^T\\[5pt] 0 & - M & \ddots & & e b^T\\ \vdots & \vdots & \ddots & \ddots \\ 0 & 0 & \cdots & - M & e b^T -A \end{pmatrix} \\[45pt] \bar c^{(r)} = \bar A^{(r)} e^{(r)} = \dfrac{1}{r} \begin{pmatrix} re-c \\ \vdots \\ e-c\\ \end{pmatrix}, \end{array}\tag{11}\] and the approximation \(y_{r,1}\) (it is in fact the last stage of the adjoint method) is given by the method \((\bar A^{(r)}, \bar b^{(r)})\) where \[\label{badjointmethod} \bar b^{(r)} \equiv ((e_s^{(r)})^T \bar A^{(r)})^T = \dfrac{1}{r} \begin{pmatrix} 0 \\ \vdots \\ b -a\\ \end{pmatrix},\tag{12}\] applied with step \(-H=-rh\). Computing the Butcher series for this method we get \[\label{seriesyr1} \begin{array}{rl} y_{r,1}= & y_{r} + \sum_{\tau} \dfrac{\alpha(\tau)}{\rho(\tau)!} \gamma(\tau) \phi_r(\tau) F(\tau)(y_{r}) (-r h)^{\rho(\tau)} \\[8pt] = & y_{r} + \sum_{\tau} \dfrac{\alpha(\tau)}{\rho(\tau)!} \gamma(\tau) \varphi_r(\tau) F(\tau)(y_{r}) h^{\rho(\tau)} \end{array}\tag{13}\] where (see e.g. [2], p. 145 or [1]) \(F(\tau)(y_{r})\) is the elementary differential associated to the tree \(\tau\) evaluated at the point \(y_{r}\), \(\alpha(\tau)\) is a symmetry factor, \(\gamma(\tau)\) is a density factor, \(\rho(\tau)\) is the order of the tree and \(\phi_r(\tau)= (-r)^{\rho(\tau)}\varphi_r(\tau)\) is the coefficient, method dependent, associated to \(\tau\). The value of \(\varphi_r(\tau)\) for the first trees \(\tau\) are given in Table ¿tbl:tableexpansion?. Note the minus signs in the coefficients of odd powers of \(h\) introduced by the negative step \(-H\).

The key point in this study is that the first coefficients of the Butcher series of \(y_{r,1}\) and \(y_{r+1,1}\) are equal.

Note that in general \(\varphi_r(\tau) \ne \varphi_{r+1}(\tau)\) for \(\rho(\tau) > r\). This is the reason why we will need to consider the composed methods with \(r=p-1\) and \(r+1=p\) to guarantee that both series have the same coefficients up to order \(p\), in order to prove the convergence of order \(p\).

Let us define the function \(z(t)\) as the Butcher series 13 of \(y_{p-1,1}\) associated to the adjoint method 11 with \(r=p-1\), evaluated at the point \(y(t)\) \[z(t)= y(t) + \sum_{\tau} \dfrac{\alpha(\tau) }{\rho(\tau)!} \gamma(\tau) \varphi_{p-1}(\tau) F(\tau)(y(t)) h^{\rho(\tau)}.\] We can state the following

We can now follow the approach in [2] (Section III.8, p. 436) for the one-step method 8 . Let us take \(y_n=y(t_n)\) and \(y_{n,1}= z(t_n)\). We will see that \(y_{n+1}-y(t_{n} +h) = \mathcal{O}(h^{p+1})\) and \(y_{n+1,1} - z(t_n +h) = \mathcal{O}(h^{p+1})\), that implies that the global error behaves as \(\mathcal{O}(h^p)\) and therefore the method is convergent with order \(p\).

First, note that by definition, \(z(t_n)\) is precisely the last stage of the adjoint method 11 with \(r=p-1\). Then, there is a starting value \(y_{n-p}\) such that the \((p-1)\)-step composed method gives \(y_n=y(t_n)\) and \(y_{n-1,s}=y_{n,1}=z(t_n)\). The next approximations \(y_{n+1}\) and \(y_{n+1,1}\) are computed with the next step, or, equivalently, with the \(p\)-step composed method starting at the point \((t_n-(p-1) h, y_{n-p+1})\). Let us denote \(y(t; t_n-(p-1)h, y_{n-p+1})\) the solution of the differential equation that at \(t_{n-p+1}= t_n -(p-1)h\) has a value \(y_{n-p+1}\) (see Figure [figura1])

Since the composed methods \((A^{(p-1)}, b^{(p-1)})\) and \((A^{(p)}, b^{(p)})\) have order \(p\), we can ensure that \[\begin{array}{l} y_n = y(t_n) = y(t_{n}; t_{n-p+1}, y_{n-p+1}) + \mathcal{O}(h^{p+1})\\[7pt] y_{n+1} = y(t_{n+1}; t_{n-p+1}, y_{n-p+1}) + \mathcal{O}(h^{p+1}). \end{array}\] Moreover, there exist a Lipschitz constant \(L\) such that \[\Vert y(t_{n+1}; t_{n-p+1}, y_{n-p+1})- y(t_{n+1}; t_n, y_n) \Vert \le L \Vert y_n - y(t_{n}; t_{n-p+1}, y_{n-p+1}) \Vert =\mathcal{O}(h^{p+1}).\] Therefore \[\begin{array}{l} y_{n+1}- y(t_{n+1})= y_{n+1} - y(t_{n+1}; t_n, y_n) \\[6pt] \quad =y_{n+1} - y(t_{n+1}; t_{n-p+1}, y_{n-p+1})+y(t_{n+1}; t_{n-p+1}, y_{n-p+1})-y(t_{n+1}) = \mathcal{O}(h^{p+1}) \end{array}\] On the other side, by definition of \(y_{n,1}=z(t_n)\) \[y_{n,1}= y_{n} + \sum_{\tau} \dfrac{1}{\rho(\tau)!} \gamma(\tau) \varphi_{p-1}(\tau) F(\tau)(y_{n+1}) h^{\rho(\tau)}\] and \[\begin{array}{l} y_{n+1,1}= y_{n+1} + \sum_{\rho(\tau)\le p} \dfrac{1}{\rho(\tau)!} \gamma(\tau) \varphi_{p}(\tau) F(\tau)(y_{n+1}) h^{\rho(\tau)} + \mathcal{O}(h^{p+1})\\ z(t_{n+1})= y(t_{n+1}) + \sum_{\rho(\tau)\le p} \dfrac{1}{\rho(\tau)!} \gamma(\tau) \varphi_{p-1}(\tau) F(\tau)(y(t_{n+1})) h^\rho(\tau) + \mathcal{O}(h^{p+1}) \end{array}\] and taking into account that from Lemma 1 \(\varphi_{p-1}(\tau)=\varphi_{p}(\tau)\) for \(\rho(\tau)\le p\), and that \(y_{n+1}-y(t_{n+1})=\mathcal{O}(h^{p+1})\), \(F(\tau)(y_{n+1})- F(\tau)(y(t_{n+1}))= \mathcal{O}(h^{p+1})\), we deduce that \[y_{n+1,1} - z(t_{n+1}) = \mathcal{O}(h^{p+1})\]

2.2 Order conditions↩︎

In order to find the conditions that the coefficients of the modified method must satisfy to have order \(p\) we will impose that the composed method has order \(p\). First let us note that the matrices \(A^{(n)}\) and \(\widehat A^{(n)}\) are related by \[\label{matrizm} \widehat A^{(n)} = A^{(n)} + \dfrac{1}{n} \begin{pmatrix} 0 & \\ M & 0 \\ 0 & M & 0 \\ \cdots & \cdots & \ddots &\ddots \\ 0 & \cdots & 0 & M & 0 \end{pmatrix} \equiv A^{(n)} + M^{(n)},\tag{14}\] with \(M=e_1 (a-b)^T\).

Let us give now some previous results

3 RK methods reusing the last stage with orders 2, 3, 4 and 5↩︎

In this section we show that, there exist standard RK methods with \(c_s=1\), order \(p= 2,3,4,5\) and minimum number of stages \(s=2,3,4,6\) whose associated RKRLS methods (\(s-1\) effective stages) have the same order \(p\).

3.1 RKRLS methods with \(s=2\) and order 2↩︎

There exists a one-parameter family of RK methods with \(s=2\) stages and order \(p=2\) given by \[\begin{array}{c|cc} 0 & \\ c_2 & c_2 \\ \hline & b_1 & b_2 \end{array} \qquad b_1 = 1 - b_2, \qquad b_2= 1/ c_2.\] With the condition \(c_2=1\) (\(b_2= 1/2\), \(b_1=1/2\)) there is a unique method. Acccording the results in Theorem 4, the associated RKRLS also has order 2.

3.2 RKRLS methods with \(s=3\) and order 3↩︎

There exists a one-parameter family of RK methods with \(s=3\) stages, \(c_3=1\) and order \(p=3\) given by \[\begin{array}{c|ccc} 0 & \\ c_2 & c_2 \\ 1 & 1-a_{32} & a_{32} \\ \hline & b_1 & b_2 & b_3 \end{array} \qquad \begin{array}{ll} b_2= \dfrac{1}{6c_2(1-c_2) }, \;& a_{32}= \dfrac{(1-c_2)}{c_2(2-3c_2)}, \\[12pt] b_3= \dfrac{2-3 c_2}{6(1-c_2)}, & b_1 = \dfrac{3c_2-1}{6c_2}. \end{array}\] None of these methods satisfies the condition \((b^TA)_1=0\) and a unique method (\(c_2=1/3\), \(a_{32}=2\), \(b_2=3/4\), \(b_3=1/4\)) satisfies the condition \(b_1=0\). Acccording the results in Theorem 4, the associated RKRLS also has order 3.

3.3 RKRLS methods with \(s=4\) and order 4↩︎

It is known that any RK method with four stages and order four must satisfy the simplifying condition \(b^T a = b^T - (b\cdot c)^T\) (usually denoted by \(D(1)\)) and therefore the last node must be \(c_4=1\). This implies that \((b^TA)_1=0\) for any method of order four and four stages. There are four families of methods of order four (see [1]).

One of the families is obtained with \(c_2=c_3\) and is given by the Butcher tableau \[\begin{array}{c|cccc} 0 & \\ \dfrac{1}{2} & \dfrac{1}{2} \\[5pt] \dfrac{1}{2} & \dfrac{a_{43}-1}{2a_{43}} & \dfrac{1}{2a_{43}} \\[8pt] 1 & 0 & 1-a_{43} & a_{43} \\[4pt] \hline \vrule height 16pt width 0pt& \dfrac{1}{6} & \dfrac{2-a_{43}}{3} & \dfrac{a_{43}}{3} & \dfrac{1}{6} \end{array}\] In this family, \(b_1 \ne 0\), but \((Ac)_4=1/2\). Therefore, according to Theorem 4 the associated RKRLS methods will have order 3. The classical RK, obtained for \(a_{43}=1\), belongs to this family and its associated RKRLS has order 3.

Another family is obtained with the condition \(c_3 (1-c_2)(c_2-c_3) \ne 0\). It is a two-parameter family depending on \(c_2\) and \(c_3\) and with coefficients \[\begin{array}{ll} b_2= \dfrac{2c_3-1}{12 c_2 (1-c_2)(c_3-c_2)}, & a_{32}= \dfrac{c_3(c_3-c_2)}{2c_2(1-2c_2)}, \\[12pt] b_3= \dfrac{1-2c_2}{12 c_3 (1-c_2)(c_3-c_2)}, & a_{42}= \dfrac{(1-c_2) (c_2+5c_3-4c_3^2-2)}{ 2 c_2 (c_3-c_2)(3+6c_2 c_3 -4c_2 -4c_3)} , \\[10pt] b_4= \dfrac{3+6c_2 c_3 -4c_2 -4c_3)}{12(1-c_2)(1-c_3)}, & a_{43} = \dfrac{(1-c_2) (1-c_3)(1-2c_2)}{ c_3(c_3-c_2)(3+6c_2 c_3 -4c_2 -4c_3)} . \end{array}\] The \(3/8\) RK method belongs to this family. For this method \(b_1=1/8 \ne 0\) and \((Ac)_4 = 1/3 \ne 1/2\). Theorem 4 guarantees order 2 for the associated RKRLS.

There is a subfamily, depending on one parameter \(c_2\), with \(b_1=0\) given by \[\begin{array}{ll} c_3=\dfrac{2 c_2-1}{6 c_2-2}, & a_{32}= \dfrac{1 - 4 c_2 + 6 c_2^2}{8 c_2(1 - 3 c_2)^2}, \\[12pt] b_2= \dfrac{1}{6 - 30 c_2 + 60 c_2^2 - 36 c_2^3}, & a_{42}= \dfrac{-(-1 + c_2)^2 (1 - 8 c_2 + 18 c_2^2)}{ 2 c_2 (1 - 10 c_2 + 36 c_2^2 - 60 c_2^3 + 36 c_2^4)}, \\[12pt] b_3= \dfrac{2 (-1 + 3 c_2)^3}{3 (-1 + 4 c_2) (1 - 4 c_2 + 6 c_2^2)}, & a_{43} = \dfrac{2 (1 - 3 c_2)^2 (-1 + c_2) (-1 + 4 c_2)}{1 - 10 c_2 + 36 c_2^2 - 60 c2^3 + 36 c2^4}, \\[10pt] b_4= \dfrac{1 - 6 c_2 + 6 c_2^2}{6 - 30 c_2 + 24 c_2^2}. \end{array}\] Acccording the results in Theorem 4, the associated RKRLS methods also have order 4.

3.4 RKRLS methods with \(s=6\) and order 5↩︎

Let us consider the three-parameter family of methods, depending on the coefficients \(c_3, c_4\) and \(c_5\), with \(s=6\) stages and order \(p=5\) obtained in [12]. This family was obtained by imposing \(c_6=1\) and the simplifying conditions \(D(1)\) \[b^T A = b^T - (b\cdot c)^T,\] and the modified \(C(2)\) \[Ac = c^2/2- (c_2^2/2) e_1.\] Consequently, all the methods satisfy \((b^T A)_1=0\) and \((Ac)_6=1/2\). The associated RKRLS methods will have order 3, In particular for the DOPRI method of order 5, for which \(b_1= 35/384 \ne 0\), the associated RKRLS has order 3. However the subfamily of methods with \[c_5 = \dfrac{3 - 5 c_3 - 5 c_4 + 10 c_3 c_4}{5 - 10 c_3 - 10 c_4 + 30 c_3 c_4},\] satisfies \(b_1=0\) which gives a two-parameter family of RKRLS methods with order 5.

4 Numerical experiments↩︎

In this section we present some numerical experiments to verify the results obtained in sections 2 and 3. In the experiments, we will use the following methods

New2

The only method with \(s=2\), \(p=2\) and \(c_2=1\). The associated RKRLS has order 2

New3

The only method with \(s=3\), \(p=3\), \(c_3=1\) and \(b_1=0\). The associated RKRLS has order 3

RK-3/8

The \(3/8\) RK method with \(s=4\), \(p=4\), \(c_2= 1/3\) and \(c_3= 2/3\). The associated RKRLS has order 2

RKClassic

The classical RK with \(s=4\), \(p=4\), \(c_2=c_3=1/2\). The associated RKRLS has order 3

New4

The method with with \(s=4\), \(p=4\), \(b_1=0\) and \(c_2= 1/6\). The associated RKRLS has order 4

DOPRI54

The Dormand and Prince method with with \(s=6\), \(p=5\), The RKRLS has order 3

New5

The method with with \(s=6\), \(p=4\), \(b_1=0\), and \[\begin{array}{l} c_3= 0.1574989977372333627197954851966028754675, \\ c_4= 0.5649477718721229393448029991747476923492, \\ c_5= 0.62386437635858023903237445638618758645. \end{array}\] The associated RKRLS has order 5.

With all these methods we have solved the Initial Value Problems

1. The two-body problem with excentricity \(e=0.1\) (problem D1 from the DETEST package [13], [14]) integrated along the interval \([0,20]\) \[\begin{array}{lll} x'' = -\dfrac{x}{\sqrt{(x^2+y^2)^3}},& x(0)= 1-e, & x'(0)=0, \\[10pt] y'' = -\dfrac{y}{\sqrt{(x^2+y^2)^3}}, & y(0)= 0, & y'(0)=\dfrac{\sqrt{1+e}}{\sqrt{1-e}}. \end{array}\]

2. The problem A3 from the DETEST package [13], [14] integrated along the interval \([0, 20]\) \[y' =y\cos t, \qquad y(0)=1.\]

Using step sizes \(h= 0.2/2^j, j= 0, \ldots, 5\) we have computed, for each method and problem, the global errors GE obtained at the end point of the integration interval. In Figures 2 and 3 we display the global error GE against the step size used in double logarithmic scale. It is clearly observed that in all cases, the slope of the curves coincides with the order predicted by Theorem 4.

When you reduce the computational cost of a given method by modifying it, it is not unusual to pay some price, either in stability or else in accuracy. In order to test this for RKRLS methods we have computed also the golbal errors given by the underlying methods without reusing the last stage. In Tables ¿tbl:table1? and ¿tbl:table2? we collect these values together with the ones given by the RKRLS method for problem 1 and in Tables ¿tbl:table3? and ¿tbl:table4? the results for problem 2. It is clear from the results that the accuracy of RKRLS methods give very similar global errors than their counterpart RK methods. It is remarkable the case of problem 1 integrated with the method RKRLS of order 4 with time step \(h=0.2\). The global error is very large, which can be due to stability. The RKRLS method has probably a smaller stability region than the associated RK method. The effects of reusing the last stage on stability of the method is currently being studied object.

5 Conclusions↩︎

We have studied the order of RK method in which the last function evaluation in a step is reused (RKRLS), substituting the first evaluation for the next step, by means of the composition of RK methods and the adjoint (also named reverse) of an RK method.

We have proved that if the composed RKRLS methods have order \(p\), the the RKRLS method has also order \(p\). We have obtained the conditions that the coefficients of an explicit RK method must satisfy to get an RKRLS method with the same order. With a detailed study of the families of \(s\)–stage explicit RK methods with the highest order \(p=s\) we have proved that with some simple additional conditions on the coefficients the associated RKRLS methods retain the order of the original RK method. This is relevant because we have proved that it is possible to maintain the order reducing the computational cost from \(s\) to \((s-1)\) evaluations of the vector field per step.

The results of some numerical experiments have confirmed the theoretical results.

References↩︎