A nonlinear second derivative method with a variable step-size based on continued fractions for singular initial value problems

We propose a nonlinear one-step second derivative method (NSDM) with a variable step-size implementation based on continued fractions for the numerical solution of singular initial value problems (IVPs). The singular IVPs typically originate from models in Mathematical Biology which are represented by reaction diffusion partial differential equations (PDEs). These PDEs are then discretized in space through finite difference methods to yield a system of ordinary differential equations with singular behavior associated with finite time blow-up. We are motivated by the fact that numerical methods based on polynomial approximation perform poorly in the proximity of a singular point when applied to singular IVPs. Therefore, a continuous method based on a finite continued fraction approximation is derived and used to obtain the NSDM. The NSDM is implemented in a self-starting mode and in a predictor-corrector mode, equipped with a variable step-size implementation. The superiority of the NSDM over standard methods in the literature is established numerically. Subjects: Advanced Mathematics; Applied Mathematics; Applied Physics; Computational Physics


PUBLIC INTEREST STATEMENT
Standard methods such as Runge-Kutta and linear multistep methods perform poorly in the proximity of a singular point when applied to singular initial value problems. Hence, we have proposed an accurate nonlinear one-step second derivative method with a variable step-size implementation based on continued fractions for the numerical solution of singular initial value problems. The singular initial value problems typically originate from models in Mathematical Biology which are represented by reaction diffusion partial differential equations discretized in space through finite difference or finite element methods to yield a system of ordinary differential equations with singular behavior associated with finite time blow-up.

Introduction
We consider the IVP Several physical processes in science and engineering can be modeled in the form (1) which could be non-stiff, stiff or singular in nature. Specifically, we consider singular IVPs that typically originate from models in Mathematical Biology triggered by a random walk of organisms, such as bacterial in the presence of a chemical resulting in interesting spatial patterns (see Tyson, Lubkin, & Murray, 1999;Tyson, Stern, & LeVeque, 2000). This kind of motion is characterized as chemotaxis (see Tyson et al., 2000) and are traditionally modeled through a reaction diffusion system of partial differential equations (PDEs). Moreover, physical processes such as flows in porous media, heat transfer, and biochemical kinetics are modeled through PDEs (see Abia, López-Marcos, & Martinez, 1998). These PDEs are then discretized in space through finite difference or finite element methods to yield a system of ODEs with singular behavior associated with finite time blow-up (see Abia et al., 1998). We note that the solutions of the singular ODEs may become unbounded in finite time and such a phenomenon is often called blow-up, while the finite time is called the blow-up time (see Cho, 2013;Groisman, 2006;Guo & Yang, 2015;Takayasu, Matsue, Sasaki, Tanaka, Mizuguchi, & Oishi, 2017). It is also reported in Takayasu et al. (2017) that blow-up solutions describe the combustion of solid fuels via the nonlinear heat equation.
In applications, most equations in the form (1) do not possess exact solutions, hence numerical methods are imperative for solving (1). Several numerical methods based on Runge-Kutta methods, linear multistep methods, Boundary value methods, and block methods have been extensively discussed in the literature for solving non-stiff and stiff IVPs (see Akinfenwa & Jator, 2013;Brugnano & Trigiante, 1998, Jator, 2010Cash, 1981;Hairer & Wanner, 1996;Lambert, 1991;Onumanyi, Awoyemi, Jator, & Sirisena, 1994;Onumanyi, Sirisena, & Jator, 1999). It turns out that these methods do not perform well when applied to solve singular problems of the form (1) since they are based on polynomial approximation (see Hairer & Wanner, 1996;Lambert, 1991). In particular, methods constructed using polynomial basis functions perform poorly in the proximity of a singular point when applied to singular IVPs. Hence, rational interpolants have been proposed for solving singular IVPs and were first introduced by Lambert and Shaw (1965). Other techniques based on rational functions have been proposed to cope with singular IVPs (see Fatunla, 1982;Hasan, Ahamed, Alam, & Hossain, 2013;Ikhile, 2004;Odekunle, Oye, Adee, & Ademiluyi, 2004;Van Niekerk, 1988). Some authors have proposed more efficient techniques based on a special rational interpolation function expressed in the form of a finite continued fraction (see Van Niekerk, 1988). Nonstandard methods for solving singular IVPs have also been proposed (see Ramos & Vigo-Aguiar, 2008), which a onestep method based on rational logarithmic basis functions has also been developed (see Garwood & Jator, 2016).
In terms of implementation, is generally implemented in a step-by-step fashion in which on the partition H N , an approximation y n is obtained at x n only after an approximation at x n−1 has been computed, where for some constant step-size h = x N −x 0 N and integer N > 0, In this paper, this approach is used as a predictor-corrector method with a fixed step-size and variable step-size. We have also used another approach that is self-starting and implemented without predictors (Akinfenwa & Jator, 2013;Jator, 2010). We propose a new method based on a continuous finite continued fraction approximation (see Van Niekerk, 1988), which is used to obtained the NSDM. The NSDM is implemented in a SSM and PCM equipped with a variable step-size implementation. The superiority of the NSDM over standard methods, such as the Trapezoidal Method (TM) and the fourth order Runge-Kutta Method (RKM) is established numerically. We also show that the method can be (1) used to solve evolution equations and parabolic problems with blow-up solutions (see Cho, 2013;Groisman, 2006;Guo & Yang, 2015;Takayasu et al., 2017).
The paper is organized as follows. In Section 2, we derive the NSDM and examine its implementation in Section 3. In Section 4, we perform numerical experiments and give the conclusion of the paper in Section 5.

