A Reliable Analytical Method for Solving Higher-Order Initial Value Problems

Abstract

In this article, a new analytical method has been devised to solve higher-order initial value problems for ordinary differential equations. This method was implemented to construct a series solution for higher-order initial value problems in the form of a rapidly convergent series with easily computable components using symbolic computation software. The proposed method is based on the Taylor series expansion which constructs an analytical solution in the form of a polynomial and reproduces the exact solution when the solution is polynomial. This technique is applied to a few test examples to illustrate the accuracy, efficiency, and applicability of the method. The results reveal that the method is very effective, straightforward, and simple.

1. Introduction

The study of nonlinear problems is of crucial importance not only in mathematics but also in physics, engineering, economic, and other disciplines, since most phenomena in our world are essentially nonlinear and are described by nonlinear equations. It is very difficult to solve nonlinear problems and, in general, it is often more difficult to get an analytic approximation than a numerical one for a given nonlinear problem. In this paper, we focus on finding approximate solution to higher-order initial value problems (IVPs), which are a combination of higher-order ordinary differential equations subject to given initial conditions. In fact, accurate and fast numerical solution of higher-order IVPs is of great importance due to its wide application in scientific and engineering research.

In the fields of engineering and science, we come across physical and natural phenomena which, when represented by mathematical models, happen to be differential equations. For example, simple harmonic motion, equation of motion, deflection of a beam, and so forth are represented by differential equations. Hence, the solution of differential equations is a necessity in such studies. There are a number of differential equations which we studied in calculus to get closed form solutions. But all differential equations do not possess closed form of finite form solutions. Even if they possess closed form solutions, we do not know the method of getting it. In such situations, depending on the need of the hour, we go in for numerical solutions of differential equations. In researches, especially after the advent of computer, the numerical solutions of differential equations have become easy for manipulations.

In this work, we introduce a new analytical method; we call it residual power series (RPS) [1, 2] to find out series solution to strongly linear and nonlinear higher-order IVPs. The RPS method is effective and easy to use for solving higher-order IVPs without linearization, perturbation, or discretization. This method constructs an analytical approximate solution in the form of polynomial. The RPS method is different from the traditional higher-order Taylor series approach; the Taylor series approach is computationally expensive for large orders. The RPS method is an alternative procedure for obtaining analytic Taylor series solution of higher-order IVPs. By using residual error concept, we get a series solution, in practice a truncated series solution. As we will see later, the exact solution is available when the solution is polynomial. Moreover, the solution and all its derivatives are applicable for each arbitrary point in the given interval. On the other hand, RPS does not require any converting while switching from the first-order to the higher-order; as a result, the method can be applied directly to the given problem by choosing appropriate values for the initial guess approximation.

The purpose of this paper is to explain the idea of the RPS method and test the efficiency and applicability of the method. To reach our goal, we will construct symbolic approximate RPS solution for higher-order ordinary differential equation of the following form: 𝑥^(𝑛)(𝑡)=𝑓𝑡,𝑥(𝑡),𝑥^′(𝑡),…,𝑥^(𝑛−1)(𝑡), (1) subject to the initial conditions 𝑥𝑡₀=𝑥₁,𝑥^′𝑡₀=𝑥₂,…,𝑥^(𝑛−1)𝑡₀=𝑥_𝑛, (2) where 𝑡∈[𝑡₀,𝑡₀+𝑎],  𝑓∶[𝑡₀,𝑡₀+𝑎]×ℝ^𝑛→ℝ^𝑛 is nonlinear function of 𝑡,𝑥(𝑡),𝑥^′(𝑡),…,𝑥^(𝑛−1)(𝑡), 𝑥(𝑡) is an unknown function of independent variable 𝑡 to be determined, and 𝑡₀,𝑎,𝑥_𝑖,𝑖=1,2,…,𝑛, are real finite constants with 𝑎>0. Throughout this paper, we assume that 𝑓 and 𝑥 are analytic functions on the given interval.

In the recent years, there are several methods used to find approximate solutions to higher-order IVPs, such as homotopy analysis method [3], differential transformation method [4], Adomian decomposition method [5], Jacobi-Gauss collocation method [6], Schauder bases method [7], or other techniques [8–12]. On the other aspect as well, the applications of other versions of series solutions to linear and nonlinear problems can be found in [13–18] and for numerical solvability of different categories of differential equations, one can consult the references [19, 20].

The outline of the paper is as follows. In the next section, we present the basic idea of RPS method. In Section 3, convergence theorem and error analysis are discussed in order to capture the behavior of solution. In Section 4, numerical examples are given to illustrate the capability of proposed method. This paper ends in Section 5 with some concluding remarks.

2. The Idea of the RPS Method

In this section, we employ our algorithm of RPS to find out series solution for higher-order ordinary differential equation subject to given initial conditions. We first formulate and analyze the RPS method for solving such equations.

The RPS method consists in expressing the solution of (1) and (2) as a power series expansion about the initial point 𝑡₀ [1, 2]. To achieve our goal, we suppose that this solution takes the form 𝑥(𝑡)=^∞_𝑚=0𝑥_𝑚(𝑡), (3) where 𝑥_𝑚 are terms of approximations and are given as 𝑥_𝑚(𝑡)=𝑐_𝑚(𝑡−𝑡₀)^𝑚,𝑚=0,1,2,….

