# A New Numerical Technique for Index-3 DAEs Arising from Constrained Multibody Mechanical Systems

Authors: Brahim Benhammouda
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

Index-3 differential-algebraic equations (DAEs) are mathematical models for the dynamics of constrained multibody mechanical systems that arise in many applications. These DAEs are known to pose a challenge to numerical methods. The purpose of this paper is to propose a new approach to efficiently solve these equations. This approach relies on an effective combination of the power series method (PSM) with the Adomian polynomials. Here, the PSM is directly applied to these DAEs without using the usual index reduction techniques, which are costly and often lead to nonphysical solutions. We expand the nonlinear terms in a series form using the Adomian polynomials to overcome the limitation of the PSM in collecting the coefficients of the power series solution. This technique has led to a simple and efficient algorithm. The domain of convergence of the power series solution is expanded by developing a multistage PSM (MSPSM). To demonstrate the efficiency of the MSPSM and show its applicability, an index-3 nonlinear DAE problem describing a two-link planar robot arm is solved. The numerical results show that the MSPSM is a powerful tool for solving the index-3 DAEs arising from constrained multibody mechanical systems.


---

## Body

## 1. Introduction

Constrained multibody systems [1] have applications in many important applications like rocket-towed systems [2], aerospace engineering [3], vehicle industry [4, 5], and robotics [6]. Constrained mechanical multibody systems are collections of rigid and/or flexible bodies that are interconnected by kinematic joints. The role of these kinematic joints is to impose some restrictions on the relative motion of the bodies of the system. As examples of constrained mechanical multibody systems, one can mention vehicles, cranes, slider cranks, and constrained robot arms. To express restrictions on the motion of the bodies, algebraic equations are used. Due to these equations, some constraint forces are present in the system. The forces that act on the multibody system components can come from springs, dampers, and actuator or external forces. The dynamics of these complex mechanical systems is described by nonlinear index-3 differential-algebraic equations (DAEs) systems of the form(1)u′=v,Muv′=fu,v,t−GTu,tλ,0=gu,t,where t∈0,T is time, the dash ′ represents the time derivative, and ut∈ℝnu and vt∈ℝnu specify the positions and the orientations of all bodies and their velocities, respectively. The vector λ∈ℝnλ is the vector of Lagrange multipliers. The matrix Mu∈ℝnu×nu is a mass matrix. The mapping f:ℝnu×ℝnu×0,T⟶ℝnu defines the applied and internal forces (other than the constraint forces), whereas g:ℝnu×0,T⟶ℝnλ, nλ≤nu defines the constraints. The term Gu,tTλ represents the constraint forces, where Gu,t=∂g/∂u∈Rnλ×nu denotes the Jacobian of gu,t. To solve DAE (1) uniquely, some initial conditions(2)u0=α,v0=β,are required. The given vectors α and β, which specify the initial configuration and initial velocity, must be chosen consistent, that is, to satisfy the algebraic system(3)gα,0=0,Gα,0β=0.Note here that no initial condition is required for the unknownλt as λ0 is already predetermined by DAE (1) and the initial conditions (2).Throughout this paper, we assume thatMu, fu,v,t, and gu,t are analytical. A standard assumption on the Jacobian Gu,t is the full row rank condition(4)rankGu,t=nλ,which means that the constraint equations defined by gu,t are linearly independent.In addition, the matrixMu is assumed to be symmetric and positive definite on KerGu,t; that is,(5)ξTMuξ>0,for allξ∈KerGu,t.If assumptions (4)-(5) hold, then the matrix(6)MuGTu,tGu,t0,is nonsingular and therefore the DAE (1) is an index-3 DAE. If initial conditions (2) are consistent, then the DAEs system (1) has a unique solution [7]. The index used in this paper is the differential index. It is defined as the minimum number of time differentiation needed to transform the DAE into an ordinary differential equation [8]. Those DAEs with an index greater than one, called also higher index DAEs, are widely known to be difficult to solve. The reason is that these DAEs are in fact ill-posed problems in a certain sense and direct discretizations do not lead to satisfactory results in general [8–11]. There are several methods in the literature for solving DAEs, which include Runge–Kutta methods [8, 12], projected explicit and implicit Taylor series methods [13], hybrid block methods [14], variational methods [15], reduction methods [16, 17], transformation methods [18], Adomian decomposition method [19], and the differential transform method [20]. One can find several techniques for solving DAEs of the form (1) [21–28]. The common way to treat these DAEs is to transform them into lower-index DAEs using successive time differentiation before using a numerical method. However, one major drawback of this technique is that the numerical solution of the transformed DAE may not satisfy the algebraic constraints of original DAE (1). This violation of constraints leads to nonphysical solutions. To control the error in the constraints, other techniques like stabilization or augmented Lagrangian formulation were constructed [22–24, 29]. The most popular one is the Baumgarte stabilization method [22]. However, the major drawback of this method is the way of selecting the feedback parameters. The great challenge is then to construct efficient approximation methods for the DAE (1), which satisfy all the constraints of the DAE.The power series method (PSM) [30, 31] is a well-known method for solving ordinary differential equations. Although this technique is effective for solving many problems like mechanical systems [27, 32], nonlinear Schrödinger equations [33], partial differential equations [34, 35], nonlinear systems of algebraic equations [36], and epidemiology [37], it was not frequently used for nonlinear problems. The reason is that the method may converge only over smaller intervals and the collection of the coefficients of the power series solution can be a difficult task. We can overcome these drawbacks using the Adomian polynomials and the multistage technique.In this paper, we present a new technique to solve general strongly nonlinear index-3 DAEs (1) arising from constrained mechanical systems efficiently. This technique is based on an effective combination of the power series method (PSM) with the Adomian polynomials [38], where there are several algorithms for computing these polynomials [39–43]. The novelty in our approach is in the direct application of the approximation method to these DAEs without the use of the usual index reduction techniques which are complex and costly and often lead to nonphysical solutions. The proposed technique has led to a general solution procedure with a simple and efficient algorithm that has three main features. The first one is that the algorithm includes all the DAE algebraic constraints. This ensures that the numerical solution remains on the DAE constraint manifolds. The second feature is that our algorithm requires only the solution of small algebraic systems for the algebraic variables with positive definite coefficient matrices. The third feature is that our technique can cope with any type of nonlinearity in the DAE using the Adomian polynomials. Finally, to enlarge the domain of convergence of the numerical solution, we developed a multistage PSM (MSPSM). This allows us to simulate constrained mechanical multibody systems over large time intervals. The proposed technique has led to a simple and efficient algorithm that can be easily implemented to solve general nonlinear index-3 DAEs. Our technique can be used in designing software for industrial applications and provides a powerful tool for mechanical engineers in particular. To demonstrate the effectiveness and accuracy of the multistage MSPSM, an index-3 DAE problem that models a two-link planar robot arm is solved by implementing the proposed technique in Maple. Numerical simulations show that the multistage MSPSM is a powerful tool to solve the DAE initial-value problem (1)-(2).This paper is organized as follows: in Section2, we review the Adomian polynomials. Next, in Section 3, we derive the proposed method for solving the DAE initial-value problem (1)-(2). Then, in Section 4, we apply the developed technique to solve a nonlinear index-3 DAE that models a two-link planar robot arm and give a discussion of numerical results. Finally, a conclusion is given in Section 5.


