Two-point block variable order step size multistep method for solving higher order ordinary differential equations directly

1 Introduction

Problems in science and engineering are often modeled in the form of ordinary, partial and fractional differential equations depending on the type of problem. These problems may range from various field such as the interaction between chemicals as found in chemical kinetics to the mechanics of fluids which are often modelled in medicine and engineering. For recent analytical and numerical methods established for solving these differential applications, readers may refer to literature by Al-Jawary et al. (2020), Khataybeh et al. (2019), Mohd Ijam et al. (2020), Zhang et al. (2020) and Rasedee et al. (2021) for ordinary differential equations, Colbrook et al. (2019), Fu et al. (2019) and Lehrenfeld and Olshanskii (2019) for partial differential equations and Goswami et al. (2019), Singh et al. (2019), Srivastava et al. (2020) and Veeresha et al. (2020) for fractional differential equations. The algorithm development provided in this research is to specifically cater a numerical method for solving higher order ordinary differential equations directly using variable order step size formulations. We begin the study with the evolution of the variable order step size algorithm based on the multistep method.

Previous numerical methods used to solve these higher order ordinary differential equations (ODEs) by reducing them to systems of first order ODEs were so efficient such that methods for solving ODEs directly were considered as robust. These methods were almost overlooked until (Krogh, 1968) revived the field of study with a modified version of the divided difference approach. He proposed that the back values of any point of the derivative were to be interpolated. In a $d^{\mathit{th}}$ order ODE, the method used to interpolate the highest derivative value was referred as the Direct Integration (DI). Krogh (1968) provide a comparison for two second order problems which was similar to Gear (1967) whom used a Nordsieck version of the multistep method. Suleiman (1979) extended works of Krogh by adapting the DI for estimating higher order derivatives ( $d^t$ h order) which corresponds with nonstiff problems and a generalized backward differentiation method (GBDF) method for lower derivatives which corresponds with stiff problems. Omar and Suleiman (2005) adapted Suleiman’s DI algorithm with a two-point explicit block algorithm. By doing so, Omar was able to reduce computational cost by providing a more efficient divided difference code. In efforts of solving higher order ODEs directly, Rasedee (2009) proposed a direct method for solving these higher order derivatives using a backward difference formulation. The backward difference method suggested by Rasedee was in effort to overcome the drawback of the DI method, as it requires calculating the integration coefficients at every step change which proves to be costly. The method proposed by Rasedee (2009) requires calculating the integration coefficients only once in the beginning of the algorithm and once more at the end, if deemed necessary.

This article is part of a research series influenced by the works of Suleiman (1979). He initiated a series of research based on Adams type difference formulation for solving ordinary differential equations (ODEs). In his research, Suleiman established the Direct Integration method (DI) for solving both, stiff and nonstiff initial value higher order ODEs directly. This method was formulated using a multistep method in divided difference form combined with variable order stepsize (VOS) capability. Using works by Suleiman as their foundation, authors such as Omar and Suleiman (2005), Majid and Suleiman (2006, 2007), Suleiman et al. (2011), Mohd Ijam et al. (2014) and Waeleh and Abdul Majid (2016) and Rasedee et al. (2016) have produced an extended variation for solving various types of differential problems.

The current research proposes a two-point block formulation in backward difference form (2PVOSBD) for solving both artificial and real life higher order ODE problems. Equipped with a VOS algorithm similar to Rasedee et al. (2014), the proposed method calculates the integration coefficients only once at the beginning and once at the end (if required) in contrast to the DI which requires calculating integration coefficients at every step change. The 2PVOSBD only requires calculating the integration coefficients a second time if $H⩾T_{\mathit{end}}-T_b$ , where H as the current stepsize, $T_{\mathit{end}}$ as the last point of the interval, and $T_b$ as the current point.

When implementing a variable order stepsize strategy, the acceptance criteria is key. An unsuitable criteria may result in an increase in the total of steps or the loss of accuracy. The selection criteria is discussed in the upcoming sections. Following the appropriate selection criteria, strategies suggested by Lambert (1973) are used and integrated with the variable step size code proposed by Krogh (1973).