Obviously, when 𝑚=0,1,2,…,𝑛−1, since the terms 𝑥₀(𝑡),𝑥^′₁(𝑡),𝑥^′′₂(𝑡),…,𝑥^(𝑛−1)_𝑛−1(𝑡) satisfy the initial conditions (2) as 𝑥𝑡₀=𝑐₀=𝑥₀𝑡₀,𝑥^′𝑡₀=𝑐₁=𝑥^′₁𝑡₀,𝑥^′′𝑡₀=2𝑐₂=𝑥^′′₂𝑡₀,⋮𝑥^(𝑛−1)𝑡₀=(𝑛−1)!𝑐_𝑛−1=𝑥^(𝑛−1)_𝑛−1𝑡₀, (4) we have the initial guess approximation of 𝑥(𝑡) which is as follows: 𝑥_initial(𝑡)=𝑥𝑡₀+1𝑥1!^′𝑡₀𝑡−𝑡₀+1𝑥2!^′′𝑡₀𝑡−𝑡₀²1+⋯+()!𝑥𝑛−1^(𝑛−1)𝑡₀𝑡−𝑡₀^𝑛−1. (5)

On the other aspect as well, if we choose 𝑥_initial(𝑡) as initial guess approximation of 𝑥(𝑡), then we can calculate 𝑥_𝑚(𝑡) for 𝑚=𝑛,𝑛+1,𝑛+2…and approximate the solution 𝑥(𝑡) of (1) and (2) by the 𝑘th-truncated series 𝑥^𝑘(𝑡)=𝑥_initial(𝑡)+^𝑘_𝑚=𝑛𝑐_𝑚𝑡−𝑡₀^𝑚. (6)

Prior to applying RPS method for solving (1) and (2), we substituting the 𝑘th-truncated series 𝑥^𝑘(𝑡) into (1) to obtain the following definition for the 𝑘th residual function: Res^𝑘(𝑡)=^𝑘_𝑚=𝑛𝑚(𝑚−1)⋯(𝑚−𝑛+1)𝑐_𝑚𝑡−𝑡₀^𝑚−𝑛−𝑓𝑡,^𝑘_𝑚=0𝑐_𝑚𝑡−𝑡₀^𝑚,^𝑘_𝑚=1𝑚𝑐_𝑚𝑡−𝑡₀^𝑚−1,_𝑘𝑚=2𝑚(𝑚−1)𝑐_𝑚𝑡−𝑡₀^𝑚−2,…,_{𝑘𝑚=𝑛−1}𝑚()⋯()×𝑚−1𝑚−𝑛+2𝑡−𝑡₀^{𝑚−𝑛+1} (7) and the following definition of ∞th residual function: Res^∞(𝑡)=lim_𝑘→∞Res^𝑘(𝑡).

It is clear that Res^∞(𝑡)=0 for each 𝑡∈[𝑡₀,𝑡₀+𝑎]. This shows that the residual function Res^∞(𝑡) is infinitely many times differentiable at 𝑡=𝑡₀. On the other hand, it easy to see that (𝑑^𝑠/𝑑𝑡^𝑠)Res^∞(𝑡₀)=(𝑑^𝑠/𝑑𝑡^𝑠)Res^𝑘(𝑡₀)=0, for each 𝑠=1,2,…,𝑘. In fact, this relation is a fundamental rule in RPS method and its applications [1, 2].

To obtain the 𝑛th approximate solution, we put 𝑘=𝑛 and 𝑥^𝑛(𝑡)=∑^𝑛_𝑚=0𝑐_𝑚(𝑡−𝑡₀)^𝑚. On the other hand, we substitute 𝑡=𝑡₀ into (7) to get Res^𝑛𝑡₀=𝑛!𝑐_𝑛−𝑓𝑡₀,𝑐₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1. (8)

By considering the fact that Res^𝑛(𝑡₀)=0, yield 𝑐_𝑛=(1/𝑛!)𝑓(𝑡₀,𝑐₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1). Thus, using 𝑛th-truncated series (6), the 𝑛th approximation of (1) and (2) can be written as 𝑥^𝑛(𝑡)=𝑥_initial(𝑡)+1𝑓𝑡𝑛!⁰,𝑐₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1𝑡−𝑡₀^𝑛. (9)

Similarly, to find the (𝑛+1)th approximate solution, we put 𝑘=𝑛+1 and 𝑥^𝑛+1(𝑡)=∑^𝑛+1_𝑚=0𝑐_𝑚(𝑡−𝑡₀)^𝑚. On the other direction, we differentiate both sides of (7) with respect to 𝑡 and substitute 𝑡=𝑡₀ to get 𝑑𝑑𝑡Res^𝑛+1𝑡₀=(𝑛+1)!𝑐_𝑛+1−𝑑𝑓𝑡𝑑𝑡⁰,𝑐₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1−𝑐₁𝑑𝑑𝑐₀𝑓𝑡₀,𝑧₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1−2𝑐₂𝑑𝑑𝑧₁𝑓𝑡₀,𝑐₀,𝑧₁,…,(𝑛−1)!𝑐_𝑛−1−⋯−𝑛!𝑐_𝑛𝑑𝑑𝑧_𝑛−1𝑓𝑡₀,𝑐₀,𝑐₁,…,𝑧_𝑛−1, (10) where 𝑧_𝑖−1=(𝑖−1)!𝑐_𝑖−1,𝑖=1,2,…,𝑛. Thus, using the fact that (𝑑/𝑑𝑡)Res^𝑛+1(𝑡₀)=0, yield the following value of 𝑐_𝑛+1: 𝑐_𝑛+1=1()!𝑑𝑛+1𝑓𝑡𝑑𝑡⁰,𝑐₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1+𝑐₁𝑑𝑑𝑐₀𝑓𝑡₀,𝑧₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1+2𝑐₂𝑑𝑑𝑧₁𝑓𝑡₀,𝑐₀,𝑧₁,…,(𝑛−1)!𝑐_𝑛−1+⋯+𝑛!𝑐_𝑛𝑑𝑑𝑧_𝑛−1𝑓𝑡₀,𝑐₀,𝑐₁,…,𝑧_𝑛−1, (11) where 𝑧_𝑖−1=(𝑖−1)!𝑐_𝑖−1,𝑖=1,2,…,𝑛. Hence, the (𝑛+1)th approximation of (1) and (2) can be written as 𝑥^𝑛+1(𝑡)=𝑥_initial(𝑡)+1𝑓𝑡𝑛!⁰,𝑐₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1×𝑡−𝑡₀^𝑛+1()!×𝑑𝑛+1𝑓𝑡𝑑𝑡⁰,𝑐₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1+𝑐₁𝑑𝑑𝑐₀𝑓𝑡₀,𝑧₀,𝑐₁,…,(𝑛−1)!𝑐_𝑛−1+2𝑐₂𝑑𝑑𝑧₁𝑓𝑡₀,𝑐₀,𝑧₁,…,(𝑛−1)!𝑐_𝑛−1+⋯+𝑛!𝑐_𝑛𝑑𝑑𝑧_𝑛−1𝑓𝑡₀,𝑐₀,𝑐₁,…,𝑧_𝑛−1×𝑡−𝑡₀^𝑛+1, (12) where 𝑧_𝑖−1=(𝑖−1)!𝑐_𝑖−1,𝑖=1,2,…,𝑛.

