Two-step General Linear Methods for Retarded Functional Differential Equations

TL;DR

This paper introduces a new class of two-step general linear methods specifically designed for efficiently solving retarded functional differential equations, including explicit methods up to fifth order, with considerations for stiff problems.

Contribution

The paper develops a novel class of two-step general linear methods tailored for retarded functional differential equations, achieving high order and addressing stiffness issues.

Findings

Explicit methods up to order five constructed

Methods designed to prevent order reduction in mildly stiff problems

Stage order close to uniform order for stability

Abstract

This paper presents a class of Two-Step General Linear Methods for the numerical solution of Retarded Functional Differential Equations. Explicit methods up to order five are constructed. To avoid order reduction for mildly stiff problems the uniform stage order of the methods is chosen to be close to uniform order.

Equations

$$\begin{array}{rl}y{}^{\prime }(t)& =f(t,y),t\in [t_0,T],\\ y(t_0)& =y_0,\end{array}$$

$$Y_i^{[n]}$$

$$y_i^{[n]}$$

$$F_i^{[n]}$$

$$y_1^{[n-1]}$$

$$Y_i^{[n]}$$

$$y_n$$

$$F_i^{[n]}$$

$$\begin{array}{r}Vu=u,\\ Uu=1.\end{array}$$

$$u_{i2}$$

$$v_{12}$$

$$K_i^{[n]}$$

$$b_j$$

$$y_n$$

$$K_i^{[n]}$$

$$Y_i^{[n]}$$

$$\displaystyle\Bigl{(}\forall\;z=(z_{1},\dots,z_{s+2})^{T},\quad z_{j}\in{\mathbb{R}^{d}},\;j=1,\dots,s+2\Bigr{)}:$$

$$\displaystyle\Bigl{(}\forall\;\varphi=(\varphi_{1},\dots,\varphi_{s+2})^{T},\quad\varphi_{j}\in\mathcal{C},\;j=1,\dots,s+2\Bigr{)}:$$

$$y{}^{\prime }(t)=f(t,y_t),$$

$$\begin{array}{rl}y{}^{\prime }(t)& =f(t,y_t),t\in [t_0,T],\\ y_{t_0}(\theta )& ={\varphi }_0(\theta ),\theta \in [-r,0].\end{array}$$

$$\displaystyle\Bigl{(}\forall\;(t,\varphi_{\,t}),(t,\psi_{\,t})\in\Omega\Bigr{)}:\quad\bigl{\|}f(t,\varphi_{\,t})-f(t,\psi_{\,t})\bigr{\|}_{\mathbb{R}^{d}}\leq L\bigl{\|}\varphi_{\,t}-\psi_{\,t}\bigr{\|}_{\mathcal{C}},$$

$${\eta }^{[n]}(\alpha h)$$

$$\alpha \in [0,1],$$

$$K_i^{[n]}$$

$$Y^{[n]i}(\alpha h)$$

$$\alpha \in [0,c_i],$$

$${\eta }^{[n]}(\theta )$$

$$Y^{[n]i}(\theta )$$

$${\eta }^{[n-1]}:[-r,h]\to {\mathbb{R}}^d,K_i^{[n-1]}\in {\mathbb{R}}^d$$

$$\text{computed in the step }n-1,$$

$$Y^{[n]i}$$

$$K_i^{[n]}$$

$$\displaystyle u_{i}(\cdot),\ \widetilde{a}_{ij}(\cdot),a_{ij}(\cdot),\quad v(\cdot),\ \widetilde{b_{j}}(\cdot),\ b_{j}(\cdot)\$$

$$a_{ij}(\cdot ),{\tilde{a}}_{ij}(\cdot ),u_i(\cdot ),$$

$$b_j(\cdot ),{\tilde{b}}_j(\cdot ),v(\cdot ),$$

$$c_i\in \mathbb{R},c_i\ge 0,$$

$$a_{ij}(0)={\tilde{a}}_{ij}(0)=0,$$

$$b_j(0)={\tilde{b}}_j(0)=0,$$

$$u_i(\cdot )$$

$$v(\cdot )$$

$${\overline{u}}_i(\cdot )=u_i(\cdot +c_i),{\overline{\tilde{a}}}_{ij}(\cdot )={\tilde{a}}_{ij}(\cdot +c_i),{\overline{a}}_{ij}(\cdot )=a_{ij}(\cdot +c_i),i,j=1,\dots ,s.$$

