Numerical Treatment of Initial Value Problems of Nonlinear Ordinary Differential Equations by Duan-Rach-Wazwaz Modified Adomian Decomposition Method ()
1. Introduction
Phenomena as diverse as the oscillations of a suspension bridge, the spread of a disease, and the motion of the planets are governed by nonlinear differential equations. Most of these nonlinear equations do not have analytical solutions, so approximation and numerical techniques must be used. The Adomian decomposition method (ADM), introduced by Adomian [1] [2] [3] [4] provides immediate and visible symbolic terms of analytic solutions as well as numerical approximate solutions to both linear and nonlinear problems without unphysical restrictive assumptions such as required by linearization, perturbation or discretization [1] [2] [3] [4] . It provides the solution in a rapidly convergent series with easily computable components if the equation has a unique solution. The technique uses a decomposition of the nonlinear operator as a series of Adomian functions. Each term of the series is a generalized polynomial, called the Adomian polynomial. The ADM has been successfully applied to a wide class of problems arising in applied sciences and engineering [1] - [14] over three decades.
Adomian decomposition method has led to a number of modifications made by various researchers for different purposes such as to improve the accuracy, or increase the speed of convergence, or expand the application of the original method. Adomian and Rach [11] introduced modified Adomian polynomials which converge slightly faster than the original polynomials and are convenient for computer generation. Wazwaz [12] [13] used Padé approximants to the solution obtained using a modified decomposition method and found that not only does this improve the results, but also that the error decreases with the increase of the degree of the Padé approximants. The later modifications of ADM were proposed by Wazwaz [14] , Wazwaz and El-Sayed [15] , Duan [16] [17] [18] [19] [20] , Duan and Rach [21] [22] [23] , Duan, Rach and Wazwaz [24] [25] .
In this paper, we consider the applications of the Duan-Rach-Wazwaz modification of ADM to the initial value problems (IVPs) for the systems of nonlinear ordinary differential equations (ODEs). In 2013, Duan, Rach and Wazwaz [25] presented a reliable modification of the ADM which bases on the previous modification schemes [14] - [24] , and computes the solutions of variable coefficients higher-order nonlinear initial value problems (IVPs) and solutions of systems of coupled nonlinear IVPs. To implement these algorithms they also designed multistage decomposition and numeric algorithms, and presented MATHEMATICA routines PSSOL and NSOL.
The text is organized as follows. The basic principles of ADM are given in Section 2. For the numerical solutions of the IVPs for the systems of nonlinear differential ODEs, the frameworks of the Duan-Rach-Wazwaz modification are presented in Section 3. In Section 4, numerical treatments of the nonlinear IVPs using the modified technique and MATHEMATICA numerical solution are performed. The solutions of some problems are also computed by using fourth-order Runge Kutta method (RK4) and the comparisons of the results are presented. A brief conclusion is given in Section 5. All computations are carried out in MATHEMATICA.
2. Basic Principles of the Adomian Decomposition Method
Consider the general nonlinear ODE in the Adomian’s operator-theoretic form
(1)
where g is a given analytic function and u is the unknown solution, and L is the linear operator to be inverted, R is the linear remainder operator, and N is an analytic nonlinear operator. We remark that the choice of the linear operator is designed to yield an easily invertible operator with resulting trivial integrations.
This means that the choice is not unique. Generally we choose
for
n-th order ODEs, then its inverse
follows as the n-fold definite integration operator from
to t. Hence, we have
, where
is determined using the initial conditions.
Application of
to each side of Equation (1) yields
(2)
where
.
The ADM decomposes the solution into a series
(3)
and then decomposes the nonlinear term into a series of Adomian polynomials
(4)
where
are called the Adomian polynomials and generated by the definitional formula
(5)
where
is a grouping parameter of convenience. The formulas of the first four Adomian polynomials for the one-variable simple analytic nonlinearity
are listed as
In the Duan-Rach-Wazwaz modification, by using Duan’s Corollary 3 algorithm [18] the one variable Adomian polynomials are written as
(6)
where the coefficients
are defined recursively [18] as
(7)
The formulae in (7) does not involve the differentiation operator for the coefficients
[20] [21] [22] , but requires only addition and multiplication. So, it is more convenient for computer algebra systems.
The definitional formula of the Adomian polynomials for decomposing multivariable nonlinear functions occurring in either single nonlinear nth-order ODEs
or in systems of coupled nonlinear ODEs with multivariable nonlinearities are published by Adomian and Rach in [6] . By assuming f is an m-ary analytic function
, where the
, for
are the unknown functions to be determined, the solutions
,
and the nonlinear function
are decomposed as
(8)
where the multivariable Adomian polynomials
depend on the
solution components
and are defined by [6]
(9)
where
is a grouping parameter of convenience. The first m-variable Adomian polynomial
is
, where
.
Substitution of the Adomian decomposition series for the solution
and the series of Adomian polynomials for the nonlinearity Nu, Equations (3) and (4) respectively, into Equation (2) yields
(10)
The solution components
may be determined by one of the several advantageous recursion schemes, which differ from another by the choice of the initial solution component
, beginning with the classic Adomian recursion scheme
(11)
where Adomian has chosen
as the initial solution. All the solution components
of the solution
can be determined using Equation (11) and hence, the solution series follows immediately [25] . We remark that the convergence of the Adomian series has already been proved by several investigators [26] [27] .
3. Duan-Rach-Wazwaz Modification of the Adomian Decomposition Method
We illustrate the general frameworks of the Duan-Rach-Wazwaz modification of the Adomian decomposition method [25] for solving the first-order differential equations and the systems of coupled nonlinear differential equations numerically. Throughout the section we assume the equations are in canonical forms.
3.1. IVP of the First-Order Nonlinear ODE
We consider the following first-order nonlinear, nonhomogeneous differential equation subject to a bounded initial condition
(12)
(13)
We assume that the nonhomogeneous term g and system coefficients
and
are analytic functions.
In Adomian’s operator-theoretic form, the Equation (12) can be written as
(14)
where L is the linear operator, R is the linear remainder operator and N is the nonlinear operator such that
We note that, accordance by ( [27] , page 105)
and g are continuously differentiable at any order and they can be expandable in entire series in the neighbourhood of
, therefore, sum function
is real analytic and indicated series is absolutely, uniformly convergent. Therefore, the presentation formula for function
is meaningful.
From the assumption on the analyticity of the functions
and
we write the respective Taylor expansion series as
Application of the Adomian decomposition series and Adomian polynomials series result
(15)
where
is the simple nonlinearity term and can be any analytic function in u and the corresponding one-variable Adomian polynomials
have the standard formula [6] [16] [17] [18]
.
By calculating the Cauchy products
and
, respectively we write
(16)
(17)
and hence
(18)
We, next solve Equation (14) for
and apply the one-fold definite integral operator
to each side of the resulting equation to get
(19)
since
.
Equation (19) is the equivalent nonlinear Volterra integral equation for the solution
.
Evaluating the integrals
and
, we write
(20)
(21)
By substituting Equations (15), (20) and (21) into Equation (19) we get the Adomian decomposition series as
Therefore, the modified recursion scheme is written as
(22)
for
, where the one-variable Adomian polynomials
are
(23)
As a result, the (m + 1)th-stage solution approximant is given by
, for
, in the limit, it yields the exact solution, that is,
.
By calculating the first several solution components using Equations (22) and (23), we derive the following sequence
By using induction, we find for
that
(24)
where the one-variable Adomian polynomials
depend solely on the solution coefficients
, for
, and are determined as
(25)
instead of the solution components
for
.
Therefore, we have derived the desired Taylor expansion series for the solution
as
. By inspection, from the Equation (25), the solution coefficients
are obtained as the nonlinear recurrence relation
(26)
where the one-variable Adomian polynomials
are the same as shown in Equation (26). So, the rule of recursion for the solution coefficients of the first order canonical nonhomogeneous nonlinear IVP with a variable input and variable system coefficients is obtained as
, for
.
3.2. IVP of the System of Coupled Nonlinear DEs
We consider the following n-th order system of m-coupled nk-th order nonhomogeneous nonlinear IVPs
(27)
(28)
where
and
for
,
,
are variable system coefficients and
are variable inputs, and
. We assume that the system coefficients and the system inputs are analytic functions. We also assume that the problem is subject to appropriate (
) bounded initial conditions, i.e.,
bounded initial conditions for each nkth-order nonlinear DE,
.
In Adomian’s operator-theoretic form Equation (27) can written as
(29)
where
are the linear operators,
are the linear remainder operators, i.e., generally sequential-order differential operators, and
are the nonlinear operators such that
For a particular nkth-order nonlinear DE in the system represented by Equation (27) or Equation (29), we choose the corresponding solution
as the primary solution and the solutions
, for
, as the secondary solutions with respect to this same nkth-order DE. We assume that
and
are analytic, and hence have the relating Taylor expansion series
The linear differential operators
are invertible, and their inverse operators
are given by the
-fold integral
for the case of a system of m-coupled nkth-order IVPs, where the initial conditions are all specified at the origin.
Application of the Adomian decomposition series and the series of the Adomian polynomials, yields
(30)
where the multi-order differential nonlinearity
can be any analytic function in
and the relating
-variable Adomian polynomials
have the standard formula [6] [16] [17] [18] .
or equivalently,
(31)
The relating Cauchy products are
where
where
,
hence it is also determined that
Next we solve Equation (29) for
as
(32)
Applying the nk-fold integral operator
to each side of Equation (32), we obtain
(33)
By integrating left side of Equation (33) and substituting the values specified in Equation (28) we obtain
(34)
Substituting this on the left side of Equation (33), we obtain
(35)
Formula (35) is the equivalent system of m-coupled nonlinear Volterra integral equations.
Evaluating the relating integrals, we get
(36)
(37)
where
(38)
where
,
(39)
Substitution of the Equations (30), (36), (37), (38) and (39) into Equation (35) yields the following system of m-coupled modified recursion schemes
for
.
Therefore, the (s + 1)th-stage solution approximants
are given by
for
.
From the calculation of the first several solution components, we deduce the following sequence
(40)
where the (
)-variable Adomian polynomials
are now constants and depend merely on the solution coefficients
for
. They are determined by induction as
(41)
or equivalently,
(42)
instead of the solution components
and solution derivative components
for
and
. Thus we have derived the desired Taylor expansion series for each of the m solutions
as
where the solution coefficients
are given by the system of m-coupled nonlinear recurrence relations, obtained from inspection of Equation (39), as
(43)
and
, where the (
)-variable Adomian polynomials
depend on the solution coefficients
, for
, as in Equation (41). Consequently, the rule of recursion for the solution coefficients of the canonical nth-order system of m-coupled nkth order nonhomogeneous nonlinear IVPs with variable inputs and variable system coefficients are given as
for
.
4. Examples
In this section, we consider several examples of IVPs for the systems of nonlinear ODEs, which have either quadratic or cubic nonlinearities but, exhibit rather complex behavior. The modified numeric solutions of the problems are obtained by using MATHEMATICA routines PSSOL and NSOL [25] . To compare the results, we have calculated the MATHEMATICA numeric solutions for the systems of differential equations by using the command “NDSolve”. We also compute numerical solutions using RK4 in examples 1, 2 and 3. Moreover, we use diagonal Padé approximants [28] [29] [30] to improve the modified results and compute errors in the approximations in the examples 2 and 3 since the exact solutions of these problems are known.
Example 1. Consider the Abel differential equation of the first kind in canonical form. It is a first order, nonhomogeneous differential equation with a cubic nonlinearity [31] .
(44)
over the interval
. This nonlinear IVP does not have an exact solution but, a detailed qualitative analysis can be found in [31] .
Running PSSOL by taking
to output 5th-degree or equivalently 6-term approximation to the solution as
We note that the order of approximation is
.
Running NSOL for
and step size
to output the numeric solution
of order 5 which is depicted with red line as the curve of 5th order approximation
and parametric plot on the left in Figure 1.
As the comparison, MATHEMATICA numeric solution and RK4 solution are found and the curves and the parametric plots of the results are sketched with blue and black lines in the middle and on the right, respectively in Figure 1.
From Figure 1, we conclude that this problem has a limit cycle.
Example 2. Consider the first-order nonhomogeneous nonlinear differential equation with a quadratic nonlinearity [32]
(45)
Figure 1. The curves and parametric plots of the 5th-degree MADM approximate solution (red), MATHEMATICA numeric solution (blue) and RK4 solution (black) using the step-size
over the interval
.
on the interval
. It has the exact solution
.
Running PSSOL routine for
and 8 to output the 13-term, 15-term and 17-term approximants of the solution, respectively, as
Indeed, these are the first 12, 14 and 16 terms of the Taylor series of the
function
about the point
, respectively. So, if it is
possible to compute all terms of the series we shall see that the Adomian series for this problem is simply that Taylor series. All terms of the series are positive so, absolute convergence is simply the convergence of the series.
Since
on the interval
, and
for any n, accordance by ( [26] , page 105) the
respective truncation errors are
,
and
.
We note that in these computations approximation orders are
and
, respectively.
The curves of the computed approximants and the exact solution are plotted in Figure 2(a).
The MATHEMATICA command
for
and 8 output the [6/6], [7/7] and [8/8] diagonal Padé-approximants of the 13-term, 15-term and 17-term approximants, generated by the routine PSSOL respectively. The curves of the Padé approximants and the exact solution are plotted in Figure 2(b).
Running NSOL routine for
, and 15 and step size
to generate numeric solutions on the interval
. In Figure 2(c), the curves of these NSOL numeric solutions and the exact solution are depicted.
As a comparison RK4 solution is computed and depicted with exact solution in Figure 2(d).
Figure 2. PSSOL outputs and exact solution in (a), Padé-approximants and exact solution in (b), outputs of the NSOL routine and exact solution in (c) and RK4 results and exact solution in (d) for n = 6, 7 and 8.
We denote
and consider the absolute error function
for
on the interval
and the maximal error parameter
for
, and 8.
In Figure 3, the curves of the absolute error function
for
, and 8 are given.
From Figure 3 we can conclude that the ADM can be combined with the diagonal Padé approximants to estimate the blow-up time [28] . Since
is the blow-up time for this problem, this can be seen from the figure.
The maximal error parameters
for
and 8 are given in Table 1.
Table 1 shows that the maximal errors for the exact solution decrease approximately at an exponential rate.
Example 3. Consider the 2-dimensional system of nonlinear differential equations with quadratic nonlinearity [32]
(46)
over the interval
. The exact solutions are
.
Running PSSOL for
to generate 11-term approximants
and
of the solutions
and
, respectively as
We note that the respective orders of the approximation are
and
Figure 3. The absolute error function for n = 6 (Red line), n = 7 (Blue line) and n = 8 (Purple line).
Table 1. The maximal error parameters MEn for n = 6, 7 and 8.
.
The MATHEMATICA command
and the MATHEMATICA command
for
output [5/5] Padé approximants for
and
as
The outputs of the routines PSSOL, NSOL and RK4 for
and 14 and correspondingly the outputs of the MATHEMATICA command
for
and 7 together with the exact solutions x and y are depicted in Figures 4(a)-(c).
In Figure 5, the absolute errors for the exact solutions
and
are sketched in (a) and (b), respectively.
We list the maximal error parameters
for
and 7 in Table 2.
Figure 4. The curves of PSSOL, NSOL, RK4 outputs and exact solution for n = 5, 6 and 7 in (a), (b) and (c), respectively.
Figure 5. (a) The absolute error for the x-component of the solution,
; (b) The absolute error for the y-component of the solution,
.
Table 2. The maximal error parameters MEn for n = 5, 6 and 7.
From Table 2 we can conclude that the maximal error parameters for both components of the exact solution decrease approximately at an exponential rate.
Example 4. Consider a three-dimensional system of autonomous nonlinear DE with quadratic nonlinearities [33]
(47)
where
, and the state variables, and
and
are constant parameters of the system. This system has many interesting complex behaviors and exhibits chaotic behavior over a wide range of parameters. It can show two coexisting one-wing, a single two-wing, three-and four-wing when its parameters are chosen appropriately.
When we set
and
the system can display two existing one-wing chaotic attractors with different initial conditions as shown in Figure 6. The right (colored in blue) and left (colored in red) attractors are constructed with initial values
and
, respectively by running NSOL for
and step-size
(left) and MATHEMATICA (right) on the interval
.
For the values of the parameters
and
, the system generates two existing two-wing chaotic attractors as shown in Figure 7. Each wing of the attractor has a helical form. Running NSOL by taking
and step-size
outputs the 11th-order numeric solutions for the system on the interval
. Two existing two-wing chaotic attractor generated by modified ADM is sketched on the left and the MATHEMATICA result is sketched on the right of Figure 7. For each attractor we use the above initial conditions.
With the choice of the parameters
and
the system generates a three-wing attractor for the initial values as
as seen in Figure 8. On the left 11th-order modified ADM numeric solution for
and step-size
, generated by routine NSOL and on the right MATHEMATICA numeric solution on the
Figure 6. Phase portraits of the two coexisting one-wing chaotic attractors for
blue line:
, red line:
generated by modified ADM (left) and MATHEMATICA (right).
Figure 7. Phase portraits of two existing two-wing chaotic attractors for
generated by modified ADM (left) and MATHEMATICA (right).
Figure 8. Phase portraits of a single three-wing chaotic attractor for
and
generated by modified ADM (left) and MATHEMATICA (right).
interval
are displayed.
It can also display a four-wing chaotic attractor with parameters
and
and initial values
as illustrated in Figure 9.
Another four-wing chaotic attractor exists for
and for
as shown in Figure 10.
Example 5. Consider 4-dimensional system of nonlinear DE with quadratic nonlinearities [34]
(48)
where a, b, c, and d are real parameters.
The system in (48) exhibits hypertoroidal behavior when the parameters are chosen as
and the initial values are taken as
as shown in Figure 11.
Running NSOL for
and step-size
outputs the 5th-order numeric solutions for the system on the interval
. The 3-dimensional x-y-z (top) and y-z-w (bottom) projections of the modified results
Figure 9. Phase portraits of a four-wing chaotic attractor for
. 11th-order modified ADM numeric solution (left) and MATHEMATICA numeric solution (right) on the interval
.
Figure 10. Phase portraits of a four-wing chaotic attractor for
. 11th-order modified ADM numeric solution (left) and MATHEMATICA numeric solution (right) on the interval
.
are depicted on the left of Figure 11. To compare the results, we also calculated MATHEMATICA numeric solutions of the system are also depicted on the right of Figure 11.
Example 6. Consider the following seven-dimensional third-order hyperchaotic system [35] with cubic nonlinearity in each equation
Figure 11. Hyperchaotic behavior of system (48) two different 3-dimensional x-y-z and y-z-w projections of the system are shown in A and B, respectively. The parameter values are
and the initial conditions are
. On the left NSOL results and on the right MATHEMATICA results are depicted.
(49)
Since it has a chaotic attractor when we set
,
,
,
,
,
,
,
and the initial conditions are taken as
,
,
,
,
,
,
we consider the numeric solution of the system (49) for those values over the interval
.
Running NSOL by taking
and step-size
generates 5th-order numeric solutions on the interval for the system. The 2- and 3-dimensional projections of the modified results are plotted on the left of Figure 12 and Figure 13. To compare the results, we also calculated MATHEMATICA numeric solution of the system and plotted on the right of Figure 12 and Figure 13.
Figure 12. Two-dimensional x1-x6, x2-x4 and x2-x6 projections of system (49) generated by NSOL (on the left) and MATHEMATICA (on the right).
Figure 13. 3-dimensional x1-x2-x6 and x2-x4-x6 projections of system (49) generated by NSOL (left) and MATHEMATICA (right).
5. Conclusion
In this study, we use the Duan-Rach-Wazwaz modified Adomian decomposition method for solving nonlinear IVPs of the first order nonlinear ODEs and two, three, four and seven dimensional systems of nonlinear ODEs. To show the computational accuracy of the technique we consider homogeneous and nonhomogeneous equations with variable and constant systems coefficients. In each example, the solution of the modified technique is compared with that from MATHEMATICA solution and with the exact solution if it is known. In addition, in examples 1, 2 and 3 the numeric solutions are also computed by RK4. We have seen that the modification results closely agree with MATHEMATICA and RK4 solutions and also with exact solutions, if available. However, we have obtained modified solutions over a bigger time step than MATHEMATICA and RK4 solutions. Moreover, different problems have been solved in order to confirm the robustness of the modification over a wide variety of ODEs. Therefore, it may be concluded that the method has the ability of applying all types of nonlinear ODEs provided uniqueness of the solutions.