• Views 910
• Citations 10
• ePub 22
• PDF 455
`Journal of Applied MathematicsVolume 2013 (2013), Article ID 987905, 8 pageshttp://dx.doi.org/10.1155/2013/987905`
Research Article

## Solving Nonlinear Differential Algebraic Equations by an Implicit Lie-Group Method

Department of Civil Engineering, National Taiwan University, Taipei, Taiwan

Received 26 February 2013; Accepted 28 June 2013

Copyright © 2013 Chein-Shan Liu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We derive an implicit 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.

#### 1. Introduction

In this paper, we propose a novel method to solve nonlinear differential algebraic equations (DAEs), which govern the evolution of variables and , with nonlinear ordinary differential equations (ODEs) and nonlinear algebraic equations (NAEs): Usually, is larger than . When , the DAEs reduce to the ODEs. There are many numerical methods used to solve ODEs, but only a few is used to solve DAEs [15]. 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. The 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 of the linear space of all nonsingular matrices. Thus, is an -dimensional manifold. The group composition is given by the matrix multiplication.

Here we give a new form of (1) from the Lie-group structure. The vector field on the right-hand side of (1) can be written as where is the coefficient matrix. The symbol in denotes the dyadic operation of and , that is, .

Because the coefficient matrix is well defined, the Lie-group element generated from the above dynamical system (3) with satisfies , such that .

#### 3. An Implicit Lie-Group Scheme

Equation (3) is a new starting point for the development of the Lie-group algorithm. In order to develop a numerical scheme from (3) and (4), we suppose that the coefficient matrix is constant with being two constant vectors, which can be obtained by taking the values of and at a suitable mid-point of , where and is a small time stepsize. The variable 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 .

Inserting (11) for into (8) and integrating the resultant equation, we can obtain where is the initial value of at an initial time , and Let be the coefficient matrix before in (12), that is, which is one sort of elementary matrices. According to [10, 11], one can prove which means that is a Lie-group element of .

Within a small time step we can suppose that the variables , are constant in the interval of . As a consequence, we can develop the following implicit scheme for solving the ODEs (1) where at the th time step, denoted by , is viewed as a parameter.(i)Give .(ii)Give an initial at an initial time and a time stepsize .(iii)For , we repeat the following computations to a terminal time : where . With the above generated from an Euler step as an initial guess, we can iteratively solve the new by If converges according to a given stopping criterion, such that then go to (iii) to the next time step; otherwise, let and go to the computations in (17) again. In all the computations given below we will use .

#### 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 , are constant in that interval of . We give an initial guess of , and insert them into (1). Then, we apply the above implicit scheme to find the next , supposing that is already obtained in the previous time step. When are available we insert them into (2) and then apply the Newton iterative scheme to solve by until the following convergence criterion is satisfied: The component of the Jacobian matrix is given by . 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-group , 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 is known as an elastic modulus, the constant is a yield strength of material.
We can view the above equations in the plastic state, that is, , 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 , which can automatically satisfy (24).
We apply the implicit scheme to solve 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 , for example, .(ii)Give an initial condition at an initial time and a time stepsize .(iii)For , we repeat the following computations to a specified terminal time : With the above generated from an Euler step as an initial guess, we then iteratively solve the new by If converges according to a given stopping criterion, such that then go to (iv); otherwise, let and go to (28).(iv)For , 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 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 and . Here we suppose that , and the initial stresses are on the yield surface with and . As shown in Liu [13, 14], the responses of and have the following closed-form solutions: where in which .
In Figure 1, we plot the response errors of and and the error of the consistency condition, which is defined by , in a time range of . We fix 200,000 MPa,  MPa, , , and , and the time stepsize used is . Under the convergence criteria for inner iterations and 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].

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

Example 2. We consider an index two nonlinear Hessenberg DAEs [7, 8]: with , where The exact solutions of this problem are
Under the convergence criteria for inner iterations and for outer iterations, we apply the LGDAE as that in (27)–(32) to solve this problem, where In the above, and denote, respectively, the differentials of and with respect to , and denotes the th component of .
We use , and the problem is solved in a range of . In Figure 2(a), we show the numerical errors of , and , of which we can see that the numerical solutions are very accurate. In Figure 2(b), we show the error of , which is almost zero with the order . 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).

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

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 , the exact solution is and . 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: where we let , , and . At the same time the above ODEs become
Taking the time differential of (40), we can derive due to and . 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: Inserting (41) and using (40), we can derive Then, we finally come to a system of ODEs: and the Lagrange multiplier is calculated by Finally, taking the differential of the above equation and inserting (45) for and we can obtain the differential equation for : 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 : 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 and ; for this DAEs system the preservations of these two invariants are utmost important.
We apply the LGDAE to solve 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 and : In the above, and are, respectively, the th and th components of the vectors and .
Under the convergence criteria for inner iterations and for outer iterations, we apply the above numerical method with to solve (39). In Figure 3(a) we show the numerical errors of , and , of which we can see that the numerical solutions are very accurate. In Figure 3(b), we show the errors of and , which are almost zero with the order . It can be seen that the accuracy of and is in the order of , and the numbers of iterations are few with three for inner iterations and two for outer iterations as shown in Figure 3(c).

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

