A New Piecewise-Spectral Homotopy Analysis Method for Solving Chaotic Systems of Initial Value Problems

An accurate algorithm for solving initial value problems (IVPs) which are highly oscillatory is proposed. The proposed method is based on a novel technique of extending the standard spectral homotopy analysis method (SHAM) and adapting it to a sequence of multiple intervals. In this new application the method is referred to as the piecewise spectral homotopy analysis method (PSHAM). The applicability of the proposed method is examined on the differential equation system modeling HIV infection of CD4 T cells and the Genesio-Tesi system which is known to be chaotic and highly oscillatory. Also, for the first time, we present here a convergence proof for SHAM. We treat in detail Legendre collocation and Chebyshev collocation. The method is compared to MATLAB’s ode45 inbuilt solver as a measure of accuracy and efficiency.


Introduction
This paper introduces a new method for solving highly oscillatory and chaotic initial value problems (IVPs).The new method extends, for the first time, the application of the spectral homotopy analysis method [1,2] to IVPs.The spectral homotopy analysis method was developed by Motsa et al. [1,2] for solving nonlinear boundary value problems (BVPs) over finite intervals.It has been successfully been applied to other BVPs arising mainly in fluid mechanicsrelated problems [3][4][5][6].In the previously cited applications the SHAM method was applied to problems which possess smooth solutions over small regions.For rapidly oscillating chaotic systems over very large regions, the SHAM may not give accurate results.The current work seeks to develop a new method that will be valid for rapidly changing solutions over all regions, small, medium-, and large sized.A simple way of ensuring the validity of the approximations for large intervals and for all functions is to determine the solution in a sequence of equal intervals, which are subject to continuity conditions at the end points of each interval.
Recently, in an effort to increase the radius of convergence of some analytical methods of approximations, multistage or piecewise approximations have been developed for solving IVPs over general intervals.This multistage approach seeks to implement the standard approximation method on sequences of subintervals whose union makes up the domain of the underlying problem.The effect of this piece-wise (multistep) approach is to accelerate the convergence of the approximate solution over a large region and to improve the accuracy of the parent approximate method of solution, particularly over rapidly oscillating regions.Examples include the, multistage homotopy analysis method [7], piecewise homotopy perturbation methods [8], multistage Adomian decomposition method [9,10], multistage differential transformation method, [11,12], and multistage variational iteration method [13].Because, these methods attempt to obtain analytical solutions at each interval they involve time-consuming and tedious computational operations and if too many small intervals are considered, as may be the case when dealing with highly oscillatory systems, the analytical integration process will be too much to handle even with the use of symbolic scientific software.For this reason, focus is now shifted towards multistage methods which use numerical integration techniques such as the piecewise iteration method [14] which uses spectral collocation methods.
This work examines the applicability of a method that blends the piecewise implementation technique with the spectral homotopy analysis method.For this reason, the proposed method is called the piecewise homotopy analysis method (PSHAM).The validity of the method is tested against two initial value problems of practical importance, namely, the mathematical models of HIV infection of CD4 + T cells and the Genesio-Tesi system.Chaotic dynamics of the Genesio-Tesi system were introduced in [15].Since then, the system has been used as one of the benchmarks for testing the validity of new and existing methods of solving IVPs.Approximate analytical or numerical methods that have been used to find solutions of the Genesio-Tesi system include the differential transform method [16], the piecewisespectral parametric iteration method [14], variational iteration method (VIM), and the multistage VIM [17].
Mathematical models play a significant role in studying epidemic and viral dynamics.Viral models are valuable to help us improve the understanding of both diseases and various drug therapy strategies against them.In 1993, Perelson et al. [18] proposed an ODE model of cell-free viral spread of human immunodeficiency virus (HIV) in a well-mixed compartment such as the bloodstream.This model consists of four components: uninfected healthy CD4 + T cells, latently infected CD4 + T cells, actively infected CD4 + T cells, and free virus.This model has been important in the field of mathematical modelling of HIV infection, and many other models have been proposed, which take the model of Perelson et al. [18] as their inspiration.For numerically solving a model for HIV infection of CD4 + T cells, Ongun [19] has applied the Laplace Adomian decomposition method, Merdan has used the homotopy perturbation method [20], and Merdan et al. [21] have applied the Padé approximation and the modified variational iteration method [22].More recently Yuzbasi [23] used the Bessel collocation method for solving the same problem.It is worth noting that in all these studies results were given for solutions over small time intervals.Thus, it may be concluded that these methods only offer a good approximation to the true solution in a very small regions.Liao [24] introduced the homotopy analysis method (HAM) which has convergence controlling parameters to improve accuracy of series solutions and increase the region of convergence.The HAM was successfully used to solve several nonlinear IVPs [24][25][26][27][28].However, even the standard HAM approach may not be suitable to resolve solutions over very large intervals and for rapidly oscillating systems which show chaotic behaviour.For practical purposes, results obtained using the previously mentioned analytical methods of approximation may not be very useful since, in addition to the short-term behaviour of some systems, long-term behaviour of the solutions may be required and insights into chaotic systems may be sought.The multistage approach, such as the one discussed in this paper, presents an important strategy of addressing the limitations of some of the previously cited methods of solution.
The rest of this paper is organized as follows: In Section 2 the basic description of the SHAM is presented.Section 3 presents existence and uniqueness of solution of SHAM that give a guarantee of convergence of SHAM.Section 4 deals with the description and formulation of the piecewise-SHAM.In Section 5, the numerical implementation of the PSHAM to the Genesio-Tesi system and HIV infection of CD4 + T cells models is considered.In Section 6 we present the main results and discussion of the numerical simulations, and finally, the conclusion is given in Section 7.

