# On the Numerical Solution of Differential-Algebraic Equations with Hessenberg Index-3

Authors: Melike Karta; Ercan Çelik
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

Numerical solution of differential-algebraic equations with Hessenberg index-3 is considered by variational iteration method. We applied this method to two examples, and solutions have been compared with those obtained by exact solutions.


---

## Body

## 1. Introduction

Many important mathematical models can be expressed in terms of differential-algebraic equations (DAEs). Many physical problems are most easily initially modeled as a system of differential-algebraic equations (DAEs) [1]. Some numerical methods have been developed, using both BDF [1–3] and implicit Runge-Kutta methods [1], Padé and Chebysev approximations method [4–6]. These methods are only directly suitable for low-index problems and often require that the problem, have special structure. Although many important applications can be solved by these methods, there is a need for more general approaches. There are many new publication in the field of analytical sueveys such as [7–10]. The variational iteration method (VIM) was developed by He in [11]. The method is used by many researchers in a variety of scientific fields. The method has been proved by many authors [12–16] to be reliable and efficient for a variety of scientific applications, linear and nonlinear as well.The most general form of a DAE is given by(1.1)F(t,x,x′)=0,
where ∂F/∂x′ may be singular. The rank and structure of this Jacobian matrix may depend, in general, on the solution x(t), and for simplicity we will always assume that it is independent of t. The important special case is of a semiexplicit DAE or an ODE with constraints:(1.2a)x′=f(t,x,z),(1.2b)0=g(t,x,z).This is a special case of (1.1). The index is 1 if ∂g/∂z is nonsingular, because then one differentiation of (1.2b) yields z′ in principle. For the semi-explicit index-1 DAE we can distinguish between differential variables x(t) and algebraic variables z(t) [1]. The algebraic variables may be less smooth than the differential variables by one derivative. In the general case, each component of x may contain a mix of differential and algebraic components, which makes the numerical solution of such high-index problems much harder and riskier.


## 2. Special Differential-Algebraic Equations (DAEs) Forms

Most of the higher-index problems encountered in practice can be expressed as a combination of more restrictive structures of ODEs coupled with constraints. In such systems the algebraic and differential variables are explicitly identified for higher-index DAEs as well, and the algebraic variables may all be eliminated using the same number of differentiations. These are called Hessenberg forms of the DAE and are given below. In this paper, the variational iteration method has been proposed for solving differential-algebraic equations with Hessenberg index-3.


### 2.1. Hessenberg Index-1

One has(2.1)x′=f(t,x,z),0=g(t,x,z).
Here the Jacobian matrix function gz is assumed to be nonsingular for all t. This is also often referred to as a semi-explicit index-1 system. Semi-explicit index-1 DAEs are very closely related to implicit ODEs.


### 2.2. Hessenberg Index-2

One has(2.2a)x′=f(t,x,z),(2.2b)0=g(t,x).Here the product of Jacobians gxfz is nonsingular for all t. Note the absence of the algebraic variables z from the constraints (2.2b). This is a pure index-2 DAE, and all algebraic variables play the role of index-2 variables.


### 2.3. Hessenberg Index-3

One has(2.3)x′=f(t,x,y,z),y′=g(t,x,y),0=h(t,y).
Here the product of three matrix functions hygxfz is nonsingular.The index of a Hessenberg DAE is found, as in the general case, by differentiation. However, here only algebraic constraints must be differentiated.


## 2.1. Hessenberg Index-1

One has(2.1)x′=f(t,x,z),0=g(t,x,z).
Here the Jacobian matrix function gz is assumed to be nonsingular for all t. This is also often referred to as a semi-explicit index-1 system. Semi-explicit index-1 DAEs are very closely related to implicit ODEs.


## 2.2. Hessenberg Index-2

One has(2.2a)x′=f(t,x,z),(2.2b)0=g(t,x).Here the product of Jacobians gxfz is nonsingular for all t. Note the absence of the algebraic variables z from the constraints (2.2b). This is a pure index-2 DAE, and all algebraic variables play the role of index-2 variables.


## 2.3. Hessenberg Index-3

One has(2.3)x′=f(t,x,y,z),y′=g(t,x,y),0=h(t,y).
Here the product of three matrix functions hygxfz is nonsingular.The index of a Hessenberg DAE is found, as in the general case, by differentiation. However, here only algebraic constraints must be differentiated.


## 3. He’s Variational Iteration Method (VIM)

