Numerical Methods For Solving Differential Algebraic

10Chapter 2Theory of DAEs2.1 IntroductionTheory of DAEs has been developed to meet the needs of certainsolution to the systems discussed in section 1.4. In this chapter, maindefinitions and theorems are stated in order to facilitate dealing withDAEs, and to understand the differences between DAEs and ODEs.2.2 SolvabilityOne of the important definitions is solvability which means, insimple words, the existence and uniqueness of solution on the intervalof interest.Definitionis a connected open subset of R 2m 1 ,Let I be an open subinterval of R ,and F is a differentiable function fromon I into R m . Then the DAE (1.1) is solvableif there is an r - dimensional family of solutions (t,c) defined on aconnected open set I ,R r , such that :1) (t,c) is defined on all of I for each c2) ( t, (t,c), (t,c))for (t,c)I

113) Ifthen(t) is any other solution w ith ( t, (t, c),(t, c)),(t) (t, c)4) T he graph ofas a function of (t, c) is an r 1 - dim ensionalm anifold.The definition says that locally there is an r-dimensional familyof solutions. At any time t 0I , the initial conditions form an r-dimensional manifold (t 0 ,c ) and r is independent of t 0 . The solutionsare continuous function of the initial conditions on this manifold (orequivalently of c). Since the solutions exist for all tI and (3) holds,there are no bifurcations of these solutions. [6]2.3 IndexAs mentioned earlier in Chapter One, DAEs are distinguished bydifficulty. One way to simplify dealing with DAEs is to convert it as asystem of ODEs, so the measure that will play a key role in theclassification and behavior of DAEs is called an index. In order tomotivate the definition of index, let us consider the special case of a semiexplicit non-linear DAE (1.9)x 1 (t ) f 1 (t , x 1 (t ), x 2 (t ))0 f 2 (t , x 1 (t ), x 2 (t ))If we differentiate the constraint equation f 2 (t , x 1 (t ), x 2 (t )) 0with respect to t , we get an implicit ODE of the form

12f 2x 1 (t , x 1 (t ), x 2 (t ))x 1f 2x 2 (t , x 1 (t ), x 2 (t ))x 2f 2 t (t , x 1 (t ), x 2 (t )). (2.1)Substitute x 1 (t ) f 1 (t , x 1 (t ), x 2 (t )) in (2.1), we havef 2x 1 (t , x 1 (t ), x 2 (t ))f 1 (t , x 1 (t ), x 2 (t )) f 2x 2 (t , x 1 (t ), x 2 (t ))x 2f 2 t (t , x 1 (t ), x 2 (t )).that isf 2x 2 (t , x 1 (t ), x 2 (t ))x 2f 2 t (t , x 1 (t ), x 2 (t )) f 2x 1 (t , x 1 (t ), x 2 (t ))f 1 (t , x 1 (t ), x 2 (t )).If f 2 is nonsingular, then the following systemx2x 1 (t ) f 1 (t , x 1 (t ), x 2 (t ))f 2x 2 (t , x 1 (t ), x 2 (t ))x 2f 2 t (t , x 1 (t ), x 2 (t )) f 2x 1 (t , x 1 (t ), x 2 (t ))f 1 (t , x 1 (t ), x 2 (t )).is an implicit ODEs and we say that it has index one. If this is not thecase, suppose that with algebraic manipulation and coordinate changeswe can rewrite it in the form of (1.9) but with different x ' s . Again, wedifferentiate the constraint equation. If an implicit ODE results, we saythat the original problem has index two. If the new system is not animplicit ODE, we repeat the process. The number of differentiationsteps required in this procedure is the index.[6]Definition: The minimum number of times that all or part of (1.1)must be differentiated with respect to t in order to determine x as acontinuous function of x, t, is the index of DAE.