This procedure can be repeated till the arbitrary order coefficients of RPS solution for (1) and (2) are obtained. Moreover, higher accuracy can be achieved by evaluating more components of the solution. In other words, choose large 𝑘 in the truncation series (6).

3. Convergence Theorem and Error Analysis

Numerical technique is widely used by scientists and engineers to solve their problems. A major advantage for numerical technique is that a numerical answer can be obtained even when a problem has no analytical solution. However, result from numerical analysis is an approximation, in general, which can be made as accurate as desired. The reliability of the numerical result will depend on an error estimate or bound; therefore, the analysis of error and the sources of error in numerical methods is also a critically important part of the study of numerical technique. In this section, we study the convergence of RPS method.

Taylor’s theorem allows us to represent exactly and fairly general functions in terms of polynomials with a known, specified, and boundable error. Next theorem will guarantee convergence to exact analytic solution of (1) and (2) using RPS method.

Theorem 1. Suppose that 𝑥(𝑡) is the exact solution of (1) and (2). Then, the approximate solution obtained using RPS method is just the Taylor expansion of 𝑥(𝑡) about 𝑡=𝑡₀.

Proof. Assume that the approximate solution of (1) and (2) is 𝑥(𝑡)=𝑐₀+𝑐₁𝑡−𝑡₀+𝑐₂𝑡−𝑡₀²+⋯+𝑐_𝑛−1𝑡−𝑡₀^𝑛−1+𝑐_𝑛𝑡−𝑡₀^𝑛+𝑐_𝑛+1𝑡−𝑡₀^𝑛+1+⋯. (13)
In order to prove the theorem, it is enough to show that the coefficients 𝑐_𝑚 in (13) take the form 𝑐_𝑚=1𝑥𝑚!^(𝑚)𝑡₀,𝑚=0,1,2,…, (14) where 𝑥(𝑡) is the exact solution of (1) and (2). It is clear that for 𝑚=0,1,2,…,𝑛−1, (4) gives the following values for 𝑐_𝑚: 𝑐_𝑚=1𝑥𝑚!^(𝑚)𝑡₀,𝑚=0,1,2,…,𝑛−1. (15)
Moreover, for 𝑚=𝑛, we substitute 𝑡=𝑡₀ into (1); we obtain 𝑥^(𝑛)𝑡₀=𝑓𝑡,𝑥𝑡₀,𝑥^′𝑡₀,𝑥^′′𝑡₀,…,𝑥^(𝑛−1)𝑡₀. (16)
On the other aspect as well, from (13) and (15), we can write 𝑥(𝑡)=𝑥𝑡₀+1𝑥1!^′𝑡₀𝑡−𝑡₀+1𝑥2!^′′𝑡₀𝑡−𝑡₀²1+⋯+()!𝑥𝑛−1^(𝑛−1)𝑡₀𝑡−𝑡₀^𝑛−1+𝑐_𝑛𝑡−𝑡₀^𝑛+𝑐_𝑛+1𝑡−𝑡₀^𝑛+1+⋯. (17)
Also, by substituting (17) into (1) and then setting 𝑡=𝑡₀, we get 𝑛!𝑐_𝑛=𝑓𝑡,𝑥𝑡₀,𝑥^′𝑡₀,𝑥^′′𝑡₀,…,𝑥^(𝑛−1)𝑡₀=𝑥^(𝑛)𝑡₀. (18)
Further, for 𝑚=𝑛+1, differentiating both sides of (1) with respect to 𝑡 and then substituting 𝑡=𝑡₀, we obtain 𝑥^(𝑛+1)(𝑡)=𝑑𝑓(𝑡)𝑑𝑡𝑡,𝑥,𝑥^′(𝑡),…,𝑥^(𝑛−1)(𝑡)+𝑥^′(𝑡)𝑑𝑑𝑤₀𝑓𝑡,𝑥(𝑡),𝑥^′(𝑡),…,𝑥^(𝑛−1)(𝑡)+𝑥^′′(𝑡)𝑑𝑑𝑤₁𝑓𝑡,𝑥(𝑡),𝑥^′(𝑡),…,𝑥^(𝑛−1)(𝑡)+⋯+𝑥^(𝑛)(𝑡)𝑑𝑑𝑤_𝑛−1𝑓𝑡,𝑥(𝑡),𝑥^′(𝑡),…,𝑥^(𝑛−1)(𝑡), (19) where 𝑤_𝑖=𝑥^(𝑖),𝑖=0,1,2,…,𝑛−1. By substituting 𝑡=𝑡₀ into (19) and using (14) and (18), we can conclude that 𝑥^(𝑛+1)𝑡₀=𝑑𝑓𝑡𝑑𝑡⁰,𝑥𝑡₀,𝑥^′𝑡₀,…,𝑥^(𝑛−1)𝑡₀+𝑐₁𝑑𝑑𝑤₀𝑓𝑡₀,𝑥𝑡₀,𝑥^′𝑡₀,…,𝑥^(𝑛−1)𝑡₀+2𝑐₂𝑑𝑑𝑤₁𝑓𝑡₀,𝑥𝑡₀,𝑥^′𝑡₀,…,𝑥^(𝑛−1)𝑡₀+⋯+𝑛!𝑐_𝑛𝑑𝑑𝑤_𝑛−1𝑓𝑡₀,𝑥𝑡₀,𝑥^′𝑡₀,…,𝑥^(𝑛−1)𝑡₀, (20) where 𝑤_𝑖=𝑥^(𝑖),𝑖=0,1,2,…,𝑛−1. According to (17) and (18), we can write the approximation of (1) and (2) as follows: 𝑥(𝑡)=𝑥𝑡₀+1𝑥1!^′𝑡₀𝑡−𝑡₀+1𝑥2!^′′𝑡₀𝑡−𝑡₀²1+⋯+()!𝑥𝑛−1^(𝑛−1)𝑡₀𝑡−𝑡₀^𝑛−1+1𝑥𝑛!^(𝑛)𝑡₀𝑡−𝑡₀^𝑛+𝑐_𝑛+1𝑡−𝑡₀^𝑛+1+⋯. (21)
On the other hand, by substituting (21) into (19) and setting 𝑡=𝑡₀, we obtain (𝑛+1)!𝑐_𝑛+1=𝑑𝑓𝑡𝑑𝑡⁰,𝑥𝑡₀,𝑥^′𝑡₀,…,𝑥^(𝑛−1)𝑡₀+𝑐₁𝑑𝑑𝑤₀𝑓𝑡₀,𝑥𝑡₀,𝑥^′𝑡₀,…,𝑥^(𝑛−1)𝑡₀+2𝑐₂𝑑𝑑𝑤₁𝑓𝑡₀,𝑥𝑡₀,𝑥^′𝑡₀,…,𝑥^(𝑛−1)𝑡₀+⋯+𝑛!𝑐_𝑛𝑑𝑑𝑤_𝑛−1𝑓𝑡₀,𝑥𝑡₀,𝑥^′𝑡₀,…,𝑥^(𝑛−1)𝑡₀, (22) where 𝑤_𝑖=𝑥^(𝑖),𝑖=0,1,2,…,𝑛−1. Finally, by comparing (20) and (22), we can conclude that 𝑐_𝑛+1=(1/(𝑛+1)!)𝑥^(𝑛+1)(𝑡₀). By continuing the above procedure, we can easily prove (14) for 𝑚=𝑛+2,𝑛+3,𝑛+4,…. So, the proof of the theorem is complete.