Full text

Keywords: Functional differential equations; two-step Runge-Kutta methods; A-methods; order conditions.

AMS Subject Classification: 65Q05

**Two-step Runge-Kutta methods for

Retarded Functional Differential Equations**

Anton Tuzov †††Department of Control systems, Siberian State Aerospace University, Krasnoyarsk, Russia, e-mail: tuzov@sibsau.ru

1 Two-step Runge-Kutta methods for

Ordinary Differential Equations

For the numerical approximation of the solution $y(t)$ of a system of Ordinary Differential Equations

[TABLE]

where $f:\ \mathbb{R}\times{\mathbb{R}^{d}}\rightarrow{\mathbb{R}^{d}},\quad y_{0}\in{\mathbb{R}^{d}}$,

we consider the class of General Linear Methods[5, 9]

[TABLE]

where $y^{[n-1]}_{1},\dots,y^{[n-1]}_{r}$ — input vectors, available at step number $n$, $Y^{[n]}_{1},\dots,Y^{[n]}_{s}$ — stage values, $F^{[n]}_{1},\dots,F^{[n]}_{s}$ — derivative values, $a_{ij},\ u_{ij},\ b_{ij},\ v_{ij},\ c_{i}$ — coefficients of the method, $h$ — integration stepsize.

Let us restrict ourselves to two-step methods and choose $r=s+2$

[TABLE]

and (1) takes the form

[TABLE]

In the construction of GLMs it is assumed that $y^{[n-1]}_{i}=u_{i}\,y(t_{n-1})+v_{i}\,hy\,^{\prime}(t_{n-1})+O(h^{2})$ and 'preconsistensy conditions' holds

[TABLE]

For (1) we have $u_{1}=1,\ v_{1}=0,\quad u_{2}=1,\ v_{2}=-1,\quad u_{2+i}=0,\ v_{2+i}=1,\;i=1,\dots,s$. It follows from (4) that

[TABLE]

Let us denote

[TABLE]

then the method (1) satisfying 'preconsistensy conditions' (4) takes the form

[TABLE]

2 Two-step Runge-Kutta methods for

Retarded Functional Differential Equations

We begin with notations introduced in [10] (see also [3]).

Let $r\in[0,+\infty)$, and $\mathcal{C}$ be the space of continuous functions $[-r,0]\rightarrow{\mathbb{R}^{d}}$, equipped with the uniform norm \bigl{\|}\varphi\bigr{\|}_{\mathcal{C}}=\max\limits_{\theta\in[-r,0]}\bigl{\|}\varphi(\theta)\bigr{\|}_{\mathbb{R}^{d}}\ ,\quad\varphi\in\mathcal{C},   where \bigl{\|}\cdot\bigr{\|}_{\mathbb{R}^{d}} is an arbitrary norm on ${\mathbb{R}^{d}}$ . Let $u$ be a continuous function $[a-r,b)\rightarrow{\mathbb{R}^{d}}$, where $a<b$. Then $\forall\,t\in[a,b)$ shift function is defined by $u_{t}(\theta)=u(t+\theta),\quad\theta\in[-r,0]$ and $u_{t}\in\mathcal{C}$.

Let $s$ be stage number of the method, $\|\cdot\|$ be the norm on ${\mathbb{R}^{(s+2)d}}$ defined as the maximum of the norms of $s+2$ subvectors in ${\mathbb{R}^{d}}$ :

[TABLE]

and ${\overline{\mathcal{C}}}$ be the space of continuous functions $[-r,0]\rightarrow{\mathbb{R}^{(s+2)d}}$, equipped with the uniform norm:

[TABLE]

Consider a system of Retarded Functional Differential Equations (RFDEs)

[TABLE]

where $f:\ \Omega\rightarrow{\mathbb{R}^{d}},\quad\Omega=\mathbb{R}\times\mathcal{C}$.

It is assumed that the function $f$ is continuous and has the first derivative with respect to the second argument $f\,^{\prime}:\quad\Omega\rightarrow{\mathcal{L}(\mathcal{C},\mathbb{R}^{d})}$, which is bounded and continuous with respect to the second argument. Then for each $(\sigma,\varphi)\in\Omega$ there exists an unique non-continuable solution $y=y(\sigma,\varphi):\quad[\sigma-r,\overline{t}\,)\rightarrow{\mathbb{R}^{d}}$ of (6) through $(\sigma,\varphi)$, where $\overline{t}=\overline{t}(\sigma,\varphi)$, $\overline{t}\in(\sigma,\infty]$, i.e. $y(t)$ satisfies (6) for $t\in[\sigma,\overline{t}\,)$ and $y_{\sigma}(\theta)=\varphi(\theta)$ for $\theta\in[-r,0]$.