13Example 1 : take the following DAEsy20y 1 f 1 (t )y2f 2 (t )The first equation is ODE, but the seond is a constraint, so differentiateit will yieldy2f 2 (t )using ODE we havey 1 f 1 (t )y1f 2 (t )f 2 (t ) f 1 (t )differentiate this equation will givey1f 2 (t ) f 1 (t )so, after two times of differentiation we have the following ODE systemy2f 2 (t )y1f 2 (t ) f 1 (t ).so, the original DAEs is of index two.

14y20y 1 f 1 (t )y2y2y 1 f 1 (t )[7]f 2 (t )y1DAEsf 2 (t ) f 1 (t )ODEsExample 2 : Let us discuss the simple pendulum problem:xxyyx20gy2L2This system can be transformed as a system of first order as follows:x1x,x2x ,x3y,x4y ,x1x 2,x2x3x 1,x 4,x4x 12x 32x3g,L20.If we want to determine the index for this system, we mustdifferentiate the constraint until we have system of ODEs by taking thefollowing steps:

15Step 1Takex 12x 32L20,differentiate implicitly2x 1x 12x 3x 30,that isx 1x 2x 3x 40,(1)Step 2Differentiate equation (1) implicitly, yieldsx1 x 2x 1x 2substitute x 1x 22x 12x3 x4x3 x4x 2, x 20,x 1, x 3x 42x 32gx 30,(x 12x 32 ) gx 30,x 4, x 4rearrangex 22x 42Step 3Differentiate equation (2)(2)x3g , yield

162x 2 x 22x 4 x 4(2x 1x 12x 3x 3 )(x 12x 32 ) gx 30,substitute2 x 2 x 1 2 x 4 x 3 2 gx 42 x 1x 22 x 3x 4Lgx 40,rearrange3gx 44 x 3x 4L. [7]We have transform DAEs for simple pendulum problem of the formxxyy0x2gy2L2into the following ODE systemx1x2x3x4x2x1x 4,x3g,3gx 44 x 3x 4LThis ODE system resulted after three differentiation times, soDAEs is of index three. The original system consists of three equations,but the new ODE system consists of five equations with two newvariables.

17Although, we can sometimes transform DAEs into ODEs, but there aredifficulties that may appear as disadvantages for this process:1)The exact solution for ODEs is not necessarily a solution for theoriginal DAEs.2)It is not trivial to choose the right initial values that satisfy theoriginal system. [2]The concept of index has been introduced in order to quantifythe level of difficulty that is involved in solving a given DAEs. So,higher index, greater than one, means that we have more difficulty tosolve DAEs. By this concept, ODEs can be classified as 0- index, whichagrees with the simplicity of numerical methods that deal directly withODEs.2.4 Theory of linear DAEsLinear DAEs in both types, constant or varying, are bestunderstood class of DAEs systems. In this section, main definitions andtheorems are needed to understand the numerical methods in thefollowing chapter.2.4.1 Matrix pencilMatrix pencil is magnitude coupled with linear DAEs of the formA n n .x n1B n n .x n1C n 1.Definition:MatrixAB,pencilisthemagnitudebe a complex parameter. If determinant of this magnitude is

18not equal to zero then the pencil is said to be regular. The importanceof matrix pencil comes from the existence of solution for linear DAEs,this is justified in the following derivation:Take the linear DAE A n n .x nlet x replaced byA.xn1xnxnhB .x nhxn11B n n .x n1Cn1, then we haveC1multiply by h to getA .x n1hB .x nfactorize x n1(A1hB ) x nA .x n1hCmultibly by ( Axn1(AhB )A .x nhB )1n nhC1n n.( hCThis means that AA x n )n1hB must be invertible, i.e. must be nonsingular.The following theorem will relate solvability with matrix pencil:Theorem: The linear DAE is solvable if and only ifpencil.Examples :1)LetAB is a regular