Consider the differential equation(3.1)Lu+Nu=g(x),
where L and N are linear and nonlinear operators, respectively, and g(x) is the source inhomogeneous term. In [11], He proposed the variational iteration method where a correction functional for (3.1) can be written as(3.2)un+1(x)=un(x)+∫0xλ(t)(Lun(t)+Nũn(t)-g(t))dt,
where λ is a general Lagrange’s multiplier, which can be identified optimally via the variational theory and ũn as a restricted variation which means δũn=0. It is to be noted that the Lagrange multiplier λ can be a constant or a function.The variational iteration method should be employed by following two essential steps. It is required first to determine the Lagrange multiplierλ that can be identified optimally via integration by parts and by using a restricted variation. Having λ determined, an iteration formula, without restricted variation, should be used for the determination of the successive approximations un+1(x), n≥0, of the solution u(x). The zeroth approximation u0 can be any selective function. However, using the initial values u(0), u′(0), and u′′(0) are preferably used for the selective zeroth approximation u0 as will be seen later. Consequently, the solution is given by(3.3)un(x)=lim⁡n→∞un(x).


### 3.1. First-Order ODEs

We first start our analysis by studying the first-order linear ODE of a standard form(3.4)u′+p(x)u=q(x),u(0)=α.
The VIM admits the use of the correction functional for this equation by(3.5)un+1(x)=un(x)+∫0xλ(t)(un′(t)+p(t)ũn(t)-g(t))dt,
where λ is the lagrange multiplier, that in this method may be a constant or a function, and ũn is a restricted value where δũn=0.Taking the variation of both sides of (3.5) with respect to the independent variable un we have(3.6)δun+1(x)=δun(x)+δ(∫0xλ(t)(un′(t)+p(t)ũn(t)-q(t)dt))
that gives(3.7)δun+1(x)=δun(t)+δ(∫0xλ(t)un′(t)dt)
obtained upon using δũn=0 and δq(t)=0. Integrating the integral of (3.6) by parts we obtain(3.8)δun+1=δun+δλun(x)-δ∫0xλ′undt
or equivalently(3.9)δun+1=δ(1+λ|t=x)un-δ∫0xλ′undt.
The extremum condition of un+1 requires that un+1=0. This means that the left hand side of (3.9) is 0, and as a result the right hand side should be 0 as well. This yields the stationary conditions(3.10)1+λt=x=0,λ′|t=x=0.
This in turn gives(3.11)λ=-1.
Substituting this value of the lagrange multiplier into the functional (3.5) gives the iteration formula
(3.12)un+1(x)=un(x)-∫0x(un′(t)+p(t)un(t)-q(t))dt,
obtained upon deleting the restriction on un that was used for the determination of λ. Considering the given condition u(0)=α, we can select the zeroth approximation uo=α. Using the selection into (3.12) we obtain the following successive approximations:(3.13)u0(t)=α,u1(x)=α-∫0x(u0′(t)+p(t)u0(t)-q(t))dt,u2(x)=u1(x)-∫0x(u1′(t)+p(t)u1(t)-q(t))dt,u3(x)=u2(x)-∫0x(u2′(t)+p(t)u2(t)-q(t))dt,⋮un+1=un(x)-∫0x(un′(t)+p(t)un(t)-q(t))dt.
Recall that(3.14)u(x)=lim⁡n→∞un+1(x),
that may give the exact solution if a closed form solution exists, or we can use the (n+1)th approximation for numerical purposes.


## 3.1. First-Order ODEs

We first start our analysis by studying the first-order linear ODE of a standard form(3.4)u′+p(x)u=q(x),u(0)=α.
The VIM admits the use of the correction functional for this equation by(3.5)un+1(x)=un(x)+∫0xλ(t)(un′(t)+p(t)ũn(t)-g(t))dt,
where λ is the lagrange multiplier, that in this method may be a constant or a function, and ũn is a restricted value where δũn=0.Taking the variation of both sides of (3.5) with respect to the independent variable un we have(3.6)δun+1(x)=δun(x)+δ(∫0xλ(t)(un′(t)+p(t)ũn(t)-q(t)dt))
that gives(3.7)δun+1(x)=δun(t)+δ(∫0xλ(t)un′(t)dt)
obtained upon using δũn=0 and δq(t)=0. Integrating the integral of (3.6) by parts we obtain(3.8)δun+1=δun+δλun(x)-δ∫0xλ′undt
or equivalently(3.9)δun+1=δ(1+λ|t=x)un-δ∫0xλ′undt.
The extremum condition of un+1 requires that un+1=0. This means that the left hand side of (3.9) is 0, and as a result the right hand side should be 0 as well. This yields the stationary conditions(3.10)1+λt=x=0,λ′|t=x=0.
This in turn gives(3.11)λ=-1.
Substituting this value of the lagrange multiplier into the functional (3.5) gives the iteration formula
(3.12)un+1(x)=un(x)-∫0x(un′(t)+p(t)un(t)-q(t))dt,
obtained upon deleting the restriction on un that was used for the determination of λ. Considering the given condition u(0)=α, we can select the zeroth approximation uo=α. Using the selection into (3.12) we obtain the following successive approximations:(3.13)u0(t)=α,u1(x)=α-∫0x(u0′(t)+p(t)u0(t)-q(t))dt,u2(x)=u1(x)-∫0x(u1′(t)+p(t)u1(t)-q(t))dt,u3(x)=u2(x)-∫0x(u2′(t)+p(t)u2(t)-q(t))dt,⋮un+1=un(x)-∫0x(un′(t)+p(t)un(t)-q(t))dt.
Recall that(3.14)u(x)=lim⁡n→∞un+1(x),
that may give the exact solution if a closed form solution exists, or we can use the (n+1)th approximation for numerical purposes.


