# Solving Nonlinear Differential Algebraic Equations by an ImplicitGL(n,ℝ) Lie-Group Method

Authors: Chein-Shan Liu
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

We derive an implicitGL(n,ℝ) Lie-group algorithm together with the Newton iterative scheme to solve nonlinear differential algebraic equations. Four numerical examples are given to evaluate the efficiency and accuracy of the new method when comparing the computational results with the closed-form solutions.


---

## Body

## 1. Introduction

In this paper, we propose a novel method to solve nonlinear differential algebraic equations (DAEs), which govern the evolution ofn+m variables xi,i=1,…,n and yj, j=1,…,m with n nonlinear ordinary differential equations (ODEs) and m nonlinear algebraic equations (NAEs):
(1)x˙=f(x,y,t),x(0)=x0,t∈ℝ,x∈ℝn,y∈ℝm,(2)F(x,y,t)=0,F∈ℝm.
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].


## 2. TheGL(n,ℝ) 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 subsetGL(n,ℝ):={G∈ℝn×n∣detG≠0} of the linear space of all n×n nonsingular matrices. Thus, GL(n,ℝ) is an n×n-dimensional manifold. The group composition is given by the matrix multiplication.Here we give a new form of (1) from the GL(n,ℝ) Lie-group structure. The vector field f on the right-hand side of (1) can be written as
(3)x˙=Ax,
where
(4)A=f∥x∥⊗x∥x∥
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 matrixA is well defined, the Lie-group element G generated from the above dynamical system (3) with G˙=AG satisfies detG(t)≠0, such that G∈GL(n,ℝ).


## 3. An ImplicitGL(n,ℝ) Lie-Group Scheme

Equation (3) is a new starting point for the development of the Lie-group GL(n,ℝ) algorithm. In order to develop a numerical scheme from (3) and (4), we suppose that the coefficient matrix A is constant with
(5)a=f-∥x-∥,b=x-∥x-∥
being two constant vectors, which can be obtained by taking the values of f and x at a suitable mid-point of t-∈[t0=0,t], where t≤t0+h and h 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
(6)x˙=b·xa.
Let
(7)w=b·x,
and (6) becomes
(8)x˙=wa.
At the same time, from the above two equations we can derive
(9)w˙=cw,
where
(10)c=a·b
is viewed as a constant scalar. Thus, we have
(11)w(t)=w0exp(ct),
where w0=b·x0.Inserting (11) for w(t) into (8) and integrating the resultant equation, we can obtain
(12)x(t)=[In+η(t)abT]x0,
where x0 is the initial value of x at an initial time t=t0=0, and
(13)η(t)=ect-1c.
Let G(t) be the coefficient matrix before x0 in (12), that is,
(14)G=In+η(t)abT,
which is one sort of elementary matrices. According to [10, 11], one can prove
(15)detG=ect>0,
which means that G is a Lie-group element of GL(n,ℝ).Within a small time step we can suppose that the variablesyj, j=1,…,m are constant in the interval of tk<t<tk+1. As a consequence, we can develop the following implicit scheme for solving the ODEs (1) where y at the kth time step, denoted by yk, is viewed as a parameter.(i)
Give0≤θ≤1.(ii)
Give an initialx0 at an initial time t=t0 and a time stepsize h.(iii)
Fork=0,1,…, we repeat the following computations to a terminal time t=tf:
(16)xk+1=xk+hfk,where fk:=f(xk,yk,tk). With the above xk+1 generated from an Euler step as an initial guess, we can iteratively solve the new xk+1 by
(17)t-k=tk+θh,x-k=(1-θ)xk+θxk+1,f-k=f(x-k,yk,t-k),ak=f-k∥x-k∥,bk=x-k∥x-k∥,ck=ak·bk,dk=xk·bk,ηk=exp(ckh)-1ck,zk+1=xk+ηkdkak.
If zk+1 converges according to a given stopping criterion, such that
(18)∥zk+1-xk+1∥<ε2,
then go to (iii) to the next time step; otherwise, let xk+1=zk+1 and go to the computations in (17) again. In all the computations given below we will use θ=1/2.


## 4. 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 yj, j=1,…,m are constant in that interval of tk<t<tk+1. We give an initial guess of yj, j=1,…,m and insert them into (1). Then, we apply the above implicit scheme to find the next xk+1, supposing that xk is already obtained in the previous time step. When xk+1 are available we insert them into (2) and then apply the Newton iterative scheme to solve yk+1 by
(19)yk+1ℓ+1=yk+1ℓ-B-1F(xk+1,yk+1ℓ,tk+1),
until the following convergence criterion is satisfied:
(20)∥yk+1ℓ+1-yk+1ℓ∥<ε1.
The component Bij of the Jacobian matrix B is given by ∂Fi/∂yj. Below we use some examples to demonstrate the numerical processes in Sections 3 and 4.


## 5. Numerical Examples of DAEs