190 1 0 01 0 0 00 0 1 00 1 0 0x (t) x (t) f (t)0 0 0 10 0 1 00 0 0 00 0 0 1This system has a regular matrix pencil because det( AB ) 1.3 ) let1 01 1x (t) x (t ) f (t )0 00 0has singular matrix pencil because det( AB ) 0.2.4.2 Nilpotency :Other measure of DAE's difficulty which is nilpotency.Definition:Let M be a square matrix, k23, if M, M ,M , , Mk 10,and M k 0 then M is said to be nilpotent, k is said to be the nilpotency.Example:Let M M20 0,-1 00 0.0 0So we can say that M is nilpotent with nilpotency k 2Theorem: Suppose that A B is a regular pencil, A and B are square

20matrices of order n. Then there exists a nonsingular matrices P, Qsuch thatPAQ I00 N, PBQ C00Iwhere N is a matrix of nilpotency k and I is the identity matrix.The degree of nilpotency of N in this theorem, namely the integer kisthe index of pencil AB . It is also the index of the DAE. [5]2.4.3 Special DAE formsDefinition :The DAEs (1.1) is in Hessenberg form of size r if it canbe written asI00I0x1I00B 11B 1,r1B 1rB 210B 2,r10B r ,r10xrx1f1xrfrwhere x i are vectors, B ij are matrices, and the product B r ,r 1B r1, r 2B 1rare nonsingular.The Hessenberg form of size two and three are the mostcommon. The Hessenberg form of size two isx1BB 12 110 B 210 0 x21 0det(B ) 0.x1x2f1f2

21The Hessenberg form of size three is1 0 0 x10 1 0 x20 0 0 x3B 11 B 12B 21 B 220 B 32B 1300x1x2x3f1f2f3det(B ) 0.2.5 Non-linear Hessenberg formThe general DAE system (1.1) can include problems which arenot well-defined in mathematical sense, as well as problems which willresult in failure for any direct discretization method. Fortunately,many of the higher-index problems can be expressed as a combinationof more restrictive structures of ODEs coupled with contraints. One ofthe more important classes of systems are the Hessenberg forms.Definition The DAEs (1.1 ) is in non-linear Hessenberg form of size rif it is writtenx1F1 (x 1 , x 2 ,., x r ,t )x2F2 (x 1 , x 2 ,., x r 1 ,t )xiFi (x i 1 , x i ,., x r 1 ,t ), 3ir-10 Fr (x r 1 ,t )and ( Fr / x r 1 )( Fr 1 / x r 2 )( F2 / x 1 )( F1 / x r ) is nonsingular.The Hessenberg form of size two is

22x 1 F1 (x 1, x 2 ,t )0F2 (x1,t)with ( F2 / x1 )( F1 / x 2 ) nonsingular.The Hessenberg form of size 3 isx 1 F1 (x 1, x 2 , x 3 ,t )x 2 F2 (x 1, x 2 ,t )0F3 (x 2 ,t)with ( F3 / x 2 ) ( F2 / x1 )( F1 / x3 ) nonsingular.2.5.1 Main Properties for these forms1)Assuming the Fi is sufficiently differentiable, the Hessenbergform of size r is solvable and has index r .2)Many applications appear as versions of the Hessenberg formssuch as trajectory control problem of size two and three andsome beam3)deflection problems are of size four[6].They are easy to manipulate numerically and analytically.2.6 The importance of DAEsThere are several reasons to consider DAEs directly, thesereasons can be summarized as follows:

231)Each variable in DAEs is significant that represents an elementof physical problem simulated. So, changing the original modelto ODEs may produce less meaningful variables.2)Sometimes it is impossible or time consuming to obtain anexplicit or normal ODEs.3) It is easier for the scientist and engineer to explore the effect ofmodeling changes and parameter variation from original DAEs.[6]These reasons were manifested in the simple pendulum problemfrom the previous chapter. There were two new variables added to theproblem in order to get ODE system. Beside this, the original DAEsconsists of three equations, whereas, the derived ODE system consistsof five equations. In addition, the algebraic manipulations used in thependulum problem, are stuffy and complicated.2.7 Review of literatureSeveral numerical methods have been developed for solvingDAEs. However the pioneering numerical methods in the field ofDAEs could be presented here:The code DASSL (Differential Algebraic System SoLver) inFortran designed by Petzold[6] based on backward difference formula(BDF) to solve general index zero and one DAEs. It will be discussedin detail in the next chapter.Radau collocation which based on implicit Runge-Kutta method.It solves DAEs of the form Mxf (t , x ) , where M is a constant,

