A novel technique to solve nonlinear higher-index Hessenberg differential–algebraic equations by Adomian decomposition method

Since 1980, the Adomian decomposition method (ADM) has been extensively used as a simple powerful tool that applies directly to solve different kinds of nonlinear equations including functional, differential, integro-differential and algebraic equations. However, for differential–algebraic equations (DAEs) the ADM is applied only in four earlier works. There, the DAEs are first pre-processed by some transformations like index reductions before applying the ADM. The drawback of such transformations is that they can involve complex algorithms, can be computationally expensive and may lead to non-physical solutions. The purpose of this paper is to propose a novel technique that applies the ADM directly to solve a class of nonlinear higher-index Hessenberg DAEs systems efficiently. The main advantage of this technique is that; firstly it avoids complex transformations like index reductions and leads to a simple general algorithm. Secondly, it reduces the computational work by solving only linear algebraic systems with a constant coefficient matrix at each iteration, except for the first iteration where the algebraic system is nonlinear (if the DAE is nonlinear with respect to the algebraic variable). To demonstrate the effectiveness of the proposed technique, we apply it to a nonlinear index-three Hessenberg DAEs system with nonlinear algebraic constraints. This technique is straightforward and can be programmed in Maple or Mathematica to simulate real application problems.

naturally in many important application problems like constrained multibody systems (Benhammouda and Vazquez-Leal 2015;Simeon 1996;Simeon et al. 1994), vehicle system dynamics (Simeon et al. 1991), space shuttle simulation (Brenan 1983) and incompressible fluids dynamics. Unfortunately, higher-index DAEs are known to be difficult to solve. Therefore, they are usually transformed to ordinary differential systems (indexzero) or index-one DAEs before solving them. This transformation, called index-reduction, can be computationally very expensive and may also change the properties of the solution. Therefore, new techniques are required to solve higher-index DAEs efficiently.
However, for the application of the ADM to DAEs, one finds only four pieces of work in the literature (Duan and Sun 2014;Çelika et al. 2006;Hosseini 2006a, b). In Duan and Sun (2014), an index-one DAE is transformed to a second order ordinary differential equation before applying the ADM to it. In Çelika et al. (2006), the ADM is applied to simple semi-explicit index-one DAEs, where the DAE is first transformed to a system of ordinary differential equations before applying the ADM. In Hosseini (2006a), the ADM is applied to linear higher-index Hessenberg DAEs after transforming them to indexone DAEs. In Hosseini (2006b), index-one and index-two DAEs with linear constraints are solved where these DAEs are pre-processed by a transformation that relies much on the special forms they have. Therefore, in all these previous works, the ADM is not applied directly to the DAEs but rather to the transformed equations. The drawback of such transformations is that they can involve complex algorithms, can be computationally expensive and may lead to non-physical solutions.
In this work, we present a new procedure for solving a class of nonlinear higher-index Hessenberg DAEs based on the ADM. The ADM is first applied directly to the DAE where the nonlinear terms are expanded using the Adomian polynomials (Duan 2010a, b, 2011, Rach 2008Wazwaz 2000). Based on the index condition, a nonsingular algebraic recursion system is derived for the expansion components of the solution. Also, it is important to note that unlike previous works (Duan and Sun 2014;Çelika et al. 2006;Hosseini 2006a, b), our procedure does not make transformations to the DAEs before applying the ADM to them. To demonstrate the effectiveness of the proposed technique, we solve a nonlinear index-three Hessenberg DAEs system with nonlinear algebraic constraints. Further, our technique is based on a simple algorithm that can be programmed in Maple or Mathematica to simulate real application problems.
This paper is organized as follows: in "Review of the Adomian decomposition method" section, we review the ADM for solving ordinary differential equations. Next, in "The proposed method" section, we present our method for the solution of nonlinear higherindex Hessenberg DAEs systems. Then in "Application to a nonlinear index-3 DAEs system" section, we apply the developed technique to solve a nonlinear index-three Hessenberg DAEs system with nonlinear algebraic constraints. Finally, a discussion and a conclusion are given in "Discussion" and "Conclusion" sections, respectively.