Example 4. Finally, we consider a more complex pendulum problem [16]: The exact solution of the above problem is not available, which has two Lagrange multipliers and and two constraints.
The initial values are and . Under the convergence criteria for inner iterations and for outer iterations, we apply the LGDAE with to solve (50). In Figure 4(a), we show the numerical solutions. In Figure 4(b), we show the errors of and , which are almost to be zero with the order for and the order for . The numbers of iterations are few with three for inner iterations and two for outer iterations as shown in Figure 4(c).

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

#### 6. Conclusions

By recasting the nonlinear ODEs: into a quasi-linear Lie-form: , where , we have derived an implicit 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.

#### Acknowledgments

The project NSC-102-2221-E-002-125-MY3 and the 2011 Outstanding Research Award from Taiwan’s National Science Council and the 2011 Taiwan Research Front Award from Thomson Reuters, granted to the author, are highly appreciated. Also the author is a Lifetime Distinguished Professor of National Taiwan University, since 2013.

#### References

1. C. Arévalo, S. L. Campbell, and M. Selva, “Unitary partitioning in general constraint preserving DAE integrators,” Mathematical and Computer Modelling, vol. 40, no. 11-12, pp. 1273–1284, 2004.
2. U. M. Ascher and L. R. Petzold, “Projected implicit Runge-Kutta methods for differential-algebraic equations,” SIAM Journal on Numerical Analysis, vol. 28, no. 4, pp. 1097–1120, 1991.
3. R. P. K. Chan, P. Chartier, and A. Murua, “Post-projected Runge-Kutta methods for index-2 differential-algebraic equations,” Applied Numerical Mathematics, vol. 42, no. 1-3, pp. 77–94, 2002.
4. J. Sand, “On implicit Euler for high-order high-index DAEs,” Applied Numerical Mathematics, vol. 42, no. 1–3, pp. 411–424, 2002.
5. C.-S. Liu, “Preserving constraints of differential equations by numerical methods based on integrating factors,” CMES: Computer Modeling in Engineering & Sciences, vol. 12, no. 2, pp. 83–107, 2006.
6. M. M. Hosseini, “Adomian decomposition method for solution of differential-algebraic equations,” Journal of Computational and Applied Mathematics, vol. 197, no. 2, pp. 495–501, 2006.
7. M. M. Hosseini, “Adomian decomposition method for solution of nonlinear differential algebraic equations,” Applied Mathematics and Computation, vol. 181, no. 2, pp. 1737–1744, 2006.
8. F. Soltanian, S. M. Karbassi, and M. M. Hosseini, “Application of He's variational iteration method for solution of differential-algebraic equations,” Chaos, Solitons & Fractals, vol. 41, no. 1, pp. 436–445, 2009.
9. M. Saravi, E. Babolian, R. England, and M. Bromilow, “System of linear ordinary differential and differential-algebraic equations and pseudo-spectral method,” Computers & Mathematics with Applications, vol. 59, no. 4, pp. 1524–1531, 2010.
10. C.-S. Liu, “A method of Lie-symmetry $GL\left(n,ℝ\right)$ for solving non-linear dynamical systems,” International Journal of Non-Linear Mechanics, vol. 52, pp. 85–95, 2013.
11. C.-S. Liu, “A state feedback controller used to solve an ill-posed linear system by a $GL\left(n,ℝ\right)$ iterative algorithm,” Communications in Numerical Analysis, vol. 2013, Article ID cna-00181, 22 pages, 2013.
12. H.-K. Hong and C.-S. Liu, “Internal symmetry in the constitutive model of perfect elastoplasticity,” International Journal of Non-Linear Mechanics, vol. 35, no. 3, pp. 447–466, 2000.
13. C.-S. Liu, “Intermittent transition to quasiperiodicity demonstrated via a circular differential equation,” International Journal of Non-Linear Mechanics, vol. 35, no. 5, pp. 931–946, 2000.
14. C.-S. Liu, “The steady-state behavior of the Prandtl-Reuss material bifurcated under nonproportional circular strain paths,” Journal of the Chinese Institute of Engineers, vol. 26, no. 2, pp. 173–190, 2003.
15. C.-S. Liu, “A Lie-group $DSO\left(n\right)$ method for nonlinear dynamical systems,” Applied Mathematics and Letters, vol. 26, pp. 710–717, 2013.
16. K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, vol. 14 of Classics in Applied Mathematics, SIAM, Philadelphia, Pa, USA, 1996.