Derivation of Three-Derivative Two-Step Runge–Kutta Methods

Abstract

In this paper, we develop explicit three-derivative two-step Runge–Kutta (ThDTSRK) schemes, and propose a simpler general form of the order accuracy conditions ($p\le 7$) by Albrecht’s approach, compared to the order conditions in terms of rooted trees. The parameters of the general high-order ThDTSRK methods are determined by utilizing the order conditions. We establish a theory for the A-stability property of ThDTSRK methods and identify optimal stability coefficients. Moreover, ThDTSRK methods can achieve the intended order of convergence using fewer stages than other schemes, making them cost-effective for solving the ordinary differential equations.

1. Introduction

The Runge–Kutta (RK) methods [1,2,3] are the most widely employed numerical schemes for solving ordinary differential equations (ODEs) and partial differential equations (PDEs). With the increasing demand for high-order accuracy methods, the classical high-order Runge–Kutta methods have exhibited some unavoidable shortcomings [4], such as the increased number of stages and computational time. Additionally, the numerical stability of the discretization method is also a primary focus.

In recent decades, there has been a concentrated effort to develop explicit stability-preserving discretization methods, while also achieving high-order convergence and superior computational efficiency. The discretization methods [5,6,7] that incorporate a combination of multistage, multistep, and multiderivative schemes offer an efficient pathway toward reaching these objectives. The main focus lies in establishing order conditions and ensuring numerical stability. Butcher [8,9] and Albrecht [10] proposed different devices to build the order conditions for the time-stepping schemes. Concerning numerical stability, temporal methods depend on the characteristics of A-stability or L-stability [11,12] to ensure numerical stability when solving ODEs. The two-derivative Runge–Kutta (TDRK) methods [13] are developed by rooted trees theory, and the higher-order TDRK schemes require fewer stages compared to the RK methods. The three-derivative Runge–Kutta (ThDRK) methods [14] are presented and stability theory is constructed by the A-stability condition. Moreover, the multiderivative RK methods combining multistep schemes provide additional options for building higher-order discretization schemes. Several successful methods have been developed by following this approach, such as the two-step Runge–Kutta (TSRK) schemes [15], two-derivative two-step Runge–Kutta (TDTSRK) methods [16], multistep RK (MSRK) schemes [17], and two-derivative multistep (TDMS) methods [18].

In this study, we derive the three-derivative two-step Runge-Kutta (ThDTSRK) schemes, which combine the benefits of multiderivative RK methods and multistep RK methods. We propose a simpler general form of the order accuracy conditions ($p\le 7$) by Albrecht’s approach [10], compared to the order conditions in terms of rooted trees. As a result, the ThDTSRK methods can achieve the desired order of accuracy by using fewer stages than other schemes, making them more efficient for solving ODEs. The organization of this paper is as follows. Section 2 proposes the structure of ThDTSRK Methods. Section 3 presents the order conditions that ensure accuracy preservation for ThDTSRK methods. Section 4 introduces the A-stability property of ThDTSRK methods and exhibits the optimal coefficients for the high-order ThDTSRK methods. In Section 5, we test the order accuracy of ThDTSRK methods using several numerical experiments. Finally, some conclusions and future work plans are given in Section 6.

2. Structure of ThDTSRK Methods

We consider the ordinary differential equation as

$$y^{\prime }=F\left(y\right),$$

(1)

and the second and third derivatives are defined by

$$y^{″}=G\left(y\right),y^{‴}=H\left(y\right).$$

(2)

The structure of s-stage ThDTSRK methods is depicted in Figure 1, and the fixed time step size $\Delta t$ is used to discretize the time interval. The explicit s-stage ThDTSRK schemes are given as follows:

$$\begin{array}{cc}\hfill Y_i^{n-1}=& y^{n-1}+\Delta t\sum _{j=1}^{i-1}{a}_{ij}F\left(Y_j^{n-1}\right)+\Delta t^2\sum _{j=1}^{i-1}{\widehat{a}}_{ij}G\left(Y_j^{n-1}\right)+\Delta t^3\sum _{j=1}^{i-1}{\overline{a}}_{ij}H\left(Y_j^{n-1}\right),\hfill \\ \hfill Y_i^n=& y^n+\Delta t\sum _{j=1}^{i-1}{a}_{ij}F\left(Y_j^n\right)+\Delta t^2\sum _{j=1}^{i-1}{\widehat{a}}_{ij}G\left(Y_j^n\right)+\Delta t^3\sum _{j=1}^{i-1}{\overline{a}}_{ij}H\left(Y_j^n\right),\hfill \\ \hfill y^{n+1}=& (1-\theta )y^n+\theta y^{n-1}+\Delta t\sum _{i=1}^s\left(v_{i}F\left(Y_i^n\right)+w_{i}F\left(Y_i^{n-1}\right)\right)\hfill \\ \hfill & +\Delta t^2\sum _{i=1}^s\left({\widehat{v}}_{i}G\left(Y_i^n\right)+{\widehat{w}}_{i}G\left(Y_i^{n-1}\right)\right)+\Delta t^3\sum _{i=1}^s\left({\overline{v}}_{i}H\left(Y_i^n\right)+{\overline{w}}_{i}H\left(Y_i^{n-1}\right)\right),\hfill \end{array}$$

(3)

where $y^{n-1}$, $y^n$, and $y^{n+1}$ are the solution values at time steps $t^{n-1}$, $t^n$, $t^{n+1}$, and $Y_i^{n-1}$, and $Y_i^n$ ($2\le i\le s$) denote the intermediate stage values at time steps $t^{n-1}$, $t^n$. We also define first-stage values equal to the time step values $Y_1^{n-1}=y^{n-1}$, $Y_1^n=y^n$. The ThDTSRK methods can be rewritten in the matrix form

$$\begin{array}{cc}\hfill {\mathit{Y}}^{n-1}& =\mathit{e}{y}^{n-1}+\Delta t\mathit{A}{\mathit{F}}^{n-1}+\Delta t^2\widehat{\mathit{A}}{\mathit{G}}^{n-1}+\Delta t^3\overline{\mathit{A}}{\mathit{H}}^{n-1},\hfill \\ \hfill {\mathit{Y}}^n& =\mathit{e}{y}^n+\Delta t\mathit{A}{\mathit{F}}^n+\Delta t^2\widehat{\mathit{A}}{\mathit{G}}^n+\Delta t^3\overline{\mathit{A}}{\mathit{H}}^n,\hfill \\ \hfill y^{n+1}& =(1-\theta )y^n+\theta y^{n-1}+\Delta t({\mathit{v}}^T{\mathit{F}}^n+{\mathit{w}}^T{\mathit{F}}^{n-1})\hfill \\ \hfill & +\Delta t^2({\widehat{\mathit{v}}}^T{\mathit{G}}^n+{\widehat{\mathit{w}}}^T{\mathit{G}}^{n-1})+\Delta t^3({\overline{\mathit{v}}}^T{\mathit{H}}^n+{\overline{\mathit{w}}}^T{\mathit{H}}^{n-1}),\hfill \end{array}$$

(4)

where $\mathit{e}$ is a vector of one, $\mathit{A}$ is the RK matrix, $\widehat{\mathit{A}}$ is the TDRK matrix, and $\overline{\mathit{A}}$ is the ThDRK matrix. $\mathit{v}$, $\widehat{\mathit{v}}$, $\overline{\mathit{v}}$, $\mathit{w}$, $\widehat{\mathit{w}}$, and $\overline{\mathit{w}}$ are weight vectors, and ${\mathit{Y}}^{n-1}$ and ${\mathit{Y}}^n$ are vectors of time stage values, which are defined by

$$\begin{array}{cc}\hfill \mathit{A}=& \left[\begin{array}{cccc}0& 0& \cdots & 0\\ a_{21}& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮\\ a_{s1}& a_{s2}& \cdots & 0\end{array}\right],\widehat{\mathit{A}}=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ {\widehat{a}}_{21}& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮\\ {\widehat{a}}_{s1}& {\widehat{a}}_{s2}& \cdots & 0\end{array}\right],\overline{\mathit{A}}=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ {\overline{a}}_{21}& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮\\ {\overline{a}}_{s1}& {\overline{a}}_{s2}& \cdots & 0\end{array}\right],\hfill \\ \hfill & {\mathit{Y}}^{n-1}=\left[\begin{array}{c}{Y}_1^{n-1}\\ Y_2^{n-1}\\ ⋮\\ Y_s^{n-1}\end{array}\right],{\mathit{Y}}^n=\left[\begin{array}{c}{Y}_1^n\\ Y_2^n\\ ⋮\\ Y_s^n\end{array}\right],\mathit{v}=\left[\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ v_s\end{array}\right],\widehat{\mathit{v}}=\left[\begin{array}{c}{\widehat{v}}_1\\ {\widehat{v}}_2\\ ⋮\\ {\widehat{v}}_s\end{array}\right],\hfill \\ \hfill & \overline{\mathit{v}}=\left[\begin{array}{c}{\overline{v}}_1\\ {\overline{v}}_2\\ ⋮\\ {\overline{v}}_s\end{array}\right],\mathit{w}=\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_s\end{array}\right],\widehat{\mathit{w}}=\left[\begin{array}{c}{\widehat{w}}_1\\ {\widehat{w}}_2\\ ⋮\\ {\widehat{w}}_s\end{array}\right],\overline{\mathit{w}}=\left[\begin{array}{c}{\overline{w}}_1\\ {\overline{w}}_2\\ ⋮\\ {\overline{w}}_s\end{array}\right].\hfill \end{array}$$