Corollary 2. If 𝑥(𝑡) is polynomial, then the RPS method will be obtained the exact solution.

It will be convenient to have a notation for the error in the approximation 𝑥(𝑡)≈𝑥^𝑘(𝑡). Accordingly, we will let Rem^𝑘(𝑡) denote the difference between 𝑥(𝑡) and its 𝑘th Taylor polynomial obtained from the RPS method; that is, Rem^𝑘(𝑡)=𝑥(𝑡)−𝑥^𝑘(𝑡)=^∞_𝑚=𝑘+1𝑥^(𝑚)𝑡₀𝑚!𝑡−𝑡₀^𝑚. (23)

The functions Rem^𝑘(𝑡) are called the 𝑘th remainder for the RPS approximation of 𝑥(𝑡). In fact, it often happens that the remainders Rem^𝑘(𝑡) become smaller and smaller, approaching zero, as 𝑘 gets larger.

4. Numerical Results and Discussion

Numerical methods tend to emphasize the implementation of algorithms. The aim of numerical methods is therefore to provide systematic methods for solving problems in a numerical form. The process of solving problems generally involves starting from an initial data, using high precision digital computers, following the steps in the algorithms, and finally obtaining the results. Often the numerical data and the methods used are approximate ones. In this section, we consider six examples to demonstrate the performance and efficiency of the present technique. Throughout this paper, all the symbolic and numerical computations were performed by using Maple 13 software package.

First of all, error functions are introduced to study the accuracy and efficiency of the method. Actually, continuous approximations to the solution will be obtained. To show the accuracy of the present method for our problems, we report four types of error. The first one is the residual error, Res, defined as Res^𝑘(𝑡)𝑑∶=^𝑛𝑑𝑡^𝑛𝑥^𝑘(𝑡)−𝑓𝑡,𝑥^𝑘(𝑡),𝑑𝑥𝑑𝑡^𝑘(𝑡),𝑑²𝑑𝑡²𝑥^𝑘(𝑡)𝑑,…,^𝑛−1𝑑𝑡^𝑛−1𝑥^𝑘(𝑡), (24) while the exact, Ext, relative, Rel, and consecutive, Con, errors are defined, respectively, by Ext^𝑘(𝑡)∶=𝑥_exact(𝑡)−𝑥^𝑘(𝑡),Rel^𝑘(𝑡)𝑥∶=^exact(𝑡)−𝑥^𝑘(𝑡)𝑥_exact(𝑡,)Con^𝑘(𝑡)∶=𝑥^𝑘+1(𝑡)−𝑥^𝑘(𝑡), (25) for 𝑖=1,2,…,𝑛, where 𝑡∈[𝑡₀,𝑡₀+𝑎], 𝑥^𝑘 is the 𝑘th-order approximation of 𝑥(𝑡) obtained by RPS method.

The proposed method provides an analytical approximate solution in terms of an infinite power series. Then, there is a practical need to evaluate this solution and to obtain numerical values from the infinite power series. The consequent series truncation and the practical procedure are conducted to accomplish this task and transform the otherwise analytical results into an exact solution, which is evaluated to a finite degree of accuracy.

In most real-life situations, the differential equation that models the problem is too complicated to solve exactly, and there is a practical need to approximate the solution. In the next example, the exact solution cannot be found analytically.