We next consider the system (6) through $(t_{0},\phi_{\,0})\in\Omega$ on the $[t_{0},T]$, where $T\in(t_{0},\overline{t}(t_{0},\phi_{\,0}))$

[TABLE]

The first derivative of $f$ with respect to the second argument is bounded, hence $f$ satisfies the Lipschitz condition with respect to the second argument

[TABLE]

where L=\sup\limits_{(t,\,\varphi_{t})\,\in\,\Omega}\bigl{\|}f\,^{\prime}(t,\varphi_{\,t})\bigr{\|}_{\mathcal{L}(\mathcal{C},\mathbb{R}^{d})}.

We introduce the class of two-step Runge-Kutta (TSRK) methods for RFDEs on the basis of approach proposed in [3, 10]. We can reformulate the method (1) for RFDEs (7) with the equispaced mesh $\Delta_{N}=\left\{t_{n}:\quad t_{n}=t_{0}+nh,\;n=0,\dots,N\right\}$, $\quad h=\dfrac{T-t_{0}}{N}\;$ as follows

[TABLE]

[TABLE]

Thus $\eta^{[n]}(\alpha h)\approx y(t_{n-1}+\alpha h),\;\alpha\in[0,1]$, $\quad\eta^{[n]}(\theta)\approx y(t_{n-1}+\theta),\;\theta\in[-r,0]$, hence $\eta^{[n]}_{\>h}\approx y_{t_{n}}$ on $[-r-h,0]$, where $t_{n}=t_{n-1}+h$.

It is assumed that coefficients (a_{ij}(\cdot),b_{j}(\cdot),c_{i},\widetilde{a}_{ij}(\cdot),\widetilde{b}_{j}(\cdot),u_{i}(\cdot),v(\cdot))_{\,i,\,j=1,\dots,s}\ of TSRK methods satisfy the following conditions:

[TABLE]

The last two conditions correspondingly guarantee continuity of the stage functions $Y^{[n]\;i}_{\ c_{i}h}\in\mathcal{C}$ and the approximate solution $\eta^{[n]}_{\>h}\in\mathcal{C}$ provided that approximate solution computed in the previous step is continuous function $\eta^{[n-1]}_{\>h}\in\mathcal{C}.$

Remark 2.1**.**

If the conditions

[TABLE]

hold, the two-step method (9) becomes the one-step RK method for RFDEs introduced in [3, 10], where initial function $\ \varphi:=\eta^{[n-1]}_{\>h}$.

Definition 2.2**.**

The TSRK method with coefficients (a_{ij}(\cdot),b_{j}(\cdot),c_{i},\widetilde{a}_{ij}(\cdot),\widetilde{b}_{j}(\cdot),u_{i}(\cdot),v(\cdot))_{i,j=1,\dots,s}\ is called explicit if a_{ij}(\cdot)=0\text{ for all j such that}j\geq i,\quad i,j=1,\dots,s.

Make the change of the independent variables: $\alpha h=\omega+c_{i}h,\;\theta=\omega+c_{i}h$ and introduce the shifted coefficient functions:

[TABLE]

Then (9c) can be rewritten in the form:

[TABLE]

Let us introduce the following notation:

[TABLE]

This notation allows us to rewrite (9a) like (2):

[TABLE]

Denote $z^{[n]}(\omega)=\left(Y^{[n]\;1}_{\ c_{1}h}(\omega),\dots,Y^{[n]\;s}_{\ c_{s}h}(\omega),\eta^{[n]}_{\>h}(\omega),\eta^{[n]}(0)\right)^{T},\qquad n=1,2,\dots\;.$

Then the method (9) can be reformulated as follows:

[TABLE]

Denote {\overline{\Omega}}_{h}=\Bigl{\{}(\sigma,\varphi)\in\mathbb{R}\times{\overline{\mathcal{C}}}\;\Big{|}\bigl{(}\sigma+(c_{j}-1)h,\varphi_{j}\bigr{)}\in\Omega,\quad j=1,\dots,s+1\Bigr{\}}.

