A New Piecewise-Spectral 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 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.

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 mechanics-related 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 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 piecewise-spectral 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–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.

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 𝐲̇(𝑡)=𝐟(𝑡,𝐲(𝑡)),𝐲𝑡₀=𝐲⁰,𝐲∶ℝ⟶ℝ^𝑛,𝐟∶ℝ×ℝ^𝑛⟶ℝ^𝑛, (1)(2) 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 ℒ_𝑟𝑦_𝑟(𝑡)+ℱ_𝑟𝑦_𝑟(𝑡)=𝜙_𝑟(𝑡), (3) where ℒ_𝑟 is a linear operator which is derived from the entire linear part of (1) and ℱ_𝑟 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 1−𝑞ℒ_𝑟𝑌_𝑟𝑡;𝑞−𝑦_𝑟,0(𝑡)=𝑞ℏ_𝑟𝐻_𝑟(𝑡)𝒩_𝑟𝑌_𝑟𝑡;𝑞−𝜙_𝑟(𝑡), (4) where 𝑞∈[0,1] is an embedding parameter, 𝑌_𝑟(𝑡;𝑞) are unknown functions, and ℏ_𝑟 and 𝐻_𝑟(𝑡) are convergence controlling parameters and functions, respectively. The nonlinear operator 𝒩_𝑟 is defined as 𝒩_𝑟𝑌_𝑟𝑡;𝑞=ℒ_𝑟𝑌_𝑟𝑡;𝑞+ℱ_𝑟𝑌_𝑟𝑡;𝑞. (5)