Example 1. Consider the Genesio equation: 𝑥^′′′(𝑡)+𝑎𝑥^′′(𝑡)+𝑏𝑥^′(𝑡)+𝑐𝑥(𝑡)−(𝑥(𝑡))²=0,𝑡≥0, (26) subject to the initial conditions 𝑥(0)=𝐺₀,𝑥^′(0)=𝐺₁,𝑥^′′(0)=𝐺₂, (27) where 𝐺₀, 𝐺₁, and 𝐺₂ are positive real numbers, satisfying 𝐺₀𝐺₁<𝐺₂.

The Genesio equation was represented by authors in [21], as a third-order IVP that includes a simple square part and depends on three positive real parameters. The Genesio equation is one of paradigms of chaos since it captures many features of chaotic equation. The reader is kindly requested to go through [21–25] in order to know more details about applications of Genesio equation and method of solutions.

As we mentioned earlier, if we select the initial guess approximation as 𝑥_initial(𝑡)=𝑥²(𝑡)=𝐺₀+𝐺₁𝑡+(1/2)𝐺₂𝑡², then the 𝑘th-truncated series takes the form 𝑥^𝑘(𝑡)=𝑥_initial(𝑡)+^𝑘_𝑚=3𝑐_𝑚𝑡−𝑡₀^𝑚=𝐺₀+𝐺₁1𝑡+2𝐺₂𝑡²+𝑐₃𝑡³+𝑐₄𝑡⁴+⋯+𝑐_𝑘𝑡^𝑘. (28)

Here, it is easy to see that 𝑐₀=𝐺₀,  𝑐₁=𝐺₁, and  𝑐₂=(1/2)𝐺₂. Anyway, in order to find out the value of coefficients 𝑐_𝑚,𝑚=3,4,… of (28), we employ our algorithm RPS method. Therefore, we construct the residual function as follows: Res^𝑘(𝑡)=^𝑘_𝑚=3𝑚(𝑚−1)(𝑚−2)𝑐_𝑚𝑡−𝑡₀^𝑚−3+𝑎^𝑘_𝑚=2𝑚(𝑚−1)𝑐_𝑚𝑡−𝑡₀^𝑚−2+𝑏^𝑘_𝑚=1𝑚𝑐_𝑚𝑡−𝑡₀^𝑚−1𝑐^𝑘_𝑚=0𝑐_𝑚𝑡−𝑡₀^𝑚−^𝑘_𝑚=0𝑐_𝑚𝑡−𝑡₀^𝑚2. (29)

Now, to find the 3rd approximate solution, we put 𝑘=3 through (29) and use the fact that Res³(0)=0 to conclude 𝑐₃=16𝐺²₀−𝑐𝐺₀−𝑏𝐺₁−𝑎𝐺₂. (30)

Thus, we can write the 3rd approximation of the RPS solution for (26) and (27) as 𝑥³(𝑡)=𝐺₀+𝐺₁1𝑡+2𝐺₂𝑡²+16𝐺²₀−𝑐𝐺₀−𝑏𝐺₁−𝑎𝐺₂𝑡³. (31)

By continuing with the similar fashion, the 4th approximation of the RPS solution for (26) and (27) takes the form 𝑥⁴(𝑡)=𝐺₀+𝐺₁1𝑡+2𝐺₂𝑡²+16𝐺²₀−𝑐𝐺₀−𝑏𝐺₁−𝑎𝐺₂𝑡³+𝑐₄𝑡⁴. (32)

Again, to find out the value of a coefficient 𝑐₄ in (32), we put 𝑘=4 through (29) and use the fact that (𝑑/𝑑𝑡)Res⁴(0)=0, to obtain 𝑐₄1=−𝑎𝐺24²⁰−𝑐𝐺₀−𝑏𝐺₁−𝑎𝐺₂+𝑏𝐺₂+𝑐𝐺₁−2𝐺₀𝐺₁. (33)

Consequently, the 4th approximation of the RPS solution for (26) and (27) according to residual function (29) will be 𝑥⁴(𝑡)=𝐺₀+𝐺₁1𝑡+2𝐺₂𝑡²+16𝐺²₀−𝑐𝐺₀−𝑏𝐺₁−𝑎𝐺₂𝑡³−1𝑎𝐺24²⁰−𝑐𝐺₀−𝑏𝐺₁−𝑎𝐺₂+𝑏𝐺₂+𝑐𝐺₁−2𝐺₀𝐺₁𝑡⁴. (34)

For numerical results, the following values, for parameters, are considered [21]: 𝐺₀=0.2,  𝐺₁=−0.3,  𝐺₂=0.1,  𝑎=1.2,  𝑏=2.92, and 𝑐=6. By continuing with the similar fashion discussed above, the 10th-order RPS solution of 𝑥(𝑡) leads to the following truncated series: 𝑥¹⁰(𝑡)=0.2−0.3𝑡+0.05𝑡²−0.0673333333𝑡³+0.0780333333𝑡⁴−0.012064𝑡⁵−0.0022902222𝑡⁶−6.4525841270×10⁻⁴𝑡⁷+2.5788923810×10⁻⁴𝑡⁸+5.6070795062×10⁻⁵𝑡⁹−2.4439052416×10⁻⁵𝑡¹⁰. (35)

Numerical comparisons are studied next. Figure 1 shows a comparison between the numerical solution of (26) and (27) for 40th-order RPS solution together with Runge-Kutta method (RKM) of order four and Predictor-Corrector method (PCM) of order four. Throughout this figure, the step size for the RKM and PCM is fixed at 0.01, while the starting values of the PCM obtained from the 4th-order RKM. It is demonstrated that the RPS solution agrees very well with the solutions obtained by the RKM and PCM.

(a)

(b)

Figure 1 

Plots of RPS solution versus RKM and PCM solutions for (26) and (27) versus time: (a) (—): 40th terms RPS approximation; (- -): RKM solution; (b) (—): 40th term RPS approximation; (-·-): PCM solution.

