Abstract
Differentialalgebraic equations (DAEs) are important tools to model complex problems in various application fields easily. Those DAEs with an index3, even the linear ones, are known to cause problems when solving them numerically. The present article proposes a new algorithm together with its multistage form to efficiently solve a class of nonlinear implicit Hessenberg index3 DAEs. This algorithm is based on the idea of applying the differential transform method (DTM) directly to the DAE without applying the traditional index reduction methods, which can be complex and often result in violations of the DAE constraints. Also, to deal with the nonlinear terms in the DAE, we approximate them using the Adomian polynomials. This new idea has given us a simple and efficient algorithm, which involves the solution of linear algebraic systems except for the initial recursion terms. This algorithm is easy to implement in Maple or Mathematica. Furthermore, to enlarge the interval of convergence of the power series solution obtained from the DTM, an algorithm for the multistage DTM is also given. Both algorithms are applied to solve two examples of highly nonlinear implicit index3 Hessenberg DAEs. Numerical results show that the DTM can determine the exact solution in convergent power series, while the multistage DTM can compute accurate numerical solutions over large intervals.
1. Introduction
Differentialalgebraic equations (DAEs) are systems formed by differential and algebraic equations. Simulations of models from various fields such as mechanical systems [1–3], vehicle industry [4, 5], biomechanics [6], machinery [7], aircraft landing gears [8], towed rockets [9], and unmanned aerial vehicles (UAVs) [10, 11] often lead to DAEs. DAEs are classified by a natural number called the index which has several definitions. In this article we use the differential index which is defined as the minimum number of differentiations required to transform a DAE to an ordinary differential equation [12]. Roughly speaking, the index of a DAE measures the difficulty this DAE can pose. If the index is greater or equal to two, we say that the DAE is a higher index DAE. The solution techniques for higher index DAEs, even the linear ones, are known to face numerical and analytical difficulties that are not encountered when solving ordinary differential equations [13–16]. As a consequence, many numerical methods were developed to solve higher index DAEs such as rangeKutta [13, 17], projected Taylor series methods [18], hybrid block algorithms [19], stabilization methods [20–22], augmented Lagrangian method [23], sequential regularization methods [24], and the differential transform method (DTM) [25, 26]. For the solution of DAEs, one can also find in the literature the power series method combined with the Adomian polynomials [27, 28] and other methods [29–33]. A very popular approach to treat higher index DAEs is first to reduce the index by differentiating the constraints one or more times with respect to time to obtain an ordinary differential system or an index1 DAE. Then apply numerical integration methods to the index reduced DAE. However, the main difficulty with this index reduction technique is that the computed numerical solution of the resulting DAE may no longer fulfill the constraints of the original DAE. Solving this index reduced system will result in constraints violation in general and therefore gives a solution with no physical meaning. To resolve this problem when solving higher index DAEs in Lagrangian formulation form, some methods such as stabilization or augmented Lagrangian formulations were proposed to control the constraint violations during the computation of the numerical solution [20–22]. The most used method of stabilization is that of Baumgarte [20], but its drawback is the way of choosing its feedback parameters. The augmented Lagrangian formulation proposed in [23] has the same problem of the parameter selection. The challenge is then to construct efficient numerical methods that provide accurate solutions for index3 DAEs while preserving all the constraints of the original DAE.
In this article, we propose a new efficient algorithm to solve the following class of nonlinear implicit Hessenberg index3 DAEs:where the dash denotes the time derivative, , , and . The vector function is such that . The vector function , () defines the apparent constraints of the DAE where denotes the Jacobian of .
The following initial conditions must be provided to solve DAE (1) uniquely:
It should be noted that the initial conditions and must be chosen consistent, that is to fulfill the following algebraic system:
Note here that we do not need to provide an initial condition for the unknown vector as the initial value is already fixed by the initial condition (2) and the hidden constraints in the DAE (1). Throughout this paper, we assume the functions and to be analytical. In addition, the Jacobian matrix is assumed to be nonsingular on . We assume further that the square matrix is nonsingular, where , which means that DAE (1) has index3. To see this, we differentiate (3) with respect to time to get equation (4). Then, we differentiate this equation once with respect to time to get
The set of equations formed by the second equation of (1) and (5) is uniquely solvable for the variables and since its Jacobianis nonsingular due to the assumption that the matrix is nonsingular. Therefore,where and are some functions. Now, differentiating (8) once with respect to time, we get an ordinary differential system
This can be written in the form
Hence, differentiation of three times has led to an ordinary differential system for the algebraic variable and thus DAE (1) has an index3.
In this manuscript, we present a new algorithm together with its multistage form to efficiently solve the class of implicit Hessenberg index3 DAEs (1). This algorithm is based on an effective combination of the DTM and the Adomian polynomials [28, 34–39]. The main idea in this paper is to apply the DTM directly to DAE (1) without using index reduction methods which are complex and computationally expensive. These methods also often result in constraints violations. Furthermore, to approximate the nonlinear terms and , we expand them in power series using the Adomian polynomials. This hybrid technique provides the exact solution in convergent power series. To extend the interval of convergence of the DTM solution, a multistage DTM (MsDTM) algorithm is also given. The MsDTM has a simple algorithm, which can be easily coded in Maple or Mathematica. To demonstrate the effectiveness and accuracy of the DTM algorithm and that of its multistage form MsDTM, two highly nonlinear implicit index3 Hessenberg DAEs are solved by coding these algorithms in Maple 15. The numerical results show that the DTM provides the solution in convergent power series while the MsDTM gives accurate approximate solutions over large intervals.
The organisation of this manuscript is as follows: in Section 2, we review the Adomian polynomials and the DTM to solve ordinary differential equations. Next, in Section 3, we present the new proposed algorithms to solve the implicit Hessenberg index3 DAE initialvalue problem (1) and (2). Then, in Section 4, two examples of highly nonlinear implicit index3 Hessenberg DAEs are solved to illustrate the accuracy of this new method. Finally, a conclusion is given in Section 5.
2. Adomian Polynomials and the Differential Transform Method
In this section, we give a brief review to the Adomian polynomials [34–39], which are useful in the expansion of the nonlinear terms and of the DAE (1), then we review the differential transform method. Usually, a nonlinear term in a differential equation is decomposed as follows:where the Adomian polynomials are computed, for all nonlinearities, from
Here, , are the terms used in the expansion
Making use of (12), we can compute the first Adomian polynomials:where
For a function with two variables, we have the following first few Adomian polynomials:where
The differential transform of a function is defined byand the inverse differential transform of is given by
An approximate solution is given bywhere is the number of terms in the approximation.
If is expanded as (17), then the nonlinear term can be expanded using the Adomian polynomials as [40]where now are the coefficients of expansion (19). One important property of the Adomian polynomials which we will make use of to develop our algorithms is the fact that the Adomian polynomials are affine with respect to .
3. The New Algorithm
The numerical solution of index3 DAEs, including the linear ones, is known to be hard to obtain. In this section, we present a new algorithm for solving the class (1) of nonlinear implicit Hessenberg index3 DAEs. The proposed algorithm is based on an effective combination of the differential transform method (DTM) with the Adomian polynomials [34–39]. The main idea of our technique is to apply the DTM directly to this class of DAEs, where the nonlinear terms are expanded in power series using the Adomian polynomials [40]. Then, by using the fact that the Adomian polynomial is affine with respect to , a linear algebraic recursion system for the differential transforms of the solution is obtained. Next, since the DAE (1) is an index3 DAE, then the derived algebraic system is shown to have a unique solution. The main advantage of our technique is that it does not need to reduce the index of the DAE (1) before applying the DTM. This has given us a simple and efficient algorithm that can be easily coded in Maple or Mathematica. The following new theorems are important for derivation of our technique.
Theorem 1. Consider the fully implicit ordinary differential system , with the initial condition where the Jacobian is nonsingular. Assume that the function is analytical. Then, the differential transform of is given by the recursion formula , with computed from . Here is the vector of th two variable Adomian polynomials of the components of the vector function .
Proof. Assume that the solution can be expanded as follows:Then, we expand the nonlinear term in terms of the Adomian polynomials as follows:This givesOr,This leads to the following recursion formula to compute the differential transform of the derivative of :Note here that the right hand side of the above formula does not depend on
In a similar manner, one can show the following theorem.
Theorem 2. Consider the algebraic system where we assume that the function is analytical. Then, the differential transform of satisfies the following recursion formula Here, is the vector of th Adomian polynomials of the components of the vector function .
To derive our algorithm, we start by writing the second equation of DAE (1) as follows:and
Then, we expand , and in terms of their differential transforms as follows:where the unknown vectors and , will be determined later by our algorithm. Then, we expand the nonlinear term of (27) in terms of the Adomian polynomials as follows:where denotes the vector of the Adomian polynomials of the components of . From the first equation of DAE (1) and equation (26), we obtain the following recursions for the differential transforms of the solution components , and :from which, we have
Now, since , then using Theorem 1, we obtainwhere , , and . Note here that depends only on the recursion terms and . System (33) represents a set of linear algebraic equations with variables. To solve it, we need to complete it by another linear algebraic equations. To do this, we approximate the nonlinear term of equation (1) using the Adomian polynomials as follows:where , denotes the vector of the Adomian polynomials of the components of the vector .
Next substituting the expansion of from (34) in the last equation of DAE (1), we deduce thatwhich leads, using Theorem 2, to the following algebraic linear system for the unknown vector :where and . It should be noted here that depends only on the recursion terms . Before solving the linear algebraic system formed by (33) and (36), we first calculate the initial recursion terms and from the initial conditions and . Furthermore, setting in (30), we obtain . Now, we need to determine the consistent initial value . Using the second equation of DAE (1) and equations (30) and (34) for , we obtain the algebraic system for the unknown vectors and :
This algebraic system is nonlinear when is nonlinear with respect to the unknown vectors and/or . Since the DAE (1) is index3, then the square matrix is nonsingular. Therefore, the Jacobian (38) of system (37) is nonsingular, and hence, system (37) is uniquely solvable for and :
Once the unknown is determined, we can then find the values of and by setting in equations (30) and (31), respectively.
Next, we need to determine the unknown vectors and for . Using (32), system (33) leads to
Finally, we can determine the unknown vectors and by solving the linear algebraic system formed by (36) and (39), for . This linear algebraic system is uniquely solvable since its coefficient matrixis nonsingular due to the fact that the matrix is nonsingular. To solve the system formed by the (36) and (39), we multiply system (39) from the left by the matrix and then use system (35) to obtain the following linear algebraic system for the unknown vector :
From system (41), we determine the unknown vector uniquely for . Substituting the expression of obtained from (41) into system (39), we can calculate the unknown vector , for as follows:
Once we have computed the unknown , we use system (28) to find the unknown vector , for Finally, we obtain an approximate solution for the implicit Hessenberg index3 DAE initialvalue problems (1) and (2) as follows:where is the order of approximation. This completes the solution process.
The above process of the DTM for solving the DAE initialvalue problems (1) and (2) is given by Algorithm 1.