The current research also provides the discussion on the order, stability and convergence of the proposed method. For the variable order step size techniques and selection criterion implemented in this research, authors adapted techniques suggested by Gear (1967), Krogh (1968), Lambert (1973), Hall and Watt (1976) and Suleiman (1979). And in more recent works for solving higher order ODEs found in Zainuddin et al. (2016), Rasedee et al. (2017b, 2018a), Ibrahim et al. (2019) and Asnor et al. (2019) and others. The following are preliminary definitions that are used to determine order, stability and convergence of the 2PVOSBD.

2 Preliminaries

This section entails definitions used to determine the order, stability and convergence of the method.

Choose a function $f\left(t\right)$ to be differentiable as often as needed, and expand $f(t+\mathit{ih})$ and $f^{\prime }(t+\mathit{ih})$ with respect to t, and arrange as $$L\left[f\right(t\left)\text{;}h\right]≔C_{0}f\left(t\right)+C_1{\mathit{hf}}^{\prime }\left(t\right)+\cdots +C_q{h}^q{z}^{\left(q\right)}\left(t\right)+\cdots$$ where $C_i,i=0,1.\dots ,q,\dots$ are constants.

3 The higher order ordinary differential equation

Differential equations have been widely used to model various types of real-life applications, from natural phenomenon to man-made machinery. Examples of these higher order ODEs can found in classical mechanics of the two-body problem, chaotic motions in periodic-self exited oscillators to the elasticity of thin clamp plates.

The focus of this research is to establish a numerical method for solving higher order initial value ordinary differential equations. For the purpose of this research, consider the following higher order ODE

Direct methods versus reduction to first order systems have been extensively discussed by Rutishauser (1960), Collatz (1960), Henrici (1962) and Gear (1967). These opinions may vary based on individual preferences. In the case of Henrici (1962), the author favours reduction methods due to the large global error of direct methods even though lower order derivatives have small local truncation error, whereas Collatz (1960) and Gear (1967) prefer direct methods. Rutishauser (1960) provided proof that the nth linear order ODE was solved using a direct method (Euler extended method) had a global error of order 1, which was equal to the global error of the reduction method (Euler method). A N order direct method algorithm requires only storing and updating one set of back values, where a N reduction method needs N set of back values. Due to the fact that the accuracy of each derivative is of the same order thus, the direct method requires controlling the local error of only ${(N-1)}^{\mathit{th}}$ derivative because it offers the lowest accuracy $O\left(h\right)$ opposed to controlling the local error of each derivative as required by the reduction method.

The method proposed in this research approximates (3.1) directly with minimal loss of accuracy, if not any. The succeeding section will provide detailed derivation of the proposed method.

4 The two-point block method

This section explains derivation of both first and second point of the two-point block method. To formulate a two-point predictor–corrector (PeCe) block method with a variable order step size backward difference algorithm, elements such as explicit–implicit integration coefficients and order stepsize strategy are necessary.

Firstly, note that the 2PVOSBD method is formulated in predictor–corrector mode similar to the Adam-Basthforth-Moulton formulation. The derivation begins with integrating the higher order ODE as defined in (3.1). By integrating (3.1) once, $f^{\left(n\right)}$ becomes

Predictor:

Corrector:

Explicit: $${\gamma }_{b,j,i}^p={(-1)}^i{\int }_0^b\frac{{(b-s)}^{d-1}}{(d-1)!}\left(\begin{array}{c}\hfill -s\hfill \\ \hfill i\hfill \end{array}\right)\mathit{ds},$$

Implicit: $${\gamma }_{b,j,i}^c={(-1)}^i{\int }_{-b}^0\frac{{(-s)}^{d-1}}{(d-1)!}\left(\begin{array}{c}\hfill -s\hfill \\ \hfill i\hfill \end{array}\right)\mathit{ds}\text{.}$$

5 Integration coefficients

The current section provides derivation for explicit and implicit coefficients. Firstly, let the first order explicit, $G_{b,1}^p\left(t\right)$ and implicit $G_{b,1}^c\left(t\right)$ generating functions as $$G_{b,1}^p\left(t\right)={\displaystyle \sum _{i=0}^{\infty }}{\gamma }_{1,d,i}^p{t}^i$$ and $$G_{b,1}^c\left(t\right)={\displaystyle \sum _{i=0}^{\infty }}{\gamma }_{1,d,i}^c{t}^i\text{.}$$ Next, by substituting ${\gamma }_{b,d,i}^p$ into $G_{\left(d\right)}^p\left(t\right)$ , the explicit generating function can be rewritten as $$G_{b,1}^p\left(t\right)={\displaystyle \sum _{i=0}^{\infty }}{(-t)}^i{\int }_0^b(b-s)\left(\begin{array}{c}-s\hfill \\ i\hfill \end{array}\right)\mathit{ds}$$ which can be simplified as