Example 2. Consider the Lienard equation: 𝑥^′′(𝑡)+𝑙𝑥(𝑡)+𝑚𝑥³(𝑡)+𝑛𝑥⁵(𝑡)=0,𝑡≥0, (36) subject to the initial conditions 𝑥(0)=−2𝑙𝑚,𝑥^′(0)√=−−𝑙𝑙(),−2𝑙/𝑚𝑚 (37) where 𝑙,  𝑚, and 𝑛 are real coefficients.

The Lienard equation in its general form 𝑥^′′(𝑡)+𝑓(𝑥(𝑡))𝑥^′(𝑡)+𝑔(𝑥(𝑡))=ℎ(𝑡),  can be regarded not only as a generalization of the damped pendulum equation or a damped spring-mass system, but also used as nonlinear models in many physically significant fields. For example, the choices 𝑓(𝑥(𝑡))=𝜖(𝑥²(𝑡)−1),  𝑔(𝑥(𝑡))=𝑥(𝑡), and ℎ(𝑡)=0 lead to the Van der Pol equation served as a nonlinear model of electronic oscillation [26, 27]. In the general case, it is commonly believed that it is very difficult to find out the exact solution of the Lienard equation by usual ways [28]. However, authors in [29, 30] studied the special case (36) and (37) of the Lienard equation. The reader is asked to refer to [26–33] in order to know more details about the Lienard equation, including its history and kinds and method of solutions.

Our next goal is to approximate the solution of (36) and (37). According to RPS method, the initial guess approximation can be written as 𝑥_initial(𝑡)=𝑥¹(𝑡)=−2𝑙𝑚−√−𝑙𝑙−2𝑙/𝑚𝑚𝑡. (38)

Hence, as in the previous example, the first few truncated series approximations of the RPS solution for (36) and (37) are as follows: 𝑥²(𝑡)=−2𝑙𝑚−√−𝑙𝑙𝑙−2𝑙/𝑚𝑚𝑡−²4𝑡−2𝑙/𝑚𝑚²,𝑥³(𝑡)=−2𝑙𝑚−√−𝑙𝑙𝑡−𝑙−2𝑙/𝑚𝑚₂4𝑡−2𝑙/𝑚𝑚²−5√−𝑙𝑙²24𝑡−2𝑙/𝑚𝑚³, (39) and so on; for numerical results, the following values, for parameters, are considered [33]: 𝑙=−1, 𝑚=4, and 𝑛=−3. If we collect and extend the above results according to these parameters, then the 10th-order RPS solution of 𝑥(𝑡) can be truncated as follows: 𝑥¹⁰(𝑡)=√22+√24√𝑡−2𝑡16²−5√2𝑡96³+17√2𝑡768⁴+121√2𝑡7680⁵−721√2𝑡92160⁶−1369√2𝑡258048⁷+58337√2𝑡20643840⁸+698161√2𝑡371589120⁹−7734241√2𝑡7431782400¹⁰. (40)

In this manner, the components of the RPS solution are obtained as far as we like. In fact, this series is exact to the last term, as one can verify, of the Taylor series of the exact closed form solution as 𝑥(𝑡)=121tanh𝑡+2. (41)

Let us now carry out the error analysis of the RPS method for this example. Figure 2 shows the exact solution 𝑥_exact(𝑡), the initial solution 𝑥_initial(𝑡), and the four iterates of RPS approximation 𝑥^𝑘(𝑡) for 𝑘=2,3,4,10. These graphs exhibit the convergence of the approximate solutions to the exact solution with respect to the order of the approximation.

Figure 2 

Plots of four iterates of RPS approximation for (49) and (52) when 𝑘=2,3,4,10 together with the initial guess approximation and the exact solution; blue: 𝑥initial(𝑡); brown: 𝑥2(𝑡); green: 𝑥3(𝑡); red: 𝑥4(𝑡); purple: 𝑥10(𝑡); (-·-): 𝑥exact(𝑡).

In Figure 3(a), we plot the error functions Ext^𝑘(𝑡) for 𝑘=5,10,15,20 which are approaching the axis 𝑦=0 as the numbers of iterations increase. These graphs show that the exact errors are getting smaller as the order of the solution is increasing, in other words, as we progress through more iterations. On the other hand, Figure 3(b) shows the residual error functions Res^𝑘(𝑡) for 𝑘=5,10,15,20 for the obtained approximations. These error indicators confirm the convergence of the RPS method with respect to the order of the approximation.

(a)

(b)

Figure 3 

Plots of error functions for (36) and (37); blue: 𝑘=5; brown: 𝑘=10; green: 𝑘=15; red: 𝑘=20: (a) Ext𝑘(𝑡); (b) Res𝑘(𝑡).

Example 3. Consider Bratu’s equation: 𝑥^′′(𝑡)−2𝑒^𝑥(𝑡)=0,𝑡≥0, (42) subject to the initial conditions 𝑥(0)=0,𝑥^′(0)=0. (43)

The standard Bratu problem was used to model a combustion problem in a numerical slab. The Bratu model appears in a number of applications such as the fuel ignition of the thermal combustion theory, the model of the thermal reaction process, the Chandrasekhar model of the expansion of the universe, questions in geometry and relativity concerning the Chandrasekhar model, chemical reaction theory, and radiative heat transfer and nanotechnology [34–38]. The reader is kindly requested to go through [39–44] in order to know more details about applications of Bratu problem, its history and kinds, its method of solutions, and so forth.

As we mentioned earlier, if we select the initial guess approximation as 𝑥_initial(𝑡)=𝑥¹(𝑡)=0, then the first few terms approximations of the RPS solution for (42) and (43) are 𝑥₂(𝑡)=𝑡²,𝑥₃(𝑡)=0𝑡³,𝑥₄(𝑡)=16𝑡⁴,𝑥₅(𝑡)=0𝑡⁵,𝑥₆(𝑡)=2𝑡45⁶,…. (44)