Basic Idea behind the Spectral Homotopy Analysis Method (SHAM)
For convenience of the interested reader, we will first present a brief description of the basic idea behind the standard SHAM [1,2].This will be followed by a description of the piecewise version of the SHAM algorithm which is suitable for solving initial value problems.To this end, we consider the initial value problem (IVP) of dimension  given as where the dot denotes differentiation with respect to .We make the usual assumption that f is sufficiently smooth for linearization techniques to be valid.If y = ( 1 ,  2 , . . .,   ) we can apply the SHAM by rewriting (1) as where L  is a linear operator which is derived from the entire linear part of (1) and F  is the remaining nonlinear component for  = 1, 2, . . ., .
The SHAM approach imports the conventional ideas of the standard homotopy analysis method (HAM) by defining the following zeroth-order deformation equations where  ∈ [0,1] is an embedding parameter,   (; ) are unknown functions, and ℏ  and   () are convergence controlling parameters and functions, respectively.The nonlinear operator N  is defined as Using the ideas of the standard HAM approach (see e.g., [24][25][26][27][28]), we differentiate the zeroth-order equations (4)  times with respect to  and then set  = 0 and finally divide the resulting equations by ! to obtain the following equations, which are referred to as the th order (or higher order) deformation equations, where After obtaining solutions for (6), the approximate solution for each   () is determined as the series solution A HAM solution is said to be of order  if the previous series is truncated at  = , that is, if It should be emphasized that the HAM and SHAM approach will be equivalent if the same linear operator L is used.In the HAM implementation the linear operator is usually chosen in such a way that the higher order deformation equations ( 6) can be solved analytically and the resulting solution conforms to the so-called rule of coefficient ergodicity and a predefined rule of solution expression [24].If these requirements are adhered to then equations of the form (6) constitute a system of decoupled differential equations which can easily be integrated, especially with the use of symbolic scientific software such as MATLAB, MAPLE, and Mathematica.If the HAM solution does not converge to the true solution of the underlying problem the convergence controlling parameters ℏ  and   () may be adjusted to improve convergence.In some cases even this intervention will not be enough to guarantee convergence.The SHAM approach seeks to remove all of the restriction of the HAM by making it more flexible.In the SHAM application, the governing differential equations are just separated into the linear and nonlinear components as in (3).The linear operator is just chosen using the linear part of the equation.The penalty for relaxing the restrictions is that the resulting higher order deformation equations will not be solvable analytically.That is, where the "S" in the SHAM comes in.Spectral collocation methods are used to solve the higher order deformation equations.A suitable initial guess to start off the SHAM algorithm is obtained by solving the linear part of (3) subject to the given initial conditions; that is, we solve If ( 10) cannot be solved exactly, the spectral collocation method is used as a means of solution.The solution  ,0 () of ( 10) is then fed to (6) which is iteratively solved for  , () (for  = 1, 2, 3, . . ., ).

