Solving Nonlinear Differential Algebraic Equations by an Implicit GL ( n , R ) Lie-Group Method

Usually, n is larger than m. When m = 0, the DAEs reduce to the ODEs. There are many numerical methods used to solve ODEs, but only a few is used to solve DAEs [1–5]. A lot of engineering problems are modelled as a combination of ODEs and NAEs, which are abbreviated as differential algebraic equations (DAEs). The DAEs are both numerically and analytically difficult than the ODEs. Recently, there were some new methods to solve DAEs, for example, Adomian decomposition method [6, 7], variational iterative method [8], and pseudospectral method [9].


The 𝐺𝐿(𝑛,R) Structure of Differential Equations System
The Lie-group is a differentiable manifold, which is endowed with a group structure that is compatible with the underlying topology of the manifold.The Lie-group method can provide a better algorithm that retains the orbit generated from numerical solution on the manifold which is associated with the Lie-group.The general linear group is a Lie group, whose manifold is an open subset (, R) of the linear space of all  ×  nonsingular matrices.Thus, (, R) is an  × -dimensional manifold.The group composition is given by the matrix multiplication.
Here we give a new form of (1) from the (, R) Liegroup structure.The vector field f on the right-hand side of (1) can be written as where is the coefficient matrix.The symbol ⊗ in u ⊗ y denotes the dyadic operation of u and y, that is, (u ⊗ y)z = y ⋅ zu.
Because the coefficient matrix A is well defined, the Liegroup element G generated from the above dynamical system (3) with Ġ = AG satisfies det G() ̸ = 0, such that G ∈ (, R).

An Implicit 𝐺𝐿(𝑛,R) Lie-Group Scheme
Equation ( 3) is a new starting point for the development of the Lie-group (, R) algorithm.In order to develop a numerical scheme from (3) and ( 4), we suppose that the coefficient matrix A is constant with being two constant vectors, which can be obtained by taking the values of f and x at a suitable mid-point of  ∈ [ 0 = 0, ], where  ≤  0 + ℎ and ℎ is a small time stepsize.The variable y is supposed to be a constant vector in this small time interval.Thus, from (3) and (4) we have Let and ( 6) becomes At the same time, from the above two equations we can derive where is viewed as a constant scalar.Thus, we have where  0 = b ⋅ x 0 .Inserting (11) for () into (8) and integrating the resultant equation, we can obtain where x 0 is the initial value of x at an initial time  =  0 = 0, and Let G() be the coefficient matrix before x 0 in (12), that is, which is one sort of elementary matrices.According to [10,11], one can prove which means that G is a Lie-group element of (, R).Within a small time step we can suppose that the variables   ,  = 1, . . .,  are constant in the interval of   <  <  +1 .As a consequence, we can develop the following implicit scheme for solving the ODEs (1) where y at the th time step, denoted by y  , is viewed as a parameter.
(ii) Give an initial x 0 at an initial time  =  0 and a time stepsize ℎ.
(iii) For  = 0, 1, . .., we repeat the following computations to a terminal time  =   : where f  := f(x  , y  ,   ).With the above x +1 generated from an Euler step as an initial guess, we can iteratively solve the new x +1 by If z +1 converges according to a given stopping criterion, such that then go to (iii) to the next time step; otherwise, let x +1 = z +1 and go to the computations in (17) again.In all the computations given below we will use  = 1/2.

Newton Iterative Scheme for DAEs
Now, we turn our attention to the DAEs defined in ( 1) and (2).Within a small time step we can suppose that the variables   ,  = 1, . . .,  are constant in that interval of   <  <  +1 .We give an initial guess of   ,  = 1, . . .,  and insert them into (1).Then, we apply the above implicit scheme to find the next x +1 , supposing that x  is already obtained in the previous time step.When x +1 are available we insert them into (2) and then apply the Newton iterative scheme to solve y +1 by until the following convergence criterion is satisfied: The component   of the Jacobian matrix B is given by   /  .Below we use some examples to demonstrate the numerical processes in Sections 3 and 4.