If we calculate more terms and then collect the above results, then the 10th-truncated RPS solution of 𝑥(𝑡) for (42) and (43) is 𝑥¹⁰(𝑡)−1=−22𝑡²−1𝑡12⁴−1𝑡45⁶−17𝑡2520⁸−31𝑡14175¹⁰, (45) which is the first tenth term in the expansion of −2ln(cos𝑡). Thus, the exact solution of (42) and (43) has the general form which are coinciding with the exact solution 𝑥(𝑡)−1=−22𝑡²−1𝑡12⁴−1𝑡45⁶−17𝑡2520⁸−31𝑡14175¹⁰−⋯=−2ln(cos𝑡). (46)

We mention here that this RPS solution is the same as Adomian decomposition solution obtained in [41], homotopy perturbation solution obtained in [42], variational iteration solution obtained in [43], and perturbation-iteration solution obtained in [44].

Example 4. Consider the following second-order nonlinear equation: 𝑥^′′(𝑡)+𝑡𝑥^′(𝑡)²−𝑥(𝑡)𝑒^𝑥(𝑡)+𝑓(𝑡)=0,𝑡≥0, (47) subject to the initial conditions 𝑥(0)=181,𝑥^′(0)1=−9, (48) where 𝑓(𝑡)=𝑡(𝑡²−(2/3)𝑡+1/9)²−2𝑡+(2/3)−(1/81−(1/9)𝑡+(1/3)𝑡²−(1/3)𝑡³)𝑒^{1/81−(1/9)𝑡+(1/3)𝑡2−(1/3)𝑡3}.

As in the previous examples, if we select the initial guess approximation as 𝑥_initial(𝑡)=𝑥¹(𝑡)=1/81−(1/9)𝑡, then the RPS expansion of the solution takes the form 𝑥(𝑡)=1−1819𝑡+𝑐₂𝑡²+𝑐₃𝑡³+⋯. (49)

Consequently, the 3rd-order RPS solution for (47) and (48) according to this initial guess is as follows: 𝑥³(𝑡)=1−18191𝑡+3𝑡²−13𝑡³. (50)

It easy to discover that each of the coefficients 𝑐_𝑚 for 𝑚>4 in the expansions (49) is vanished. In other words, we have 𝑥³(𝑡)=³_𝑚=0𝑐_𝑚𝑡^𝑚=^∞_𝑚=0𝑐_𝑚𝑡^𝑚=𝑥(𝑡), (51) with full agreement with Corollary 2. Thus, the analytic approximate solution of (47) and (48) agrees well with the exact solution 𝑥(𝑡)=1/81−(1/9)𝑡+(1/3)𝑡²−(1/3)𝑡³.

A numerical experiment is carried out to verify the mathematical results, and the theoretical statement for the solution is supported by the results of numerical experiments. Next, we present more complicated IVPs in order to show the efficiency and applicability of the present method.

Example 5. Consider the following fourth-order nonlinear equation: 𝑥^′′′′(𝑡)−𝑥^′′′(𝑡)³𝑥^′′(𝑡)⁴−sin𝑥^′(𝑡)⁴cos𝑥^′′′(𝑡)⁵−𝑥²(𝑡)𝑒^{1−(𝑥′(𝑡))2}−𝑥(𝑡)+𝑓(𝑡)=0,𝑡≥0, (52) subject to the initial conditions 𝑥(0)=0,𝑥^′(0)=1,𝑥^′′(0)=0,𝑥^′′′(0)=−1, (53) where 𝑓(𝑡)=sin⁴𝑡cos³𝑡+sin²𝑡𝑒^sin2𝑡−sin(cos⁴𝑡)cos⁵(cos𝑡).

Now, we apply the RPS method to solve (52) and (53) numerically. Therefore, we construct the residual function (7) as follows: Res^𝑘(𝑡)=𝑥^𝑘_𝑚=4𝑚(𝑚−1)(𝑚−2)(𝑚−3)𝑐_𝑚𝑡−𝑡₀^𝑚−4−^𝑘_𝑚=3𝑚(𝑚−1)(𝑚−2)𝑐_𝑚𝑡−𝑡₀^𝑚−33×^𝑘_𝑚=2𝑚(𝑚−1)𝑐_𝑚𝑡−𝑡₀^𝑚−24−sin^𝑘_𝑚=1𝑚𝑐_𝑚𝑡−𝑡₀^𝑚−14×cos^𝑘_𝑚=3𝑚(𝑚−1)(𝑚−2)𝑐_𝑚𝑡−𝑡₀^𝑚−35−^𝑘_𝑚=0𝑐_𝑚𝑡−𝑡₀^𝑚2𝑒^{1−(∑𝑘𝑚=1𝑚𝑐𝑚(𝑡−𝑡0)𝑚−1)2}−^𝑘_𝑚=0𝑐_𝑚𝑡−𝑡₀^𝑚+𝑓(𝑡), (54) where ∑^𝑘_𝑚=0𝑐_𝑚(𝑡−𝑡₀)^𝑚 is 𝑘th-truncated series that approximate the solution 𝑥(𝑡).

Assuming that the initial guess approximation has the form 𝑥_initial(𝑡)=𝑥³(𝑡)=𝑡−(1/6)𝑡³, then, the 10th-truncated RPS solution of 𝑥(𝑡) for (52) and (53) is 𝑥¹⁰(𝑡)1=𝑡−6𝑡³+1𝑡120⁵−1𝑡5040⁷+1𝑡362880⁹=⁴_𝑗=0(−1)^𝑗𝑡^2𝑗+1!.2𝑗+1 (55)

In this manner, the components of the RPS solution are obtained as far as we like. Thus, the exact solution of (52) and (53) has the general form which are coinciding with the exact solution 𝑥(𝑡)=^∞_𝑗=0(−1)^𝑗𝑡^2𝑗+12𝑗+1!=sin𝑡. (56)