Derivation of method
Although (1) is a system of ODEs, in our derivation, we initially assume the scalar form of (1) for simplicity, since the corresponding system can also be solved by obvious notational modifications. We are motivated by the fact that numerical methods based on polynomial approximation perform poorly in the proximity of a singular point when applied to singular IVPs. Therefore, a continuous method based on a rational function is proposed since rational are considered better approximants for singular IVPs (see Fatunla, 1982;Ikhile, 2004;Odekunle et al., 2004). The derivation of the method is based on a special kind of rational interpolation function expressed as a continued fraction (see Van Niekerk, 1988). In order to solve the IVP (1) on the sequence of points {x n }, defined by where a 0 , … , a 3 , are coefficients to be uniquely determined. We further assume that y n+j = u(x n + jh) and y(x n+j ) ≈ y n+j , f n+j = U � (x n + jh) is an approximation to y � (x n+j ), g n+j = U �� (x n + jh) is an approximation to y �� (x n+j ), j = 0, 1. In order to determine the coefficients (a 0 , a 1 , a 2 , a 3 ) in (2), we demand that the following conditions must be satisfied and lead to a system of four equations and four unknowns, which is solved to obtain the coefficients (a 0 , a 1 , a 2 , a 3 ). These coefficients are now known and substituted into in (2) to yield the continuous form u(x) which is then evaluated at x n+1 to yield the NSDM given by We note that y n+j is the numerical approximation to the analytical solution y(x n+j ), f n+j = f (x n+j , y n+j ),

The Local Truncation Error (LTE) of the method (3) is obtained by carry out a Taylor Series expan-
sion about the point x n to get hence the method is of order 3.