Using the ideas of the standard HAM approach (see e.g., [24–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, ℒ_𝑟𝑦_𝑟,𝑚(𝑡)−𝜒_𝑚𝑦_{𝑟,𝑚−1}(𝑡)=ℏ_𝑟𝐻_𝑟(𝑡)𝑅_{𝑟,𝑚−1}, (6) where 𝑅_{𝑟,𝑚−1}=1()!𝜕𝑚−1^𝑚−1𝑌_𝑟𝑡;𝑞−𝜙_𝑟(𝑡)𝜕𝑞^𝑚−1_𝑞=0,𝜒_𝑚=0,𝑚⩽1,1,𝑚>1. (7)

After obtaining solutions for (6), the approximate solution for each 𝑦_𝑟(𝑡) is determined as the series solution 𝑦_𝑟(𝑡)=𝑦_𝑟,0(𝑡)+𝑦_𝑟,1(𝑡)+𝑦_𝑟,2(𝑡)+⋯. (8)

A HAM solution is said to be of order 𝑀 if the previous series is truncated at 𝑚=𝑀, that is, if 𝑦_𝑟(𝑡)=^𝑀_𝑚=0𝑦_𝑟,𝑚(𝑡). (9)

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 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 ℒ_𝑟𝑦_𝑟,0(𝑡)=𝜙_𝑟(𝑡),𝑦_𝑟,0(0)=𝑦⁰_𝑟. (10)

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,…,𝑀).

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 well-known Legendre polynomials are defined on the interval (−1,1) and can be determined with the aid of the following recurrence formula [30]: 𝒯₀(𝑥)=1,𝒯₁(𝑥)𝒯=𝑥,_𝑗+1(𝑥)=2𝑗+1𝑗+1𝑥𝒯_𝑗(𝑥)−𝑗𝒯𝑗+1^𝑗−1(𝑥),𝑗≥1. (11) In order to use these polynomials on the interval 𝑥∈(0,𝑇) we defined the so-called shifted Legendre polynomials by introducing the change of variable 𝑥=(2𝑡/𝑇)−1. Let the shifted Legendre polynomials 𝒯_𝑗((2𝑡/𝑇)−1) be denoted by 𝒯_𝑇,𝑗(𝑡). Then 𝒯_𝑇,𝑗(𝑡) can be generated by using the following recurrence relation: 𝒯_𝑇,𝑗+1(𝑡)=2𝑗+1𝑗+12𝑡𝑇−1𝒯_𝑇,𝑗(𝑡)−𝑗𝒯𝑗+1^{𝑇,𝑗−1}(𝑡),𝑗=1,2,…, (12) where 𝒯_𝑇,0(𝑡)=1 and 𝒯_𝑇,1(𝑡)=(2𝑡/𝑇)−1. The orthogonality condition is ^𝑇₀𝒯_𝑇,𝑖(𝑡)𝒯_𝑇,𝑘(𝑡)𝑇𝑑𝑡=𝛿2𝑖+1^𝑖,𝑘, (13) where 𝛿_𝑖,𝑘 is the Kronecker function. Any function 𝑦(𝑡), square integrable in (0,𝑇), may be expressed in terms of shifted Legendre polynomials as 𝑦(𝑡)=^∞_𝑗=0𝑎_𝑗𝒯_𝑇,𝑗(𝑡), (14) where the coefficients 𝑎_𝑗 are given by 𝑎_𝑗=2𝑙+1𝑇^𝑇₀𝑦(𝑡)𝒯_𝑇,𝑗(𝑡)𝑑𝑡,𝑗=0,1,2,…. (15) In practice, only the first (𝑁+1) terms of shifted Legendre polynomials are considered. Hence we can write 𝑦(𝑡)=^𝑁_𝑗=0𝑎_𝑗𝒯_𝑇,𝑗(𝑡). (16) Now, we turn to Legendre-Gauss interpolation. We denote by 𝑡^𝑁_𝑗,  0≤𝑗≤𝑁, the nodes of the standard Legendre-Gauss interpolation on the interval (−1,1). The corresponding Christoffel numbers are 𝜔^𝑁_𝑗,0≤𝑗≤𝑁. The nodes of the shifted Legendre-Gauss interpolation on the interval (0,𝑇) are the zeros of 𝒯_𝑇,𝑁+1(𝑡), which are denoted by 𝑡^𝑁_𝑇,𝑗,0≤𝑗≤𝑁. Clearly  𝑡^𝑁_𝑇,𝑗=(𝑇/2)(𝑡^𝑁_𝑗+1). The corresponding Christoffel numbers are 𝜔^𝑁_𝑇,𝑗,=(𝑇/2)𝜔^𝑁_𝑗. Let 𝒫_𝑁(0,𝑇) be the set of all polynomials of degree at most 𝑁. Due to the property of the standard Legendre-Gauss quadrature, it follows that for any Φ∈𝒫_2𝑁+1(0,𝑇)^𝑇₀Φ(𝑡)𝑇𝑑𝑡=2¹₋₁Φ𝑇2()=𝑇𝑡+1𝑑𝑡2^𝑁_𝑗=0𝜔^𝑁_𝑗Φ𝑇2𝑡^𝑁_𝑗=+1_𝑁𝑗=0𝜔^𝑁_𝑇,𝑗Φ𝑡^𝑁_𝑇,𝑗. (17)

Definition 1. Let (𝑢,𝑣)_𝑇 and ||𝑣||_𝑇 be the inner product and the norm of space 𝐿²(0,𝑇), respectively. We introduce the following discrete inner product and norm: (𝑢,𝑣)_𝑇,𝑁=^𝑁_𝑗=0𝑢𝑡^𝑁_𝑇,𝑗𝑣𝑡^𝑁_𝑇,𝑗𝜔^𝑁_𝑇,𝑗,‖𝑣‖_𝑇,𝑁=(𝑣,𝑣)^1/2_𝑇,𝑁. (18) From (17), for any Φ𝜓∈𝒫_2𝑁+𝜆(0,𝑇), Φ,𝜓_𝑇=Φ,𝜓_𝑇,𝑁, (19) where 𝜆=−1,0,1 for the Legendre-Gauss interpolation, the Legendre Gauss-Radau interpolation, and the Legendre Gauss-Lobatto integration, respectively.
Moreover, for the Legendre-Gauss integration and the Legendre-Gauss-Radau integration, 𝜙_𝑇=𝜙_𝑇,𝑁,𝜙∈𝒫_𝑁(0,𝑇). (20) For the Legendre-Gauss-Lobatto integration, ||𝜙||_𝑇≠||𝜙||_𝑇,𝑁 usually. But for mostly used orthogonal systems in [0,𝑇], they are equivalent, namely, for certain positive constants 𝑐_𝑙 and 𝑐₂, 𝑐₁𝜙_𝑇≤𝜙_𝑇,𝑁≤𝑐₂𝜙_𝑇. (21)

3.2. Existence and Uniqueness of the SHAM Solution

We consider the initial value problem (IVP) of dimension 𝑛 given as 𝐲̇(𝑡)=𝐟(𝑡,𝐲(𝑡)),𝐲𝑡₀=𝐲⁰,𝐲∶ℝ⟶ℝ^𝑛,𝐟∶ℝ×ℝ^𝑛⟶ℝ^𝑛, (22) 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 ℒ𝑦(𝑡)+𝒩𝑦(𝑡)=𝜙(𝑡), (23) 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 {𝑌_𝑛}^∞_𝑛=0 as 𝒩𝑦(𝑡)=^∞_𝑘=0𝑁_𝑘𝑦₀,𝑦₁,…,𝑦_𝑘,𝑌₀=𝑦₀,𝑌₁=𝑦₀+𝑦₁,⋮𝑌_𝑛=𝑦₀+𝑦₁+𝑦₂+⋯+𝑦_𝑛. (24)(25) 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: ℒ𝑦_𝑚(𝑡)−𝜒_𝑚𝑦_𝑚−1(𝑡)=ℏ𝐻(𝑡)𝑅_𝑚𝑦_𝑚−1(𝑡), (26) subject to the initial condition 𝑦_𝑚(0)=0, (27) where 𝑅_𝑚𝑦_𝑚−1=ℒ𝑦_𝑚−1+𝒩_𝑚−1𝑦₀,𝑦₁,…,𝑦_𝑚−1−1−𝜒_𝑚𝜙(𝑡). (28) Therefore, ℒ𝑦₁(𝑡)=ℏ𝐻(𝑡)ℒ𝑦₀+𝒩₀(𝑡,ℒ𝑦−𝜙)₂(𝑡)−𝑦₁(𝑡)=ℏ𝐻(𝑡)ℒ𝑦₁+𝒩₁,ℒ𝑦₃(𝑡)−𝑦₂(𝑡)=ℏ𝐻(𝑡)ℒ𝑦₂+𝒩₂,⋮ℒ𝑦_𝑚(𝑡)−𝑦_𝑚−1(𝑡)=ℏ𝐻(𝑡)ℒ𝑦_𝑚−1+𝒩_𝑚−1. (29) Upon summing these equations, we have ℒ𝑦_𝑚(𝑡)=ℏ𝐻(𝑡)^𝑚−1_𝑘=0ℒ𝑦_𝑘+^𝑚−1_𝑘=0𝒩_𝑘−𝜙(𝑡), (30) and from (25) we obtain ℒ𝑌_𝑚(𝑡)−𝑌_𝑚−1(𝑡(𝑡ℒ𝑌)=ℏ𝐻)_𝑚−1+𝒩𝑌_𝑚−1−𝜙(𝑡), (31) subject to the initial condition 𝑌_𝑚(0)=0. (32) Consequently, the collocation method is based on a solution 𝑌^𝑁(𝑡)∈𝒫_𝑁+1(0,𝑇), for (31) such that ℒ𝑌^𝑁_𝑚𝑡^𝑁_𝑇,𝑘−𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘=ℏ𝐻^𝑁𝑡^𝑁_𝑇,𝑘ℒ𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘+𝒩𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,k−𝜙^𝑁𝑡^𝑁_𝑇,𝑘, (33)(34) subject to the initial condition 𝑌^𝑁_𝑚(0)=0. (35)

Definition 2. A mapping 𝑓 of space 𝐿²(0,𝑇) into itself is said to satisfy a Lipschitz condition with Lipschitz constant 𝛾 if for any 𝑧 and 𝑧^∗, 𝑓(𝑧,𝑡)−𝑓𝑧^∗,𝑡≤𝛾𝑧−𝑧^∗. (36) 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 𝑌=𝑇(𝑌).(ii)For any 𝑌₀∈𝐾, the sequence {𝑌_𝑛}^∞_𝑛=0⊂𝐾 defined by 𝑌_𝑛+1=𝑇(𝑌_𝑛), 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 ℒ𝑌^𝑁_𝑚𝑡^𝑁_𝑇,𝑘=1+ℏ𝐻^𝑁𝑡^𝑁_𝑇,𝑘ℒ𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘+ℏ𝐻^𝑁𝑡^𝑁_𝑇,𝑘𝒩𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘−𝜙^𝑁𝑡^𝑁_𝑇,𝑘,𝑌0≤𝑘≤𝑁,𝑚≥1,_𝑁𝑚(0)=0,𝑚≥0. (37) Since 𝑓(𝑡,𝑧) satisfies the Lipschitz-continuous condition, then there exists a constant 𝛾≥0 such that 𝑓(𝑡,𝑧)−𝑓𝑡,𝑧^∗≤𝛾𝑧−𝑧^∗, (38) for all 𝑡∈[0,𝑇], and all 𝑧 and 𝑧^∗.
Now, for the problem (23), we choose 𝐿[𝑌(𝑡=𝑑)](𝑡)𝑁[𝑌(𝑡(𝑡)()𝑑𝑡𝑌+𝛼𝑌,)]=−𝛼𝑌−𝑓𝑡,𝑌 (39) and 𝜙(𝑡)≡0 where 𝛼(𝑡) is an arbitrary analytic function.
Let 𝑌^∼𝑁_𝑚(𝑡)=𝑌^𝑁_𝑚(𝑡)−𝑌^𝑁_𝑚−1(𝑡), then we have from (37) that ℒ𝑌^∼𝑁_𝑚𝑡^𝑁_𝑇,𝑘=1+ℏ𝐻𝑡^𝑁_𝑇,𝑘×ℒ𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘−𝑌^𝑁_𝑚−2𝑡^𝑁_𝑇,𝑘+ℏ𝐻𝑡^𝑁_𝑇,𝑘×𝒩𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘−𝒩𝑌^𝑁_𝑚−2𝑡^𝑁_𝑇,𝑘0≤𝑘≤𝑁,𝑚≥1 (40) or according to the definitions of 𝐿[𝑌(𝑡)] and 𝑁[𝑌(𝑡)], 𝑑𝑑𝑡^∼𝑌^𝑁_𝑚𝑡^𝑁_𝑇,𝑘+𝛼𝑡𝑌^𝑁_𝑇,𝑘^∼𝑁_𝑚=1+ℏ𝐻𝑡^𝑁_𝑇,𝑘𝑑𝑑𝑡^∼𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘+𝛼𝑡𝑌^𝑁_𝑇,𝑘^∼𝑁_𝑚−1−ℏ𝐻𝑡^𝑁_𝑇,𝑘𝑓𝑡^𝑁_𝑇,𝑘,𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘+𝛼𝑡^𝑁_𝑇,𝑘𝑌^𝑁_𝑚−1−𝑓𝑡^𝑁_𝑇,𝑘,𝑌^𝑁_𝑚−2𝑡^𝑁_𝑇,𝑘−𝛼𝑡^𝑁_𝑇,𝑘𝑌^𝑁_𝑚−2,0≤𝑘≤𝑁,𝑚≥1, (41) from where 𝑑𝑑𝑡^∼𝑌^𝑁_𝑚𝑡^𝑁_𝑇,𝑘=1+ℏ𝐻𝑡^𝑁_𝑇,𝑘𝑑𝑑𝑡^∼𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘+𝛼𝑡𝑌^𝑁_𝑇,𝑘^∼𝑁_𝑚−1−𝛼𝑡𝑌^𝑁_𝑇,𝑘^∼𝑁_𝑚−ℏ𝐻𝑡^𝑁_𝑇,𝑘×𝑓𝑡^𝑁_𝑇,𝑘,𝑌^𝑁_𝑚−1𝑡^𝑁_𝑇,𝑘−𝑓𝑡^𝑁_𝑇,𝑘,𝑌^𝑁_𝑚−2𝑡^𝑁_𝑇,𝑘0≤𝑘≤𝑁,𝑚≥1. (42) It is obvious that 𝑌^∼𝑁_𝑚(0)=0. clearly, 𝑌^∼𝑁_𝑚(𝑇)²=2𝑌^∼𝑁_𝑚,𝑑𝑑𝑡^∼𝑌^𝑁_𝑚𝑇≤2𝑌^∼𝑁_𝑚𝑇𝑑𝑑𝑡^∼𝑌^𝑁_𝑚𝑇. (43) Furthermore, for any 𝑡∈[0,𝑇], 𝑌^∼𝑁_𝑚(𝑡)²=𝑌^∼𝑁_𝑚(𝑇)²−^𝑇_𝑡𝑑𝑑𝑥^∼𝑌^𝑁_𝑚(𝑥)²≤𝑑𝑥_∼𝑌^𝑁_𝑚(𝑇)²+2𝑌^∼𝑁_𝑚𝑇𝑑𝑑𝑥^∼𝑌^𝑁_𝑚𝑇. (44) Integrating the previously mentioned with respect to 𝑡 yields that 𝑌^∼𝑁_𝑚²_𝑇≤𝑇𝑌^∼𝑁_𝑚(𝑇)²+2𝑇𝑌^∼𝑁_𝑚𝑇𝑑𝑑𝑡^∼𝑌^𝑁_𝑚𝑇, (45) from where 𝑌^∼𝑁_𝑚(𝑇)²≥1𝑇^∼𝑌^𝑁_𝑚²_𝑇−2𝑌^∼𝑁_𝑚𝑇𝑑𝑑𝑡^∼𝑌^𝑁_𝑚𝑇. (46) Using (46) and (43) we have 𝑌^∼𝑁_𝑚𝑇𝑑≤4𝑇𝑑𝑡^∼𝑌^𝑁_𝑚𝑇,𝑚≥0. (47) Let 𝛼=max_{𝑡∈[0,𝑇]}|𝛼(𝑡)| and 𝐻=max_{𝑡∈[0,𝑇]}|𝐻(𝑡)|. Therefore, a combination of (21), (36), and (43) results in 𝑑𝑑𝑡^∼𝑌^𝑁_𝑚𝑇≤𝑑𝑑𝑡^∼𝑌^𝑁_{𝑚𝑇,𝑁}≤(|ℏ|𝐻𝑑1+)𝑑𝑡^∼𝑌^𝑁_{𝑚−1𝑇,𝑁}+𝑐𝛾|ℏ|𝐻𝑌^∼𝑁_{𝑚−1𝑇,𝑁}+𝛼𝑌^∼𝑁_{𝑚𝑇,𝑁}+𝛼𝑌^∼𝑁_{𝑚−1𝑇,𝑁}≤𝑐₂(|ℏ|𝐻𝑑1+)𝑑𝑡^∼𝑌^𝑁_𝑚−1𝑇+𝑐₂𝐶𝛾|ℏ|𝐻𝑌^∼𝑁_𝑚−1𝑇+𝑐₂𝛼𝑌^∼𝑁_𝑚𝑇+𝑐₂𝛼𝑌^∼𝑁_𝑚−1𝑇. (48) Using (47) and (48) results in 𝑌^∼𝑁_𝑚𝑇≤4𝑐₂𝑇(1+|ℏ|𝐻)1−4𝑐₂d𝛼𝑇d𝑡^∼𝑌^𝑁_𝑚−1𝑇+4𝑐₂𝑇𝛼+𝑐𝛾|ℏ|𝐻1−4𝑐₂𝛼𝑇^∼𝑌^𝑁_𝑚−1𝑇, (49) where 𝑐 is a positive constant. Then, by (47) with 𝑚−1, instead of 𝑚 and multiplying the resulting inequality by (𝑐₂(1+|ℏ|𝐻))/(1−4𝑐₂𝛼𝑇)≥0, we have 𝑐₂(1+|ℏ|𝐻)1−4𝑐₂𝛼𝑇^∼𝑌^𝑁_𝑚−1𝑇≤4𝑐₂𝑇(1+|ℏ|𝐻)1−4𝑐₂d𝛼𝑇d𝑡^∼𝑌^𝑁_𝑚−1𝑇. (50) Subtracting (50) from (49), after simplifying, we obtain 𝑌^∼𝑁_𝑚𝑇≤𝑐₂(1+|ℏ|𝐻)+4𝑐₂𝑇𝛼+𝑐𝛾|ℏ|𝐻1−4𝑐₂𝛼𝑇^∼𝑌^𝑁_𝑚−1𝑇. (51) Therefore, if (𝑐₂(1+|ℏ|𝐻)+4𝑐₂𝑇(𝛼+𝑐𝛾|ℏ|𝐻))/(1−4𝑐₂𝛼𝑇)≤𝛽<1, then ||𝑌^∼𝑁_𝑚||_𝑇→0 as 𝑚→∞. According to Theorem 3, it implies the existence and uniqueness of solution of (33).

4. 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 Λ=[𝑡⁰,𝑡^𝐹], we divide the interval Λ into subintervals Λ_𝑖=[𝑡^𝑖−1,𝑡^𝑖], 𝑖=1,2,3,…,𝐹 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 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 ℒ_𝑟𝑦^𝑖_𝑟,𝑚(𝑡)−𝜒_𝑚𝑦^𝑖_{𝑟,𝑚−1}(𝑡)=ℏ_𝑟𝐻_𝑟(𝑡)𝑅^𝑖_{𝑟,𝑚−1},𝑡∈𝑡^𝑖−1,𝑡^𝑖 (52) subject to the initial condition 𝑦^𝑖_𝑟,𝑚𝑡^𝑖−1=0. (53) The initial approximations for solving (8) are obtained as solutions of the system ℒ𝑦^𝑖_𝑟,0(𝑡)=𝜙_𝑟(𝑡),𝑡∈𝑡^𝑖−1,𝑡^𝑖 (54) subject to the initial condition 𝑦^𝑖_𝑟,0𝑡^𝑖−1=𝑦_𝑟^𝑖−1𝑡^𝑖−1. (55) 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 𝑡𝑡=^𝑖−𝑡^𝑖−1𝑥2+𝑡^𝑖+𝑡^𝑖−12 (56) 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 𝑥^𝑖₀=−1, 𝑥^𝑖_𝑁=1 and for 1≤𝑗≤𝑁−1, 𝑥^𝑖_𝑗 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 𝑑𝑦^𝑖_𝑟,𝑚=𝑑𝑡^𝑁_𝑘=0𝐃_𝑗𝑘𝑦^𝑖_𝑟,𝑚𝑥^𝑖_𝑘=𝐃𝐘^𝑖_𝑟,𝑚,𝑗=0,1,…,𝑁, (57) where 𝐃=2𝐷/(𝑡^𝑖−𝑡ⁱ⁻¹) and 𝐘^𝑖_𝑟,𝑚=[𝑦^𝑖_𝑟,𝑚(𝑥^𝑖₀),𝑦^𝑖_𝑟,𝑚(𝑥^𝑖₂),…,𝑦^𝑖_𝑟,𝑚(𝑥^𝑖_𝑁)]^𝑇 is the vector function at the collocation points 𝑥^𝑖_𝑗. The matrix 𝐷 is of size (𝑁+1)×(𝑁+1) and its entries are defined [31, 32] as 𝐷_𝑗𝑘=−𝑁(𝑁+1)4𝑁(),if𝑗=𝑘=0,𝑁+14𝐿,if𝑗=𝑘=𝑁,_𝑁𝑡_𝑗𝐿_𝑁𝑡_𝑘1𝑡_𝑗−𝑡_𝑘,if𝑗≠𝑘,0,otherwise. (58)

Applying the the Legendre spectral collocation method in (52)–(55) gives 𝐀𝐘^𝑖_𝑟,0=Φ_𝑟,𝐘^𝑖_𝑟,0𝑥_𝑁^𝑖−1=𝐘_𝑟^𝑖−1𝑥_𝑁^𝑖−1,𝐀𝐘^𝑖_𝑟,𝑚−𝜒_𝑚𝐘^𝑖_{𝑟,𝑚−1}=ℏ_𝑟𝐑^𝑖_{𝑟,𝑚−1},𝐘^𝑖_𝑟,𝑚𝑥_𝑁^𝑖−1=0, (59)(60) where 𝐀 is an (𝑁+1)×(𝑁+1) matrix that comes from transforming the linear operator ℒ using the derivative matrix 𝐷, Φ_𝑟, and 𝐑^𝑖_{𝑟,𝑚−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 𝐘^𝑖_𝑟,0 which is the solution of (59), the recurrence formula 𝐘^𝑖_𝑟,𝑚=𝜒_𝑚𝐘^𝑖_{𝑟,𝑚−1}+ℏ_𝑟𝐀⁻¹𝐑^𝑖_{𝑟,𝑚−1} (61) can be used to obtain the solution 𝑦^𝑖_𝑟(𝑡) using (9) in the interval [𝑡^𝑖,𝑡^𝑖−1]. The solution approximating 𝑦_𝑟(𝑡) in the entire interval [𝑡⁰,𝑡^𝐹] is given by 𝑦_𝑟(𝑡)=𝑦¹_𝑟(𝑡),𝑡∈𝑡⁰,𝑡¹𝑦²_𝑟(𝑡),𝑡∈𝑡¹,𝑡²⋮𝑦^𝐹_𝑟(𝑡),𝑡∈𝑡^𝐹−1,𝑡^𝐹. (62) It must be noted that when 𝐹=1, the proposed piecewise-spectral homotopy analysis method (PSHAM) becomes equivalent to the original SHAM algorithm.

5. Numerical Experiments

To demonstrate the applicability of the proposed piecewise-spectral 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.

5.1. Gensio-Tesi System

The Genesio-Tesi system [15] considered in this paper is given as follows: 𝑑𝑦₁𝑑𝑡=𝑦₂,𝑦₁(0)=𝑦⁰₁,𝑑𝑦₂𝑑𝑡=𝑦₃,𝑦₂(0)=𝑦⁰₂,𝑑𝑦₃𝑑𝑡=−𝑐𝑦₁−𝑏𝑦₂−𝑎𝑦₃+𝑦²₁,𝑦₃(0)=𝑦⁰₃. (63) In this example, the parameters used in the SHAM and PSHAM algorithms are ℒ_𝑟=𝑑0𝑑𝑑𝑡−10𝑑𝑑𝑡−1𝑐𝑏,00𝑑𝑡+𝑎,𝐀=𝐃−𝐼𝐎𝐎𝐃−𝐼𝑐𝐼𝑏𝐼𝐃+𝑎𝐼ℱ=−𝑦²₁000,𝑅,𝜙=_{𝑖𝑟,𝑚−1}=ℒ_𝑟𝑦^𝑖_{𝑟,𝑚−1}+𝑄^𝑖_{𝑟,𝑚−1},𝑄^𝑖_{𝑟,𝑚−1}=00−∑^𝑚−1_𝑗=0𝑦^𝑖_1,𝑗𝑦^𝑖_{1,𝑚−1−𝑗}. (64) With these definitions, the PSHAM algorithm gives 𝐀𝐘^𝑖_𝑟,𝑚=𝜒_𝑚+ℏ_𝑟𝐀𝐘^𝑖_{𝑟,𝑚−1}+ℏ_𝑟𝐐^𝑖_{𝑟,𝑚−1}. (65) Because the right-hand 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: 𝑑𝑦₁𝑑𝑡=𝑠−𝛼𝑦₁+𝑟^∗𝑦₁𝑦1−¹+𝑦₂𝑇_max−𝑘𝑦₁𝑦₃,𝑑𝑦₂𝑑𝑡=𝑘𝑦₁𝑦₃−𝛽𝑦₂,𝑑𝑦₃𝑑𝑡=𝑁₁𝛽𝑦₂−𝛾𝑦₃, (66) with the initial conditions 𝑦₁(0)=𝑟₁,𝑦₂(0)=𝑟₂,𝑦₃(0)=𝑟₃. (67)

In this model 𝑦₁, 𝑦₂, and 𝑦₃ 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−((𝑦₁+𝑦₂)/(𝑇_max))), and proliferation of infected CD4⁺ T cells is neglected. The term 𝑘𝑦₂𝑦₃ 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 𝑁₁ virus particles during its life time, including any of its daughter cells [21, 33]. Throughout this paper, we set 𝑠=0.1,𝛼=0.02,𝛽=0.3,𝑟^∗=3,𝛾=2.4,𝑘=0.0027,𝑇_max=1500,𝑁₁=10. (68) In this example, the parameters used in the SHAM and PSHAM algorithms are ℒ_𝑟=𝑑−𝑟𝑑𝑡^∗0𝑑−𝛼00𝑑𝑡+𝛽00−𝑁₁𝛽𝑑,𝑟𝑑𝑡+𝛾𝐀=𝐃−𝐼^∗−𝛼𝐎𝐎𝐎𝐃+𝛽𝐼𝐎𝐎−𝑁₁,𝛽𝐃+𝛾𝐼ℱ=𝑟𝑦₁𝑦₁+𝑦₂𝑇_max+𝑘𝑦₁𝑦₃−𝑘𝑦₁𝑦₃0𝑠00,𝑅,𝜙=_{𝑖𝑟,𝑚−1}=ℒ_𝑟𝑦^𝑖_{𝑟,𝑚−1}+𝑄^𝑖_{𝑟,𝑚−1},𝑄^𝑖_{𝑟,𝑚−1}=𝑟∑^∗𝑚−1_𝑗=0𝑦^𝑖_1,𝑗𝑦^𝑖_{1,𝑚−1−𝑗}+𝑦^𝑖_{2,𝑚−1−𝑗}+𝑘∑^𝑚−1_𝑗=0𝑦^𝑖_1,𝑗𝑦^𝑖_{3,𝑚−1−𝑗}−𝑘∑^𝑚−1_𝑗=0𝑦^𝑖_1,𝑗𝑦^𝑖_{3,𝑚−1−𝑗}0. (69) With these definitions, the PSHAM algorithm gives 𝐀𝐘^𝑖_𝑟,𝑚=𝜒_𝑚+ℏ_𝑟𝐀𝐘^𝑖_{𝑟,𝑚−1}+ℏ_𝑟𝐐^𝑖_{𝑟,𝑚−1}−1−𝜒_𝑚Φ. (70) Because the right-hand 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 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 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 ℏ_𝑟=−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.

Table 1 

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

Figure 1 gives the graphical profiles of the classes 𝑦₁(𝑡),𝑦₂(𝑡),𝑦₃(𝑡) in the time interval [0,40]. The initial conditions used to generate the graphs are 𝑦₁(0)=0.2,𝑦₂(0)=−0.3,  and  𝑦₃(0)=0.1. 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)

Figure 1 

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

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.

(a)

(b)

(c)

Figure 2 

Phase portraits of the Genesio-Tesi system using the PSHAM approach for 𝑐=6, 𝑏=2.92, and 𝑎=1.2.

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 𝑦₁(0)=0.1,𝑦₂(0)=0, and 𝑦₃(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.

Table 2 

Numerical comparison between PSHAM and ode45 for the HIV infection of CD4⁺ T cells model.

Figure 3 gives a comparison between the PSHAM and the ode45 results (small dots) for the variation of the classes 𝑦₁(𝑡),𝑦₂(𝑡),𝑦₃(𝑡) 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.

(a)

(b)

(c)

Figure 3 

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

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)

Figure 4 

Phase portraits of the HIV model using the PSHAM approach.

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 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.