Convergence Theorem of the Spectral Homotopy Analysis Method
In this section, we give some definitions and theorems that explain the convergence properties of the SHAM from a functional analysis framework [29].We remark that the Newton's iteration approach was used in [22] to address the nonlinearities, whereas in this work we use the homotopy analysis method approach to decompose the nonlinear terms.
In addition, our proposed method is applied in a series of subintervals as opposed to the entire interval.We begin by giving a brief description of the properties of the Legendre polynomials, followed by the existence and uniqueness properties of the solution of the SHAM and it's convergence properties.

Existence and Uniqueness of the SHAM Solution.
We consider the initial value problem (IVP) of dimension  given as where the dot denotes differentiation with respect to .We make the usual assumption that f is sufficiently smooth for linearization techniques to be valid, rewriting (22) as where L is a linear operator which is derived from the entire linear part of (22) and N is the remaining nonlinear component.
Let us define the nonlinear operator N and the sequence . . .
Using the ideas of the HAM approach [24], we obtain the following equation, which is referred to as the th order (or higher order) deformation equation: subject to the initial condition where Therefore, Upon summing these equations, we have and from (25) we obtain subject to the initial condition Consequently, the collocation method is based on a solution   () ∈ P +1 (0, ), for (31) such that subject to the initial condition If this condition is satisfied with a Lipschitz constant  such that 0 ≤  < 1 then  is called a contraction mapping.
Theorem 3 (Banach's fixed point theorem).Assume that K is a nonempty closed set in a Banach space V, and further, that  :  →  is a contractive mapping with contractivity constant , 0 ≤  < 1.Then the following results hold.
(i) There exists a unique  ∈  such that  = ().
Proof.From (33) we have Since (, ) satisfies the Lipschitz-continuous condition, then there exists a constant  ≥ 0 such that

Piecewise-Spectral Homotopy Analysis Method
It is worth noting that the SHAM method described previously is ideally suited for boundary value problems whose solutions do not rapidly change in behaviour or oscillate over small regions of the domain of the governing problem.The SHAM solution can thus be considered to be local in nature and may not be suitable for initial value problems at very large values of the independent variable .A simple way of ensuring the validity of the approximations for large  is to determine the solution in a sequence of equal intervals, which are subject to continuity conditions at the end points of each interval.To extend this solution over the interval Λ = [ 0 ,   ], we divide the interval Λ into subintervals Λ  = [ −1 ,   ],  = 1, 2, 3, . . .,  where  0 ≤  1 ≤ ⋅ ⋅ ⋅ ≤   .We solve (6) in each subinterval Λ  .Let  1  () be the solution of ( 6) in the first subinterval [ 0 ,  1 ] and let    () be the solutions in the subintervals Λ  for 2 ≤  ≤ .The initial conditions used in obtaining the solutions in the subinterval Λ  (2 ≤  ≤ ) are obtained from the initial conditions of the subinterval Λ −1 .Thus, we solve subject to the initial condition The initial approximations for solving (8) are obtained as solutions of the system subject to the initial condition The Legendre spectral collocation method is then applied to solve (52)-(54) on each interval [ −1 ,   ].Before applying the spectral method, it is convenient to transform the region [ −1 ,   ] to the interval [−1, 1] on which the spectral method is defined.This can be achieved by using the linear transformation in each interval [ −1 ,   ] (for  = 1, . . ., ).After the transformation, the interval [ −1 ,   ] is discretized using the Legendre-Gauss-Lobatto (LGL) nodes.These points,    ,  = 0, 1, . . ., , are unevenly distributed on [−1, 1] and are defined by   0 = −1,    = 1 and for 1 ≤  ≤  − 1,    are the zeros of L  , the derivative of the Legendre polynomial of degree ,   .
The Legendre spectral differentiation matrix  is used to approximate the derivatives of the unknown variables   , () at the collocation points as the matrix vector product where D = 2/(  −     .The matrix  is of size ( + 1) × ( + 1) and its entries are defined [31,32] as Applying the the Legendre spectral collocation method in (52)-( 55) gives where A is an ( + 1) × ( + 1) matrix that comes from transforming the linear operator L using the derivative matrix , Φ  , and R  ,−1 are (+1)×1 vectors corresponding to the functions   () and   ,−1 , respectively, when evaluated at the collocation points.We remark that, throughout this paper, for simplicity, we have set the auxiliary function   () appearing in (52) to be 1.Thus, starting from the initial approximation Y  ,0 which is the solution of (59), the recurrence formula can be used to obtain the solution    () using ( 9) in the interval [  ,  −1 ].The solution approximating   () in the entire interval [ 0 ,   ] is given by It must be noted that when  = 1, the proposed piecewisespectral homotopy analysis method (PSHAM) becomes equivalent to the original SHAM algorithm.

Numerical Experiments
To demonstrate the applicability of the proposed piecewisespectral homotopy analysis method (PSHAM) algorithm as an appropriate tool for solving nonlinear IVPs, we apply the proposed algorithm to the Genesio-Tesi system [15] and a nonlinear initial value problem that models dynamics of HIV Infection of CD4 + T Cells [18].Both these examples exhibit solutions which are chaotic and highly oscillatory over small regions.Thus the standard SHAM approach would not be suitable as a method of solving these problems.The results obtained by the piecewise SLM are compared with the results of the standard SLM approach and Runge-Kutta-generated results.

Gensio-Tesi
System.The Genesio-Tesi system [15] considered in this paper is given as follows: In this example, the parameters used in the SHAM and PSHAM algorithms are = [ [ With these definitions, the PSHAM algorithm gives Because the right-hand side of (65) is known, the solution can easily be obtained by using methods for solving linear system of equations.[18] model.Of particular interest to us is the model in [18], which is given by the following system of differential equations:

HIV Infection of 𝐶𝐷4 + T Cells
with the initial conditions In this model  1 ,  2 , and  3 denote the concentration of susceptible CD4 + T cells, infected CD4 + T cells, and free HIV virus particles in the blood, respectively.To explain   the parameters, we note that  is the source of CD4 + T cells from precursors, , , and  are natural turnover rates of uninfected T cells, infected T cells, and virus particles, respectively, and  max is the maximum level of CD4 + T cells concentration in the body.Because of the viral burden on the HIV infected T cells, we assume that  ≤ .The logistic growth of the healthy CD4 + T cells is now described by  *  1 (1 − (( 1 +  2 )/( max ))), and proliferation of infected CD4 + T cells is neglected.The term  2  3 describes the incidence of HIV infection of health CD4 + T cells, where   > 0 is the infection rate.Each of the infected CD4 + T cell is assumed to produce  1 virus particles during its life time, including any of its daughter cells [21,33].Throughout this paper, we set In this example, the parameters used in the SHAM and PSHAM algorithms are

A = [ [
With these definitions, the PSHAM algorithm gives Because the right-hand side of (70) is known, the solution can easily be obtained by using methods for solving linear system of equations.

Results and Discussion
In this section we present the results of the implementation of the piecewise spectral homotopy analysis method (PSHAM) described in the previous section on the Genesio-Tesi and the HIV Infection of CD4 + T cells model.Unless otherwise specified, the results presented in this paper were generated using  = 10 order of the PSHAM approximation with  = 20 collocation points in each [ −1 ,   ] interval.The width of each interval was taken to be 0.1.We remark that the values of  and  were chosen through numerical experimentation by testing a few values that give converging results for the governing parameters of the problems under consideration.The convergence controlling auxiliary parameter used in generating all the results was chosen to be ℏ  = −1.We remark that the accuracy of the homotopy analysis derivative methods, such as the one discussed in this paper, can be improved by seeking an optimal value of ℏ  using different techniques.For instance, in the traditional homotopy analysis method, Liao [24] proposed that the optimal value of ℏ was the one located in the flat region of the so-called ℏ-curves.More recently, Marinca and Heris ¸anu [34] suggested the use of the minimum of the square of the residuals to obtain the optimal values of ℏ.A detailed review of the various methods for obtaining the optimal ℏ can be found in [35].Since the PSHAM approach requires the implementation of the SHAM in numerous intervals, it may not be practical to seek an optimal value of ℏ in every interval.Thus, for illustration purposes ℏ  = −1 was used in this paper, and it was found that this value gives accurate results for all the parameters considered in this paper.In order to assess the accuracy and performance of the proposed PSHAM approach, the present results are compared with those obtained with the MATLAB inbuilt solver ode45.The ode45 solver integrates a system of ordinary differential equations using explicit 4ℎ&5ℎ Runge-Kutta (4,5) formula, the Dormand-Prince pair [36].
Table 1 gives a comparison between the present PSHAM approximate results and the numerically generated ode45 for all the governing dependent variables of the Genesio-Tesi system at selected values of time .It can be seen from the table that there is good agreement between the two results.
In Figure 2, we give the various phase portraits of the Genesio-Tesi problem using the PSHAM method.The phase portraits are qualitatively similar to those obtained in [12] wherein the multistage differential transform method was used to solve the same problem.Numerical simulations using ode45 yield exactly the same results.
In Table 2 we give a comparison between the present PSHAM approximate results and the ode45 for all the governing dependent variables of the HIV infection model at selected values of time .This was done with the standard parameter values given in (68) and initial values  1 (0) = 0.1,  2 (0) = 0, and  3 (0) = 0.1.For this problem, in order to get accurate results the order of the PSHAM used was  = 25.We observe, again, that there is good agreement between the two results.This confirms the validity of the proposed PSHAM approach as a potential method for solving initial value problems including those with chaotic behaviour.
Figure 3 gives a comparison between the PSHAM and the ode45 results (small dots) for the variation of the classes  1 (),  2 (),  3 () in the time interval [0, 100].Again, excellent agreement between the two results is observed.The graphic illustrations depicted in Figure 3 indicate that our results are qualitatively similar to other results from the literature (see e.g., [12,33]).In particular, we observe that for the selected parameters the numerical simulations show the existence of periodic solutions.Because of this periodic nature of the solutions, it would be difficult to resolve the solutions for very large time intervals using the standard SHAM approach.
Figure 4 is an illustration of the phase portraits of the HIV infection model problem generated using the PSHAM method.We observe that the graphical results in Figure 4 are qualitatively similar with the multistage differential transform method results of [12].

Conclusion
In this paper we presented a new application of the spectral homotopy analysis method in solving a class of nonlinear differential equations whose solutions show chaotic behaviour.The proposed method (referred to as the piecewise-spectral homotopy analysis method (PSHAM)) of solution extends, for the first time, the application of the SHAM to initial value problems.The PSHAM accelerates the convergence of the SHAM solution over a large region and improves the accuracy of the method.The approximate PSHAM numerical results were compared with results generated using the MATLAB ode45 solver and excellent agreement was obtained.This confirms the validity of the proposed PSHAM approach as a suitable method for solving a wide variety of initial value problems in practical applications including those with chaotic behaviour.

Figure 1 :
Figure 1: Comparison between the PSHAM (solid line) and ode45 results for the Genesio-Tesi system.

Figure 3 :
Figure 3: Comparison between the PSHAM (solid line) and ode45 results for the HIV infection model.

Figure 4 :
Figure 4: Phase portraits of the HIV model using the PSHAM approach.

Table 1 :
Numerical comparison between PSHAM and ode45 for the Genesio-Tesi system.

Table 2 :
Numerical comparison between PSHAM and ode45 for the HIV infection of CD4 + T cells model.