When one step of the method (2) is applied with stepsize $h$ to (7) for the computation of the solution through $(t_{n-1},z^{[n-1]})\in{\overline{\Omega}}_{h}$ on $[t_{n-1},t_{n}]$, it yields a continuous function $z^{[n]}=z^{[n]}(h,t_{n-1},z^{[n-1]})$ in ${\overline{\mathcal{C}}}$ , where $z^{[n]\;i}\approx y_{t_{n-1}+c_{i}h},\;i=1,\dots,{s+1},$ $z^{[n],\;{s+2}}\approx y(t_{n-1})$. Note that $z^{[n],\;{s+1}}=\eta^{[n]}_{h}\approx y_{t_{n}}$, since $c_{s+1}=1$. The approximation $z^{[n]}$ is defined for stepsizes $h\in(0,\overline{h}\,)$, where $\overline{h}=\overline{h}(t_{n-1},z^{[n-1]})$, moreover $(t_{n},z^{[n]})\in{\overline{\Omega}}_{h}$.

Let $z^{[0]}$ be computed by some starting method $S:\;z^{[0]}=S(h,t_{0},\phi_{\,0})$. When several steps of the method (2) are applied to (7) for the computation of the solution through $(t_{0},z^{[0]})\in{\overline{\Omega}}_{h}$ on $[t_{0},T]$, they yield the finite sequence of continuous functions $\{z^{[n]}\}_{n=1,\dots,N\;}$ in ${\overline{\mathcal{C}}}$ . These sequence is defined for stepsizes $h\in(0,\overline{h}_{N}\,)$, where $\overline{h}_{N}=\min\limits_{n=1,\dots,N\;}\overline{h}(t_{n-1},z^{[n-1]})$, moreover $(t_{n},z^{[n]})\in{\overline{\Omega}}_{h}$ for all $n=1,\dots,N\;$.

Rewrite the method (2) in matrix-vector form. For this purpose we define

[TABLE]

where e=\Bigl{(}1,\dots,1\Bigr{)}^{T},\quad\overline{u\vphantom{\widetilde{a}}}\,\Bigl{(}\dfrac{\omega}{h}\Bigr{)}=\Bigl{(}\overline{u\vphantom{\widetilde{a}}}_{1}\Bigl{(}\dfrac{\omega}{h}\Bigr{)},\dots,\overline{u\vphantom{\widetilde{a}}}_{s+1}\Bigl{(}\dfrac{\omega}{h}\Bigr{)}\Bigr{)}^{T},\quad

O is $(s+1)\times s$ - dimensional zero matrix.

Then the method (2) can be represented in the form of $A$-method [1, 2]:

[TABLE]

[TABLE]

Denote the exact value function of the method (17) by $Z^{[n]}(\omega)$ and the global error function by $q^{[n]}(\omega):$

[TABLE]

Definition 2.3**.**

The method (17) has uniform order of convergence $q$ if

[TABLE]

and discrete order of convergence $p$ if

[TABLE]

where $q$ and $p$ are positive integers.

Denote by $hd^{[n]}$ the local discretization error function obtained as the residual function by replacing $z^{[n-1]},\;z^{[n]}$ in (17) by the exact value functions $Z^{[n-1]},\;Z^{[n]}$ correspondingly:

[TABLE]

Definition 2.4**.**

The method (17) has uniform order of consistency $q$ if

[TABLE]

and discrete order of consistency $p$ if

[TABLE]

It follows from (18b) that $d^{[n]\;i}(\omega)=0$ for $\omega\in[-r,-c_{i}h]$. Subtracting (17) from (18) we have

[TABLE]

where $\Delta\varPhi^{[n]}=\varPhi(t_{n},Z^{[n]},h)-\varPhi(t_{n},z^{[n]},h)$.

We further extend the shifted coefficient functions \overline{u\vphantom{\widetilde{a}}}_{i}\Bigl{(}\dfrac{\omega}{h}\Bigr{)},\;\overline{\widetilde{a}\vphantom{\widetilde{a}}}_{ij}\Bigl{(}\dfrac{\omega}{h}\Bigr{)},\;\overline{a\vphantom{\widetilde{a}}}_{ij}\Bigl{(}\dfrac{\omega}{h}\Bigr{)} to $\omega\in[-r,-c_{i}h)$ by

[TABLE]