In order to assess the performance of the newly developed scheme based on the Lie-groupGL(n,ℝ), 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]:
(21)Q˙=keq˙-keλQ0Q,
where λ is subjected to
(22)λ=1Q0Q·q˙>0if∥Q∥=Q0,Q·q˙>0,λ=0if∥Q∥<Q0orQ·q˙≤0.
While ke>0 is known as an elastic modulus, the constant Q0>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:
(23)Q˙=keq˙-keλQ0Q=:f,(24)∥Q∥2=Q02.
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
(25)λ=1Q0Q·q˙.
Inserting it into (23) we obtain a nonlinear ODEs system:
(26)Q˙=keq˙-keQ02Q·q˙Q.
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 SOo(n,1), which can automatically satisfy (24).
We apply the implicitGL(n,ℝ) 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.(i)
Give an initial guess ofλ0, for example, λ0=0.(ii)
Give an initial conditionQ0 at an initial time t=t0 and a time stepsize h.(iii)
Fork=0,1,…, we repeat the following computations to a specified terminal time t=tf:
(27)Qk+1=Qk+hfk.
With the above Qk+1 generated from an Euler step as an initial guess, we then iteratively solve the new Qk+1 by
(28)Q-k=(1-θ)Qk+θQk+1,ak=f-k∥Q-k∥,bk=Q-k∥Q-k∥,ck=ak·bk,dk=Qk·bk,ηk=exp(ckh)-1ck,zk+1=Qk+ηkdkak.
If zk+1 converges according to a given stopping criterion, such that
(29)∥zk+1-Qk+1∥<ε2,
then go to (iv); otherwise, let Qk+1=zk+1 and go to (28).(iv)
Forj=0,1,…, we repeat the following computations:
(30)λj+1=λj-FjFj′,
where the prime denotes the differential with respect to λ, and
(31)ak′=f-k′∥Q-k∥,ck′=bk·ak′,ηk′=ck′[(hck-1)exp(ckh)+1]ck2,Qk+1′=ηk′dkak+ηkdkak′,Fj=∥Qk+1∥2-Q02,Fj′=2Qk+1·Qk+1′.
If λj converges according to
(32)|λj+1-λj|<ε1,
then go to (iii) with λj as an initial guess of λ for the next time step; otherwise, let λj=λj+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 beingq1=e0cos(ωt) and q2=e0sin(ωt). Here we suppose that β=kee0/Q0>1, and the initial stresses are on the yield surface with Q1=Q0cosθ0 and Q2=Q0sinθ0. As shown in Liu [13, 14], the responses of Q1 and Q2 have the following closed-form solutions:
(33)Q1Q0=z+βcosθ0-1βzcos(ωt)-z˙βωzsin(ωt),Q2Q0=z˙βωzcos(ωt)+z+βcosθ0-1βzsin(ωt),
where
(34)z(t)=1+β2-βcosθ0m2[cosh(mωt)-1]+βsinθ0msinh(mωt),
in which m=β2-1.
In Figure1, we plot the response errors of Q1 and Q2 and the error of the consistency condition, which is defined by |Q12+Q22-Q0|, in a time range of t∈[0,10]. We fix ke= 200,000 MPa, Q0=200 MPa, e0=0.002, ω=1, and θ0=0, and the time stepsize used is h=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 SO(n) scheme [15].Example1 of a plastic equation to compute the stresses, showing the numerical errors of (a) Q1, (b) Q2, and (c) consistency condition.
(a)
(b)
(c)Example 2.
We consider an index two nonlinear Hessenberg DAEs [7, 8]:
(35)x˙1=tx22+λ+g1(t),x˙2=texp(x1)+tλ+g2(t),x1+tx2+g3(t)=0,
with x1(0)=x2(0)=λ(0)=0, where
(36)g1(t)=1-t2-t3(1+t)2,g2(t)=1-t-4t2-4t3-t4(1+t)2,g3(t)=-ln(1+t)-t21+t.
The exact solutions of this problem are
(37)x1(t)=ln(1+t),x2(t)=λ(t)=t1+t.
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
(38)Fj=x1k+1+tk+1x2k+1+g3(tk+1),(ak1)′=1∥x-k∥,(ak2)′=t-k∥x-k∥,(x1k+1)′=ak1dkηk′+(ak1)′dkηk,(x2k+1)′=ak2dkηk′+(ak2)′dkηk,Fj′=(x1k+1)′+tk+1(x2k+1)′.
In the above, (x1k+1)′ and (x2k+1)′ denote, respectively, the differentials of x1k+1 and x2k+1 with respect to λ, and aki denotes the ith component of ak.
We useh=10-3, and the problem is solved in a range of t≤1. In Figure 2(a), we show the numerical errors of x1, x2 and λ, of which we can see that the numerical solutions are very accurate. In Figure 2(b), we show the error of |x1(t)+tx2(t)+g3(t)|, 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).Example2 showing (a) the numerical errors of solutions, (b) the error of invariant, and (c) the numbers of iterations.
(a)
(b)
(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:
(39)u¨1=2u2+λu1,u¨2=-2u1+λu2,u12+u22=1.
For (u1(0),u2(0))=(0,1),λ(0)=0, the exact solution is u1(t)=sint2,u2(t)=cost2 and λ(t)=-4t2. The above problem can be viewed as a mechanical control problem to select a suitable controller λ(t) 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:(40)I1=x12+x32-1=u12+u22-1=0,
where we let x1=u1, x2=u˙1, x3=u2 and x4=u˙2. At the same time the above ODEs become
(41)x˙1=x2,x˙2=2x3+λx1,x˙3=x4,x˙4=-2x1+λx3.
Taking the time differential of (40), we can derive
(42)x1x˙1+x3x˙3=x1x2+x3x4=0,
due to x2=x˙1 and x4=x˙3. The above is a first differential form of the constraint in (40). However, it is still not available for the determination of the Lagrange multiplier λ. Thus, the second differential form of (42) is required:
(43)x1x˙2+x˙1x2+x3x˙4+x˙3x4=0.
Inserting (41) and using (40), we can derive
(44)I2=λ+x22+x42=0.
Then, we finally come to a system of ODEs:
(45)x˙1=x2,x˙2=2x3-x1(x22+x42),x˙3=x4,x˙4=-2x1-x3(x22+x42),
and the Lagrange multiplier λ is calculated by
(46)λ=-x22-x42.
Finally, taking the differential of the above equation and inserting (45) for x˙2 and x˙4 we can obtain the differential equation for λ:
(47)λ˙=2x2[x1(x22+x42)-2x3]+2x4[2x1+x3(x22+x42)].
In the above, we have transformed the DAEs in (39) into a set of ODEs system (45) and (47) through three times differentiations. This DAEs system is thus said to be of index three.
In general, we may hope the solution of the system (45) can automatically satisfy the following constraint on (x1,x3):
(48)x12+x32=1,
which is the last equation in (39). However, we should remind that there exists no general numerical integrator which can automatically satisfy (48). For this purpose, the numerical integrator must be particularly designed to meet that requirement. In the above, we have introduced two invariants I1 and I2; for this DAEs system the preservations of these two invariants are utmost important.
We apply the LGDAE to solvexi,i=1,…,4 by using (41) and then iteratively solve the weak form constraint in (42) by the Newton method to find λ. The processes are the same as that in (27)–(32), but now we have different Fj and Fj′:
(49)Fj=x1k+1x2k+1+x3k+1x4k+1,ck′=x-1bk2+x-3bk4∥x-k∥,ηk′=ck′[(hck-1)exp(ckh)+1]ck2,(x1k+1)′=ak1dkηk′,(x2k+1)′=ak2dkηk′+dkηkx-1,(x3k+1)′=ak3dkηk′,(x4k+1)′=ak4dkηk′+dkηkx-3,Fj′=x1k+1(x2k+1)′+(x1k+1)′x2k+1+(x3k+1)′x4k+1+x3k+1(x4k+1)′.
In the above, aki and bkj are, respectively, the ith and jth components of the vectors ak and bk.
Under the convergence criteriaε2=10-15 for inner iterations and ε1=10-6 for outer iterations, we apply the above numerical method with h=10-4 to solve (39). In Figure 3(a) we show the numerical errors of u1,u2, and λ, of which we can see that the numerical solutions are very accurate. In Figure 3(b), we show the errors of |I1| and |I2|, which are almost zero with the order 10-10. It can be seen that the accuracy of u1 and u2 is in the order of h2.5, and the numbers of iterations are few with three for inner iterations and two for outer iterations as shown in Figure 3(c).Example3 showing (a) the numerical errors of solutions, (b) the errors of invariants, and (c) the numbers of iterations.
(a)
(b)
(c)Example 4.
Finally, we consider a more complex pendulum problem [16]:
(50)x˙1=x3-λ2x1,x˙2=x4-λ2x2,x˙3=-λ1x1,x˙4=-λ1x2-1,I1=x12+x22-1=0,I2=x1x3+x2x4=0.
The exact solution of the above problem is not available, which has two Lagrange multipliers λ1 and λ2 and two constraints.
The initial values arex1(0)=1 and x2(0)=x3(0)=x4(0)=0. Under the convergence criteria ε2=10-15 for inner iterations and ε1=10-10 for outer iterations, we apply the LGDAE with h=10-4 to solve (50). In Figure 4(a), we show the numerical solutions. In Figure 4(b), we show the errors of |I1| and |I2|, which are almost to be zero with the order 10-13 for |I1| and the order 10-17 for |I2|. The numbers of iterations are few with three for inner iterations and two for outer iterations as shown in Figure 4(c).Example4 showing (a) the numerical solutions, (b) the errors of invariants, and (c) the numbers of iterations.
(a)
(b)
(c)


## 6.  Conclusions

By recasting the nonlinear ODEs:x˙=f into a quasi-linear Lie-form: x˙=Ax, where A=f/∥x∥⊗x/∥x∥, we have derived an implicit GL(n,ℝ) 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.



---

*Source: 987905-2013-07-16.xml*