The derivatives ${\mathit{F}}^n$, ${\mathit{F}}^{n-1}$, ${\mathit{G}}^n$, ${\mathit{G}}^{n-1}$${\mathit{H}}^n$, and ${\mathit{H}}^{n-1}$ of the time stage vectors are defined by

$$\begin{array}{c}\hfill {\mathit{F}}^n=\left[\begin{array}{c}F\left(Y_1^n\right)\\ F\left(Y_2^n\right)\\ ⋮\\ F\left(Y_s^n\right)\end{array}\right],{\mathit{F}}^{n-1}=\left[\begin{array}{c}F\left(Y_1^{n-1}\right)\\ F\left(Y_2^{n-1}\right)\\ ⋮\\ F\left(Y_s^{n-1}\right)\end{array}\right],{\mathit{G}}^n=\left[\begin{array}{c}G\left(Y_1^n\right)\\ G\left(Y_2^n\right)\\ ⋮\\ G\left(Y_s^n\right)\end{array}\right],\\ \hfill {\mathit{G}}^{n-1}=\left[\begin{array}{c}G\left(Y_1^{n-1}\right)\\ G\left(Y_2^{n-1}\right)\\ ⋮\\ G\left(Y_s^{n-1}\right)\end{array}\right],{\mathit{H}}^n=\left[\begin{array}{c}H\left(Y_1^n\right)\\ H\left(Y_2^n\right)\\ ⋮\\ H\left(Y_s^n\right)\end{array}\right],{\mathit{H}}^{n-1}=\left[\begin{array}{c}H\left(Y_1^{n-1}\right)\\ H\left(Y_2^{n-1}\right)\\ ⋮\\ H\left(Y_s^{n-1}\right)\end{array}\right].\end{array}$$

The size of the matrices is $s\times s$ and the size of the vectors is $s\times 1$. The time-levels for discretization methods are referred to as the abscissas $\mathit{c}$, which can be computed by $\mathit{c}=\mathit{A}\mathit{e}$. The ThDTSRK methods in Equation (3) can be represented by the extended Butcher tableau:

$$\begin{array}{cccc}\mathit{c}\hfill & \mathit{A}& \widehat{\mathit{A}}& \overline{\mathit{A}}\\ \hfill & {\mathit{v}}^T& {\widehat{\mathit{v}}}^T& {\overline{\mathit{v}}}^T\\ \hfill & {\mathit{w}}^T& {\widehat{\mathit{w}}}^T& {\overline{\mathit{w}}}^T\end{array}$$

(5)

3. Order Conditions of ThDTSRK Methods

In this section, we develop a general form of the order accuracy conditions following Albrecht’s approach [10]; the detailed proof of order conditions for ThDTSRK methods is provided below.

First, the abscissas $\mathit{c}$ are utilized to characterize the time stage vectors ${\mathit{t}}^n$ and ${\mathit{t}}^{n-1}$, and the expressions are given as ${\mathit{t}}^n=t^n\mathit{e}+\mathit{c}\Delta t$ and ${\mathit{t}}^{n-1}=t^{n-1}\mathit{e}+\mathit{c}\Delta t$. By considering the truncation error ${\tau }^{n+1}$ and stage truncation errors ${\mathit{\tau }}^n$ and ${\mathit{\tau }}^{n-1}$, the exact solution of ThDTSRK methods in Equation (4) can be expressed as follows:

$$\begin{array}{cc}\hfill {\tilde{\mathit{Y}}}^{n-1}& =\mathit{e}{\tilde{y}}^{n-1}+\Delta t\mathit{A}{\tilde{\mathit{F}}}^{n-1}+\Delta t^2\widehat{\mathit{A}}{\tilde{\mathit{G}}}^{n-1}+\Delta t^3\overline{\mathit{A}}{\tilde{\mathit{H}}}^{n-1}+{\mathit{\tau }}^{n-1},\hfill \\ \hfill {\tilde{\mathit{Y}}}^n& =\mathit{e}{\tilde{y}}^n+\Delta t\mathit{A}{\tilde{\mathit{F}}}^n+\Delta t^2\widehat{\mathit{A}}{\tilde{\mathit{G}}}^n+\Delta t^3\overline{\mathit{A}}{\tilde{\mathit{H}}}^n+{\mathit{\tau }}^n,\hfill \\ \hfill {\tilde{y}}^{n+1}& =(1-\theta ){\tilde{y}}^n+\theta {\tilde{y}}^{n-1}+\Delta t({\mathit{v}}^T{\tilde{\mathit{F}}}^n+{\mathit{w}}^T{\tilde{\mathit{F}}}^{n-1})\hfill \\ \hfill & +\Delta t^2({\widehat{\mathit{v}}}^T{\tilde{\mathit{G}}}^n+{\widehat{\mathit{w}}}^T{\tilde{\mathit{G}}}^{n-1})+\Delta t^3({\overline{\mathit{v}}}^T{\tilde{\mathit{H}}}^n+{\overline{\mathit{w}}}^T{\tilde{\mathit{H}}}^{n-1})+{\tau }^{n+1},\hfill \end{array}$$

(6)

where the tilde superscript denotes the exact values of the variables. Then, we expand truncation errors and exact stage values to the Taylor series expansions at the time step $t^n$

$$\begin{array}{cc}\hfill {\mathit{\tau }}^{n-1}=& \sum _{k=1}^{\infty }{\mathit{\tau }}_k^{n-1}\Delta t^k{\tilde{y}}^{n\left(k\right)},{\mathit{\tau }}^n=\sum _{k=1}^{\infty }{\mathit{\tau }}_k^n\Delta t^k{\tilde{y}}^{n\left(k\right)},{\tau }^{n+1}=\sum _{k=1}^{\infty }{\tau }_k^{n+1}\Delta t^k{\tilde{y}}^{n\left(k\right)},\hfill \\ \hfill {\tilde{Y}}_i^{n-1}=& \sum _{k=0}^{\infty }\frac{{\left(\Delta t(c_i-1)\right)}^k}{k!}{\tilde{y}}^{n\left(k\right)},{\tilde{Y}}_i^n=\sum _{k=0}^{\infty }\frac{\Delta t^k{c}_i^k}{k!}{\tilde{y}}^{n\left(k\right)},{\tilde{Y}}^{n+1}=\sum _{k=0}^{\infty }\frac{\Delta t^k}{k!}{\tilde{y}}^{n\left(k\right)},\hfill \\ \hfill F\left({\tilde{Y}}_i^{n-1}\right)=& \sum _{k=1}^{\infty }\frac{{\left(\Delta t(c_i-1)\right)}^{k-1}}{(k-1)!}{\tilde{y}}^{n\left(k\right)},F\left({\tilde{Y}}_i^n\right)=\sum _{k=1}^{\infty }\frac{\Delta t^{k-1}{c}_i^{k-1}}{(k-1)!}{\tilde{y}}^{n\left(k\right)},\hfill \\ \hfill G\left({\tilde{Y}}_i^{n-1}\right)=& \sum _{k=2}^{\infty }\frac{{\left(\Delta t(c_i-1)\right)}^{k-2}}{(k-2)!}{\tilde{y}}^{n\left(k\right)},G\left({\tilde{Y}}_i^n\right)=\sum _{k=2}^{\infty }\frac{\Delta t^{k-2}{c}_i^{k-2}}{(k-2)!}{\tilde{y}}^{n\left(k\right)},\hfill \\ \hfill H\left({\tilde{Y}}_i^{n-1}\right)=& \sum _{k=3}^{\infty }\frac{{\left(\Delta t(c_i-1)\right)}^{k-3}}{(k-3)!}{\tilde{y}}^{n\left(k\right)},H\left({\tilde{Y}}_i^n\right)=\sum _{k=3}^{\infty }\frac{\Delta t^{k-3}{c}_i^{k-3}}{(k-3)!}{\tilde{y}}^{n\left(k\right)}.\hfill \end{array}$$

(7)

Taking Equation (7) into Equation (6), we obtain

