An efficient variable step-size rational Falkner-type method for solving the special second-order IVP

Abstract

In this paper, firstly a rational one-parameter family of Falkner-type explicit methods is presented for directly solving numerically special second order initial value problems in ordinary differential equations. The proposed family of methods has second algebraic order of convergence. Imposing that the principal term of the local truncation error of the proposed family vanishes, we get an expression for the free parameter at the grid point (x_n, y_n). By substituting this value of the free parameter in the family, a new rational third order method is obtained. Further, by combining the third order method with any member of the second order family, their variable step-size formulation as an embedded pair is considered. Some numerical experiments are given to illustrate the performance and efficiency of the proposed methods.

Introduction

The present paper is concerned with the following so-called special second-order differential equation

$$y^{\prime \prime }\left(x\right)=f(x,y\left(x\right)),y\left(a\right)=y_0,y^{\prime }\left(a\right)=y_0^{\prime },$$

where $x\in [a,b],y:[a,b]\to \mathbb{R}$  and $f:[a,b]\times \mathbb{R}\to \mathbb{R}$ are sufficiently differentiable functions. It is possible that (1) can be integrated by reformulating it as a system of two first order ODEs and applying one of the methods available for those systems (see for example [7], [13], [14], [17]). Nevertheless, it seems more natural to provide numerical methods to integrate (1) directly without transforming it in a first order system. These approaches result to be more efficient. For instance, it is well-known that a Runge–Kutta–Nyström method for solving (1) has a real improvement as compared to standard Runge–Kutta methods. In case of a linear k-step method for first order ODEs, it becomes a 2k-step method for (1), thus increasing the computational work.

For the general second order initial value problem of the form

$$y^{\prime \prime }\left(x\right)=f(x,y\left(x\right),y^{\prime }\left(x\right)),y\left(a\right)=y_0,y^{\prime }\left(a\right)=y_0^{\prime }$$

one of the methods for solving it directly is the explicit method due to Falkner [1] which can be written in the form

$$y_{n+1}=y_n+h{y}_n^{\prime }+h^2\sum _{j=0}^{k-1}{\beta }_j{\nabla }^j{f}_n,$$

$$y_{n+1}^{\prime }=y_n^{\prime }+h\sum _{j=0}^{k-1}{\gamma }_j{\nabla }^j{f}_n,$$

where h is the step-size, y_n and $y_n^{\prime }$ are numerical approximations to the theoretical solution and its derivative respectively, at the grid point $x_n=a+nh;n=0,1,2,3,\dots N,h=\frac{(b-a)}{N},f_n=f(x_n,y_n,y_n^{\prime })$ and ∇^jf_n is the standard notation for the backward differences.

There exist similar implicit Falkner formulas [2] given by

$$y_{n+1}=y_n+h{y}_n^{\prime }+h^2\sum _{j=0}^k{\beta }_j^{*}{\nabla }^j{f}_{n+1},$$

$$y_{n+1}^{\prime }=y_n^{\prime }+h\sum _{j=0}^k{\gamma }_j^{*}{\nabla }^j{f}_{n+1}.$$

Note here that the formulas given in (3) and (5) are the Adams–Bashforth and Adams–Moulton methods respectively for solving the problem

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

which are used to obtain the approximate values of the first derivative.

The usual and unusual implementation of these methods have been considered in the literature. For instance, in molecular dynamics, when the acceleration at time only depends on position and not on velocity, the direct integration methods are usually implemented in a semi-implicit formulation. This is the case for the well-known Velocity Verlet algorithm [3]. This method uses the one-step explicit method in (2) to compute the positions:

$$y_{n+1}=y_n+h{y}_n^{\prime }+\frac{h^2}{2}{y}_n^{\prime \prime }$$

and the one-step implicit method in (5) to update the velocities:

$$y_{n+1}^{\prime }=y_n^{\prime }+\frac{h}{2}(y_n^{\prime \prime }+y_{n+1}^{\prime \prime }).$$