This extension is continuous, since \overline{u\vphantom{\widetilde{a}}}_{i}\Bigl{(}\dfrac{\omega}{h}\Bigr{)}\Bigr{|}_{\omega=-c_{i}h}=u_{i}(0)=1,\quad\overline{\widetilde{a}\vphantom{\widetilde{a}}}_{ij}\Bigl{(}\dfrac{\omega}{h}\Bigr{)}\Bigr{|}_{\omega=-c_{i}h}=\widetilde{a}_{ij}(0)=0,\quad \overline{a\vphantom{\widetilde{a}}}_{ij}\Bigl{(}\dfrac{\omega}{h}\Bigr{)}\Bigr{|}_{\omega=-c_{i}h}=a_{ij}(0)=0,\qquad i=1,\dots,s+1,\ j=1,\dots,s. Thus, the elements of $M(\cdot)$, $\widetilde{\mathcal{A}}(\cdot)$, $\mathcal{A}(\cdot)$ are continuous functions $[-r,0]\rightarrow{\mathbb{R}}$. From (19a) we have

[TABLE]

It follows from (19a) for $\omega=0$ that

[TABLE]

Definition 2.5**.**

The TSRK method (17) is called zero-stable if $M(0)$ is power bounded, i.e.

[TABLE]

Lemma 2.6**.**

The method (17) is zero-stable iff   $0\leq\overline{u\vphantom{\widetilde{a}}}_{s+1}(0)<2$.

Proof.

The condition (22) holds iff [5]

1. i)

the minimal polynomial of $M(0)$ has all its zeros in the closed unit disc and 2. ii)

all its multiple zeros in the open unit disc.

Denote by $\psi(\lambda)$ and $\Delta(\lambda)$ the minimal and characteristic polynomials of $M(0)$ correspondingly. Then $\psi(\lambda)=\dfrac{\Delta(\lambda)}{D_{s+1}(\lambda)}\,,\;$ where $D_{s+1}(\lambda)$ is the greatest common divisor of $(s+1)$ - t order minors of the characteristic matrix $\lambda I-M(0)$.

[TABLE]

Hence $\psi(\lambda)=\lambda(\lambda-1)(\lambda-(v-1))$ and zeros of the minimal polynomial $\psi(\lambda)$ are

[TABLE]

The condition (i) holds iff $0\leq v\leq 2$; (ii) holds iff $v\neq 2$. Thus, (22) holds iff $0\leq v<2$. ∎

Theorem 2.7**.**

Assume that the method (17) is zero-stable and the starting procedure $S$ for it, which specifies the starting value $z^{[0]}=S(h,t_{0},\phi_{\,0})$, such that \bigl{\vvvert}q^{[0]}\bigr{\vvvert}=O(h^{p-1}), \bigl{\|}q^{[0]}(0)\bigr{\|}=O(h^{p}) as $h\rightarrow 0$.

If the method has uniform order of consistency p-1:\;\max\limits_{n=1,\dots,N}\bigl{\vvvert}d^{[n]}\bigr{\vvvert}=O(h^{p-1}) and discrete order of consistency p:\;\max\limits_{n=1,\dots,N-1}\bigl{\|}d^{[n]}(0)\bigr{\|}=O(h^{p}) as $h\rightarrow 0$, then it has uniform order of convergence $p:$ \max\limits_{n=1,\dots,N}\bigl{\vvvert}q^{[n]}\bigr{\vvvert}=O(h^{p}),\quad h\rightarrow 0.

Proof.

Since the function $f$ satisfies the Lipschitz condition (8) with respect to the second argument with constant $L$, we have

[TABLE]

Hence

[TABLE]

The method (17) is zero-stable, i.e. (22) holds, where it can be assumed without loss of generality that $C\geq 1$. Denote $K=CL$. It follows from (21), (22) and (23) that

[TABLE]

where we assume that the sum is zero if the lower summation index exceeds the upper one.

By (20) and (23) it follows that

[TABLE]

Consider (2) for $n=1,\dots,N$. Since h\,\sum\limits_{l=1}^{n-1}\bigl{\|}d^{[l]}(0)\bigr{\|}\leq h\,(n-1)\max\limits_{l=1,\dots,N-1}\bigl{\|}d^{[l]}(0)\bigr{\|} and $h\,(n-1)\leq h\,N=T-t_{0}$ we obtain h\sum\limits_{l=1}^{n-1}\bigl{\|}d^{[l]}(0)\bigr{\|}=O(h^{p}) as $h\rightarrow 0,\quad n=1,\dots,N$.