## 2. Adomian Polynomials

In this section, we briefly review the Adomian polynomials [38]. These polynomials will be used to expand the nonlinear terms in the differential equation. A nonlinear term Au in a differential equation is often decomposed as(7)Au=∑k=0∞Aku0,u1,…,uk,where the Adomian polynomials Ak are generated, for all types of nonlinearities, from(8)Ak≔Aku0,u1,…,uk=1k!dkdξkA∑i=0∞ξiuiξ=0,k≥0.Hereuk, k=0,1,2,…, denote the components used in the expansion(9)u=∑k=0∞uk.Using (8), one can easily compute the first few Adomian polynomials(10)A0=Au0,A1=u1A1,A2=u2A1+u122!A2,A3=u3A1+u1u2A2+u133!A3,A4=u4A1+u222!+u1u3A2+u12u22!A3+u144!A4,where Ak:=dkAu0/du0k.In a similar way, one can generate any desired Adomian polynomial from (8). For a two-variable function Au,v, the first few Adomian polynomials are(11)A1=u1A1,0+v1A0,1,A2=u2A1,0+v2A0,1+u122A2,0+v122A0,2+u1v1A1,1,A3=u3A1,0+v3A0,1+u1u2A2,0+v1v2A0,2+u1v2+u2v1A1,1+u136A3,0+v136A0,3+u12v12A2,1+u1v122A1,2,where Ak,l≔∂k+lAu0,v0/∂u0k∂v0l.Whenu is expanded in a power series u=∑k=0∞uktk, the nonlinear term Au can be expanded using the Adomian polynomials as follows [44]:(12)A∑k=0∞uktk=∑k=0∞Aku0,u1,…,uktk.One silent property of the Adomian polynomials which we will exploit when developing our algorithm is the fact thatAk is affine with respect to uk for all k≥1.