Derivation of the predictor:
Since the method is also implemented in the PCM, the predictor is derived similarly. Thus, the predictor is constructed by assuming that on the interval [x n , x n+1 ], the exact solution y(x), is approximated by the finite continued fraction u p ; i.e. y(x) ≈ u p (x), thus, where b 0 , … , b 3 , are coefficients that are uniquely determined. We further assume that y n = u(x n ) and y(x n ) ≈ y n , f n = U � (x n ) is an approximation to y � (x n ), g n = U �� (x n ) is an approximation to y �� (x n ), and z n+j = U ��� (x n ) is an approximation to y ��� (x n ). In order to determine the coefficients (5), we demand that the following conditions must be satisfied and lead to a system of four equations and four unknowns, which is solved to obtain the coefficients . These coefficients are now known and substituted into in (5) to yield the continuous form of the predictor u p (x) which is then evaluated at x n+1 to yield The LTE for the predictor is given by hence the method is also of order 3.
Remark 2.1 While the conditions imposed to construct (2) and (5) are based on interpolation and collocation, Hermite interpolation generally involves the use of a table of divided differences. We also note that g(x, y(

Implementation of method
PCM-fixed step size: The method is implemented as a predictor-corrector using (3) and (6) in a stepby-step fashion on the partition H N . The approximation at x n is used to obtain the approximation at x n+1 , such that is the constant step-size, N is a positive integer, and n is the grid index. The PCM-fixed step-size implementation is given in Algorithm 1.
y p n+1 = −3h 2 g 2 n − 6g n y n + 2hy n z n + 2hf n −3g n + hz n −6g n + 2hz n (7) The method is implemented in the SSM as follows: Using method (3), n = 0, solve for the values of y 1 with the aid of Newton's Method on the sub-interval [x 0 , x 1 ], as y 0 is known from the IVP (1). Next, for n = 1, the values of y 2 is obtained over the sub-interval [x 1 , x 2 ], as y 1 is known from the previous sub-interval. The process is continued for n = 2, … , N − 1 to obtain the numerical solution to (1) on the sub-intervals [x 2 , The details of the implementation is given in Algorithm 2.

Variable-step implementation of method:
The method is also implemented in a PCM in a variable-step fashion on the partition H N , in which case, we need an error estimate ( n+1 ) for y n+1 . In order to obtain the error estimate, we use (7) to estimate the error in the predictor (y p n+1 ), which is then used to estimate the error in the corrector (y n+1 ) given by (4). Thus, using and In order to proceed with the variable step-size implementation, we use n+1 as an estimate for the LTE in y p n+1 and control the step-size of the integration on the basis of this estimate (see Cash, 1981).
If a local tolerance (Tol) is imposed at each step, then, the new step-size h new and the old step-size h old are related as described in Algorithm 3.

Numerical examples
In this section, we give numerical examples to illustrate the accuracy of the method in the neighborhood of a singularity. Let y(x n ) be the exact solution and y n the approximate solution on the partition H N , we find the error of the approximate solution as | SSM and PCM equipped with a variable step-size implementation. The superiority of the NSDM over the TM and the fourth order RKM is established numerically.
Example 4.1 We consider the singular IVP that is also solved in Garwood and Jator (2016), whose exact solution is This example is solved using the NSDM in the PCM and SSM as well as solved using the TM and RKM and the absolute errors in the neighborhood of the singularity (x = ∕4 ≈ 0.7853981633974483) are compared. The results displayed in Table 1 show the superiority of the NSDM over the TM and RKM. Figure 1 also shows that the NSDM is in good agreement with the exact solution in the neighborhood of the singularity.

whose exact solution is
This example is solved using the NSDM in the PCM and SSM as well as solved using the TM and RKM and the absolute errors in the neighborhood of the singularity (x = 1) are compared. The results displayed in Table 2 show the superiority of the SSM over TM and RKM. We note that the PCM does not perform well in the proximity of the singularity. Nevertheless, it performs excellently when applied adaptively as the results in Tables 3 and 4 indicate.

Example 4.3
We consider the following nonlinear ODE whose blow-up time T = √ 6 that was also solved in Cho (2013). whose exact solution is This example is solved using the NSDM in the PCM and SSM and the results given in Figure 2 show that the NSDM is in good agreement with the exact solution in the neighborhood of the singularity.

Parabolic equations with blow-up solutions
In this subsection, we SSM to solve parabolic equations with blow-up solutions, since the SSM has the advantage of being self-starting and used without predictors (Figure 3).

Example 4.4
We consider the PDE given (Guo & Yang, 2015). and the exact solution is given by u(x, t) = e −t sin( x).
We obtain the following first order system after discretizing the spatial variable x using the second order finite difference scheme  which is written in the form u � = f (t, u), subject to the initial conditions u(t 0 ) = u 0 where f (t, u) = Au + q, q is the matrix arising from the central difference approximations and q is a vector of constants. We then define u m The problem is now a system of ODEs which is solved by the SSM. In Figure 4, we show that the numerical method (SSM) is in good agreement with the exact solution at T = 0.1. We note that this consistent is also observed at x = 1 where a numerical method could perform (see Guo & Yang, 2015). In Figure 5, we that the SSM performs well when applied to parabolic PDEs as the numerical  Example 4.5 We consider the nonlinear PDE given (Guo & Yang, 2015).
we obtain the following first order system after discretizing the spatial variable x using second order finite difference schemes which is written in the form u � = f (t, u), subject to the initial conditions u(t 0 ) = u 0 where f (t, u) = Au + q, A is the matrix arising from the central difference approximations and q is a vector of constants. We then define u m (t) ≈ u(x m , t), u = [u 1 (t), u 2 (t), … , u M−1 (t)] T , T is the transpose.
The problem is now a system of ordinary differential equations which is solved by the SSM. In Figure 6, we show that the SSM provides results for the reference blow-up time of T = 8.24371 × 10 −2 for N = 80 and N = 160 that are consistent with the results given in Guo and Yang (2015). It observed that the blow-up is located at only one point x = 0.5 (see Guo & Yang, 2015).   Note: The variable step version of the PCM performs better that the fixed step version as well as the SSM as expected. Note: The variable step version of the PCM performs better that the fixed step version as well as the SSM as expected. Note: The variable step version of the PCM performs better that the fixed step version as well as the SSM as expected.