Moreover,

[TABLE]

where $\widetilde{C}_{0}>0,\;D>0$.

By hypothesis, we have that \bigl{\vvvert}q^{[0]}\bigr{\vvvert}\leq\widetilde{\widetilde{C}}_{0}\,h^{p} for some $\widetilde{\widetilde{C}}_{0}>0$ and sufficiently small $h$. Denote {C}_{0}=\max\Bigl{\{}\widetilde{C}_{0},\widetilde{\widetilde{C}}_{0}\Bigr{\}}. Then

[TABLE]

Let $K_{n}$ be defined by

[TABLE]

Using (25), (26) and (19b) it is easy to prove by induction that

[TABLE]

Finally, it follows from (26) that $K_{n}=(1-Dh)^{-1}K_{n-1}=\dots=(1-Dh)^{-n}K_{0}$, where

(1-Dh)^{-n}=\Bigl{(}1+\dfrac{Dh}{1-Dh}\Bigr{)}^{n}\leq\exp\Bigl{(}\dfrac{Dnh}{1-Dh}\Bigr{)}\leq\exp\Bigl{(}\dfrac{D(T-t_{0})}{1-D\overline{h}}\Bigr{)}. Hence

[TABLE]

which concludes the proof. ∎

3 Order conditions

Assume that $f$ is of class $C^{l}$ with respect to the second argument for a sufficiently large $l$ and solution $y(t)$ of (7) is of piecewise class $C^{m}$ for a sufficiently large $m$. Let $c_{1}^{*},\dots,c_{s^{*}}^{*}$ such that $c_{1}^{*}<c_{2}^{*}<\dots<c_{s^{*}}^{*}$ and $\{c_{1}^{*},\dots,c_{s^{*}}^{*}\}=\{c_{1},\dots,c_{s}\}$, i.e. $c_{i}^{*}$ are distinct $c_{i}$ in increasing order.

Observe that convergence and consistency of the method (17) in Definitions 2.3 and 2.4 means stage convergence and consistency of the corresponding method (9). Now we consider weaker definitions.

Definition 3.1**.**

The method (9) has uniform order of consistency $q$ if

[TABLE]

and discrete order of consistency $p$ if

[TABLE]

Definition 3.2**.**

The method (9) has uniform order of convergence $q$ if for the corresponding method (17) the following condition holds:

[TABLE]

and discrete order of convergence $p$ if

[TABLE]

Lemma 3.3**.**

Let $p$ be a positive integer. If $y$ is of piecewise class $C^{p}$ then the local discretization error function $hd^{[n]}$ satisfy

[TABLE]

[TABLE]

Proof.

Consider $hd^{[n]\;i}$ defined by (18a) for $\omega\in[-c_{i}h,0]$, make the change of the independent variable: $\omega=(\alpha-c_{i})h$ and use (12). The proof follows by Taylor series expansion about $t_{n-1}$. ∎

For convenience we denote $\Gamma_{k}:=\Gamma_{s+1,\;k}$. Using (14) and (30) we obtain

[TABLE]

Remark 3.4**.**

If the conditions (11) hold the $\Gamma_{ik},\ \Gamma_{k}$

are the same as for the one-step RK method [3, 10].

In the following we assume that the TSRK method satisfies the conditions $\Gamma_{1}=0$ and $\Gamma_{i1}=0,\;i=1,\dots,s,$ that is

[TABLE]

The above condition is an equivalent form of uniform stage order one condition.

The proofs of the Theorems 3.5, 3.6, 3.7 are not difficult but rather technical. We omit them for the sake of brevity.

Theorem 3.5**.**

The TSRK method satisfying (3) has uniform order two iff $\ \Gamma_{2}=0$.

Theorem 3.6**.**

Let the TSRK method satisfy (3) and have uniform order two.

If \ \Gamma_{3}=0\ and $\sum\limits_{\begin{subarray}{c}i=1\\ c_{i}=c_{m}^{*}\end{subarray}}^{s}{b_{i}(\alpha)\Gamma_{i\,2}(\beta)}=0,\quad\alpha\in[0,1],\quad\beta\in[0,c_{m}^{*}],\quad m=1,\dots,s^{*}$

then the method has uniform order three.

Theorem 3.7**.**

Let the TSRK method satisfy (3) and have uniform order three.

[TABLE]

then the method has uniform order four.

