\(\S \) This work has been supported in part by grant no. PN II RU TE 2011-3-0013

\(^\ast \)Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, P.O. Box 68–1, \(400110\), Cluj-Napoca, Romania, e-mail: fpatrulescu@ictp.acad.ro.

1 Introduction

Consider a scalar initial value problem (IVP):

where: \(I\subseteq \mathbb {R}\), \(y_0\in \mathbb {R}\), \(f:I\times \mathbb {R}\rightarrow \mathbb {R}\) and \(x_0\in I\). We assume that \(I=[x_0,x_0+T]\), \(0{\lt}T{\lt}\infty \) and the function \(f\) satisfies all requirements necessary to insure the existence of a unique solution \(y\) on the bounded interval \(I\), see [ 1 ] , [ 2 ] for details.

In this paper we present a numerical method to approximate the solution of IVP (3) using a particular case of a nonlinear interpolation formula analyzed in [ 5 ] and [ 7 ] . Another particular case was utilized in [ 4 ] .

For a \(5\)-times differentiable function \(h:I\rightarrow \mathbb {R}\) consider the function \(g:I\rightarrow \mathbb {R}\) given by

where \(\alpha _1=\tfrac {3-\sqrt{3}}{6}\) and \(\alpha _2=\tfrac {3+\sqrt{3}}{6}\).

It is showed in [ 5 ] that the function \(g\) verifies

where \(\displaystyle M_5=\sup _{x\in I}|h^{(5)}(x)|\). Also, we observe that the coefficients \(\alpha _1\), \(\alpha _2\) satisfy the equalities

In the next sections, for simplicity, we consider only the autonomous case, i.e. \(f=f(y)\). The general case can be treated using similar arguments.

The paper is structured as follows. In Section \(2\) we derive the numerical method and in Section \(3\) we analyze its versions corresponding to the particular cases of the parameter. In Section \(4\) numerical examples are provided to exemplify the results obtained in the previous sections.

2 The numerical method

The interval \(I=[x_0,x_0+T]\) is partitioned by the point set \(\{ x_n|n=\overline{0,N}, N\in \mathbb {N}\} \), where \(x_n\) are given by a prescribed rule, and we denote by \(y_n\) a numerical approximation to the exact solution \(y\) of (3) at \(x_n\).

We suppose that the exact solution \(y\) of initial value problem (3) is \(5\)-times differentiable at least. Conditions for regularity of exact solution of an initial value problem can be found in [ 2 ] .

Using the results presented above we deduce that there exists an approximation \(\tilde{y}\) of \(y\) given by

for all \(x\in I\). From (??) we deduce that this approximation verifies

where \(\displaystyle M_5=\sup _{x\in I}|y^{(5)}(x)|\).

The unknown quantities \(y(x_0+\alpha _i(x-x_0))\), \(i=1,2\), in (??) can be approximated in the same manner and we continue this approximation procedure for the next unknown values \(y(x_0+\alpha ^p_1\alpha ^q_2(x-x_0))\), \(p,\, q=1,2,\ldots \quad \).

Using the notation

the algorithm described above can be written for the firsts \(p\) steps in the following compact form

and we can define \(\tilde{y}(x)=u_{00}(x)\), \(x\in I\).

Taking into account that

we can choose \(p=p_0\) such that \(y(x_0)\) approximates \(u_{p_00}\), \(u_{p_0-10},\ldots , u_{0p_0}\) with any given accuracy, since \(y\) is a continuous function. From (??) we obtain

We use this recurrence relation to construct a numerical method for the scalar initial value problem (3). Replacing \(x\) by \(x_1\) in (??) we obtain an approximation \(\tilde{u}_{00}(x_1)\) for the exact value \(y(x_1)\). We denote this approximation by \(y_1\) and we apply again the algorithm (??) for \(x_2\), but, instead of \(y(x_0)=y_0\), we consider the value \(y_1\) previously computed. We repeat this procedure for each \(x_n\), \(n=\overline{1,N}\).

Using the notation

the above algorithm can be written in the following way

where \(h_n=x_{n+1}-x_n\), \(n=\overline{0,\, N-1}\).

For \(p_0=1\) we obtain the Euler forward method given by

We have the following equivalence result.

the approximation for the new point \(x_{n+1}=x_n+h_n\) is given by

where

It easy to see that (14) provides for autonomous case and for the above Butcher array a rule of the form (12), which concludes the proof.

3 Method (12) when \(p_0=2\), \(3\) and \(4\)

For \(p_0=2\) we obtain

The method (17) can be written in a more suitable form

and can be identified as a \(3\)–stages explicit Runge-Kutta method with the Butcher array given by

As in [ 6 ] , we suppose that

where \(\| f\| =\rm {sup}\{ |f(\it {t})|: \it {t}\in I\} \) and \(M\), \(L\) are positive real numbers.

The convergence order of the method (17) is provided in the following result.