Beeman [4] proposed the modification of the Verlet family of methods for the calculation of velocities. The method used to compute the positions at time $x_n+h$ is the following two-step explicit method given in (2)

$$y_{n+1}=y_n+h{y}_n^{\prime }+\frac{h^2}{6}(4y_n^{\prime \prime }-y_{n-1}^{\prime \prime })$$

while the formula to update the velocities is the following two-step method (which we note it is not the two-step formula in (5))

$$y_{n+1}^{\prime }=y_n^{\prime }+\frac{h}{6}(2y_{n+1}^{\prime \prime }+5y_n^{\prime \prime }-y_{n-1}^{\prime \prime }).$$

Two important aspects must be considered in this paper in order to obtain the methods that will be proposed: the way to implement the formulas, and the variable stepsize character of the numerical scheme. We briefly outline some comments on them.

The implementation of the explicit Falkner method for the special second order IVP is considered in usual and unusual way in the literature, which are described below (see [5]).

The usual implementation of the explicit Falkner method on each step of integration for solving (1) consists in the following steps

• Evaluate $y_{n+1}$ using the formula given in (2).

• Evaluate $y_{n+1}^{\prime }$ using the formula given in (3).

• Evaluate $f_{n+1}=f(x_{n+1},y_{n+1}).$

The unusual implementation of the explicit Falkner method (also called reformed Falkner method in [15], [16], at each step of integration for solving (1) consists in the following steps

• Evaluate $y_{n+1}$ using the formula given in (2).

• Evaluate $f_{n+1}=f(x_{n+1},y_{n+1}).$

• Evaluate $y_{n+1}^{\prime }$ using the formula given in (5),

which can be accomplished due to the absence of the derivative on the function f. Therefore, having obtained the value $y_{n+1}$ it is a straightforward task to calculate $f_{n+1}$ to be used in (5). In this way, the implicit method given in (5) is no longer implicit, resulting in an explicit formulation of the method. For instance, the two-step procedure provided by Falkner formulas is given by

$$\begin{array}{ccc}\hfill y_{n+1}& =& y_n+h{y}_n^{\prime }+\frac{h^2}{6}(4y_n^{\prime \prime }-y_{n-1}^{\prime \prime })\hfill \\ \hfill y_{n+1}^{\prime }& =& y_n^{\prime }+\frac{h}{12}(5y_{n+1}^{\prime \prime }+8y_n^{\prime \prime }-y_{n-1}^{\prime \prime }).\hfill \end{array}$$

In virtually all modern codes for ODEs, the step-size is selected automatically to achieve reliability and efficiency [8]. Practically speaking, any discretization integrator with constant step-size performs poorly if the solution varies rapidly in some parts of the integration interval and slowly in other parts of the integration interval. As some authors have remarked, to be efficient, an integrator based on a particular formula must be suitable for a variable step-size formulation [10], [11]. There are mainly two ways to do that, using a formula where the coefficients depend on the ratios of the step-sizes (as in the variable-coefficient linear multistep codes) or using a second method (as in the embedding Runge–Kutta pairs). In all the cases the goal is to adjust the step-sizes so as to keep the estimated local errors smaller than a given tolerance and, at the same time, to solve the problem as efficiently as possible [12]. So, a reliable estimate of the local error will be needed.

This paper aims at developing some variable step-size rational Falkner-type pairs for directly solving numerically the problem in (1). In the present paper some rational 3(2) pairs of Falkner type are proposed and numerical results are given in comparison with some of the existing classical methods. The paper is organized as follows: Section 2 is concerned with the derivation of the methods. In Section 3, the error analysis of the methods are carried out. In Section 4 some rational Falkner methods are listed. In Section 5 the linear stability analysis of the methods is given. In Section 6 a variable step-size formulation of the proposed methods is considered. Some numerical experiments are presented in Section 7 to validate the performance of the proposed method, and some conclusions put at the end to the paper.