Theorem 3.8**.**

*The TSRK method has uniform stage order $\widetilde{q}\;$ iff

$\Gamma_{ik}=0,\ \Gamma_{k}=0,\quad i=1,\dots,s,\quad k=1,\dots,\widetilde{q}.$*

Proof.

It follows from (29). ∎

The following results can be obtained as corollaries of Theorems 3.8 and 3.6, 3.7.

Corollary 3.9**.**

Let the TSRK method have uniform stage order two.

If \ \Gamma_{3}=0\ then the method has uniform order three.

Corollary 3.10**.**

Let the TSRK method have uniform stage order three.

If $\ \Gamma_{4}=0$ then the method has uniform order four.

The results of Corollary 3.9 and 3.10 can be easily generalized as follows.

Theorem 3.11**.**

Let the TSRK method have uniform stage order $\widetilde{q}$. It has uniform order $q=\widetilde{q}+1$ iff $\ \Gamma_{\widetilde{q}+1}=0$.

4 Construction of explicit TSRK methods of

uniform stage order four and five

Consider a two-stage explicit TSRK method satisfying (3). Its Butcher tableau is

[TABLE]

A natural choice will be to space out the abscissae $c_{i},\;i=1,\dots,s$ uniformly in the interval $[0,1]$ so that [6] $c_{1}=0,\ c_{2}=\dfrac{1}{s-1},\dots,\ c_{s-1}=\dfrac{s-2}{s-1},\ c_{s}=1$. In the case of $s=2$ we have $c_{1}=0,\ c_{2}=1$.

Since $c_{1}=0,$ conditions $\Gamma_{1k}(\alpha)=0,\quad\alpha\in[0,c_{1}],\;k=1,2,\dots$ reduce to $\Gamma_{1k}(0)=0,$ $\ k=1,2,\dots$ that follows from (10c). It also follows that $u_{1}(\cdot)=1,\ \widetilde{a}_{11}(\cdot)=0,\ \widetilde{a}_{12}(\cdot)=0$.

For the sake of brevity we omit the argument $\alpha$ of the method coefficient functions. By Theorem 3.11, the method has uniform order four and uniform stage order three if

$\Gamma_{k}=0,\;k=1,2,3,4$ and $\Gamma_{2\,k}=0,\;k=1,2,3,\;$ that is

[TABLE]

where $\alpha\in[0,1]$.

The coefficients are defined by

[TABLE]

where $\widetilde{a}_{2,2},\ \widetilde{b}_{2}$ remain free. The relation $\Gamma_{5}(1)=\dfrac{4}{15}\neq 0$ implies that it is impossible to attain discrete order five.

The uniform order and the uniform stage order can be increased by finding a suitable value for $c_{2}$. Assume that $c_{1}=0,\ c_{2}\neq 0$ (in general case $c_{2}\neq 1$). By theorem (3.11), the method has uniform order five and uniform stage order four if $\Gamma_{k}=0,\;k=1,2,3,4,5$ and $\Gamma_{2\,k}=0,\;k=1,2,3,4,\;$ that is

[TABLE]

where $\alpha\in[0,1]$ in the first five equations (4) and $\alpha\in[0,c_{2}]$ in other ones.

The coefficients are defined by

[TABLE]

To attain the discrete stage order five, we determine $c_{2}$ from $\Gamma_{2\,5}(1)=0$. We have

[TABLE]

The relation \Gamma_{6}(1)=-{\dfrac{16\ \ \bigl{(}17-2\sqrt{41}\,\bigr{)}}{75\,\bigl{(}71-11\sqrt{41}\,\bigr{)}}\neq 0} implies that it is impossible to attain discrete order six.

So we construct the explicit TSRK method of uniform order five, uniform stage order four and discrete stage order five.

Remark 4.1**.**

There is not a method of uniform stage order two in a class of explicit one-step RK methods for RFDEs.

Indeed, for explicit one-step RK methods $c_{1}=0,\ c_{2}\neq 0$ and $a_{2,j}=0,\ j=2,\dots,s$, hence $\Gamma_{2\,k}=-\dfrac{\alpha^{k}}{k!}\neq 0,\quad\alpha\in(0,c_{2}],\ k=2,3,\dots\ \ .$

It is known [8] that methods with low stage order suffer from the order reduction phenomenon when applied to stiff ODEs. Hence, the explicit TSRK methods may be more appropriate for some mildly stiff RFDEs.