$$\begin{array}{cc}\hfill {\mathit{\tau }}_k^n& =\frac{{\mathit{c}}^k}{k!}-\frac{\mathit{A}{\mathit{c}}^{k-1}}{(k-1)!}-\frac{\widehat{\mathit{A}}{\mathit{c}}^{k-2}}{(k-2)!}-\frac{\overline{\mathit{A}}{\mathit{c}}^{k-3}}{(k-3)!},\hfill \\ \hfill {\tau }_k^{n+1}& =\frac{1-\theta {(-1)}^k}{k!}-\frac{{\mathit{v}}^T{\mathit{c}}^{k-1}}{(k-1)!}-\frac{{\mathit{w}}^T{(\mathit{c}-\mathit{e})}^{k-1}}{(k-1)!}-\frac{{\widehat{\mathit{v}}}^T{\mathit{c}}^{k-2}}{(k-2)!}\hfill \\ \hfill & -\frac{{\widehat{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^{k-2}}{(k-2)!}-\frac{{\overline{\mathit{v}}}^T{\mathit{c}}^{k-3}}{(k-3)!}-\frac{{\overline{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^{k-3}}{(k-3)!}.\hfill \end{array}$$

(8)

The vectors ${\tilde{\mathit{Y}}}^{n-1}$, ${\tilde{\mathit{F}}}^{n-1}$, ${\tilde{\mathit{H}}}^{n-1}$ and ${\tilde{\mathit{G}}}^{n-1}$ can be expanded using Taylor series expansions at the time step $t^{n-1}$. By substituting these expansions into Equation (7), we can compute ${\mathit{\tau }}^{n-1}$ using the following expressions:

$${\mathit{\tau }}^{n-1}=\sum _{k=1}^{\infty }{\mathit{\tau }}_k^n\Delta t^k{\tilde{y}}^{n-1\left(k\right)}.$$

(9)

The derivatives ${\tilde{y}}^{n-1\left(k\right)}$ can be expanded by Taylor series expansions at the time step $t^n$, and we have ${\tilde{y}}^{n-1\left(k\right)}={\sum }_{m=0}^{\infty }\frac{{(-\Delta t)}^m}{m!}{\tilde{y}}^{n(k+m)}$, so the stage truncation errors ${\mathit{\tau }}_k^{n-1}$ can be devised by

$${\mathit{\tau }}_k^{n-1}=\sum _{j=1}^k\frac{{(-\Delta t)}^{k-j}}{(k-j)!}{\mathit{\tau }}_j^n.$$

(10)

After determining the truncation error ${\tau }^{n+1}$ and stage truncation errors ${\mathit{\tau }}^n$ and ${\mathit{\tau }}^{n-1}$, the subsequent task is to define the global error of the ThDTSRK methods by subtracting the exact value from the discretized numerical value. The specific procedure entails subtracting Equation (6) from Equation (4); we have

$$\begin{array}{cc}\hfill {\mathit{ϵ}}^{n-1}& =\mathit{e}{ϵ}^{n-1}+\Delta t\mathit{A}{\mathit{\delta }}^{n-1}+\Delta t^2\widehat{\mathit{A}}{\mathit{\sigma }}^{n-1}+\Delta t^3\overline{\mathit{A}}{\mathit{\eta }}^{n-1}-{\mathit{\tau }}^{n-1},\hfill \\ \hfill {\mathit{ϵ}}^n& =\mathit{e}{ϵ}^n+\Delta t\mathit{A}{\mathit{\delta }}^n+\Delta t^2\widehat{\mathit{A}}{\mathit{\sigma }}^n+\Delta t^3\overline{\mathit{A}}{\mathit{\eta }}^n-{\mathit{\tau }}^n,\hfill \\ \hfill {ϵ}^{n+1}& =(1-\theta ){ϵ}^n+\theta {ϵ}^{n-1}+\Delta t({\mathit{v}}^T{\mathit{\delta }}^n+{\mathit{w}}^T{\mathit{\delta }}^{n-1})\hfill \\ \hfill & +\Delta t^2({\widehat{\mathit{v}}}^T{\mathit{\sigma }}^n+{\widehat{\mathit{w}}}^T{\mathit{\sigma }}^{n-1})+\Delta t^3({\overline{\mathit{v}}}^T{\mathit{\eta }}^n+{\overline{\mathit{w}}}^T{\mathit{\eta }}^{n-1})-{\tau }^{n+1},\hfill \end{array}$$

(11)

where the global error $ϵ$, the global stage errors $\mathit{ϵ}$, and the derivatives $\mathit{\delta }$, $\mathit{\sigma }$, and $\mathit{\eta }$ can be described as

$$ϵ=y-\tilde{y},\mathit{ϵ}=\mathit{Y}-\tilde{\mathit{Y}},\mathit{\delta }=\mathit{F}-\tilde{\mathit{F}},\mathit{\sigma }=\mathit{G}-\tilde{\mathit{G}},\mathit{\eta }=\mathit{H}-\tilde{\mathit{H}}.$$

(12)

The above global errors and derivatives are power series in $\Delta t$:

$$\begin{array}{cc}\hfill {\mathit{ϵ}}^{n-1}& =\sum _{k=1}^p{\mathit{ϵ}}_k^{n-1}\Delta t^k+O(\Delta t^{p+1}),{\mathit{ϵ}}^n=\sum _{k=1}^p{\mathit{ϵ}}_k^n\Delta t^k+O(\Delta t^{p+1}),\hfill \\ \hfill {\mathit{\delta }}^{n-1}& =\sum _{k=1}^{p-1}{\mathit{\delta }}_k^{n-1}\Delta t^k+O(\Delta t^p),{\mathit{\delta }}^n=\sum _{k=1}^{p-1}{\mathit{\delta }}_k^n\Delta t^k+O(\Delta t^p),\hfill \\ \hfill {\mathit{\sigma }}^{n-1}& =\sum _{k=1}^{p-2}{\mathit{\sigma }}_k^{n-1}\Delta t^k+O(\Delta t^{p-1}),{\mathit{\sigma }}^n=\sum _{k=1}^{p-2}{\mathit{\sigma }}_k^n\Delta t^k+O(\Delta t^{p-1}),\hfill \\ \hfill {\mathit{\eta }}^{n-1}& =\sum _{k=1}^{p-3}{\mathit{\eta }}_k^{n-1}\Delta t^k+O(\Delta t^{p-2}),{\mathit{\eta }}^n=\sum _{k=1}^{p-3}{\mathit{\eta }}_k^n\Delta t^k+O(\Delta t^{p-2}).\hfill \end{array}$$

(13)

Substituting Equations (7) and (13) into Equation (11) yields

$$\begin{array}{cc}\hfill {\mathit{ϵ}}^{n-1}=& \mathit{e}{ϵ}^{n-1}+\sum _{k=1}^{p-1}\mathit{A}{\mathit{\delta }}_k^{n-1}\Delta t^{k+1}+\sum _{k=1}^{p-2}\widehat{\mathit{A}}{\mathit{\sigma }}_k^{n-1}\Delta t^{k+2}+\sum _{k=1}^{p-3}\overline{\mathit{A}}{\mathit{\eta }}_k^{n-1}\Delta t^{k+3}\hfill \\ & -\sum _{k=1}^p{\mathit{\tau }}_k^{n-1}{\tilde{y}}^{n\left(k\right)}\Delta t^k+\mathit{e}O(\Delta t^{p+1}),\hfill \\ \hfill {\mathit{ϵ}}^n=& \mathit{e}{ϵ}^n+\sum _{k=1}^{p-1}\mathit{A}{\mathit{\delta }}_k^n\Delta t^{k+1}+\sum _{k=1}^{p-2}\widehat{\mathit{A}}{\mathit{\sigma }}_k^n\Delta t^{k+2}+\sum _{k=1}^{p-3}\overline{\mathit{A}}{\mathit{\eta }}_k^n\Delta t^{k+3}\hfill \\ & -\sum _{k=1}^p{\mathit{\tau }}_k^n{\tilde{y}}^{n\left(k\right)}\Delta t^k+\mathit{e}O(\Delta t^{p+1}),\hfill \\ \hfill {ϵ}^{n+1}=& (1-\theta ){ϵ}^n+\theta {ϵ}^{n-1}+\sum _{k=1}^{p-1}({\mathit{v}}^T{\mathit{\delta }}_k^n+{\mathit{w}}^T{\mathit{\delta }}_k^{n-1})\Delta t^{k+1}+\sum _{k=1}^{p-2}({\widehat{\mathit{v}}}^T{\mathit{\sigma }}_k^n+{\widehat{\mathit{w}}}^T{\mathit{\sigma }}_k^{n-1})\Delta t^{k+2}\hfill \\ & +\sum _{k=1}^{p-3}({\overline{\mathit{v}}}^T{\mathit{\eta }}_k^n+{\overline{\mathit{w}}}^T{\mathit{\eta }}_k^{n-1})\Delta t^{k+3}-\sum _{k=1}^p{\tau }_k^{n+1}{\tilde{y}}^{n\left(k\right)}\Delta t^k+O(\Delta t^{p+1}).\hfill \end{array}$$

(14)

Then, we analyze the global error ${ϵ}^{n+1}$ for the ThDTSRK methods. Since the previous time discretization format was the p-order accuracy method, we can deduce the p-order conditions for the ThDTSRK methods:

$$\left\{\begin{array}{c}{\tau }_k^{n+1}=0,\hfill \\ {\mathit{v}}^T{\mathit{\delta }}_{k-1}^n+{\mathit{w}}^T{\mathit{\delta }}_{k-1}^{n-1}=0,\hfill \\ {\widehat{\mathit{v}}}^T{\mathit{\sigma }}_{k-2}^n+{\widehat{\mathit{w}}}^T{\mathit{\sigma }}_{k-2}^{n-1}=0,\hfill \\ {\overline{\mathit{v}}}^T{\mathit{\eta }}_{k-3}^n+{\overline{\mathit{w}}}^T{\mathit{\eta }}_{k-3}^{n-1}=0.\hfill \end{array}\right.$$

(15)

To complete the order conditions for the ThDTSRK methods, we need to establish the vectors ${\mathit{\delta }}_k^{}$, ${\mathit{\sigma }}_k^{}$, and ${\mathit{\eta }}_k^{}$. Indeed, these vectors can be recursively derived from the global stage errors $\mathit{ϵ}$. By employing the Taylor series expansions at the exact stage values $\tilde{\mathit{y}}$, the interrelation between the derivatives and the global stage errors of derivatives can be revealed as

$$\left\{\begin{array}{c}\mathit{F}=\tilde{\mathit{F}}+\sum _{j=1}^{\infty }\frac{1}{j!}{\left({\mathit{ϵ}}^{}\right)}^j·{\tilde{\mathit{F}}}^{\left(j\right)}\hfill \\ \mathit{G}=\tilde{\mathit{G}}+\sum _{j=1}^{\infty }\frac{1}{j!}{\left({\mathit{ϵ}}^{}\right)}^j·{\tilde{\mathit{G}}}^{\left(j\right)}\hfill \\ \mathit{H}=\tilde{\mathit{H}}+\sum _{j=1}^{\infty }\frac{1}{j!}{\left({\mathit{ϵ}}^{}\right)}^j·{\tilde{\mathit{H}}}^{\left(j\right)}\hfill \end{array}\right.⟹\left\{\begin{array}{c}\mathit{\delta }=\mathit{F}-\tilde{\mathit{F}}=\sum _{j=1}^{\infty }\frac{1}{j!}{\left({\mathit{ϵ}}^{}\right)}^j·{\tilde{\mathit{F}}}^{\left(j\right)}\hfill \\ \mathit{\sigma }=\mathit{G}-\tilde{\mathit{G}}=\sum _{j=1}^{\infty }\frac{1}{j!}{\left({\mathit{ϵ}}^{}\right)}^j·{\tilde{\mathit{G}}}^{\left(j\right)}\hfill \\ \mathit{\eta }=\mathit{H}-\tilde{\mathit{H}}=\sum _{j=1}^{\infty }\frac{1}{j!}{\left({\mathit{ϵ}}^{}\right)}^j·{\tilde{\mathit{H}}}^{\left(j\right)}\hfill \end{array}\right.$$

(16)

Here, the dot product represents component-wise multiplication, and the vectors of derivations are given as

$$\begin{array}{cc}\hfill {\tilde{\mathit{F}}}^{\left(j\right)}& ={\left[\frac{{\partial }^{j}F\left({\tilde{Y}}_1\right)}{{(\partial {\tilde{Y}}_1)}^j},\frac{{\partial }^{j}F\left({\tilde{Y}}_2\right)}{{(\partial {\tilde{Y}}_2)}^j},\dots ,\frac{{\partial }^{j}F\left({\tilde{Y}}_s\right)}{{(\partial {\tilde{Y}}_s)}^j}\right]}^T,\hfill \\ \hfill {\tilde{\mathit{G}}}^{\left(j\right)}& ={\left[\frac{{\partial }^{j}G\left({\tilde{Y}}_1\right)}{{(\partial {\tilde{Y}}_1)}^j},\frac{{\partial }^{j}G\left({\tilde{Y}}_2\right)}{{(\partial {\tilde{Y}}_2)}^j},\dots ,\frac{{\partial }^{j}G\left({\tilde{Y}}_s\right)}{{(\partial {\tilde{Y}}_s)}^j}\right]}^T,\hfill \\ \hfill {\tilde{\mathit{H}}}^{\left(j\right)}& ={\left[\frac{{\partial }^{j}H\left({\tilde{Y}}_1\right)}{{(\partial {\tilde{Y}}_1)}^j},\frac{{\partial }^{j}H\left({\tilde{Y}}_2\right)}{{(\partial {\tilde{Y}}_2)}^j},\dots ,\frac{{\partial }^{j}H\left({\tilde{Y}}_s\right)}{{(\partial {\tilde{Y}}_s)}^j}\right]}^T.\hfill \end{array}$$

Similarly, we derive expressions for ${\tilde{\mathit{F}}}^{\left(j\right)}$, ${\tilde{\mathit{G}}}^{\left(j\right)}$, and ${\tilde{\mathit{H}}}^{\left(j\right)}$ of Equation (16) by employing the Taylor series expansions at time step $t^n$

$$\begin{array}{cc}\hfill {\tilde{\mathit{F}}}^{n-1\left(j\right)}& =\sum _{l=0}^{\infty }\frac{\Delta t^l}{l!}{(\mathit{c}-\mathit{e})}^l·{\tilde{\mathit{F}}}^{n(j,l)},{\tilde{\mathit{F}}}^{n\left(j\right)}=\sum _{l=0}^{\infty }\frac{\Delta t^l}{l!}{\mathit{c}}^l·{\tilde{\mathit{F}}}^{n(j,l)},\hfill \\ \hfill {\tilde{\mathit{G}}}^{n-1\left(j\right)}& =\sum _{l=0}^{\infty }\frac{\Delta t^l}{l!}{(\mathit{c}-\mathit{e})}^l·{\tilde{\mathit{G}}}^{n(j,l)},{\tilde{\mathit{G}}}^{n\left(j\right)}=\sum _{l=0}^{\infty }\frac{\Delta t^l}{l!}{\mathit{c}}^l·{\tilde{\mathit{G}}}^{n(j,l)},\hfill \\ \hfill {\tilde{\mathit{H}}}^{n-1\left(j\right)}& =\sum _{l=0}^{\infty }\frac{\Delta t^l}{l!}{(\mathit{c}-\mathit{e})}^l·{\tilde{\mathit{H}}}^{n(j,l)},{\tilde{\mathit{H}}}^{n\left(j\right)}=\sum _{l=0}^{\infty }\frac{\Delta t^l}{l!}{\mathit{c}}^l·{\tilde{\mathit{H}}}^{n(j,l)},\hfill \end{array}$$

(17)

where the vector of derivation is as follows:

$$\begin{array}{c}\hfill {\tilde{\mathit{F}}}^{(j,l)}={\left[\frac{{\partial }^{j+l}F\left({\tilde{Y}}_1\right)}{{(\partial {\tilde{Y}}_1)}^j{(\partial t^n)}^l},\frac{{\partial }^{j+l}F\left({\tilde{Y}}_2\right)}{{(\partial {\tilde{Y}}_2)}^j{(\partial t^n)}^l},\dots ,\frac{{\partial }^{j+l}F\left({\tilde{Y}}_s\right)}{{(\partial {\tilde{Y}}_s)}^j{(\partial t^n)}^l}\right]}^T,\\ \hfill {\tilde{\mathit{G}}}^{(j,l)}={\left[\frac{{\partial }^{j+l}G\left({\tilde{Y}}_1\right)}{{(\partial {\tilde{Y}}_1)}^j{(\partial t^n)}^l},\frac{{\partial }^{j+l}G\left({\tilde{Y}}_2\right)}{{(\partial {\tilde{Y}}_2)}^j{(\partial t^n)}^l},\dots ,\frac{{\partial }^{j+l}G\left({\tilde{Y}}_s\right)}{{(\partial {\tilde{Y}}_s)}^j{(\partial t^n)}^l}\right]}^T,\\ \hfill {\tilde{\mathit{H}}}^{(j,l)}={\left[\frac{{\partial }^{j+l}H\left({\tilde{Y}}_1\right)}{{(\partial {\tilde{Y}}_1)}^j{(\partial t^n)}^l},\frac{{\partial }^{j+l}H\left({\tilde{Y}}_2\right)}{{(\partial {\tilde{Y}}_2)}^j{(\partial t^n)}^l},\dots ,\frac{{\partial }^{j+l}H\left({\tilde{Y}}_s\right)}{{(\partial {\tilde{Y}}_s)}^j{(\partial t^n)}^l}\right]}^T.\end{array}$$

According to Equations (16) and (17), we ultimately deduce the intended expansion:

$$\begin{array}{cc}\hfill {\mathit{\delta }}^{n-1}& =\sum _{j=1}^{\infty }\sum _{l=0}^{\infty }\frac{\Delta t^l}{j!l!}{\overline{\mathit{C}}}^l{\left({\mathit{ϵ}}^{n-1}\right)}^j·{\tilde{\mathit{F}}}^{n(j,l)},{\mathit{\delta }}^n=\sum _{j=1}^{\infty }\sum _{l=0}^{\infty }\frac{\Delta t^l}{j!l!}{\mathit{C}}^l{\left({\mathit{ϵ}}^n\right)}^j·{\tilde{\mathit{F}}}^{n(j,l)},\hfill \\ \hfill {\mathit{\sigma }}^{n-1}& =\sum _{j=1}^{\infty }\sum _{l=0}^{\infty }\frac{\Delta t^l}{j!l!}{\overline{\mathit{C}}}^l{\left({\mathit{ϵ}}^{n-1}\right)}^j·{\tilde{\mathit{G}}}^{n(j,l)},{\mathit{\sigma }}^n=\sum _{j=1}^{\infty }\sum _{l=0}^{\infty }\frac{\Delta t^l}{j!l!}{\mathit{C}}^l{\left({\mathit{ϵ}}^n\right)}^j·{\tilde{\mathit{G}}}^{n(j,l)},\hfill \\ \hfill {\mathit{\eta }}^{n-1}& =\sum _{j=1}^{\infty }\sum _{l=0}^{\infty }\frac{\Delta t^l}{j!l!}{\overline{\mathit{C}}}^l{\left({\mathit{ϵ}}^{n-1}\right)}^j·{\tilde{\mathit{H}}}^{n(j,l)},{\mathit{\eta }}^n=\sum _{j=1}^{\infty }\sum _{l=0}^{\infty }\frac{\Delta t^l}{j!l!}{\mathit{C}}^l{\left({\mathit{ϵ}}^n\right)}^j·{\tilde{\mathit{H}}}^{n(j,l)},\hfill \end{array}$$

(18)

where $\mathit{C}$ and $\overline{\mathit{C}}$ are diagonal matrices composed of vectors $\mathit{c}$ and $\mathit{c}-\mathit{e}$.

Finally, the terms of vectors $\mathit{\delta }$, $\mathit{\sigma }$, $\mathit{\eta }$, and $\mathit{ϵ}$ are provided in Appendix A. Based on Equation (15), the order conditions for the ThDTSRK schemes are presented in

Table 1. Order conditions ($p\le 7$) for the ThDTSRK schemes.

p	Conditions
1	${\mathit{v}}^T\mathit{e}+{\mathit{w}}^T\mathit{e}=1+\theta$.
2	${\mathit{v}}^T\mathit{c}+{\widehat{\mathit{v}}}^T\mathit{e}+{\mathit{w}}^T(\mathit{c}-\mathit{e})+{\widehat{\mathit{w}}}^T\mathit{e}=\frac{1-\theta }{2}$.
3	${\mathit{v}}^T{\mathit{c}}^2+2{\widehat{\mathit{v}}}^T\mathit{c}+2{\overline{\mathit{v}}}^T\mathit{e}+{\mathit{w}}^T{(\mathit{c}-\mathit{e})}^2+2{\widehat{\mathit{w}}}^T(\mathit{c}-\mathit{e})+2{\overline{\mathit{w}}}^T\mathit{e}=\frac{1+\theta }{3}$,
	$({\mathit{v}}^T+{\mathit{w}}^T)(\frac{1}{2}{\mathit{c}}^2-\mathit{A}\mathit{c}-\widehat{\mathit{A}}\mathit{e})=0$.
4	${\mathit{v}}^T{\mathit{c}}^3+3{\widehat{\mathit{v}}}^T{\mathit{c}}^2+6{\overline{\mathit{v}}}^T\mathit{c}+{\mathit{w}}^T{(\mathit{c}-\mathit{e})}^3+3{\widehat{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^2+6{\overline{\mathit{w}}}^T(\mathit{c}-\mathit{e})=\frac{1-\theta }{4}$,
	$({\mathit{v}}^T\mathit{C}+{\mathit{w}}^T\overline{\mathit{C}})(\frac{1}{2}{\mathit{c}}^2-\mathit{A}\mathit{c}-\widehat{\mathit{A}}\mathit{e})=0$,   $({\mathit{v}}^T\mathit{A}+{\mathit{w}}^T\mathit{A})(\frac{1}{2}{\mathit{c}}^2-\mathit{A}\mathit{c}-\widehat{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T+{\mathit{w}}^T)(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})-{\mathit{w}}^T(\frac{1}{2}{\mathit{c}}^2-\mathit{A}\mathit{c}-\widehat{\mathit{A}}\mathit{e})=0$.
	$({\widehat{\mathit{v}}}^T+{\widehat{\mathit{w}}}^T)(\frac{1}{2}{\mathit{c}}^2-\mathit{A}\mathit{c}-\widehat{\mathit{A}}\mathit{e})=0$,
5	${\mathit{v}}^T{\mathit{c}}^4+4{\widehat{\mathit{v}}}^T{\mathit{c}}^3+12{\overline{\mathit{v}}}^T{\mathit{c}}^2+{\mathit{w}}^T{(\mathit{c}-\mathit{e})}^4+4{\widehat{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^3+12{\overline{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^2=\frac{1+\theta }{5}$,
	$\frac{1}{2}{\mathit{c}}^2-\mathit{A}\mathit{c}-\widehat{\mathit{A}}\mathit{e}=\mathbf{0}$,   $({\mathit{v}}^T\mathit{C}+{\mathit{w}}^T\overline{\mathit{C}})(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T\mathit{A}+{\mathit{w}}^T\mathit{A})(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,  $({\widehat{\mathit{v}}}^T+{\widehat{\mathit{w}}}^T)(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T+{\mathit{w}}^T)(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})-{\mathit{w}}^T(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$.
6	${\mathit{v}}^T{\mathit{c}}^5+5{\widehat{\mathit{v}}}^T{\mathit{c}}^4+20{\overline{\mathit{v}}}^T{\mathit{c}}^3+{\mathit{w}}^T{(\mathit{c}-\mathit{e})}^5+5{\widehat{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^4+20{\overline{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^3=\frac{1-\theta }{6}$,
	$({\mathit{v}}^T{\mathit{A}}^2+{\mathit{w}}^T{\mathit{A}}^2)(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T\mathit{C}\mathit{A}+{\mathit{w}}^T\overline{\mathit{C}}\mathit{A})(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T\mathit{A}\mathit{C}+{\mathit{w}}^T\mathit{A}\overline{\mathit{C}})(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T{\mathit{C}}^2+{\mathit{w}}^T{\overline{\mathit{C}}}^2)(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T\widehat{\mathit{A}}+{\mathit{w}}^T\widehat{\mathit{A}})(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\widehat{\mathit{v}}}^T\mathit{A}+{\widehat{\mathit{w}}}^T\mathit{A})(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\widehat{\mathit{v}}}^T\mathit{C}+{\widehat{\mathit{w}}}^T\overline{\mathit{C}})(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,   $({\overline{\mathit{v}}}^T+{\overline{\mathit{w}}}^T)(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T\mathit{C}+{\mathit{w}}^T\overline{\mathit{C}})(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})$$-{\mathit{w}}^T\overline{\mathit{C}}(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T\mathit{A}+{\mathit{w}}^T\mathit{A})(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})$$-{\mathit{w}}^T\mathit{A}(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\widehat{\mathit{v}}}^T+{\widehat{\mathit{w}}}^T)(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})$$-{\widehat{\mathit{w}}}^T(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$,
	$({\mathit{v}}^T+{\mathit{w}}^T)(\frac{1}{120}{\mathit{c}}^5-\frac{1}{24}\mathit{A}{\mathit{c}}^4-\frac{1}{6}\widehat{\mathit{A}}{\mathit{c}}^3-\frac{1}{2}\overline{\mathit{A}}{\mathit{c}}^2)$$-{\mathit{w}}^T(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})$
	$+\frac{1}{2}{\mathit{w}}^T(\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e})=0$.
7	${\mathit{v}}^T{\mathit{c}}^6+6{\widehat{\mathit{v}}}^T{\mathit{c}}^5+30{\overline{\mathit{v}}}^T{\mathit{c}}^4+{\mathit{w}}^T{(\mathit{c}-\mathit{e})}^6+6{\widehat{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^5+30{\overline{\mathit{w}}}^T{(\mathit{c}-\mathit{e})}^4=\frac{1+\theta }{7}$,
	$\frac{1}{6}{\mathit{c}}^3-\frac{1}{2}\mathit{A}{\mathit{c}}^2-\widehat{\mathit{A}}\mathit{c}-\overline{\mathit{A}}\mathit{e}=\mathbf{0}$,    $({\mathit{v}}^T{\mathit{A}}^2+{\mathit{w}}^T{\mathit{A}}^2)(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\mathit{v}}^T\mathit{A}\mathit{C}+{\mathit{w}}^T\mathit{A}\overline{\mathit{C}})(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\mathit{v}}^T\mathit{C}\mathit{A}+{\mathit{w}}^T\overline{\mathit{C}}\mathit{A})(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\mathit{v}}^T{\mathit{C}}^2+{\mathit{w}}^T{\overline{\mathit{C}}}^2)(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\mathit{v}}^T\widehat{\mathit{A}}+{\mathit{w}}^T\widehat{\mathit{A}})(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\widehat{\mathit{v}}}^T\mathit{A}+{\widehat{\mathit{w}}}^T\mathit{A})(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\widehat{\mathit{v}}}^T\mathit{C}+{\widehat{\mathit{w}}}^T\overline{\mathit{C}})(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\overline{\mathit{v}}}^T+{\overline{\mathit{w}}}^T)(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\mathit{v}}^T\mathit{A}+{\mathit{w}}^T\mathit{A})(\frac{1}{120}{\mathit{c}}^5-\frac{1}{24}\mathit{A}{\mathit{c}}^4-\frac{1}{6}\widehat{\mathit{A}}{\mathit{c}}^3-\frac{1}{2}\overline{\mathit{A}}{\mathit{c}}^2)$$-{\mathit{w}}^T\mathit{A}(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\mathit{v}}^T\mathit{C}+{\mathit{w}}^T\overline{\mathit{C}})(\frac{1}{120}{\mathit{c}}^5-\frac{1}{24}\mathit{A}{\mathit{c}}^4-\frac{1}{6}\widehat{\mathit{A}}{\mathit{c}}^3-\frac{1}{2}\overline{\mathit{A}}{\mathit{c}}^2)$$-{\mathit{w}}^T\overline{\mathit{C}}(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\widehat{\mathit{v}}}^T+{\widehat{\mathit{w}}}^T)(\frac{1}{120}{\mathit{c}}^5-\frac{1}{24}\mathit{A}{\mathit{c}}^4-\frac{1}{6}\widehat{\mathit{A}}{\mathit{c}}^3-\frac{1}{2}\overline{\mathit{A}}{\mathit{c}}^2)$$-{\widehat{\mathit{w}}}^T(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$,
	$({\mathit{v}}^T+{\mathit{w}}^T)(\frac{1}{720}{\mathit{c}}^6-\frac{1}{120}\mathit{A}{\mathit{c}}^5-\frac{1}{24}\widehat{\mathit{A}}{\mathit{c}}^4-\frac{1}{6}\overline{\mathit{A}}{\mathit{c}}^3)-{\mathit{w}}^T(\frac{1}{120}{\mathit{c}}^5-\frac{1}{24}\mathit{A}{\mathit{c}}^4-\frac{1}{6}\widehat{\mathit{A}}{\mathit{c}}^3-\frac{1}{2}\overline{\mathit{A}}{\mathit{c}}^2)$
	$+\frac{1}{2}{\mathit{w}}^T(\frac{1}{24}{\mathit{c}}^4-\frac{1}{6}\mathit{A}{\mathit{c}}^3-\frac{1}{2}\widehat{\mathit{A}}{\mathit{c}}^2-\overline{\mathit{A}}\mathit{c})=0$.

. Moreover, further information can be obtained from our source code [19].

4. Stability Analysis

4.1. A-Stability Property

To investigate the A-stability property of the ThDTSRK methods, the classical test equation $y^{\prime }=\lambda y$ is utilized, where $\lambda$ is a complex parameter. The ThDTSRK methods in Equation (4) can be expressed in a compact matrix form:

$$\left[\begin{array}{c}\mathit{T}\\ {\mathit{y}}^{n+1}\end{array}\right]=\left[\begin{array}{cccc}{\mathit{A}}_1& {\mathit{B}}_1& {\mathit{C}}_1& {\mathit{D}}_1\\ {\mathit{A}}_2& {\mathit{B}}_2& {\mathit{C}}_2& {\mathit{D}}_2\end{array}\right]\left[\begin{array}{c}\Delta tF\left(\mathit{T}\right)\\ \Delta t^{2}G\left(\mathit{T}\right)\\ \Delta t^{3}H\left(\mathit{T}\right)\\ {\mathit{y}}^n\end{array}\right],$$

(19)

where the vectors $\mathit{T}$, ${\mathit{y}}^n$, and ${\mathit{y}}^{n+1}$ are specified by

$$\mathit{T}=\left[\begin{array}{c}{\mathit{Y}}^{n-1}\\ {\mathit{Y}}^n\end{array}\right],{\mathit{y}}^n=\left[\begin{array}{c}{\mathit{Y}}^{n-1}\\ y^{n-1}\\ y^n\end{array}\right],{\mathit{y}}^{n+1}=\left[\begin{array}{c}{\mathit{Y}}^n\\ y^n\\ y^{n+1}\end{array}\right],$$

(20)

and the matrices ${\mathit{A}}_1,{\mathit{A}}_2,{\mathit{B}}_1,{\mathit{B}}_2,{\mathit{C}}_1,{\mathit{C}}_2,{\mathit{D}}_1,{\mathit{D}}_2$ are given by

$$\begin{array}{cc}\hfill & {\mathit{A}}_1=\left[\begin{array}{cc}\mathit{A}& \mathbf{0}\\ \mathbf{0}& \mathit{A}\end{array}\right],{\mathit{B}}_1=\left[\begin{array}{cc}\widehat{\mathit{A}}& \mathbf{0}\\ \mathbf{0}& \widehat{\mathit{A}}\end{array}\right],{\mathit{C}}_1=\left[\begin{array}{cc}\overline{\mathit{A}}& \mathbf{0}\\ \mathbf{0}& \overline{\mathit{A}}\end{array}\right],\hfill \\ \hfill & {\mathit{D}}_1=\left[\begin{array}{ccc}\mathbf{0}& \mathit{e}& 0\\ \mathbf{0}& 0& \mathit{e}\end{array}\right],{\mathit{A}}_2=\left[\begin{array}{cc}\mathbf{0}& \mathit{A}\\ 0& 0\\ {\mathit{w}}^T& {\mathit{v}}^T\end{array}\right],{\mathit{B}}_2=\left[\begin{array}{cc}\mathbf{0}& \widehat{\mathit{A}}\\ 0& 0\\ {\widehat{\mathit{w}}}^T& {\widehat{\mathit{v}}}^T\end{array}\right],\hfill \\ \hfill & {\mathit{C}}_2=\left[\begin{array}{cc}\mathbf{0}& \overline{\mathit{A}}\\ 0& 0\\ {\overline{\mathit{w}}}^T& {\overline{\mathit{v}}}^T\end{array}\right],{\mathit{D}}_2=\left[\begin{array}{ccc}{\mathbf{0}}_{s\times s}& {\mathbf{0}}_{s\times 1}& \mathit{e}\\ {\mathbf{0}}_{1\times s}& 0& 1\\ {\mathbf{0}}_{1\times s}& \theta & 1-\theta \end{array}\right].\hfill \end{array}$$

(21)

Based on the stability theory for TSRK schemes [20], the recurrence relation can be given by ${\mathit{y}}^{n+1}=\mathit{M}\left(z\right){\mathit{y}}^n$ and the stability matrix $\mathit{M}\left(z\right)$ can be determined as

$$\mathit{M}\left(z\right)={\mathit{D}}_2+\left(z{\mathit{A}}_2+z^2{\mathit{B}}_2+z^3{\mathit{C}}_2\right){\left({\mathit{I}}_{2s}-z{\mathit{A}}_1-z^2{\mathit{B}}_1-z^3{\mathit{C}}_1\right)}^{-1}{\mathit{D}}_1,$$

(22)

where $z=\lambda \Delta t$, and the stability polynomial $f(\alpha ,z)$ of the ThDTSRK schemes Equation (19) is defined by

$$f(\alpha ,z)=det\left(\alpha {\mathit{I}}_{s+2}-\mathit{M}\left(z\right)\right).$$

(23)

Then, we compute the roots of equation $f(\alpha ,z)=0$ and obtain the expression of the roots as $\alpha =g\left(z\right)$. The region of absolute stability $\mathbb{D}$ is the definition that all roots are inside the unit circle; we have

$$\begin{array}{c}\hfill \mathbb{D}=\{z\in \mathbb{C}:\mid {\alpha }_k\mid =\mid g_k\left(z\right)\mid <1,k=1,2,\cdots ,s+2\}.\end{array}$$

(24)

4.2. Optimal ThDTSRK Schemes

By employing stability theory, A-stability preserving ThDTSRK${}_{sp}$ methods have been proposed, in which the subscript denotes the p-order scheme with s-stage. Here, we display the ThDTSRK schemes with two stages; the stability polynomial $f(\alpha ,z)$ is expressed as

$$\begin{array}{cc}\hfill & f(\alpha ,z)={\alpha }^2\left({\alpha }^2-{\phi }_1\alpha -{\phi }_2\right),\hfill \\ \hfill {\phi }_1=& 1+(v_1+v_2)z+(a_{21}{v}_2+{\widehat{v}}_1+{\widehat{v}}_2)z^2+(\frac{1}{2}{a}_{21}^2{v}_2+a_{21}{\widehat{v}}_2+{\overline{v}}_1+{\overline{v}}_2)z^3+\hfill \\ \hfill & (\frac{1}{6}{a}_{21}^3{v}_2+\frac{1}{2}{a}_{21}^2{\widehat{v}}_2+a_{21}{\overline{v}}_2)z^4+(\frac{1}{6}{a}_{21}^3{\widehat{v}}_2+\frac{1}{2}{a}_{21}^2{\overline{v}}_2)z^5+\frac{1}{6}{a}_{21}^3{\overline{v}}_2{z}^6-\theta ,\hfill \\ \hfill {\phi }_2=& (w_1+w_2)z+(a_{21}{w}_2+{\widehat{w}}_1+{\widehat{w}}_2)z^2+(\frac{1}{2}{a}_{21}^2{w}_2+a_{21}{\widehat{w}}_2+{\overline{w}}_1+{\overline{w}}_2)z^3+\hfill \\ \hfill & (\frac{1}{6}{a}_{21}^3{w}_2+\frac{1}{2}{a}_{21}^2{\widehat{w}}_2+a_{21}{\overline{w}}_2)z^4+(\frac{1}{6}{a}_{21}^3{\widehat{w}}_2+\frac{1}{2}{a}_{21}^2{\overline{w}}_2)z^5+\frac{1}{6}{a}_{21}^3{\overline{w}}_2{z}^6+\theta .\hfill \end{array}$$

(25)

Next, we calculate the roots of equation $f(\alpha ,z)=0$ and the expression of the roots is shown as

$$\left\{\begin{array}{c}{\alpha }_1=\frac{1}{2}({\phi }_1+\sqrt{{\phi }_1^2+4{\phi }_2}),\hfill \\ {\alpha }_2=\frac{1}{2}({\phi }_1-\sqrt{{\phi }_1^2+4{\phi }_2}),\hfill \\ {\alpha }_3=0,\hfill \\ {\alpha }_4=0.\hfill \end{array}\right.$$

(26)

Utilizing Equations (24) and (26), we can obtain the region of absolute stability and determine the optimal corresponding parameters for the ThDTSRK schemes. The aim of this optimization process is to construct ThDTSRK methods with large region intervals $[\mathcal{L},0]$ of absolute stability. Here, we define the intersection $\mathcal{L}$ of the absolute stability region with the negative real line. Similarly to the authors of reference [21], this will be carried out by solving the optimization problem to find the optimal ThDTSRK schemes with the largest value of $\left|\mathcal{L}\right|$.

(1) 

Two-Stage Fifth-Order ThDTSRK Scheme (ThDTSRK${}_{25}$)

The coefficients of the general two-stage fifth-order ThDTSRK method are then given as a eight-parameter system, depending on $\theta ,a_{21},w_1,w_2,{\widehat{v}}_2,{\widehat{w}}_2,{\overline{v}}_2,{\overline{w}}_2,$ which are related through

$$\begin{array}{cc}\hfill v_1=& 1-w_1+\theta ;v_2=-w_2;{\widehat{a}}_{21}=a_{21}^2/2;{\overline{a}}_{21}=a_{21}^3/6;\hfill \\ \hfill {\widehat{v}}_1=& (3-20{\widehat{v}}_2+60a_{21}^2{\widehat{v}}_2+40a_{21}^3{\widehat{v}}_2+120a_{21}{\overline{v}}_2+120a_{21}^2{\overline{v}}_2+10w_1+10w_2+\hfill \\ & 20a_{21}{w}_2-40a_{21}^3{w}_2-60a_{21}^2{\widehat{w}}_2+40a_{21}^3{\widehat{w}}_2-120a_{21}{\overline{w}}_2+120a_{21}^2{\overline{w}}_2-7\theta )/20;\hfill \\ \hfill {\widehat{w}}_1=& (1-\theta )/2+w_1+w_2-{\widehat{v}}_1-{\widehat{v}}_2-{\widehat{w}}_2;\hfill \\ \hfill {\overline{v}}_1=& (23-60a_{21}{\widehat{v}}_2-120a_{21}^2{\widehat{v}}_2-60a_{21}^3{\widehat{v}}_2-60{\overline{v}}_2-240a_{21}{\overline{v}}_2-180a_{21}^2{\overline{v}}_2-5w_1-5w_2+\hfill \\ & 30a_{21}^2{w}_2+60a_{21}^3{w}_2+60a_{21}^2{\widehat{w}}_2-60a_{21}^3{\widehat{w}}_2+120a_{21}{\overline{w}}_2-180a_{21}^2{\overline{w}}_2+3\theta )/60;\hfill \\ \hfill {\overline{w}}_1=& (8-60a_{21}^2{\widehat{v}}_2-60a_{21}^3{\widehat{v}}_2-120a_{21}{\overline{v}}_2-180a_{21}^2{\overline{v}}_2+5w_1+5w_2-30a_{21}^2{w}_2+60a_{21}^3{w}_2-\hfill \\ & 60a_{21}{\widehat{w}}_2+120a_{21}^2{\widehat{w}}_2-60a_{21}^3{\widehat{w}}_2-60{\overline{w}}_2+240a_{21}{\overline{w}}_2-180a_{21}^2{\overline{w}}_2-2\theta )/60.\hfill \end{array}$$

(27)

The coefficients matrices of the optimal ThDTSRK${}_{25}$ scheme are derived by

$$\begin{array}{cc}\hfill \theta =& 0,\mathit{c}=\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}0\\ 0.1983891070202614\end{array}\right],\hfill \\ \hfill \mathit{A}=& \left[\begin{array}{cc}0& 0\\ c_2& 0\end{array}\right],\widehat{\mathit{A}}=\left[\begin{array}{cc}0& 0\\ \frac{c_2^2}{2}& 0\end{array}\right],\overline{\mathit{A}}=\left[\begin{array}{cc}0& 0\\ \frac{c_2^3}{6}& 0\end{array}\right],\hfill \\ \hfill {\mathit{v}}^T& =\left[\begin{array}{cc}0.4988123289876567& -0.1677439748133182\end{array}\right],\hfill \\ \hfill {\mathit{w}}^T& =\left[\begin{array}{cc}0.5011876710123433& 0.1677439748133182\end{array}\right],\hfill \\ \hfill {\widehat{\mathit{v}}}^T& =\left[\begin{array}{cc}-0.0958493173039603& 0.6579633161995648\end{array}\right],\hfill \\ \hfill {\widehat{\mathit{w}}}^T& =\left[\begin{array}{cc}-0.8843764374259575& 1.4911940843560145\end{array}\right],\hfill \\ \hfill {\overline{\mathit{v}}}^T& =\left[\begin{array}{cc}-0.0202481631489146& 0.1199846505868748\end{array}\right],\hfill \\ \hfill {\overline{\mathit{w}}}^T& =\left[\begin{array}{cc}-0.1160041365433313& 0.0621952996182998\end{array}\right].\hfill \end{array}$$

(2) 

Two-Stage Sixth-Order ThDTSRK Scheme (ThDTSRK₂₆)

The coefficients of the general two-stage sixth-order ThDTSRK method are then given as a five-parameter system, depending on $\theta ,a_{21},{\widehat{w}}_2,{\overline{v}}_2,{\overline{w}}_2,$ which are related through

$$\begin{array}{cc}\hfill v_2=& 0;w_2=0;{\widehat{v}}_2=-{\widehat{w}}_2;{\widehat{a}}_{21}=a_{21}^2/2;{\overline{a}}_{21}=a_{21}^3/6;\hfill \\ \hfill v_1=& (15-120a_{21}{\overline{v}}_2-360a_{21}^2{\overline{v}}_2-240a_{21}^3{\overline{v}}_2+240a_{21}^3{\widehat{w}}_2-120a_{21}{\overline{w}}_2+\hfill \\ \hfill & 360a_{21}^2{\overline{w}}_2-240a_{21}^3{\overline{w}}_2+\theta )/2;\hfill \\ \hfill w_1=& 1+\theta -v_1;\hfill \\ \hfill {\widehat{v}}_1=& (-31+360a_{21}{\overline{v}}_2+960a_{21}^2{\overline{v}}_2+600a_{21}^3{\overline{v}}_2+10{\widehat{w}}_2-60a_{21}^2{\widehat{w}}_2-600a_{21}^3{\widehat{w}}_2+\hfill \\ \hfill & 240a_{21}{\overline{w}}_2-840a_{21}^2{\overline{w}}_2+600a_{21}^3{\overline{w}}_2-\theta )/10;\hfill \\ \hfill {\widehat{w}}_1=& (-29+240a_{21}{\overline{v}}_2+840a_{21}^2{\overline{v}}_2+600a_{21}^3{\overline{v}}_2-10{\widehat{w}}_2+60a_{21}^2{\widehat{w}}_2-600a_{21}^3{\widehat{w}}_2+\hfill \\ \hfill & 360a_{21}{\overline{w}}_2-960a_{21}^2{\overline{w}}_2+600a_{21}^3{\overline{w}}_2+\theta )/10;\hfill \\ \hfill {\overline{v}}_1=& (111-120{\overline{v}}_2-1080a_{21}{\overline{v}}_2-2160a_{21}^2{\overline{v}}_2-1200a_{21}^3{\overline{v}}_2+120a_{21}{\widehat{w}}_2+\hfill \\ \hfill & 360a_{21}^2{\widehat{w}}_2+1200a_{21}^3{\widehat{w}}_2-360a_{21}{\overline{w}}_2+1440a_{21}^2{\overline{w}}_2-1200a_{21}^3{\overline{w}}_2+\theta )/120;\hfill \\ \hfill {\overline{w}}_1=& (-49+360a_{21}{\overline{v}}_2+1440a_{21}^2{\overline{v}}_2+1200a_{21}^3{\overline{v}}_2-120a_{21}{\widehat{w}}_2+360a_{21}^2{\widehat{w}}_2-\hfill \\ \hfill & 1200a_{21}^3{\widehat{w}}_2-120{\overline{w}}_2+1080a_{21}{\overline{w}}_2-2160a_{21}^2{\overline{w}}_2+1200a_{21}^3{\overline{w}}_2+\theta )/120.\hfill \end{array}$$

(28)

The coefficients matrices of the optimal ThDTSRK${}_{26}$ scheme are given by

$$\begin{array}{cc}\hfill \theta =& 0,\mathit{c}=\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}0\\ 0.5873258965737987\end{array}\right],\hfill \\ \hfill \mathit{A}=& \left[\begin{array}{cc}0& 0\\ c_2& 0\end{array}\right],\widehat{\mathit{A}}=\left[\begin{array}{cc}0& 0\\ \frac{c_2^2}{2}& 0\end{array}\right],\overline{\mathit{A}}=\left[\begin{array}{cc}0& 0\\ \frac{c_2^3}{6}& 0\end{array}\right],\hfill \\ \hfill {\mathit{v}}^T& =\left[\begin{array}{cc}1.0471220060600115& 0\end{array}\right],{\mathit{w}}^T=\left[\begin{array}{cc}-0.0471220060600116& 0\end{array}\right],\hfill \\ \hfill {\widehat{\mathit{v}}}^T& =\left[\begin{array}{cc}0.4467995963745828& 0.1411691523070592\end{array}\right],\hfill \\ \hfill {\widehat{\mathit{w}}}^T& =\left[\begin{array}{cc}0.0060783975654054& -0.1411691523070592\end{array}\right],\hfill \\ \hfill {\overline{\mathit{v}}}^T& =\left[\begin{array}{cc}0.0482868172625281& 0.0243580486114999\end{array}\right],\hfill \\ \hfill {\overline{\mathit{w}}}^T& =\left[\begin{array}{cc}0.0052528132887524& -0.0227607642077618\end{array}\right].\hfill \end{array}$$

(3) 

Two-Stage Seventh-Order ThDTSRK Scheme (ThDTSRK₂₇)

The coefficients of the general two-stage seventh-order ThDTSRK method are then given as a two-parameter system, depending only on $\theta ,a_{21}$, which are related through

$$\begin{array}{cc}\hfill v_2=& 0;w_2=0;{\widehat{v}}_2=0;{\widehat{w}}_2=0;{\widehat{a}}_{21}=a_{21}^2/2;{\overline{a}}_{21}=a_{21}^3/6;\hfill \\ \hfill v_1=& (105-627a_{21}+1050a_{21}^2+7\theta +3a_{21}\theta +70a_{21}^2\theta )/\left(14(1+10a_{21}^2)\right);\hfill \\ \hfill w_1=& (-91+627a_{21}-910a_{21}^2+7\theta -3a_{21}\theta +70a_{21}^2\theta )/\left(14(1+10a_{21}^2)\right);\hfill \\ \hfill {\widehat{v}}_1=& (-45+627a_{21}-868a_{21}^2-3\theta -3a_{21}\theta -28a_{21}^2\theta )/\left(28(1+10a_{21}^2)\right);\hfill \\ \hfill {\widehat{w}}_1=& (-123+627a_{21}-812a_{21}^2+3\theta -3a_{21}\theta +28a_{21}^2\theta )/\left(28(1+10a_{21}^2)\right);\hfill \\ \hfill {\overline{v}}_1=& (-209+300a_{21}-6270a_{21}^2+15540a_{21}^3+\theta +20a_{21}\theta +\hfill \\ \hfill & 30a_{21}^2\theta +140a_{21}^3\theta )/\left(1680(a_{21}+10a_{21}^3)\right);\hfill \\ \hfill {\overline{v}}_2=& (209-\theta )/\left(1680(a_{21}+10a_{21}^3)\right);{\overline{w}}_2=-{\overline{v}}_2;\hfill \\ \hfill {\overline{w}}_1=& (209-1940a_{21}+6270a_{21}^2-6860a_{21}^3-\theta +20a_{21}\theta -\hfill \\ \hfill & 30a_{21}^2\theta +140a_{21}^3\theta )/\left(1680(a_{21}+10a_{21}^3)\right).\hfill \end{array}$$

(29)

The coefficient matrices of the optimal ThDTSRK${}_{27}$ scheme are obtained by

$$\begin{array}{cc}\hfill & \theta =0,\mathit{c}=\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}0\\ \frac{1}{2}\end{array}\right],\mathit{A}=\left[\begin{array}{cc}0& 0\\ c_2& 0\end{array}\right],\widehat{\mathit{A}}=\left[\begin{array}{cc}0& 0\\ \frac{c_2^2}{2}& 0\end{array}\right],\overline{\mathit{A}}=\left[\begin{array}{cc}0& 0\\ \frac{c_2^3}{6}& 0\end{array}\right],\hfill \\ \hfill & {\mathit{v}}^T=\left[\begin{array}{cc}\frac{54}{49}& 0\end{array}\right],{\mathit{w}}^T=\left[\begin{array}{cc}-\frac{5}{49}& 0\end{array}\right],{\widehat{\mathit{v}}}^T=\left[\begin{array}{cc}\frac{103}{196}& 0\end{array}\right],\hfill \\ \hfill & {\widehat{\mathit{w}}}^T=\left[\begin{array}{cc}-\frac{25}{196}& 0\end{array}\right],{\overline{\mathit{v}}}^T=\left[\begin{array}{cc}\frac{79}{735}& \frac{209}{2940}\end{array}\right],{\overline{\mathit{w}}}^T=\left[\begin{array}{cc}-\frac{17}{980}& -\frac{209}{2940}\end{array}\right].\hfill \end{array}$$