Numerical Examples of DAEs
In order to assess the performance of the newly developed scheme based on the Lie-group (, R), let us investigate the following four examples of DAEs.
Example 1.Using the on-off switching criteria, we can synthesize the flow model of perfect plasticity into a two-phase system [12]: where  is subjected to While   > 0 is known as an elastic modulus, the constant  0 > 0 is a yield strength of material.We can view the above equations in the plastic state, that is,  > 0, as a system of DAEs: Now we explain that ( 23) and ( 24) are index two DAEs.
Taking the time differential of (24) and inserting (23) into the resultant equation we can solve  by Inserting it into (23) we obtain a nonlinear ODEs system: A further differential of (25) and inserting (26) leads to a differential equation for .Usually, when one applies the general purpose numerical integration method to solve (26), it cannot guarantee that the yield condition in (24) can be automatically satisfied.Hong and Liu [12] have developed the exponential-based scheme from the Lorentz group   (, 1), which can automatically satisfy (24).We apply the implicit (, R) scheme to solve Q through (23) and then iteratively solve the unknown function  through (24) by the Newton iterative method.The numerical processes of this implicit Lie-group DAE (LGDAE) are given below.
(ii) Give an initial condition Q 0 at an initial time  =  0 and a time stepsize ℎ.
(iii) For  = 0, 1, . .., we repeat the following computations to a specified terminal time  =   : With the above Q +1 generated from an Euler step as an initial guess, we then iteratively solve the new Q +1 by If z +1 converges according to a given stopping criterion, such that then go to (iv); otherwise, let Q +1 = z +1 and go to (28).
(iv) For  = 0, 1, . .., we repeat the following computations: where the prime denotes the differential with respect to , and If   converges according to then go to (iii) with   as an initial guess of  for the next time step; otherwise, let   =  +1 and go to (28).
In order to assess the performance of the above numerical method, we consider a strain control case with the strain components being  1 =  0 cos() and  2 =  0 sin().Here we suppose that  =    0 / 0 > 1, and the initial stresses are on the yield surface with  1 =  0 cos  0 and  2 =  0 sin  0 .As shown in Liu [13,14], the responses of  1 and  2 have the following closed-form solutions: where in which  = √ 2 − 1.
In Figure 1, we plot the response errors of  1 and  2 and the error of the consistency condition, which is defined by |√ 2 1 +  2 2 −  0 |, in a time range of  ∈ [0, 10].We fix   = 200,000 MPa,  0 = 200 MPa,  0 = 0.002,  = 1, and  0 = 0, and the time stepsize used is ℎ = 0.001.Under the convergence criteria  2 = 10 −8 for inner iterations and  1 = 10 −8 for outer iterations, we apply the LGDAE to solve the above problem.From Figure 1(c), we can observe that the LGDAE can retain the consistency condition very well.As shown in Figures 1(a) and 1(b), the accuracy of LGDAE is better than that obtained by the exponential-based scheme [12] and the () scheme [15].
Example 2. We consider an index two nonlinear Hessenberg DAEs [7,8]: with  1 (0) =  2 (0) = (0) = 0, where Under the convergence criteria  2 = 10 −15 for inner iterations and  1 = 10 −10 for outer iterations, we apply the LGDAE as that in ( 27)-(32) to solve this problem, where We use ℎ = 10 −3 , and the problem is solved in a range of  ≤ 1.In Figure 2(a), we show the numerical errors of  1 ,  2 and , of which we can see that the numerical solutions are very accurate.In Figure 2(b), we show the error of | 1 () +  2 () +  3 ()|, which is almost zero with the order 10 −11 .It can be seen that the numbers of iterations are few with three to six for inner iterations and two or three for outer iterations as shown in Figure 2(c).
Example 3. We consider an index three differential algebraic equations system given by Sand [4], which describes the position of a particle on a circular track: For ( 1 (0),  2 (0)) = (0, 1), (0) = 0, the exact solution is  1 () = sin  2 ,  2 () = cos  2 and () = −4 2 .The above problem can be viewed as a mechanical control problem to select a suitable controller () changing the system's stiffness such that the orbit of the mechanical system can really trace a circle in time.
We use this example to demonstrate how to transform the above DAEs to a full system of ODEs with the following strong form constraint:

Conclusions
By recasting the nonlinear ODEs: ẋ = f into a quasi-linear Lie-form: ẋ = Ax, where A = f/‖x‖ ⊗ x/‖x‖, we have derived an implicit (, R) Lie-group algorithm together with the Newton iterative scheme, namely, the LGDAE, to solve the nonlinear differential algebraic equations.Four numerical examples were given to assess the performance of the novel method, which is easily implemented to solve the nonlinear differential algebraic equations with a high efficiency and accuracy.

4 JournalFigure 1 :
Figure 1: Example 1 of a plastic equation to compute the stresses, showing the numerical errors of (a)  1 , (b)  2 , and (c) consistency condition.

Figure 2 :
Figure 2: Example 2 showing (a) the numerical errors of solutions, (b) the error of invariant, and (c) the numbers of iterations.

Figure 3 :Figure 4 :
Figure 3: Example 3 showing (a) the numerical errors of solutions, (b) the errors of invariants, and (c) the numbers of iterations.