## 3. Derivation of the Proposed Technique

In this section, we construct our method for the solution of DAEs (1). These systems are index-3 DAEs which are known to be difficult to solve even numerically. The technique we propose here is an effective combination of the power series method (PSM) with the Adomian polynomials [38]. To make our algorithm robust and efficient, we apply the PSM directly to the DAE (1) and do not use index reduction techniques which are complex and computationally expensive. Then, we expand the nonlinear terms in power series using the Adomian polynomials [44]. To incorporate the algebraic constraints in the solution procedure, we use the fact that the Adomian polynomials Ak are affine polynomials with respect to uk for all k≥1. All this has simplified our algorithm and led to the solution of only small nλ×nλ nonsingular linear algebraic systems. This has also forced the numerical solution to satisfy all the algebraic constraints of the original DAE. Our algorithm is simple and can be easily implemented in Maple or Mathematica. The following new theorem gives some important results for the Adomian polynomials in the case of vectors and matrices.Theorem 1.
Letq=q1,…,qnT∈ℝn, h=h1,…,hnT∈ℝn, and u=u1,…,unT∈ℝn be three vectors and M:=Mui,j∈ℝn×ℝn be an invertible matrix. Let q=∑k=0∞qk and h=∑k=0∞hk with qk=q1,k,…,qn,kT and hk=h1,k,…,hn,kT, where qi,k and hi,k are the k th Adomian polynomials of the components qi and hi, i=1…,n, respectively. Let Mk=Mi,j,k, where Mi,j,k=Mi,j,ku1,0,…,un,0,…,u1,k,…,un,k denotes the k th Adomian polynomial of the entry Mi,j. Then(a)
h=Mq⇒hk=∑l=0kMk−lql(b)
q=M−1h⇒M0qk=hk−∑l=0k−1Mk−lql, k=0,1,…Proof.
(a)
Leth=Mq; then,(13)hk=Mqk=∑j=1nM1,jqjk⋮∑j=1nMn,jqjk=∑j=1n∑l=0kM1,j,k−lqj,l⋮∑j=1n∑l=0kMn,j,k−lqj,l=∑l=0kMk−lql.(b)
Letq=M−1h; then h=Mq. Using (a), we have(14)hk=Mqk=∑l=0kMk−lql=M0qk+∑l=0k−1Mk−lql.
Therefore,qk is given by(15)M0qk=hk−∑l=0k−1Mk−lql.
This completes the proof of the theorem.
We start our technique by writing first the second equation of DAE (1) as(16)v′=M−1ufu,v,t−GTu,tλ.
Now, let(17)hu,v,λ,t=fu,v,t−GTu,tλand(18)qu,v,λ,t=M−1uhu,v,λ,t.
Then, (16) can be written as(19)v′=qu,v,λ,t.
Then, we expand the solution componentsut,vt, and λt in power series as(20)ut=∑k=0∞uktk,vt=∑k=0∞vktk,λt=∑k=0∞λktk,where the unknown vectors uk∈ℝnu,vk∈ℝnu, and λk∈ℝnλ, k=0,1,2,…, will be determined later by our method.
The nonlinear termqu,v,t is expanded using the Adomian polynomial as(21)q∑k=0∞uktk,∑k=0∞vktk,t=∑n=0∞qktk,where qk denotes the vector of the Adomian polynomials of the components of q.
Substituting (20) into the first equation of the DAE (1) and (19), we obtain(22)∑k=1∞kuktk−1=∑k=0∞vktkand(23)∑k=1∞kvktk−1=∑k=0∞qktk.
Equations (22) and (23) lead to the recursions(24)vk−1=kuk,k≥1,and(25)kk−1uk=qk−2,k≥2.
Further, the nonlinear termsfu,v,t, gu,t, Mu, and Gu,t are also expanded in power series using the Adomian polynomials as(26)f∑k=0∞uktk,∑k=0∞vktk,t=∑k=0∞fktk,g∑k=0∞ukt,t=∑k=0∞gktk,and(27)M∑k=0∞ukt=∑k=0∞Mktk,G∑k=0∞ukt,t=∑k=0∞Gktk.where fk≔fku0,v0,…,uk,vk and gk≔gku0,…,uk denote the vectors of the Adomian polynomials of the components of the vectors f and g, respectively. The matrices Mku0,…,uk and Gku0,…,uk are the matrices of the Adomian polynomials of the entries of the matrices M and G.
Then, from (18) and the result (b) of the previous theorem, we have(28)qk=M0−1hk−∑l=0k−1Mk−lql.
Now, using (17) and the result (a) of the previous theorem, we obtain(29)hk=fk−∑l=0kGk−lTλl.
Substituting equation (29) into equation (28), we deduce that(30)qk=M0−1fk−∑l=0kGk−lTλl−∑l=0k−1Mk−lql,which can be written as(31)qk=−M0−1G0Tλk+M0−1fk−∑l=0k−1Gk−lTλl+Mk−lql.
Equation (31) can be written as(32)qk−2=−M0−1G0Tλk−2+rk−2,where(33)rk−2=M0−1fk−2−∑l=0k−3Gk−2−lTλl+Mk−2−lql.
Here, it is important to note thatrk−2 depends only on u0,…,uk−1, v0,…,vk−1, and λ0,…,λk−3.
Equations (25) and (32) lead to the algebraic linear system for uk and λk−2.(34)kk−1uk+M0−1G0Tλk−2=rk−2,k≥2.
Using the expansion ofg from (26) in the last equation of DAE (1), we have(35)gk=0,k=0,1,….
We splitgk as(36)gk=−sk+G0uk,k≥2.
From the definition of the Adomian polynomialgk of g, the term sk depends only on u0,…,uk−1.
Using now equation (35), we obtain the algebraic linear system for uk.(37)G0uk=sk,k≥2.
Finally, the unknown vectorsuk and λk are determined by solving the system linear algebraic system formed by (34) and (37). Then, the unknown vector vk is obtained from (24). Since the Jacobian Gu,t has full row rank and the matrix M(u) is positive definite, then the linear algebraic system formed by (34) and (37) has a unique solution. To solve this system, we multiply system (34) from the left by G0 then use (37) to obtain the following linear algebraic system for λk−2.(38)G0M0−1G0Tλk−2=G0rk−2−kk−1sk,k≥2.
Since the rank condition (4) holds and Mu is positive definite, then the matrix G0M0−1G0T is nonsingular and system (38) determines λk−2 uniquely for k≥2. Substituting the expression of λk−2 into system (34), we can compute the unknown vector uk, for k≥2. Using now (24), we find the unknown vector vk, for k≥1. The initial terms for the recursion (38), (34), and (24) are u0=u0, v0=v0.
Finally, we obtain an approximate solution for the DAE initial-value problem (1)-(2) as(39)ut=∑k=0Kuktk,vt=∑k=0K−1vktk,λt=∑k=0K−2λktk,where K is the order of approximation.
The process of the PSM for solving the DAE initial-value problem (1)-(2) is given by (Algorithm 1).
The PSM may lead to solutions of the DAE initial-value problem (1)-(2) that converge only for small intervals. Therefore, we extend the interval of convergence of the PSM by developing a multistage PSM (MSPSM). The process of the MSPSM is described as follows (Algorithm 2).Algorithm 1: PSM Algorithm.
inputreadK (order of approximation);readnu and nλ<nu, M∈ℝnu×nu, f∈ℝnu and g∈ℝnλ;readu0∈ℝnu and v0∈ℝnu (consistent initial conditions);outputthe approximate solutionut,vt,λt;initializationu0:=u0,v0:=v0,u1:=v0;while0≤k≤KdocomputeMk, Gk, fk, gk.endH:=G0M0−1G0T;while2≤k≤Kdosk:=−gk+G0uk;rk−2:=M0−1fk−2−∑l=0k−3Gk−2−lTλl+Mk−2−lql;SolveHλk−2=G0rk−2−kk−1sk for λk−2;qk−2:=rk−2−M0−1G0Tλk−2;uk:=qk−2/kk−1;vk−1:=kukendAlgorithm 2:MSPSM 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 PSM to the DAE (1) on 0,t1 using u01=u0, v01=v0 to obtain the approximate solution u1t=∑k=0Kuk1tk, and v1t=∑k=0K−1vk1tkλ1t=∑k=0K−2λk1tk,0≤t≤t1(4)
Fromi=2 until N, apply the PSM to DAE (1) on ti−1,ti using u0i=ui−1ti−1, v0i=vi−1ti−1 to obtain the approximate solution 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)
The approximate solution of DAE initial-value problem (1)-(2) on 0,T is ut=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 Simulation and Discussion