To evaluate the performance of the designed methods, a comparison was conducted between high-order ThDTSRK methods, TDRK methods [13], and ThDRK methods [14]. The absolute stability regions of high-order ThDTSRK methods are displayed in Figure 2, Figure 3 and Figure 4. We define the scaled intersection ${\mathcal{L}}^{*}$ of absolute stability, which is derived by

$${\mathcal{L}}^{*}=\frac{p}{\gamma s}\mathcal{L},$$

(30)

where $\gamma =2$ for two-derivative schemes and $\gamma =3$ for three-derivative schemes. In

Table 2. The intersection $\mathcal{L}$ and scaled intersection ${\mathcal{L}}^{*}$ of absolute stability for high-order ThDTSRK${}_{sp}$, ThDRK${}_{sp}$, and TDRK${}_{sp}$ schemes.

	ThDTSRK${}_{25}$	ThDTSRK${}_{26}$	ThDTSRK${}_{27}$
$\mathcal{L}$	$-8.181$	$-6.266$	$-3.610$
${\mathcal{L}}^{*}$	$-6.818$	$-6.266$	$-4.212$
	ThDRK${}_{25}$	ThDRK${}_{36}$	ThDRK${}_{37}$
$\mathcal{L}$	$-3.990$	$-7.263$	$-5.213$
${\mathcal{L}}^{*}$	$-3.325$	$-4.842$	$-4.055$
	TDRK${}_{35}$	TDRK${}_{46}$	TDRK${}_{57}$