Example 6. Consider the following third-order nonlinear equation: 𝑥^′′′(𝑡)−𝑡𝑥^′′(𝑡)²𝑥^′(𝑡)−𝑒^{𝑥′(𝑡)}+𝑡cosh(𝑥(𝑡))²(𝑡)𝜋+𝑓=0,𝑡≥2, (57) subject to the initial conditions 𝑥𝜋2=12,𝑥^′𝜋2=√36,𝑥^′′𝜋21=−,18 (58) where 𝑓(𝑡)=𝑒^{((1/3)cos((1/3)𝑡)−(𝑡−(𝜋/2))2)}+1127cos3𝑡𝜋+2−𝑡𝑡−2²−131cos3𝑡×12𝑡−𝜋+91sin3𝑡²1−𝑡coshsin3𝑡−13𝜋𝑡−2³². (59)

As in the construction of RPS method, we assume that the initial guess approximation takes the form 𝑥_initial(𝑡)=𝑥²(𝑡)=1/2+(√3/6)(𝑡−(𝜋/2))−(1/36)(𝑡−(𝜋/2))². Then, the 10th-truncated RPS solution of 𝑥(𝑡) for (57) and (58) is as follows: 𝑥¹⁰(𝑡)=12+√36𝜋𝑡−2−1𝜋36𝑡−2²−13𝜋𝑡−2³−√3𝜋324𝑡−2³+1𝜋3888𝑡−2⁴+√3𝜋58320𝑡−2⁵−1𝜋1049760𝑡−2⁶−√3𝜋22044960𝑡−2⁷+1𝜋529079040𝑡−2⁸+√3𝜋14285134080𝑡−2⁹−1𝜋428554022400𝑡−2¹⁰. (60)

It easy to see that the 10th-truncated RPS solution of 𝑥(𝑡) above can be separated in the form of the following series that agree well with the general form: 𝑥¹⁰(𝑡)=_𝑛=0,4,812(3^𝑛)𝜋𝑛!𝑡−2^𝑛+_𝑛=2,6,10(−1)2(3^𝑛)𝜋𝑛!𝑡−2^𝑛+_𝑛=1,5,9√32(3^𝑛)𝜋𝑛!𝑡−2^𝑛+_𝑛=3,7,9√3(−1)2(3^𝑛)𝜋𝑛!𝑡−2^𝑛−13𝜋𝑡−2³. (61)

Thus, the exact solution of (57) and (58) has the general form which are coinciding with the exact solution 𝑥(𝑡)=_{𝑛=0,4,8,…}12(3^𝑛)𝜋𝑛!𝑡−2^𝑛+_{𝑛=2,6,10,…}(−1)2(3^𝑛)𝜋𝑛!𝑡−2^𝑛+_{𝑛=1,5,9,…}√32(3^𝑛)𝜋𝑛!𝑡−2^𝑛+_{𝑛=3,7,9,…}√3(−1)2(3^𝑛)𝜋𝑛!𝑡−2^𝑛−13𝜋𝑡−2³. (62)

Clearly, the RPS method provides the approximate solution in terms of an infinite power series. Now, we collect and transform these series in order to discover the exact solution. Anyway, it is not difficult to see that the Taylor series of the exact closed form solution of (57) and (58) is 𝑥(𝑡)1=sin3𝑡−13𝜋𝑡−2³. (63)

Our next goal is to show how the value of 𝑘 in the truncation series (6) affects the approximate solutions. To determine this effect, an error analysis is performed. We substitute the RPS approximations of 𝑥^𝑘(𝑡) for various 𝑘 into the exact error formula and obtain the exact error values. Next, the maximum and average errors at 𝑘=5,10,20 of (57) and (58) have been listed in Table 1 for 𝑡_𝑖=𝜋/2+(𝜋/20)𝑖, 𝑖=0,1,2,…,10.

Table 1 

The maximum and average exact error function of 𝑥^𝑘(𝑡) at 𝑘=5,10,15,20.

In Table 2, the exact, relative, and residual errors have been calculated for various 𝑡 in [𝜋/2,𝜋] with step size 𝜋/20 to measure the extent of agreement between the 10th-order RPS solution and the exact solution. From the table, it can be seen that the RPS method provides us with the accurate approximate solution for (57) and (58). Also, we can note that the RPS solution is more accurate at the beginning values of the independent interval [𝜋/2,𝜋].

Table 2 

The values of different error functions at 𝑘=10 for different values of 𝑡.

In Table 3, the value of consecutive error function for the two consecutive RPS solutions 𝑥^𝑘+1(𝑡) and 𝑥^𝑘(𝑡) has been calculated for various 𝑡 in [𝜋/2,𝜋] with step size 𝜋/20 to measure the difference between consecutive solutions obtained from the 8th-, 9th-, 10th-, 11th-, and 12th-order RPS solutions for (57) and (58). However, the computational results below provide a numerical estimate for the convergence of the RPS method. Also, it is clear that the accuracy obtained using present method is advanced by using only few terms approximations. In addition, we can conclude that higher accuracy can be achieved by evaluating more components of the RPS solution.

Table 3 

The values of consecutive error function at 𝑘=8,9,10,11 for different values of 𝑡.

5. Conclusion

The main concern of this work has been to propose an efficient algorithm for the solution of higher-order IVPs. The goal has been achieved by constructing a new method to solve this class of ordinary differential equations. We can conclude that the method is powerful and efficient technique in finding approximate solution for IVPs of different orders and types. The proposed algorithm produced a rapidly convergent series with easily computable components using symbolic computation software. There is an important point to make here, the results obtained by the method are very effective and convenient with less computational work and time. This confirms our belief that the efficiency of our technique gives it much wider applicability for general classes of linear and nonlinear problems.