Efficient Numerical Methods for Solving Differential Algebraic Equations

Ampon Dhamacharoen

Department of Mathematics, Burapha University, Chonburi, Thailand.

DOI: 10.4236/jamp.2016.41007                

Abstract

This research aims to solve Differential Algebraic Equation (DAE) problems in their original form, wherein both the differential and algebraic equations remain. The Newton or Newton-Broyden technique along with some integrators such as the Runge-Kutta method is coupled together to solve the problems. Experiments show that the method developed in this paper is efficient, as it demonstrates that implementation of the method is not difficult, and such method is able to provide approximate solutions with ease within some desired accuracy standards.

Keywords

Differential-Algebraic Equations, Newton-Broyden Method, Index-2 Hessenberg DAE

Received 19 November 2015; accepted 8 January 2016; published 11 January 2016

1. Introduction

A differential algebraic Equation (DAE) is an equation involving an unknown function and its derivatives. A DAE in its most general form is given by the following:

(1)

where R is the set of real number, , , and. In this form, the relationship between the variables and derivatives may be implicit. In some systems, the equations may be written in the explicit form of derivatives, as follows:

(2a)

(2b)

where, This set of equations is called a semi-explicit DAE system.

Differential algebraic equations arise in the mathematical modeling of a wide variety of problems found in engineering and science, such as multi-body and flexible body mechanics, electrical circuit design, optimal control, incompressible fluids, molecular dynamics, chemical kinetics and chemical process control [1] -[5] .

DAEs can be transformed into ODE problems via differentiation. The number of differentiations needed in the transforming process is called the differentiation index. This number can describe some characteristics of the problem. In general, the higher the index of a DAE, the more difficulties one can expect in its numerical solution [1] - [3] [6] .

Although DAEs can be transformed into an explicit ODE so that it can be solved using the methods of ODE, there are still many numerical methods that can solve DAE directly. In the DAE solvers software, the numerical approaches for the solution of DAEs can be divided roughly into two classes: a) direct discretizations of the given system; and b) methods which involve a reformulation (e.g. index reduction), combined with a discretization. Direct discretizations are easier to use, but are limited in their utilitisation to essentially index-1, index-2 Hessenberg, and index -3 Hessenberg DAE systems, while a reformulation may be costly, and it may also require more input from the user and involve more user intervention [3] [6] .

In this research, we will place emphasis on constructing an algorithm for solving a semi-explicit DAE. The approach employed is to firstly solve the system of algebraic equations, and subsequently solve the differential equations using the derived information. The Newton-Broyden method plays a key role in solving algebraic equations, since it performs almost as well as the Newton method, but requires less energy, and in addition, it also outperforms other methods of the same order of convergence [7] - [9] .

2. The Newton-Broyden Method

The method, which was first proposed by Dhamacharoen in 2011 [7] , [8] , aims to solve the equation of the form:

(3)

The Newton scheme of this problem is:

(4)

Replacing by, gives

(5)

with updating:

(6)

where. Equation (6) is called the Broyden rank-1 update, and (5) with the update (6) is called the Newton-Broyden method. This method retains the good part of the Newton method, and replaces the difficult part of Newton’s with Broyden’s. With good initial guesses for z₀ and D₀, the Newton-Broyden Scheme (5) will produce a sequence that converges to a solution of (3), with q-super linear order of convergence.

Although the order of convergence of the Newton-Broyden method is equal to that of the Broyden method, and less than that of the Newton method, in practice, the Newton-Broyden performs well in the sense that a good initial guess is easily found, and it reaches the solution in a reasonable number of ierations. As shown in [7] and [8] , for the same problem and using the same initial guess, the sequences from the Newton method and from the Newton-Broyden method reach the solution, while that from the Broyden does not. In addition, the Newton-Broyden method requires less amount of work in comparison with the Newton method.

Note that if function F is linear, then Broyden’s update and Newton-Broyden’s update coincide.

3. Constructing the Method

The initial value problems:

Assume that the DAE is expressed in the form (2a), (2b) (renumbering)

(7a)

(7b)

(7c)

in which the system of equations is to be solved for x(t) and y(t), where t Î [a, b]. In solving the system numerically, the conditions on the initial values of y must be imposed sufficiently for the system to have a unique solution. If (7b) can express x(t) in terms of t and y(t) explicitly, then the system becomes a pure ODE. Therefore, we consider the case when (7b) expresses x(t) implicitly.