$\mathcal{L}$	$-5.144$	$-4.063$	$-4.134$
${\mathcal{L}}^{*}$	$-4.287$	$-3.047$	$-2.894$

, the results indicate that the ThDTSRK methods exhibit a larger scaled intersection ${\mathcal{L}}^{*}$ compared to other schemes of the same order.

5. Numerical Examples

5.1. Prothero–Robinson Problem

We test the order accuracy of ThDTSRK methods using the Prothero–Robinson problem [13], which is expressed as

$$y^{\prime }\left(x\right)=\lambda (y\left(x\right)-\sin \left(x\right))+\cos \left(x\right),\lambda =-10.$$

(31)

The exact solution is given by $y\left(x\right)=\sin \left(x\right)$ and the results are integrated to $2.8\pi$. In Figure 5, the results clearly demonstrate that the ThDTSRK methods attain the desired order accuracy. Additionally, the ThDTSRK methods are more efficient than the classical high-order Runge–Kutta methods [11,13].

5.2. Kaps Problem

We test the order accuracy of ThDTSRK methods employing the Kaps problem [13], which is expressed as

$$y^{\prime }\left(x\right)=\left[\begin{array}{c}-y_1\left(1+y_1\right)+y_2\\ \lambda \left(y_1^2-y_2\right)-2y_2\end{array}\right],y\left(0\right)=\left[\begin{array}{c}1\hfill \\ 1\hfill \end{array}\right],\lambda =10.$$

(32)

The exact result is $y\left(x\right)=[e^{(-x)},e^{(-2x)}]$, and the solutions integrated up to five are presented in Figure 6. These results also reveal that the ThDTSRK methods attain the designed order accuracy. The ThDTSRK methods are more efficient compared with the classical high-order Runge–Kutta methods in solving ordinary differential equations.

6. Conclusions

In this study, we developed the order conditions for the ThDTSRK schemes that ensure these methods can achieve theoretical order accuracy. The parameters of the general high-order ThDTSRK methods were provided by utilizing the order conditions. Additionally, we established the A-stability property theory for ThDTSRK methods to determine the optimal stability coefficients. By comparing the ThDTSRK methods with ThDRK and TDRK schemes, the results revealed that the ThDTSRK schemes can maintain the A-stability property and attain the desired order accuracy while utilizing fewer stages. This indicates that the ThDTSRK methods are computationally efficient for solving ordinary differential equations. In further studies, our proposed ThDTSRK schemes can be applied to partial differential equations which are reduced to ordinary differential equation system using various spatial discretization schemes [5,22]. The TDMSRK methods can demonstrate high computational efficiency, especially when employing high-order spatial discretization schemes.