24square matrix which may be singular. The code is applicable toproblems of index 1,2,3. The higher-index variables must be identifiedby the user. [4]DAETS (Differential Algebraic Equations by Taylor Series) thatsolves initial value problems for (DAEs) using Taylor series expansion.DAETS designed by Nedialkov and Pryce.[10]Power series method, which was suggested by Bayram[3], itdeals with the DAEs directly by constructing a polynomial toapproximate the solution. It will be discussed in detailed in the nextchapter.There are instant solvers for specified DAEs problems which canbe available on Mathwork website. For example, ode15s.m, hb1dae.mand dae4.m.

25Chapter 3Numerical MethodsIn this chapter, we want to use and explore some numericalmethods for solving differential algebraic equations (DAEs).3.1 Backward difference formula (BDF) Method:The first general technique for the numerical solution of DAEs ,proposed in 1971 by C. W. Gear. This method was initially defined forsystems of differential equations of the formF (t,x(t),x (t) ) 0G (t,x(t) ) 0x (t 0 ) x 0 , x (t 0 ) x1 .This method was soon extended to apply to any fully-implicitDAEsF (t , x , x ) 0. It is considered the most popular multistepmethod.The simplest first order BDF method is the implicit Euler method,which consists of replacing the derivative x' by a backward differenceF (t n , x n ,where hxntnx n-1) 0ht n 1.The resulting system of nonlinear equations for x n at each timestep is then usually solved by Newton's method. The k-step (constant

26stepsize ) BDF consists of replacing x' by the derivative of thepolynomial which interpolates the computed solution at k 1 timest n , t n 1 ,., t n k , evaluated at t n . This yieldsF (t n , x n ,yn) 0,hwhere y nki 0iyniandi,i1,., k are the coefficients of theBDF method.[6]3.1.1 Software for BDFThe most widely used production code for DAEs that based onBDF is the code DASSL (Differential Algebraic System Solver) whichdesigned by Petzold. It can be used for the solution of DAEs of indexzero and one.DASSL uses a variable stepsize variable order fixed leading coefficientimplementation of BDF formula to advance the solution from one timestep to the next. The fixed leading coefficient implementation is oneway of extending the fixed stepsize BDF methods to variable stepsize.There are three main approaches to extending fixed stepsize multi stepmethod to variable stepsize. These formulations are called fixedcoefficient, variable coefficient and fixed leading coefficient. The fixedcoefficient method have the property that they can be implementedvery efficiently for smooth problems, but suffer from inefficiency, orpossible instability, for problems which require stepsize adjustments.The variable coefficient methods are the most stable implementation,but have the disadvantage that they have tend to require more

27evaluation of the Jacobian matrix in intervals when the stepsize ischanging, and hence are usually considered to be less efficient than afixed coefficient implementation for most problems. The fixed leadingcoefficient formulation is a compromise between the fixed coefficientand variable coefficient approaches, offering somewhat less stability,along with usually fewer Jacobian evaluation, than the variablecoefficient formulation. [6]3.1.2 Basic formula used in DASSLSuppose we have approximation y nforiito the true solution y (t n i )0,1,., k where k is the order of the BDF method that we arecurrently planning to use. We would like to find an approximation tothe solution at time t n 1 . First, an initial guess for the solution and itsderivative at t n 1 is formed by evaluating the predictor polynomial andthe derivative of the predictor polynomial at t n 1 . The predictorpolynomialk 1 times,Pn 1Pn 1(t ) is the polynomial which interpolates y n i at the last(t n i )yn i, i0,1,., k .The predicted values for y and y' at t n 1 are obtained by evaluatingPn 1(t ) andPn 1(t ) at t n 1 , so thaty n(0)1Pn 1y n(0)1Pn 1(t n 1 )(t n 1 ).The approximation y n 1 to the solution at t n 1 which is finallyaccepted by DASSL is the solution to the corrector formula. The

