# The Differential Transform Method as an Effective Tool to Solve Implicit Hessenberg Index-3 Differential-Algebraic Equations

Authors: Brahim Benhammouda
Journal: Journal of Mathematics
 (2023)

---

## Abstract

Differential-algebraic equations (DAEs) are important tools to model complex problems in various application fields easily. Those DAEs with an index-3, 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 index-3 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 index-3 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.


---

## Body

## 1. Introduction

Differential-algebraic 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], bio-mechanics [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 range-Kutta [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 index-1 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 index-3 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 index-3 DAEs:(1)u′=v,fv′,v,u,λ=0,t∈0,T,gu=0,where the dash ′ denotes the time derivative, ut∈ℝnu, vt∈ℝnu, and λt∈ℝnλ. The vector function f is such that f:ℝnu×ℝnu×ℝnu×ℝnλ⟶ℝnu. The vector function g:ℝnu⟶ℝnλ, (nλ≤nu) defines the apparent constraints of the DAE where Gu=gu′u∈ℝnλ×nu denotes the Jacobian of g.The following initial conditions must be provided to solve DAE (1) uniquely:(2)u0=α,v0=β.It should be noted that the initial conditionsα and β must be chosen consistent, that is to fulfill the following algebraic system:(3)gu=0,(4)Guv=0.Note here that we do not need to provide an initial condition for the unknown vectorλt as the initial value λ0 is already fixed by the initial condition (2) and the hidden constraints in the DAE (1). Throughout this paper, we assume the functions fw,v,u,λ and gu to be analytical. In addition, the Jacobian matrix A≔fw′w,v,u,λ∈ℝnu×nu is assumed to be nonsingular on KerGu. We assume further that the square matrix GA−1B∈ℝnλ×nλ is nonsingular, where B≔fλ′w,v,u,λ∈ℝnu×nλ, which means that DAE (1) has index-3. 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(5)Guv′+Guuv,v=0.The set of equations formed by the second equation of (1) and (5) is uniquely solvable for the variables v′ and λ since its Jacobian(6)ABG0∈ℝnu+nλ×nu+nλ,is nonsingular due to the assumption that the matrix GA−1B is nonsingular. Therefore,(7)v′=pv,u,λ,(8)λ=qv,u,where p and q are some functions. Now, differentiating (8) once with respect to time, we get an ordinary differential system(9)λ′=qv′v,upv,u,λ+qu′v,uv.This can be written in the form(10)λ′=sv,u,λ.Hence, differentiation of three times has led to an ordinary differential system for the algebraic variableλ and thus DAE (1) has an index-3.In this manuscript, we present a new algorithm together with its multistage form to efficiently solve the class of implicit Hessenberg index-3 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 f and g, 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 index-3 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 Section2, 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 index-3 DAE initial-value problem (1) and (2). Then, in Section 4, two examples of highly nonlinear implicit index-3 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 f and g of the DAE (1), then we review the differential transform method. Usually, a nonlinear term hu in a differential equation is decomposed as follows:(11)hu=∑k=0∞hku0,u1,…,uk,where the Adomian polynomials hk are computed, for all nonlinearities, from(12)hk≔hku0,u1,…,uk=1k!dkdξkh∑i=0∞ξiuiξ=0,k=0,1,2,….Here,uk,k=0,1,2,…, are the terms used in the expansion(13)ut=∑k=0∞uk.Making use of (12), we can compute the first Adomian polynomials:(14)h0=hu0,h1=u1h1,h2=u2h1+u122h2,h3=u3h1+u1u2h2+u133!h3,h4=u4h1+u222+u1u3h2+u12u22h3+u144!h4,h5=u5h1+u2u3+u1u4h2+u1u222+u12u32h3+u13u23!h4+u155!h5,where hk≔dkhu0/du0k,k=0,1,2,…For a functionhu,v with two variables, we have the following first few Adomian polynomials:(15)h0=hu0,v0,h1=u1h1,0+v1h0,1,h2=u2h1,0+v2h0,1+u122h2,0+v122h0,2+u1v1h1,1,h3=u3h1,0+v3h0,1+u1u2h2,0+v1v2h0,2+u1v2+u2v1h1,1+u136h3,0+v136h0,3+u12v12h2,1+u1v122h1,2,h4=u4h1,0+v4h0,1+u1u3+u222h2,0+v1v3+v222h0,2+u1v3+u2v2+u3v1h1,1+u12u22h3,0+v12v22h0,3+u1u2v1+u12v22h2,1+u1v1v2+u2v122h1,2+u144!h4,0+u13v16h3,1+u12v124h2,2+u1v134h1,3+v144!h0,4,where hk,l≔∂k+lhu0,v0/∂u0k∂v0l,k=0,1,2,…,l=0,1,2,…The differential transform of a functionut is defined by(16)uk=1k!dkutdtkt=0,k=0,1,2,…and the inverse differential transform of uk,k=0,1,2,… is given by(17)ut=∑k=0∞uktk.From (16) and (17), we have(18)ut=∑k=0∞1k!dkutdtkt=0tk.An approximate solution is given by(19)ut=∑k=0Kuktk,where K is the number of terms in the approximation.Ifut is expanded as (17), then the nonlinear term hu can be expanded using the Adomian polynomials as [40](20)h∑k=0∞uktk=∑k=0∞hku0,u1,…,uktk,where now uk,k=0,1,2,… 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 hk,k=1,2,3,… are affine with respect to uk.


## 3. The New Algorithm

The numerical solution of index-3 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 index-3 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 hk is affine with respect to uk, a linear algebraic recursion system for the differential transforms of the solution is obtained. Next, since the DAE (1) is an index-3 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 systemfu′,u=0, with the initial condition u0=u0 where the Jacobian fu′u′,u is nonsingular. Assume that the function f is analytical. Then, the differential transform of u′ is given by the recursion formula k+1fu′u1,u0uk+1=−fku1,…,kuk,0,u0,…,uk,k=1,2,3,…, with u1=u′0 computed from fu′0,u0=0. Here fk=fku1,…,kuk,k+1uk+1,u0,…,uk,k=1,2,3,… is the vector of k-th two variable Adomian polynomials of the components of the vector function f.Proof.
Assume that the solutionut can be expanded as follows:(21)ut=∑k=0∞uktk.
Then, we expand the nonlinear termfu′,u in terms of the Adomian polynomials as follows:(22)f∑k=0∞k+1uk+1tk,∑k=0∞uktk=∑k=0∞fktk=0.
This gives(23)fku1,…,kuk,k+1uk+1,u0,…,uk=0,k=1,2,3,…
Or,(24)k+1fu′u1,u0uk+1+fku1,…,kuk,0,u0,…,uk,k=1,2,3,…
This leads to the following recursion formula to compute the differential transform of the derivative ofu:(25)k+1fu′u1,u0uk+1=−fku1,…,kuk,0,u0,…,uk,k=1,2,3,…
Note here that the right hand side of the above formula does not depend onuk+1,k=1,2,3,…
In a similar manner, one can show the following theorem.Theorem 2.
Consider the algebraic systemgu=0 where we assume that the function g is analytical. Then, the differential transform of u satisfies the following recursion formula guu0uk=−gku0,…,uk−1,0,k=1,2,3,… Here, gk=gku0,…,uk−1,uk,k=1,2,3,… is the vector of k-th Adomian polynomials of the components of the vector function g.To derive our algorithm, we start by writing the second equation of DAE (1) as follows:(26)v′=w,and(27)fw,v,u,λ=0.Then, we expandwt,ut,vt, and λt in terms of their differential transforms as follows:(28)wt=∑k=0∞wktk,ut=∑k=0∞uktk,vt=∑k=0∞vktk,λt=∑k=0∞λktk,where the unknown vectors wk∈ℝnu,uk∈ℝnu,vk∈ℝnu and λk∈ℝnλ,k=0,1,2,…, will be determined later by our algorithm. Then, we expand the nonlinear term fw,v,u,λ of (27) in terms of the Adomian polynomials as follows:(29)f∑k=0∞wktk,∑k=0∞vktk,∑k=0∞uktk,∑k=0∞λktk=∑k=0∞fktk,where fk=fkw0,…,wk,v0,…,vk,u0,…,uk,λ0,…,λk,k=0,1,2,… denotes the vector of the Adomian polynomials of the components of f. From the first equation of DAE (1) and equation (26), we obtain the following recursions for the differential transforms of the solution components ut,vt, and wt:(30)vk−1=kuk,k=1,2,3,…,(31)wk−1=kvk,k=1,2,3,…,from which, we have(32)kk−1uk=wk−2,k=2,3,4,…Now, sincefk=0, then using Theorem 1, we obtain(33)A0wk+B0λk=−rk,k=1,2,3,…,where rk=fkw0,…,wk−1,0,v0,…,vk,u0,…,uk,λ0,…,λk−1,0, A0=fw′w0,v0,u0,λ0∈ℝnu×nu, and B0=fλ′w0,v0,u0,λ0∈ℝnu×nλ. Note here that rk depends only on the recursion terms w0,…,wk−1,v0,…,vk,u0,…,uk and λ0,…,λk−1. System (33) represents a set of nu linear algebraic equations with nu+nλ variables. To solve it, we need to complete it by another nλ linear algebraic equations. To do this, we approximate the nonlinear term gu of equation (1) using the Adomian polynomials as follows:(34)g∑k=0∞ukt=∑k=0∞gktk,where gk≔gku0,…,uk,k=0,1,2,…, denotes the vector of the Adomian polynomials of the components of the vector g.Next substituting the expansion ofg from (34) in the last equation of DAE (1), we deduce that(35)gk=0,k=0,1,2,…,which leads, using Theorem 2, to the following algebraic linear system for the unknown vector uk:(36)G0uk=−sk,k=1,2,3,…,where sk=gku0,…,uk−1,0 and G0=gu′u0∈ℝnλ×nu. It should be noted here that sk depends only on the recursion terms u0,…,uk−1. Before solving the linear algebraic system formed by (33) and (36), we first calculate the initial recursion terms u0 and v0 from the initial conditions u0=u0 and v0=v0. Furthermore, setting k=1 in (30), we obtain u1=v0. Now, we need to determine the consistent initial value λ0=λ0. Using the second equation of DAE (1) and equations (30) and (34) for k=2, we obtain the nu+nλ×nu+nλ algebraic system for the unknown vectors w0 and λ0:(37)fw0,v0,u0,λ0=0,G0w0=−2s2.This algebraic system is nonlinear whenf is nonlinear with respect to the unknown vectors w0 and/or λ0. Since the DAE (1) is index-3, then the square matrix G0A0−1B0∈ℝnλ×nλ is nonsingular. Therefore, the Jacobian (38) of system (37) is nonsingular, and hence, system (37) is uniquely solvable for w0 and λ0:(38)A0B0G00∈ℝnu+nλ×nu+nλ.Once the unknownw0 is determined, we can then find the values of u2=0.5w0 and v1=w0 by setting k=2 in equations (30) and (31), respectively.Next, we need to determine the unknown vectorsuk and λk−2 for k≥3. Using (32), system (33) leads to(39)kk−1A0uk+B0λk−2=−rk−2,k≥3.Finally, we can determine the unknown vectorsuk and λk−2 by solving the linear algebraic system formed by (36) and (39), for k≥3. This linear algebraic system is uniquely solvable since its coefficient matrix(40)kk−1A0B0G00∈ℝnu+nλ×nu+nλ,is nonsingular due to the fact that the matrix G0A0−1B0∈ℝnλ×nλ is nonsingular. To solve the system formed by the (36) and (39), we multiply system (39) from the left by the matrix G0A0−1 and then use system (35) to obtain the following linear algebraic system for the unknown vector λk−2:(41)G0A0−1B0λk−2=−G0A0−1rk−2+kk−1sk,k≥3.From system (41), we determine the unknown vector λk−2 uniquely for k≥3. Substituting the expression of λk−2 obtained from (41) into system (39), we can calculate the unknown vector uk, for k≥3 as follows:(42)kk−1A0uk=−rk−2+B0G0A0−1B0−1G0A0−1rk−2−kk−1sk,k≥3.Once we have computed the unknownuk, we use system (28) to find the unknown vector vk, for k=2,3,4,… Finally, we obtain an approximate solution for the implicit Hessenberg index-3 DAE initial-value problems (1) and (2) as follows:(43)ut=∑k=0Ktkuk,vt=∑k=0K−1tkvk,λt=∑k=0K−2tkλk,where K is the order of approximation. This completes the solution process.The above process of the DTM for solving the DAE initial-value problems (1) and (2) is given by Algorithm 1.Algorithm 1: DTM algorithm.
inputreadK;readnu, nλ, f∈ℝnu and g∈ℝnλ;readu0∈ℝnu and v0∈ℝnuoutputthe approximate solution:ut,vt,λt;while0≤k≤Kdocomputefk, gkendinitializationu0≔u0,v0≔v0,u1≔v0;ComputeA0=fw′w0,v0,u0,λ0,B0=fλ′w0,v0,u0,λ0,G0=gu′u0;Solve the algebraic systemfw0,v0,u0,λ0=0 and G0w0=−2g2 for w0 and λ0;u2≔0.5w0, v1=2u0;H≔G0A0−1B0;while3≤k≤KdoSolveHλk−2=kk−1sk−G0A0−1rk−2 for λk−2;Solvekk−1A0uk≔−rk−2−B0λk−2 for uk;vk−1≔kuk;endut≔∑k=0Ktkuk, vt≔∑k=0K−1tkvk, and λt≔∑k=0K−2tkλk.We can extend the interval of convergence of the solution by using a multistage DTM (MsDTM). The procedure of the MsDTM is described in Algorithm2.Algorithm 2: MsDTM algorithm.
(1)
Choose the numberN of subdivisions and the order of approximation K.(2)
Subdivide the interval0,T into N equal subintervals ti−1,ti, i=1,…,N of size h=T/N where ti=ih.(3)
Leti=1, apply the DTM to the DAE (1) on 0,t1 using u01=u0 and v01=v0 to obtain the approximate solution u1t=∑k=0Kuk1tk, v1t=∑k=0K−1vk1tk, and λ1t=∑k=0K−2λk1tk, 0≤t≤t1.(4)
Fromi=2 until N, apply the DTM to DAE (1) on ti−1,ti using u0i=ui−1ti−1 and v0i=vi−1ti−1 to obtain the approximate solutions uit=∑k=0Kukit−ti−1k, vit=∑k=0K−1vkit−ti−1k, and λit=∑k=0K−2λkit−ti−1k, ti−1≤t≤ti.(5)
Repeat step 4.(6)
An approximate solution to the DAE initial-value problems (1) and (2) on 0,T is given byut=u1t,0≤t≤t1u2t,t1≤t≤t2⋮uNt,tN−1≤t≤T,vt=v1t,0≤t≤t1v2t,t1≤t≤t2⋮vNt,tN−1≤t≤T and λt=λ1t,0≤t≤t1λ2t,t1≤t≤t2⋮λNt,tN−1≤t≤T.


## 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 f given by(44)fv′,v,u,λ=1+4cosu11+cosu11+cosu11−cosu1v1′3+v1′+u13+u1v2′3+v2′+u23+u2−v13−v2−λ3−2λcosu1+cosu1+u2v13−v2−λ3−2λcosu1+u2,and where the function g defining the constraint is given by(45)gu=sinu1+sinu1+u2,where u=u1,u2T∈ℝ2, v=v1,v2T∈ℝ2, and λ∈ℝ.
The Jacobian ofg is Gu=gu′u=cosu1+cosu1+u2,cosu1+u2. We can easily check that the exact solution of the DAE of this example is ut=sint,−2sintT, vt=cost,−2costT, and λt=cost. This exact solution is obtained from the initial conditions(46)u0=0,0T,v0=1,−2T.
It is easy to check that the above initial conditions are consistent, since they fulfillgu0=0 and gu′u0v0=0. The DAE corresponding to this example has an index-3, 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 K=10, our DTM algorithm has computed the following power series approximations:(47)ut=t1−2−t316;−13+t51120−160−t715040−12520+t91362880−1181440,vt=1−2−t212−1+t4124−112−t61720−1360+t8140320−120160,λt=1−12t2+124t4−1720t6+140320t8.
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 approximationK=10 and subdivided the solution interval 0,T=0,10 into N=100 subintervals. Starting the recursion with u0=u0=0,0T and v0=v0=1,−2T, 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 4×10−12. The drift errors in the position and the velocity constraints gu=0 and gu′uv=0 are found to be 10−14 and 4×10−15, 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 K or the number N of subdivisions of the solution interval 0,T.Figure 1
The approximate solution componentu1 (a) and the absolute error eu1 (b).
(a)(b)Figure 2
The approximate solution componentu2 (a) and the absolute error eu2 (b).
(a)(b)Figure 3
The approximate solution componentv1 (a) and the absolute error ev1 (b).
(a)(b)Figure 4
The approximate solution componentv2 (a) and the absolute error ev2 (b).
(a)(b)Figure 5
The approximate solution componentλ (a) and the absolute error eλ (b).
(a)(b)Figure 6
The error in the position constraintgu=0 (a) and the error in the velocity constraint gu′uv=0 (b).
(a)(b)Example 2.
In this example, we consider the DAE (1) with the vector function f given by(48)fv′,v,u,λ=tanv1′+2u2−u1v12+v22+8u1λ+2v1′+2u2−u1v12+v22+8u1λtanv2′−2u1−u2v12+v22+8u2λ+3v2′−2u1−u2v12+v22+8u2λ,and the function g defining the constraint is given by(49)gu=u12+u22−1,where u=u1,u2T∈ℝ2, v=v1,v2T∈ℝ2, and λ∈ℝ.
The Jacobian ofg is Gu=gu′u=2u1,2u2. We can easily check that the exact solution of the DAE for this Example 2 is ut=cost2,sint2T, vt=−2tsint2,2tcost2T, and λt=t2. This exact solution is obtained from the initial conditions(50)u0=1,0T,v0=0,0T.
We can easily check that the above initial conditions are consistent, since they satisfygu0=0 and gu′u0v0=0. The DAE corresponding to this example has an index-3, 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 K=10, our DTM algorithm has computed the following power series approximations:(51)ut=10+t201−t4120−t6016+t81240+t1001120,vt=t02−t320−t501+t7130+t90112,λt=t2.
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 approximationK=12 and subdivided the solution interval 0,T=0,5 into N=300 subintervals. Starting the recursion with u0=u0=1,0T and v0=v0=0,0T, 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 10−11. The drift errors in the position and the velocity constraints gu=0 and gu′uv=0 are found to be 6×10−13 and 3×10−12, 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 K or the number N of subdivisions of the solution interval 0,T.Figure 7
The approximate solution componentu1 (a) and the absolute error eu1 (b).
(a)(b)Figure 8
The approximate solution componentu2 (a) and the absolute error eu2 (b).
(a)(b)Figure 9
The approximate solution componentv1 (a) and the absolute error ev1 (b).
(a)(b)Figure 10
The approximate solution componentv2 (a) and the absolute error ev2 (b).
(a)(b)Figure 11
The approximate solution componentλ (a) and the absolute error eλ (b).
(a)(b)Figure 12
The error in the position constraintgu=0 (a) and the error in the velocity constraint gu′uv=0 (b).
(a)(b)


## 5. Conclusion

Differential-algebraic equations (DAEs) are systems of differential and algebraic equations. Index-3 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 index-3 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 index-3 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 index-3 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.



---

*Source: 3620870-2023-01-23.xml*