Research Article  Open Access
H. Saberi Nik, Sohrab Effati, Sandile S. Motsa, Stanford Shateyi, "A New PiecewiseSpectral Homotopy Analysis Method for Solving Chaotic Systems of Initial Value Problems", Mathematical Problems in Engineering, vol. 2013, Article ID 583193, 13 pages, 2013. https://doi.org/10.1155/2013/583193
A New PiecewiseSpectral Homotopy Analysis Method for Solving Chaotic Systems of Initial Value Problems
Abstract
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 GenesioTesi 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.
1. 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–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 piecewise (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 timeconsuming 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 GenesioTesi system. Chaotic dynamics of the GenesioTesi 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 GenesioTesi 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 cellfree viral spread of human immunodeficiency virus (HIV) in a wellmixed compartment such as the bloodstream. This model consists of four components: uninfected healthy T cells, latently infected T cells, actively infected 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 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–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 shortterm behaviour of some systems, longterm 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 piecewiseSHAM. In Section 5, the numerical implementation of the PSHAM to the GenesioTesi system and HIV infection of 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.
2. 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 is sufficiently smooth for linearization techniques to be valid. If we can apply the SHAM by rewriting (1) as where is a linear operator which is derived from the entire linear part of (1) and is the remaining nonlinear component for .
The SHAM approach imports the conventional ideas of the standard homotopy analysis method (HAM) by defining the following zerothorder deformation equations where is an embedding parameter, are unknown functions, and and are convergence controlling parameters and functions, respectively. The nonlinear operator is defined as
Using the ideas of the standard HAM approach (see e.g., [24–28]), we differentiate the zerothorder equations (4) times with respect to and then set 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 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 socalled 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 of (10) is then fed to (6) which is iteratively solved for (for ).
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.
3.1. Properties of Shifted Legendre Polynomials
The wellknown Legendre polynomials are defined on the interval and can be determined with the aid of the following recurrence formula [30]: In order to use these polynomials on the interval we defined the socalled shifted Legendre polynomials by introducing the change of variable . Let the shifted Legendre polynomials be denoted by . Then can be generated by using the following recurrence relation: where and . The orthogonality condition is where is the Kronecker function. Any function , square integrable in , may be expressed in terms of shifted Legendre polynomials as where the coefficients are given by In practice, only the first terms of shifted Legendre polynomials are considered. Hence we can write Now, we turn to LegendreGauss interpolation. We denote by , , the nodes of the standard LegendreGauss interpolation on the interval . The corresponding Christoffel numbers are . The nodes of the shifted LegendreGauss interpolation on the interval are the zeros of , which are denoted by . Clearly . The corresponding Christoffel numbers are . Let be the set of all polynomials of degree at most . Due to the property of the standard LegendreGauss quadrature, it follows that for any
Definition 1. Let and be the inner product and the norm of space , respectively. We introduce the following discrete inner product and norm:
From (17), for any ,
where for the LegendreGauss interpolation, the Legendre GaussRadau interpolation, and the Legendre GaussLobatto integration, respectively.
Moreover, for the LegendreGauss integration and the LegendreGaussRadau integration,
For the LegendreGaussLobatto integration, usually. But for mostly used orthogonal systems in , they are equivalent, namely, for certain positive constants and ,
3.2. 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 is sufficiently smooth for linearization techniques to be valid, rewriting (22) as where is a linear operator which is derived from the entire linear part of (22) and is the remaining nonlinear component.
Let us define the nonlinear operator and the sequence as 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 , for (31) such that subject to the initial condition
Definition 2. A mapping of space into itself is said to satisfy a Lipschitz condition with Lipschitz constant if for any and , If this condition is satisfied with a Lipschitz constant such that 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 , . Then the following results hold.(i)There exists a unique such that .(ii)For any , the sequence defined by , converges to [29].
Theorem 4 (existence and uniqueness of solution). Assume that in the initial value problem (IVP) (22) satisfies condition of (36); then (33) has a unique solution.
Proof. From (33) we have
Since satisfies the Lipschitzcontinuous condition, then there exists a constant such that
for all , and all and .
Now, for the problem (23), we choose
and where is an arbitrary analytic function.
Let , then we have from (37) that
or according to the definitions of and ,
from where
It is obvious that . clearly,
Furthermore, for any ,
Integrating the previously mentioned with respect to yields that
from where
Using (46) and (43) we have
Let and . Therefore, a combination of (21), (36), and (43) results in
Using (47) and (48) results in
where is a positive constant. Then, by (47) with , instead of and multiplying the resulting inequality by , we have
Subtracting (50) from (49), after simplifying, we obtain
Therefore, if , then as . According to Theorem 3, it implies the existence and uniqueness of solution of (33).
4. PiecewiseSpectral 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 , we divide the interval into subintervals , where . We solve (6) in each subinterval . Let be the solution of (6) in the first subinterval and let be the solutions in the subintervals for . The initial conditions used in obtaining the solutions in the subinterval are obtained from the initial conditions of the subinterval . 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 . Before applying the spectral method, it is convenient to transform the region to the interval [−1, 1] on which the spectral method is defined. This can be achieved by using the linear transformation in each interval (for ). After the transformation, the interval is discretized using the LegendreGaussLobatto (LGL) nodes. These points, , , are unevenly distributed on and are defined by , and for , are the zeros of , 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 and is the vector function at the collocation points . The matrix is of size and its entries are defined [31, 32] as
Applying the the Legendre spectral collocation method in (52)–(55) gives where is an matrix that comes from transforming the linear operator using the derivative matrix , , and are vectors corresponding to the functions and , 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 which is the solution of (59), the recurrence formula can be used to obtain the solution using (9) in the interval . The solution approximating in the entire interval is given by It must be noted that when , the proposed piecewisespectral homotopy analysis method (PSHAM) becomes equivalent to the original SHAM algorithm.
5. 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 GenesioTesi system [15] and a nonlinear initial value problem that models dynamics of HIV Infection of 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 RungeKuttagenerated results.
5.1. GensioTesi System
The GenesioTesi 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 righthand side of (65) is known, the solution can easily be obtained by using methods for solving linear system of equations.
5.2. HIV Infection of CD4^{+} T Cells [18] model
Of particular interest to us is the model in [18], which is given by the following system of differential equations: with the initial conditions
In this model , , and denote the concentration of susceptible T cells, infected T cells, and free HIV virus particles in the blood, respectively. To explain the parameters, we note that is the source of T cells from precursors, , and are natural turnover rates of uninfected T cells, infected T cells, and virus particles, respectively, and is the maximum level of 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 T cells is now described by , and proliferation of infected T cells is neglected. The term describes the incidence of HIV infection of health T cells, where is the infection rate. Each of the infected T cell is assumed to produce 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 With these definitions, the PSHAM algorithm gives Because the righthand side of (70) is known, the solution can easily be obtained by using methods for solving linear system of equations.
6. 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 GenesioTesi and the HIV Infection of T cells model. Unless otherwise specified, the results presented in this paper were generated using order of the PSHAM approximation with collocation points in each 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 . 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 socalled curves. More recently, Marinca and Heriş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 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 RungeKutta (4,5) formula, the DormandPrince 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 GenesioTesi system at selected values of time . It can be seen from the table that there is good agreement between the two results.

Figure 1 gives the graphical profiles of the classes in the time interval . The initial conditions used to generate the graphs are and . The PSHAM results are plotted against the ode45 results (small dots). It can be seen from the figure that there is good agreement between the two results.
(a)
(b)
(c)
In Figure 2, we give the various phase portraits of the GenesioTesi 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.
(a)
(b)
(c)
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 , and . For this problem, in order to get accurate results the order of the PSHAM used was . 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 in the time interval . 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.
(a)
(b)
(c)
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].
(a)
(b)
(c)
7. 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 piecewisespectral 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.
References
 S. S. Motsa, P. Sibanda, and S. Shateyi, “A new spectralhomotopy analysis method for solving a nonlinear second order BVP,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 9, pp. 2293–2302, 2010. View at: Publisher Site 