28formula used is the fixed leading coefficient form of the k th order BDFmethod. The solution to the corrector formula is the vector y nthat the corrector polynomialCn 11such(t ) and its derivative satisfy the DAEat t n 1 , and the corrector polynomial interpolates the predictorpolynomial at k equally spaced points behind t n 1 ,Cn 1(t n 1 )Cn 1(t nF (t n 1 ,yn1ihn 1 )1Cn 1(t n 1 ),Cn 1Cn 1(t n1ihn 1 ), 1 ik,(t n 1 )) 0.The values of the predictor y n(0)1 , y n(0)1 and the corrector y n1att n 1 are defined in terms of polynomials which interpolate the solutionat the previous time steps. These polynomials are represented by inDASSL in terms of modified divided differences of y.The divided difference are defined by the recurrence[y n ] y n[ y n , y n 1 ,., y n k ][ y n , y n 1 ,., y n] [ y n 1 , y n 2 ,., y n k ].tn tn kk 1The predictor polynomial is given in terms of the divided difference byPn 1(t(t )ynt n )(t(tt n 1)t n )[ y n , y n 1 ] (t(ttnk 1t n )(t)[ y n , y n 1 ,t n 1 )[ y n , y n 1 , y n 2 ] ., y n k ].

29The corrector polynomial can be formulated as followsCn 1Pn 1(t )(t ) b (t )( y ny n(0)1),1whereb (t n1ihn 1) 0, i 1,2, , kb (t n 1) 1.by differentiating at t n 1 givess(y n1y n(0)1) hn 1( y n1y n(0)1 ) 0Solving for y n 1 , we find that the corrector iteration must solveF (t n 1 , y n 1 , y n(0)1shn(y n1y n(0)1 )) 0.1This equation must be solved at each time step using modifiedNewton's iteration by simplifying the equation as

30F (t , y , y) 0.wheresy n(0)11y n(0)1.kand/ hnsj11 jso the iterative formula will bey( m1)y ( m ) G 1F (t , y ( m ) , y ( m )),with y(0) can be taken from pridector polynomial,and GFyF.[6]y3.1.3 AlgorithmThis section presents a new algorithm that facilitates dealingwith linear DAEs, this algorithm based on BDF with fixed leadingcoefficient, fixed step size, variable order polynomial.BDF method for solving linear DAEsInput : system of ODEs with initial conditionsOutput : approximated valuesStep 1 evaluate starting values using Runge Kutta 4th order method forsystem

31Step 2 repeat step 3-5 until end of intervalStep 3 interpolate the previous k points for each yStep 4 substitute tn 1 in the polynomial of degree k-1, 2 k6Step 5 check consistencyStep 6 save new yn 1See flowchart at Appendix A3.2 Power Series MethodAssume that differential algebraic equations (DAEs) has theform (1.1) with initial values y (x 0 )y 0 , y (x 0 )y 1 which areconsistent , i.e.F ( y 0 , y 1, x 0 ) 0G (y 0,x 0) 0The solution of (1.1) can be assumed that vy0y 1xex 2 (3.1)where e is a vector function which is the same size as y 0 and y 1 .Substitute (3.1) into (1.1) will yield vector function which is equalto zero, then choose appropriate terms in order to produce linearequations of e in which can be written in the form Se T , where S isconstant square matrix and T is constant vector, solve this system fore, return to (3.1) with new value of e. Add new term to y and repeatthe above steps until achieve specified number of iterations. Theresultingpower series can be transf

methods for solving linear Dierential-algebraic equations (DAEs): Power series and Backward Difference Formula (BDF). It aims to facilitate dealing with DAEs by analyzing these two numerical methods, designing simple algorithms, and using MATLAB programs