Review of the Adomian decomposition method
In this section, we give a brief review for the Adomian decomposition method (ADM) (Adomian and Rach 1985;Almazmumy et al. 2012;Fatoorehchi et al. 2015;Hendi et al. 2012; Pue-on and Viryapong 2012; Ramana and Raghu Prasad 2014;Wazwaz 2001) to solve ordinary differential equations. For this purpose, let us consider the following nonlinear differential equation where L is an easily invertible operator (usually taken as the highest-order derivative), R is an operator grouping the remaining lower-order derivatives, N (u) is the nonlinear term and f is a given analytical function.
Solving Eq. (1) for Lu then applying the inverse operator L −1 to both sides, we obtain If Lu =u = du/dt and the initial condition u(t 0 ) = c 0 is given, then L −1 represents the integral from t 0 to t and L −1 Lu = u − c 0 . If Lu =ü = u (2) = d 2 u/dt 2 and the initial conditions u(t 0 ) = c 0 and u(t 0 ) = c 1 are given, then L −1 is the double fold integral from t 0 to t and L −1 Lu = u − c 0 − c 1 (t − t 0 ). In this case from (2), we have To apply the ADM to Eq. (3), we first assume that the solution u of (1) to have the infinite series form where the unknown solution components u n , n = 0, 1, 2, . . . are to be determined later by the method.
Second, the nonlinear term N(u) is expanded in an infinite series in terms of the Adomian polynomials N n (Duan 2010a(Duan , b, 2011Rach 2008Rach , 1984Wazwaz 2000) as Substituting (4) and (5) into (3) and choosing u 0 as we obtain Comparing the general terms on the left hand side and the right hand side of (7), we derive the following recursion scheme for the ADM Since u 0 is known from (6), recursion (8) can be used to generate as many solution components u n as one wants. Further, if series (4) converges then it gives the exact solution of (1) and an approximation of order K to solution can be obtained from To compute the Adomian polynomials N n , n = 0, 1, . . . associated with the nonlinearity N (u), one can use the following definition for all forms of nonlinearity Using this formula, we obtain the following first few Adomian polynomials where the dash ( ′ ) represents the differentiation with respect to u 0 .
In a similar manner, one can easily generate the remaining polynomials from ( 10) and expand (11). In the literature, there are several algorithms for computing the Adomian polynomials without the need for formula (10), but a more convenient algorithm for the m-variable case is recently proposed in Duan (2011) The proposed method In this section, we present our method for solving a class of nonlinear higher-index Hessenberg differential-algebraic equations (DAEs). This technique is based on the Adomian decomposition method (ADM). To solve this class of DAEs, we first apply the ADM directly to it and expand the nonlinear terms using the Adomian polynomials. Then, an algebraic recursion system for the solution expansion components is derived. Taking account of the index of the DAE, this system is shown to be uniquely solvable for the solution expansion components. Also, it is important to note that unlike previous works (Çelika et al. 2006;Sun 2014, 2006a, b), our technique does not make transformations to DAEs before applying the ADM to them. The class of nonlinear higher-index Hessenberg DAEs we consider here is This DAE is supplied with some consistent initial conditions where c n , n = 0, . . . , m − 1 are given constants. Note here that no initial conditions are prescribed to the variable v. We also assume that DAE initial-value problem (13)-(14) has a unique analytical solution.
The vectors u and v are called the differential and the algebraic variables respectively. System (13) Systems (13) arise frequently in many important applications like Navier-Stokes equations in incompressible fluids dynamics or Euler-Lagrange equations in constrained multibody systems. In what follows, we assume that system (13) (15) holds.
To solve DAE initial-value problem (13)-(14) by the ADM, we apply the operator L −1 (which is in the present case the m-fold integral from t 0 to t) to both sides of the first equation of (13) to get First, we assume the components u and v of the solution of (13) to have the infinite series form where the unknowns u n and v n , n = 0, 1, 2, . . . will be determined later by our method.
Second, we expand the nonlinear terms M(u, v) and N(u) in infinite series using the Adomian polynomials M n and N n as where M n := M n (u 0 , v 0 , . . . , u n , v n ) and N n := N n (u 0 , . . . , u n ).  (17) and (18) into equation (16), we get Choosing the first m terms u n as then comparing the general term on the left hand side with that on the right hand side of (19), we derive the following recursion system for the unknowns u n and v n−m Using (12), the Adomian polynomials M n and N n , n = 0, 1, 2, . . . can be written as and One important property of the Adomian polynomials M n and N n we exploit here is their linearity with respect to v n and u n , n ≥ 1 respectively.
For the value n = m, system (21) can be written as The unknowns in this algebraic system are u m and v 0 . System (24) is linear if the given function M(u, v) is linear with respect to v and this system is nonlinear with respect to v 0 if M(u, v) is nonlinear with respect to v. Since index condition (15) holds, the Jacobian matrix of system (24) is nonsingular. Therefore, system (24) is uniquely solvable for the unknowns u m and v 0 .
To solve system (24) for u m and v 0 , we multiply its first equation by − ∂N 0 ∂u 0 and substitute − ∂N 0 ∂u 0 u m by its expression from the second equation of (24). Then, we obtain the following equation for the unknown v 0 (20) u n = (c n /n!)(t − t 0 ) n , n = 0, . . . , m − 1, Since index condition (15) holds, Eq. (26) can be solved uniquely for v 0 . Substituting this computed value of v 0 into the first equation of (24), we can determine u m . Now for the values n ≥ m + 1, we use the first equation of (21) and the second equation of (22) to get where From the second equations of (21) and (23), we deduce that where Set of Eqs. (27) and (29) is a linear algebraic system for the unknowns u n and w n−m . Since index condition (15) holds, the coefficient matrix (25) of this system is nonsingular. Therefore, system (27) and (29) is uniquely solvable for u n and w n−m , for n ≥ m + 1.
To solve set of Eqs. (27) and (29), we multiply Eq. (27) by − ∂N 0 ∂u 0 and substitute − ∂N 0 ∂u 0 u n by its expression from (29), to obtain the following algebraic system for the unknown w n−m Since index condition (15) holds, Eq. (31) can be solved uniquely for w n−m , and we have Then, substituting the expression of w n−m into (27), we determine the unknown u n (26)

Now differentiating Eq. (32) m times with respect to t, we determine the unknown v n−m
Using (33) and (34) for n ≥ m with (20), we can determine all solution components u n and v n . An approximate solution for DAE initial-value problem (13)-(14) can be obtained using (9) as where K is the order of approximation of u(t).
It is worth noting that Eq. (21) is linear with respect to v n−m for all values of n ≥ m, except for the case n = m where (21) is nonlinear with respect to v 0 if the given function M (u, v) is nonlinear with respect to v. This linearity property of Eq. (21) has a great positive impact on the simplicity of our method and its efficiency. Note here also that many important problems arising from applications like constrained mechanical systems and the semi-discrete form of Navier-Stokes equations lead to DAEs systems of the form of (13 ), where M(u, v) is linear with respect to v. For these problems, system ( 21) is linear for all values of n ≥ m.

Application to a nonlinear index-3 DAEs system
In this section, we illustrate and demonstrate the effectiveness of our technique to solve nonlinear higher-index Hessenberg DAEs systems, which are known to be difficult to solve even numerically. Following the procedure developed in the previous section, we first apply the ADM directly to the DAEs system and expand the nonlinear terms using the Adomian polynomials. Then, taking account of the index of the DAE, we derive a nonsingular algebraic recursion system for the expansion components of the solution. Finally, by solving this algebraic system we obtain the solution of the DAE.
As a test problem, we consider the following nonlinear index-three Hessenberg DAEs system which describes the constrained motion of a particle to a circular track Benhammouda (2015) DAEs system (36) is supplied with the following (consistent) initial conditions Note that no initial condition v(0) is prescribed to the variable v(t) because v(0) is predetermined by the DAE and initial conditions (37). System (36) is index-three since three time differentiations of the algebraic equation (third equation) of (36) will lead to an ordinary differential equation for v(t). As a consequence, this DAEs system is difficult (37) u 1 (0) = 1,u 1 (0) = 0, u 2 (0) = 0,u 2 (0) = 1.
to solve numerically. Note also that the third equation of (36) is nonlinear which make this DAEs system harder to solve.
To solve DAEs initial-value problem (36)-(37) by the procedure developed in the previous section, we let Lu =ü and have L −1 u = t 0 t 0 udtdt. We assume that the solution components u 1 , u 2 , and v of (36) to have the form where the unknowns u 1,n , u 2,n and v n , n = 0, 1, 2, . . . will be determined later by our technique.
Applying L −1 to both sides of the first two equations of (36), we obtain Substituting expressions (38) into (39), we get where A k,n , B k,n , k = 1, 2 and C n , n = 0, 1, 2, . . . are the Adomian polynomials associated to the nonlinear terms u 3 k , u k v and u 2 1 + u 2 2 − 1, respectively. Choosing the first two terms in the series expansions of u 1 and u 2 as Then comparing the general terms on the left and right hand sides of (40), we derive the following recursion system Since B k,n−2 = n−2 i=0 u k,n−2−i v i , k = 1, 2 and C n = n i=0 u 1,n−i u 1,i + u 2,n−i u 2,i , for n ≥ 2, system (40) gives u 1,0 = 1, u 1,1 = 0, u 2,0 = 0, u 2,1 = t.
be used. This example shows that the direct application of the ADM is a simple powerful technique to obtain the exact or approximate solutions of nonlinear higher-index Hessenberg DAEs. In the case we want to solve a DAE with an unknown solution, one way to measure the error for the approximate solution is to use the mean square residual (MSR) as in Benhammouda and Vazquez-Leal (2015) since the convergence of the method is still not shown.

Conclusion
This work presents the analytical solution of a class of nonlinear higher-index Hessenberg DAEs using the Adomian decomposition method (ADM). A procedure for solving this class of DAEs is presented. For this class, the technique was tested on a nonlinear higher-index Hessenberg DAEs system with nonlinear algebraic constraints. The results obtained show that the method can be applied to nonlinear higher-index Hessenberg DAEs efficiently to obtain the exact or an approximate solution. On the one hand, it is important to note that these types of DAEs are difficult to solve and on the other, the direct application of the ADM was able to solve this class of nonlinear higher-index Hessenberg DAEs. Also, it is important to note that unlike previous works (Duan and Sun 2014;Çelika et al. 2006;Hosseini 2006a, b), our procedure does not make transformations to DAEs before applying the ADM to them. Our technique is based on a straightforward procedure that can be programmed in Maple or Mathematica to simulate real application problems. Finally, further work is needed to show the convergence of the proposed method, apply it with its modifications (for example a multistage ADM) to solve nonlinear higher-index Hessenberg partial differential-algebraic equations and other nonlinear higher-index DAEs.