We can extend the interval of convergence of the solution by using a multistage DTM (MsDTM). The procedure of the MsDTM is described in Algorithm 2.

4. Numerical Examples
In this section, two examples of highly nonlinear implicit DAEs (1) and (2) are solved to demonstrate the effectiveness and accuracy of the presented technique. The implementation of our algorithms was performed in Maple 15. For each example, we first apply the algorithm of the standard DTM given in Section 3 to compute the solution in power series. Then, we apply the multistage DTM algorithm given in Section 3 to compute a numerical solution over large intervals.
Example 1. In this first example, we consider the DAE (1) with the vector function given byand where the function defining the constraint is given bywhere , , and .
The Jacobian of is . We can easily check that the exact solution of the DAE of this example is , , and . This exact solution is obtained from the initial conditionsIt is easy to check that the above initial conditions are consistent, since they fulfill and . The DAE corresponding to this example has an index3, hence difficult to treat even numerically. To compute the solution in power series, we apply the algorithm of the standard DTM given in Section 3. This DTM algorithm was coded in Maple 15. By choosing the order of approximation , our DTM algorithm has computed the following power series approximations:One can easily check that the above approximations represent the first terms of the Maclauran power series of the exact solution. Note here that all our computations were performed exactly using Maple. To compute accurate numerical solutions over large internals and extend the interval of convergence of the solution, we apply the multistage DTM (MsDTM) algorithm given in Section 3. For this, we used an order of approximation and subdivided the solution interval into subintervals. Starting the recursion with and , we obtained the numerical results shown in Figures 1–6. Each of Figures 1–5 show the numerical solution (solid line), the exact solution (dashed line), and the corresponding absolute error. From the graphs (a), we can see that the numerical solution approximates well the exact solution. In each graph (b) of Figures 1–5, we can see the corresponding absolute error. These graphs show that the errors of approximation using the MsDTM are less than . The drift errors in the position and the velocity constraints and are found to be and , as shown in Figure 6. These numerical results show that DAE constraints are satisfied with a high accuracy and that the MsDTM computes accurate numerical solutions over large intervals while satisfying all the DAE constraints. To increase the accuracy in the approximate MsDTM solution, one can increase the approximation order or the number of subdivisions of the solution interval .
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
Example 2. In this example, we consider the DAE (1) with the vector function given byand the function defining the constraint is given bywhere , , and .
The Jacobian of is . We can easily check that the exact solution of the DAE for this Example 2 is , , and . This exact solution is obtained from the initial conditionsWe can easily check that the above initial conditions are consistent, since they satisfy and . The DAE corresponding to this example has an index3, hence difficult to treat even numerically. To compute the solution in power series, we apply the algorithm of the standard DTM given in Section 3. By choosing the order of approximation , our DTM algorithm has computed the following power series approximations:One can easily check that the above approximations represent the first terms of the Maclauran power series of the exact solution. Note here that all our computations were performed exactly using Maple. To compute accurate numerical solutions over large internals and extend the interval of convergence of the solution, we apply the multistage DTM (MsDTM) algorithm given in Section 3. For this, we used an order of approximation and subdivided the solution interval into subintervals. Starting the recursion with and , we obtained the numerical results shown in Figures 7–12. Each of Figures 7 to 11 show the numerical solution (solid line), the exact solution (dashed line), and the corresponding absolute error. From the graphs (a), we can see that the numerical solution approximates well the exact solution. In each graph (b) of Figures 7 to 11, we can see the corresponding absolute error. These graphs show that the errors of approximation using the MsDTM are less than . The drift errors in the position and the velocity constraints and are found to be and , as shown in Figure 12. These numerical results show that DAE constraints are satisfied with a high accuracy and that the MsDTM computes accurate numerical solutions over large intervals while satisfying all the DAE constraints. To increase the accuracy in the approximate MsDTM solution, one can increase the approximation order or the number of subdivisions of the solution interval .
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
5. Conclusion
Differentialalgebraic equations (DAEs) are systems of differential and algebraic equations. Index3 DAEs are known to be hard to treat even by using numerical methods. In this manuscript, we have proposed a new efficient algorithm that combines the differential transform method (DTM) with the Adomian polynomials to solve the nonlinear implicit Hessenberg index3 DAEs (1). The main advantage of this algorithm is that it is simple and does not require the use of the usual index reduction techniques, which are complex and often lead to DAEs constraints violations. The second advantage of this algorithm is that it can be easily coded in Maple or Mathematica. To illustrate the effectiveness and accuracy of our technique, two highly nonlinear implicit Hessenberg index3 DAEs are solved by the DTM and the MsDTM. Numerical results show that the DTM algorithm has computed the exact solutions in a power series forms. To extend the interval of convergence of the DTM solution, we have developed a multistage DTM (MsDTM) version of the DTM algorithm. This DTM multistage algorithm is useful in particular when simulating engineering problems over large intervals. The numerical results show that the MsDTM has successfully handled this class of nonlinear implicit index3 DAEs over large time intervals without violating the original DAE constraints. The MsDTM provides more accurate numerical solutions if more terms are included in the series or if more interval subdivisions are used. Our algorithms can be easily modified to solve other implicit nonlinear Hessenberg DAEs. Future work will be the design of efficient algorithms to solve more general higher index DAEs.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The author declares that there are no conflicts of interest.