Explicit Coefficients: $${\gamma }_{b,d,0}^p={\gamma }_{b,d-1,1}^p,{\gamma }_{b,d,k}^p={\gamma }_{b,d-1,k+1}^p-{\displaystyle \sum _{i=0}^{k-1}}\left(\frac{{\gamma }_{b,d,i}^p}{k-i+1}\right),k=1,2,\dots \text{.}$$

Implicit Coefficients: $${\gamma }_{b,d,0}^c={\gamma }_{b,d-1,1}^c,{\gamma }_{b,d,k}^c={\gamma }_{b,d-1,k+1}^c-{\displaystyle \sum _{i=0}^{k-1}}\left(\frac{{\gamma }_{b,d,i}^c}{k-i+1}\right),k=1,2,\dots \text{.}$$ Inline with the aim of this research to reduce computational cost, a recursive relationship between explicit and implicit coefficients is established. The relationship obtained begins in terms of the generating functions. Referring to Eqs. (5.10) and (5.11), which are the first order generating functions

6 Order and stability

Various researchers have studied conditions of zero stability and the order of a block method which can be found in works of Watts and Shampine (1972), Mohd Ijam et al. (2018), Mohd Ijam and Ibrahim (2019) and Rasedee et al. (2021a) and etc. Discussions on order and stability of the proposed method applies definitions which are similar to Ola Fatunla (1991).

7 Convergence of the backward difference method

As previously shown by Suleiman (1993), there are certain conditions that need to be met to show the convergence of the backward difference method.

In the previous section, both the predictor and corrector satisfies conditions in Definition 2.4, making the method zero stable; whereas the order method as previously shown in Definition 2.3, with the predictor $$C_0=C_1=C_2=\cdots =C_7=\left(\begin{array}{c}0\hfill \\ 0\hfill \end{array}\right),C_8\ne \left(\begin{array}{c}0\hfill \\ 0\hfill \end{array}\right)$$ and the corrector $$C_0=C_1=C_2=\cdots =C_8=\left(\begin{array}{c}0\hfill \\ 0\hfill \end{array}\right)C_9\ne \left(\begin{array}{c}0\hfill \\ 0\hfill \end{array}\right)$$ are consistent with Definition 2.5 of order of 8 and 9 respectively. Satisfying both conditions $i)$ and $\mathit{ii})$ of Theorem 7.1 shows that the method to be convergent.

8 Variable order step size strategies

The literature on the variable order stepsize codes show an incomplete theory with numerous strategies of order and stepsize selection are given. Firstly, varying the order. The variable strategies are dependent on the amount of back values stored. With the appropriated amount of back values, the order can simply be increased by retaining the values of previous steps and can be decreased by disregarding unnecessary values. Previous researches have shown that the order strategies which are unbiased to lower order is efficient in an Adams based code for solving nonstiff problems.

As mentioned earlier, the acceptance criteria in a variable order strategy is crucial. In this article, the second step evaluation is determined on whether the estimated error, $$E_k^{(d-p)}=h^d{\gamma }_{b,d-p,k}^c{\nabla }^k{\varphi }_{n+1}$$ satisfies the local accuracy requirement, $${\mathrm{\Theta }}_{n+1}\left|E_k^{(d-p)}\right|<\mathit{TOL},{\mathrm{\Theta }}_{n+1}=\frac{1}{\left(A+B+P_n\right)}$$ where A and B determine the type of error test, which is discussed later in the numerical results section. For reason of simplicity, the error estimate $E_k^{(d-p)}$ is denoted as $E_k$ , where $k=n-2,n-1,n,n+1$ . The order condition is set for $n>2,\max \{\left|E_{n-1}\right|,\left|E_{n-2}\right|\}⩽\left|E_n\right|$ and if both conditions $\left|E_{n-1}\right|⩽0.5\left|E_n\right|$ and $n=2$ are satisfied, the order is reduced by 1. If $E_{n+1}$ is available, the condition to reduce the order is when $n>1$ and $\left|E_{n-1}\right|⩽\min \{\left|E_n\right|,\left|E_{n+1}\right|\}$ is satisfied. For the latter, the algorithm must first complete $n+1$ successful steps in constant stepsize, only then for $n>1$ where $\left|E_{n+1}\right|<\left|E_n\right|<\max \{\left|E_{n-1}\right|,\left|E_{n-2}\right|\}$ , and both $\left|E_{n+1}\right|<0.5\left|E_n\right|$ and $n=1$ are satisfied will the order be increased.