## 4. Applications

Example 4.1.
We first considered the following differential-algebraic equations with Hessenberg index-3 form:(4.1)x2′+x1-1=0,xx2′+x3′+2x2-2x=0,xx2+x3-ex=0
with initial conditions
(4.2)(x1(0)x2(0)x3(0))=(0-11).
The exact solutions are
(4.3)x1(x)=ex-1,x2(x)=2x-ex,x3(x)=(1+x)ex-2x2,
where x2, x3 represent the differential variables and x1 represents the algebraic variables. After three times of differentiation of (4.1) we have the following ODE system:
(4.4)x1′=ex,x2′=2-ex,x3′=xx1-2x2+x.
Differential-algebraic equation (DAE) is a Hessenberg index-3 form.
To solve system (4.4), we can construct following correction functionals:
(4.5)x1(n+1)(x)=x1(n)(x)+∫0xλ1(t)(x1′(n)(t)-et)dt,x2(n+1)(x)=x2(n)(x)+∫0xλ2(t)(x2′(n)(t)-2+et)dt,x3(n+1)(x)=x3(n)(x)+∫0xλ3(t)(x3′(n)(t)-tx̃1(n)+2x̃2(n)-t)dt,
where λ1(t), λ2(t), and, λ3(t) are general Lagrange multipliers and x̃1(n), x̃2(n) denote restricted variations, that is, δx̃1(n)=δx̃2(n)=0.
Making the above correct functional stationary,(4.6)δx1(n+1)(x)=δx1(n)(x)+δ∫0xλ1(t)(x1′(n)(t)-et)dt,δx2(n+1)(x)=δx2(n)(x)+δ∫0xλ2(t)(x2′(n)(t)-2+et)dt,δx3(n+1)(x)=δx3(n)(x)+δ∫0xλ3(t)(x3′(n)(t)-tx̃1(n)(t)+2x̃2(n)(t)-t)dt,δx1(n+1)(x)=δx1(n)(x)+δλ1(t)x1(n)(t)|0x-∫0xλ1′(t)δx1(n)(t)dt=0,δx2(n+1)(x)=δx2(n)(x)+δλ2(t)x2(n)(t)|0x-∫0xλ2′(t)δx2(n)(t)dt=0,δx3(n+1)(x)=δx3(n)(x)+δλ3(t)x3(n)(t)|0x-∫0xλ3′(t)δx3(n)(t)dt=0.
Its stationary conditions can be obtained as follows:
(4.7)λ1′(t)=λ2′(t)=λ3′(t)=0,1+λ1(t)|t=x=1+λ2(t)|t=x=1+λ3(t)|t=x=0.
The Lagrange multipliers can be identified as follows:
(4.8)λ1(t)=λ2(t)=λ3(t)=-1,
and the following multipliers can be obtained as
(4.9)x1(n+1)(x)=x1(n)(x)-∫0x(x1′(n)(t)-et)dt,x2(n+1)(x)=x2(n)(x)-∫0x(x2′(n)(t)-2+et)dt,x3(n+1)(x)=x3(n)(x)-∫0x(x3′(n)(t)-tx̃1(n)+2x̃2(n)-t)dt.
Beginning with x1(0)(x)=0, x2(0)(x)=-1, x3(0)(x)=1 by the iteration formula (4.9), we have
(4.10)x1*(x)=x+12x2+16x3+124x4+1120x5+1720x6+15040x7+140320x8-1362880x9+⋯,x2*(x)=-1+x-12x2-16x3-124x4-1120x5-1720x6-15040x7-140320x8-1362880x9-⋯,x3*(x)=1+2x-12x2+23x3+524x4+120x5+7720x6+1630x7+14480x8+136288x9+⋯.Example 4.2.
One has(4.11)(100010000)(v1′v2′v3′)+(11xexx+100x20)(v1v2v3)=(2xx2+x+2x3)
with initial conditions
(4.12)x(0)=(101).
The exact solutions are v1(x)=e-x, v2(x)=x, v3(x)=1:
(4.13)v1′+v1+v2+xv3-2x=0,v2′+exv1+(x+1)v2-x2-x-2=0,x2v2-x3=0,
where v1, v2 represent the differential variables and v3 represents the algebraic variables. After three times of differentiation of (4.11) we have the following ODE system:
(4.14)v1′=-v1-v2-xv3+2x,v2′=-exv1-(x+1)v2+x2+x+2,v3′=0.
Differential-algebraic equation (DAE) is a Hessenberg index-3 form.
To solve system (4.14), we can construct the following correction functionals:
(4.15)v1(n+1)(x)=v1(n)(x)+∫0xλ1(t)(v1′(n)(t)+ṽ1(n)(t)+ṽ2(n)(t)+tṽ3()(t)-2t)dt,v2(n+1)(x)=v2(n)(x)+∫0xλ2(t)(v2′(n)(t)+etv1(n)(t)+(t+1)ṽ2(n)(t)-t2-t-2)dt,v3(n+1)(x)=v3(n)(x)+∫0xλ3(t)v2(n)(t)dt.
By using the basic definition of the variational iteration method can obtain that
(4.16)v1*=1-x+12x2-16x3+124x4-1120x5+1720x6-15040x7+140320x8-1362880x9+⋯,v2*=x,v3*=1.