Let the interval [a, b] be divided into n subintervals. In each interval, we will solve (7a) numerically for. In an initial value problem, the value is specified, and the value can subsequently be solved from (7b). In order to use the Runge-Kutta method, in each interval the value, and, which are needed for computing, and respectively, can each be solved from (7b) using the values of, and. Once the value is computed, can then be solved from (7b). Advance to the next interval, and repeat this procedure until we reach the last subinterval.

Suppose the term x(t) is missing from the Equation (7b). Therefore, the system becomes:

(8a)

(8b)

(8c)

in which the system is called index-2 Hessenberg. In solving the differential Equation (8a), the variable x(t) are treated as unknowns. In each interval, the value of is given as a guess, and (8a) is then solved

numerically for. Check the condition (8b). Update the value of and iterate the process until the condition (8b) is met. Subsequently, advance to the next interval using the last value of as the guessing value for. Repeat this procedure until we complete the last subinterval.

In solving the differential equations, we may use the 4^th order Runge-Kutta method, or the 4^th order Taylor method, since their local error is and the total error is, which is acceptable and also easy to im-

plement. In solving the algebraic equations associated in the problem, the Newton method is used in (7b), and the Newton-Broyden method is used in (8b).

As per the procedure described above, the formulation will be as follows:

Partition the interval [a, b] into n subintervals., and

Problem (7a), (7b):

Define the function F:,.

The problem become (3)

which can be solved using the Newton method.

P1: Newton Method.

Prescribe a small positive real number e.

Guess the value z.

A: Compute F(z).

Check if ||F(z)|| < e

If so, proceed to B. If not, carry out the next step.

Compute F'(z).

Solve the system F'(z)b = −F(z), for b.

Set z = z + b.

Go to A.

B: End.

The process for solving the DAE will be as follows:

Algorithm A:

Initial step:

1) Set the given initial condition.

2) Solve for, using P1: Set.

Main step:

For; process the following steps:

Step 1. Solve the initial value problem (7a) 1 step to obtain, by proceeding using the following steps:

1) Compute:, , Solve for. (Using P1).

2) Compute:

, Solve for. (Using P1).

3) Compute:

, Solve for. (Using P1).

4) Compute:

5)

Step 2. Solve (7b) for, using P1: Set.

Proceed to the next interval.

Problem (8a), (8b):

Define the function F:,.

The problem become (3)

which can be solved using the Newton-Broyden method.

P2: One Step Runge-Kutta Method

Given the value and b.

Compute:, ,

Compute:, ,

Compute:, ,

Compute:

Then

The process for solving the DAE will be as follows:

Algorithm B:

Initial step:

1) Set the given initial condition.

2) Guess the initial condition..

3) Guess the initial matrix D (may be = I), Guess the value b.

4) Solve the initial value problem (8a) 1 step, to obtain, using P2.

5) Compute the value

Main step:

For; process the following steps:

Step 1. Check the condition where e is a prescribed small number.

If it passes, proceed to the next interval (next i).

If it fails, carry out step 2.

Step 2.

1) Set

2) Compute the matrix

3) Solve for b.

4) Set

5) Solve the initial value problem (8a) 1 step (using P2) to obtain

6) Compute the new value

7) Update the matrix

8) Go to step 1.

End.

4. Experiments

Three examples are used to illustrate the method. The first problem is an index-1 Hessenberg DAE system, with nonlinear differential equations and a nonlinear algebraic equation, as follows:

initial conditions,

algebraic equation

Using Algorithm A, this problem is solved and has a nice result with a small error as compared to the exact solution

,

.

Some results are illustrated in Table 1. (xs, ys and zs are values from the exact solution).

The second problem is a classical example of the DAE problem which is “The pendulum problem”, as expressed in the xy co-ordinate plane as follows:

(8)

where g = 9.8. This problem is an index 2 semi-explicit DAE problem. [3] , [10] . Proceeding to solve the problem numerically, we let L = 1, and impose the initial conditions x(0) = 1, y(0) = 0. Using the new variables, we have a system of differential algebraic equations as follows:

Table 1. Resulting values for x(t), y(t) and z(t) from Algorithm A, and the errors.

initial conditions:

(9a)

algebraic equation:

(9b)

If we let, , substitute in the problem (8) and manipulate some derivatives and algebra, we will obtain the equivalent initial value problem of an ordinary differential equation:

(10)

From (9b) if we differentiate the equation twice, we will obtain the equation:

(11)

The system (9a) with (11) becomes a pure ODE, since the function l(t) is expressed in the terms and explicitly.

In this experiment, we solve (9a), (9b) using Algorithm B, and provide some of the results in Table 2. System (10) and (9a) with (11) are also solved, and they provide the same results, wherein their values are used to com-

Table 2. Result values for x (t), y (t) and l (t) from Algorithm B, and the errors.

pare to the results from Algorithm B, as shown in Table 2. The value xs, ys and ls are from problem (10). Note that the value of l(t)'s are at the mid-points t + h/2.

The third problem is an index-2 Hessenberg DAE system with nonlinear differential equations and a nonlinear algebraic equation, as follows:

initial conditions,

algebraic equation

Note that this problem is similar to the first problem, but there is no variable z(t) in the algebraic equation.

Using Algorithm B, this problem is solved, and provides a nice result as shown in Table 3. The exact solution is as follows:

Algorithm B gives the results for x(t) and y(t), with a small error. However, the solution for z(t) has some errors less than 0.005. Some of the results are illustrated in Table 3. Note that the value of z(t)’s are at the mid- points t + h/2.

Table 3. Result values for x (t), y (t) and z (t) from Algorithm A, and the errors.

5. Conclusion

The methods are constructed with the objective of solving DAE systems in their original forms. Both algorithms use the Runge-Kutta method as the integrator, and couple this with a method to solve the algebraic systems associated in the problem. The Newton method is used in Algorithm A for index-1 DAE, while in Algorithm B the Newton-Broyden method is needed for an index-2 DAE system. The methods can give approximate solutions for the problem very well, with only small errors. Experiments have also shown that high index DAEs are harder to solve than lower index DAEs.

Acknowledgements

The author would like to thank the referees for their comments. This work was financially supported by the Research Grant of Burapha University through National Research Council of Thailand (Grant No. 35/2556).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Ascher, U.M. and Petzold, L.R. (1998) Computer Method for Ordinary Differential Equations and Differential-Algebraic Equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia. http://dx.doi.org/10.1137/1.9781611971392
[2]	Brenan, K.E., Campbell, S.L. and Petzold, L.R. (1996) Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, Revised and Corrected Reprint of 1989 Original, with an Additional Chapter and Additional References. Classic in Applied Mathematics, 14. SIAM, Philadelphia.
[3]	Campbell, S.L., et al. (2008) Differential-Algebraic Equations. Scholarpedia, 3, 2849. http://dx.doi.org/10.4249/scholarpedia.2849
[4]	Rabier, P.J. and Rheinboldt, W.C. (2002) Theoretical and Numerical Analysis of Differential-Algebraic Equations. Hand Book of Numerical Analysis, 8, 183-540. http://dx.doi.org/10.1016/S1570-8659(02)08004-3
[5]	Riaza, R. (2008) Differential-Algebraic System: Analytical Espects and Circuit Applications. World Scientific, Singapore City.
[6]	Hairer, E. and Wanner, G. (1998) Solving Ordinary Differential Equation. II. Stiff and Differential-Algebraic Problems, 2nd Edition, Springer Series in Computational Mathematics, 14. Springer-Verlag, Berlin.
[7]	Chompuvised, K. and Dhamacharoen, A. (2011) Solving Boundary Value Problems of Ordinary Differential Equations with Non-Separated Boundary Conditions. Applied Mathematics and Computation, 217, 10355-10360. http://dx.doi.org/10.1016/j.amc.2011.05.044
[8]	Dhamacharoen, A. (2014) An Efficient Hybrid Method for Solving Systems of Nonlinear Equations. Journal of Computational and Applied Mathematics, 263, 59-68. http://dx.doi.org/10.1016/j.cam.2013.12.006
[9]	Dhamacharoen, A. and Chompuvised, K. (2013) An Efficient Method For Solving Multipoint Equation Boundary Value Problems. World Academy of Science, Engineering and Technology, 7, 329-333.
[10]	Differential-Algebraic Equation, Wikipedia, the Free Encyclopedia.