In this section, an example of DAEs (1)-(2) taken from [25] is solved to illustrate and demonstrate the effectiveness of our technique. To this end, a maple code is developed and we obtained the numerical results discussed in this section. The example used here describes the dynamics of a two-link planar robot arm subject to a prescribed trajectory (see Figure 1).Figure 1
A constrained two-link planar robot arm.The robot arm is composed by two rigid links. The left end of the first link (link on the left) is fixed at the origin of thexy-plane but allowed to rotate. The angle of rotation of this end is measured from the x-axis and is denoted by θ1. The right end of this first link is connected to the left end of the second link (link on the right) at x1,y1. This second link is allowed to rotate in the x, y-plane and the angle of rotation is denoted by θ2 as shown in Figure 1. The masses and lengths of the links are denoted by mi and li, i=1,2, respectively. The coordinates x1,y1 of the middle joint of the robot arm are given by(40)x1=l1cosθ1,y1=l1sinθ1,and the coordinates x2,y2 of the end-effector of the robot arm are(41)x2=x1+l2cosθ1+θ2,y2=y1+l2sinθ1+θ2.If the end-effector of the robot arm is constrained to move along thex-axis, then y2=0 and we obtain the following constraint on the angles θ1 and θ2.(42)l1sinθ1+l2sinθ1+θ2=0.For this robot arm example, the mass matrix is(43)Mθ1,θ2=m1l123+m2l12+l223+l1l2cosθ2m2l223+0.5l1l2cosθ2m2l223+0.5l1l2cosθ2m2l223.The force term is(44)fθ1,θ2,dθ1dt,dθ2dt,t=l1cosθ1+l2cosθ1+θ2dθ1/dt−3θ1l2cosθ1+θ2dθ1/dt+1−1.5cosθ2θ1.And the constraint function is(45)gθ1,θ2,t=l1sinθ1+l2sinθ1+θ2.In this numerical simulation, we takel1=l2=1 and m1=m2=3. Using the notation of the DAE system (1), with nu=2, nλ=1, u=u1,u2T=θ1,θ2T, and v=du/dt, we have(46)Mu=5+3cosu21+1.5cosu21+1.5cosu21and(47)fu,v,t=cosu1+cosu1+u2v1−3u1cosu1+u2v1+1−1.5cosu2u1.The constraint function becomes(48)gu,t=sinu1+sinu1+u2.And its Jacobian Gu,t=gu′u,t=cosu1+cosu1+u2,cosu1+u2 is full row rank nλ=1. It is easy to check that the exact solution of the DAE corresponding to this robot arm problem is ut=sint,−2sintT, vt=cost,−2costT, and λt=cost. This solution corresponds to the initial conditions(49)u0=0,0T,v0=1,−2T.These initial conditions are consistent, since they satisfygu0,0=0 and gu′u0,0v0=0. They correspond to the position where the robot arm initially lies along the x-axis and starts moving with the initial velocity v0=1,−2T. The DAE corresponding to this robot arm example is an index-3 DAE, therefore being difficult to solve even numerically. To solve it using our technique, we have implemented Algorithm 1 presented in Section 3 in Maple 15. Taking the number of approximation terms K=10, our technique was able to compute the first terms of the Maclaurin power series of the exact solution as(50)ut=t1−2−t33!1−2+t55!1−2−t77!1−2+t99!1−2,vt=1−2−t22!1−2+t44!1−2−t66!1−2+t88!1−2and(51)λt=1−t22!+t44!−t66!+t88!.This truncated power series may give accurate solutions only on small intervals.To extend the domain of convergence of the numerical solution, we implemented Algorithm2 of the multistage PSM (MSPSM) given in Section 3. For the numerical simulations, we applied the MSPSM over the time interval 0,T=0,12, with an order of approximation K=5, and a number of interval subdivisions N=240. Starting the recursion with u0=u0=0,0T, v0=v0=1,−2T and using the MSPSM, we obtain the numerical results shown in Figures 1 to 9. These figures show the solution components together with their corresponding absolute errors. In the figures on the left, the solid line represents the numerical solution and the dashed line represents the exact solution. Figure 2 on the left shows the numerical solution for the angle θ1=u1 together with the exact angle θ1 whereas the figure on the right shows the absolute error for the angle θ1. From Figure 2, we can see that there is a good agreement between the computed and exact values for the angle θ1 and the absolute error is 2×10−7. In Figure 3, the left figure shows the numerical solution for the angle θ2=u2 together with the exact angle θ2 whereas the figure on the right shows the absolute error for the angle θ2. From Figure 3, we can see that there is a good agreement between the computed and exact values for the angle θ2 and the absolute error is 5×10−7. In Figure 4, the left figure shows the numerical solution for the velocity of the angle θ1 together with the exact velocity of this angle whereas the figure on the right shows the absolute error for the velocity of this angle. From Figure 4, we can see that there is a good agreement between the computed and exact velocities for this angle θ1. The absolute error is 2×10−7. Similarly, in Figure 5, the left figure shows the numerical solution for the velocity of the angle θ2 together with the exact velocity of this angle whereas the figure on the right shows the absolute error for the velocity of this angle. From Figure 5, we can see that there is a good agreement between the computed and exact velocities for this angle θ2. The absolute error is 5×10−7. In Figure 6, the left figure shows the numerical solution for the Lagrange multiplier λ together with the exact value of λ whereas the figure on the right shows the absolute error for this algebraic variable λ. From Figure 6, we can see that there is a good agreement between the computed and exact values for λ. The absolute error for this algebraic variable is 3×10−7. In Figure 7, the left figure shows the x component x1 of the computed position x1,y1 of the middle joint of the robot arm together with the exact values of x1. The figure on the right shows the absolute error in this variable. From Figure 7, we can see that there is a good agreement between the computed and exact values of x1. The absolute error for this variable is 10−7. Similarly, in Figure 8, the left figure shows the x-component x2 of the computed position x2,y2 of the end-effector of the robot arm together with the exact values of x2. The figure on the right shows the absolute error in this variable. From Figure 8, we can see that there is a good agreement between the computed and exact values of x2. The absolute error for this variable is 2×10−7. The numerical values of the y-coordinate of the end-effector of the robot arm are found to be y2=0 with an error of 5×10−18. Therefore, the robot arm end-effector follows the prescribed trajectory (the x-axis) with an accuracy of 7×10−18. Finally, in Figure 9, the left figure shows the y-component y1 of the computed position x1,y1 of the middle joint of the robot arm together with the exact values of y1. The figure on the right shows the absolute error in this variable. From Figure 9, we can see that there is a good agreement between the computed and exact values of y1. The absolute error for this variable is 2×10−7. We have also computed the drift (error) in the position and the velocity constraints gu,t=0 and gu′u,tv=0 to test the accuracy of our new algorithm. These errors are found to be 5×10−18 and 4×10−18, respectively. These results show that DAE constraints are satisfied with a high accuracy and that the end-effector of the robot arm is following the prescribed path with high accuracy. These numerical results also show that the MSPSM has computed the solution accurately for large intervals without violating the constraints. Finally, to increase accuracy in the approximate solution, one can increase the order of approximation K or increase the number of interval subdivisions N.Figure 2
Approximate solution componentu1 (a) and absolute error eu1 (b).
(a)(b)Figure 3
Approximate solution componentu2 (a) and absolute error eu2 (b).
(a)(b)Figure 4
Approximate solution componentv1 (a) and absolute error ev1 (b).
(a)(b)Figure 5
Approximate solution componentv2 (a) and absolute error ev2 (b).
(a)(b)Figure 6
Approximate solution componentλ (a) and absolute error eλ (b).
(a)(b)Figure 7
Approximate solution componentx1 (a) and absolute error ex1 (b).
(a)(b)Figure 8
Approximate solution componentx2 (a) and absolute error ex2 (b).
(a)(b)Figure 9
Approximate solution componenty1 (a) and absolute error ey1 (b).
(a)(b)


## 5. Conclusion

The differential and algebraic equations (DAEs) with index-3 are known to be difficult to solve numerically. In this paper, we propose a novel technique to solve the index-3 DAEs (1). Our technique avoids transformation techniques like index reductions which are known to be costly and complex. Our algorithm instead applies the approximation method directly to the original DAE (1) using the power series method and the Adomian polynomials. This has led to a simple and efficient algorithm that can be easily implemented in Maple or Mathematica to solve the index-3 DAEs (1). It should be emphasized here that only small linear algebraic systems are solved in this technique. The PSM has computed the exact solution in a power series form. In order to extend the domain of convergence of the solution, we have constructed a multistage PSM (MSPSM). To illustrate the effectiveness of our technique, a nonlinear DAE problem describing a two-link planar robot arm is solved. The numerical simulations show that the MSPSM has successfully handled this nonlinear index-3 DAE over large time intervals without violating the DAE constraints. The numerical simulations also show that the robot arm end-effector follows the prescribed trajectory with an accuracy of5×10−18. The method can provide a much greater accuracy by adding few terms into the power series or more time interval subdivisions. Future research will be to generalize this technique to more implicit higher-index DAE systems.



---

*Source: 2712196-2022-02-27.xml*