## 5. Conclusion

The method has been proposed for solving differential-algebraic equations with Hessenberg index-3. Results show the advantages of the method. Tables1–4 and Figures 1–4 show that the numerical solution approximates the exact solution very well in accordance with the above method.Table 1
Numerical solution ofx1(x).
xx1(x)x1*(x)|x1(x)-x1*(x)|0.00000.10.1051709180.10517091810.1×10-90.20.2214027580.22140275820.2×10-90.30.3498588080.34985880760.4×10-90.40.4918246980.49182469770.3×10-90.50.6487212710.64872127050.5×10-90.60.8221188000.82211879870.13×10-80.71.0137527071.0137526990.8×10-80.81.2255409281.2255408970.31×10-70.91.4596031111.4596030060.105×10-61.01.7182818281.7182815270.301×10-6Table 2
Numerical solution ofx2(x).
xx2(x)x2*(x)|x2(x)-x2*(x)|0.0−1−100.1−0.905170918−0.90517091810.1×10-90.2−0.821402758−0.82140275820.2×10-90.3−0.749858808−0.74985880760.4×10-90.4−0.691824698−0.69182469770.3×10-90.5−0.648721271−0.64872127050.5×10-90.6−0.622118800−0.62211879870.13×10-80.7−0.613752707−0.61375269920.78×10-80.8−0.625540928−0.62554089670.313×10-70.9−0.659603111−0.65960300660.1044×10-61.0−0.718281828−0.71828152560.3024×10-6Table 3
Numerical solution ofx3(x).
xx3(x)x3*(x)|x3(x)-x3*(x)|0.01100.11.1956880101.19568801000.21.3856833101.3856833090.1×10-80.31.5748164501.5748164510.1×10-80.41.7685545771.7685545760.1×10-80.51.973819061.9730819030.3×10-80.62.1953900802.1953900610.19×10-70.72.4433796022.4433795110.91×10-70.82.7259736702.7259733170.353×10-60.93.0532459113.0532447510.1160×10-51.03.4365636563.4365602950.3361×10-5Table 4
Numerical solution ofv1(x).
xv1(x)v1*(x)|v1-v1*(x)|0.01100.10.90483741800.90483741810.1×10-90.20.81873075310.81873075320.1×10-90.30.74081822070.74081822060.1×10-90.40.67032004600.67032004610.1×10-90.50.60653065970.60653065950.2×10-90.60.54881163610.54881163450.16×10-80.70.49658530380.49658529660.72×10-80.80.49932896410.49932893650.276×10-70.90.40656965970.40656957100.887×10-71.00.36787944120.36787918880.2524×10-6Figure 1
Values ofx1(x) and its x1*(x) variational iteration.Figure 2
Values ofx2(x) and its x2*(x) variational iteration.Figure 3
Values ofx3(x) and its x3*(x) variational iteration.Figure 4
Values ofv1(x) and its v1*(x) variational iteration.



---

*Source: 147240-2012-02-26.xml*