As for varying the stepsize, this research adopts the halving and doubling of the stepsize algorithm suggested by Krogh (1973). Although doubling the stepsize might seem necessary because the previous back values of ${\varphi }_{n-2},{\varphi }_{n-4},\dots ,{\varphi }_{n-2k+2}$ can be used, Krogh argues that such technique will contribute to the loss of accuracy. In fact, doubling the stepsize after every evaluation may also have a toll on the accuracy. This is when halving the stepsize is required. In halving the step size algorithm, new back values are obtained. The use of back values in the form of ${\varphi }_{n-\frac{1}{2}},{\varphi }_{n-1},\dots ,{\varphi }_{n-\frac{k}{2}+\frac{1}{2}}$ has shown to be more accurate. For more details on the algorithm, reader may refer to Rasedee (2009) and Rasedee et al. (2014). Fig. 1 illustrates the general flow of the proposed 2PVOSBD algorithm.

Close

Fig. 1

The general flowchart for the 2PVOSBD algorithm.

Export to PPT

9 Numerical results and discussions

As numerical method becomes more robust, the accuracy of a method alone is no longer sufficient hence, the need for efficiency. When dealing with complex problems that require massive calculation, computational cost becomes a concern. The current research highlights the efficiency of proposed method. With little to none loss of accuracy, the 2PVOSBD offers an accuracy to total step ratio which proves to be efficient, similar as proposed in Rasedee (2009).

The 2PVOSBD method is tested with a combination of linear and nonlinear higher order initial value ODEs. The selected test problems are also inclusive of artificial and real-life, practical ODEs. The approximation obtained by the 2PVOSBD was compared with approximations of 1PVOSBD and DI methods which were reconstructed for the purpose of this research. Due to its divided difference based algorithm, the DI method was considered as a suitable counterpart of the 2PVOSBD backward difference based algorithm. The comparison is based on accuracy, total steps and efficiency of the methods. To compare approximated solution, the error solution in this research is defined by $${\mathrm{Err}}_n=\left|\frac{{\left(f_i\right)}_n-f{\left(t_i\right)}_n}{A+B\left(f{\left(t_i\right)}_n\right)}\right|$$ where ${\left(f_i\right)}_n$ is the $n^{\mathit{th}}$ component of the exact solution and $f{\left(t_i\right)}_n$ is the $n^{\mathit{th}}$ component of the approximated solution of f. There are 3 types of error test provided in the algorithm, where the most suitable estimation is selected. The choice of error tests are, absolute error when the corresponding $A=1,B=0$ , relative error when $A=0,B=1$ . and mixed error when $A=1,B=1$ . For the current research, the mixed error test was deemed best suited. On the other hand, the maximum error (MaxErr) is defined as $$\mathrm{MaxErr}={\displaystyle \underset{1⩽j⩽\mathit{Succstp}}{\max }}\left({\displaystyle \underset{1⩽j⩽B⩽N}{\max }}\left({\mathrm{Err}}_n\right)\right)$$ and the average error (AvErr) $$\mathrm{AvErr}=\frac{{\sum }_{i=1}^{\mathrm{SuccStp}}{\sum }_{i=1}^N{\sum }_{i=1}^B{\mathrm{Err}}_n}{\left(B*N*\mathrm{SuccStp}\right)}$$ with N as the number of equations and B as the number of blocks. Contrary, the efficiency of the methods is defined as the ratio of accuracy compared to the corresponding total steps. These are the abbreviations that are used in the current section:

The following are test problems that were selected for comparisons.

Tables 3 and 4 are the approximated results for Problems 1–4 using the 2PVOSBD, but with the absence of the variable order stepsize component. The result presented compare the constant stepsizes, H $={10}^{-1}$ , H $={10}^{-2}$ , H $={10}^{-3}$ , H $={10}^{-4}$ and H $={10}^{-5}$ . Starting with Euler’s method, the algorithm increase the back values (order) each step until the total of 12 back values were obtained. Then, the algorithm completed its estimation using 12 back values for the remaining steps.