where \(z\) is an arbitrary function defined on \(I\), sufficiently differentiable and \(z'(x)=f(z(x))\), for all \(x\in I\).

Expanding in Taylor series with respect to \(x\) and using (??) we obtain

Then, using the definition for the convergence order given in [ 3 ] we deduce that the method has second-order of accuracy. Also, substituting \(z\) by the exact solution \(y\), \(x\) by \(x_n\) and supposing localizing assumption \(y_i=y(x_i)\), \(i=\overline{1,\, n}\) the local truncation error can be written as

Moreover, the coefficient of principal local truncation error is given by

Next, using (??) we have for \(C_3\) the following bound

which concludes the proof.

Following [ 2 ] , for the Runge-Kutta methods in addition of the row-sum condition, \(Ae=Ce\), and consistency condition, \(b^Te=1\), the sufficient conditions for order \(2\) are given by

where \(e=[1,\ldots ,1]\) and \(C=\mbox{diag}(c)\). A simple calculus shows that the method (17) verifies conditions (??) and it represents a validation for (19).

The stability function and absolute stability region of the method (17) are provided in the following result.

and we obtain the difference equation

Denoting \(z=\lambda h_n\) and using (??) we obtain the stability function

Note that this function is the same as the stability function for \(2\) stages explicit Runge-Kutta methods of order \(2\).

The absolute stability region is given by

and it is plotted using scanning technique in Figure 1. The computations have been carried out using Matlab software.

For \(p_0=3\) we obtain the method

which can be identified as a \(6\)–stages explicit Runge-Kutta method with the Butcher array given by

The convergence order of the method (25) is provided in the following result.

Therefore we deduce that this method has convergence order \(3\). Also, we can determine the local truncation error as

and the coefficient of principal local truncation error is given by

Using (??) we obtain \(\| C_4\| {\lt}\tfrac {1}{24}ML^3\), which concludes the proof.

A simple calculus shows that the method (25) verifies the sufficient conditions for order \(3\), see [ 2 ] for details,

and it represents a validation for (??).

We observe that method (25) has the same stability function as \(3\) stages explicit Runge-Kutta methods of order \(3\). The absolute stability region defined by

is plotted in Figure 1.

For \(p_0=4\) we obtain the method

and can be identified as a \(10\)–stages explicit Runge-Kutta method with the Butcher array given by

Using (??) and similar arguments to those used above we deduce that the method (31) has convergence order \(4\). The stability function is given by

and we observe that it is the same as stability function for \(4\) stages explicit Runge-Kutta methods of order \(4\). The absolute stability region defined by

is plotted in Figure 1.

We restrict ourselves to the cases when \(p_0\) takes the values \(2\), \(3\) and \(4\). For \(p_0\geq 5\) methods obtain from (12) have a form much more complicated with a high cost of calculus and the results concerning accuracy are not so outstanding. From the next theorem we deduce that the convergence order is not grater than \(4\).

where \(z\) is an arbitrary function defined on \(I\), \(5\)-times differentiable at least and \(z'(x)=f(z(x))\), for all \(x\in I\).

Next, using Taylor series with respect to \(x\) and (??) we obtain

Then, from the definition of convergence order given in [ 3 ] we deduce that for \(p_0\geq 5\) method (12) has convergence order \(4\), which concludes the proof.

Note that the method (12) is a zero-stable method because it verifies root-condition. Also, since the convergence order is grater than \(2\) for all \(p_0\geq 2\) we conclude that it satisfies the consistency condition. It follows that our method represents a convergent method, see [ 3 ] for details.

4 Numerical examples

We consider the initial value problems.

We apply the methods (17), (25) and (31) with constant steplength to determine numerical solutions for the above examples.

As a measure of the performance we consider the errors obtained as the maximum of the absolute errors on the mesh points \(x_n=nh\), \(n=\overline{0,\, N}\)

for different values of step \(h\).

In the next tables we present results given by the proposed methods (17), (25) (31) in comparison with results given by Adams-Bashforth methods of orders \(2\), \(3\), \(4\) defined by

• \(A-B(2):\quad y_{n+1}-y_n=\tfrac {h}{2}(3f_n-f_{n-1}),\)

• \(A-B(3):\quad y_{n+1}-y_n=\tfrac {h}{12}(23f_n-16f_{n-1}+5f_{n-2}),\)

• \(A-B(4):\quad y_{n+1}-y_n=\tfrac {h}{24}(55f_n-59f_{n-1}+37f_{n-2}-9f_{n-3})\).

Also, we present results given by the explicit Runge-Kutta methods with order equal with number of stages defined by the following Butcher arrays

We analyze the results and we can see that for the above examples the proposed methods are at least comparable to the classical ones.

We observe that if length of the step decreases \(10\) times then the error magnitude for methods (17), (25) and (31) decreases \(10^2\), \(10^3\) and \(10^4\) times, respectively. These results represent a validation of the fact that the convergence orders for the methods (17), (25) and (31) are \(2\), \(3\) and \(4\), respectively.

Because the method is explicit we indicate to be used especially for non-stiff problems, where the requirements on the steplength imposed by stability are no restrictive than the requirements imposed by accuracy, see [ 1 ] and [ 3 ] for more details. Also, we expect that the method (12) gives better results in the case of a variable steplength.