Modified Extended Lie Group Method for Hessenberg Differential Algebraic Equations with Index 3

1. Introduction

Differential-Algebraic Equations (DAEs) are a class of mathematical models that arise in various fields, such as multibody dynamics [1], electrical circuits [2], fluid mechanics [3], optimal control [4] and chemical reactions [5]. DAEs are a natural extension of ordinary differential equations that also involve algebraic equations. As such, they offer a convenient way to model complex systems that involve both dynamic and algebraic components. The study of DAEs has been an active area of research for many years, and numerous numerical methods have been developed to solve these DAEs, particularly for low index systems. These methods include well known Runge-Kutta methods [6,7,8], backward differentiation formulae [9,10],and others. However, the efficient solution of DAEs remains a challenging problem due to the presence of algebraic equations that may be nonlinear, implicit, and of high index.

Hessenberg differential algebraic equations (Hessenberg-DAEs) are an important type of differential equation which are also commonly used in the modeling of mechanical systems and electrical circuits.Because of their unique structure, numerical methods for solving Hessenberg-DAEs can differ significantly from those used for other types of differential equations. Over the years, some approaches have been proposed for the numerical solution of Hessenberg-DAEs, including variational iteration method [11], adomian decomposition method for [12], exponential integrators [13]. These methods either lack rigorous theoretical analysis or are only suitable for low-index systems.

One promising approach is the Lie group method whose several variations have been developed specifically for solving (Hessenberg) DAEs with high index in recent years. Lie group methods mostly constitute one-step extensions of classical methods, including Runge-Kutta type formulas [14]. In mechanical engineering, BDF methods, which have traditionally been used, have been extended to Lie groups as well [15]. Despite their usefulness, these methods suffer from a numerical disadvantage. Various numerical techniques can be used to compute the matrix exponential, but the properties of these methods, such as stability, accuracy, and efficiency, can impact the results of the Lie-group methods. Based on the structural information of Hessenberg DAEs, C.S.Liu [16,17,18] propose a Lie-group differential algebraic equations (LGDAE) method. By converting the canonical representation of the governing equations for nonholonomic systems, the LGDAE method expresses them in the form of DAE systems, comprising of nonlinear ordinary differential equations and nonlinear algebraic equations. The method is composed of two components: the first component pertains to the Lie-group structure and concerns the numerical solution of nonlinear ordinary differential equations. The second component, the Newton iterative scheme, is used to solve the nonlinear algebraic equations. The LGDAE method has been successfully employed to solve various practical problems such as the heat source recovery problem [19,20], sliding mode control problem [21], and inverse nonlinear vibration problems [22]. However, it is noteworthy that the LGDAE method is only applicable to systems with index 2.

The main focus of this paper is to investigate Lie group methods for Hessenberg-DAEs with an index of 3 or higher. We present a modified extended Lie group method that incorporates the benefits of Liu’s original method, along with the structural characteristics of Hessenberg-DAEs. Our proposed method enables the accurate solution of high-index Hessenberg-DAEs and effectively captures the behavior of solutions to algebraic constraints. We demonstrate the effectiveness of our method through numerical experiments on some benchmark problems.

The structure of this paper is as follows. Section 2 presents a brief introduction to Hessenberg-DAEs and Lie group methods. In Section 3, we present our novel modified extended Lie group DAE method for Hessenberg-DAEs with an index of 3,and provide a theoretical justification for the proposed method. Section 4 presents the results of numerical experiments on several benchmark problems. Finally, we conclude with a discussion of our results and directions for future research in Section 5.

2. Preliminaries

In this section, we provide a concise overview of Hessenberg-DAEs and an implicit $GL(n,\mathbb{R})$ Lie group methods.

2.1. Hessenberg DAEs

A system of ordinary differential equations (ODEs) with algebraic constraints which can be written in the general form as

$$\mathbf{f}(t,\mathbf{x}\left(t\right),{\mathbf{x}}^{\prime }\left(t\right))=0$$

(1)

is called differential algebraic equation, where $\frac{\partial \mathbf{f}}{\partial \mathbf{x}}$ may be singular and the structure and rank of the of this Jacobian matrix may depend on the solution $x\left(t\right)$. An important special case of 1 is the semi-explicit DAE of the form

$$\left\{\begin{array}{ccc}\hfill {{\mathbf{x}}_1}^{\prime }\left(t\right)& =& {\mathbf{f}}_1(t,{\mathbf{x}}_1\left(t\right),{\mathbf{x}}_2\left(t\right))\hfill \\ \hfill 0& =& {\mathbf{f}}_2(t,{\mathbf{x}}_1\left(t\right),{\mathbf{x}}_2\left(t\right)),\hfill \end{array}\right.$$

(2)

whose index is 1 if the Jacobian matrix function $\frac{\partial {\mathbf{f}}_2}{\partial {\mathbf{x}}_2}$ is assumed to be nonsingular for all t. Hessenberg Index-1 is also often referred to as a semi-explicit index-1 system.

Definition 1. 

Hessenberg Index-2

$$\left\{\begin{array}{ccc}{{\mathbf{x}}_1}^{\prime }\left(t\right)\hfill & =& {\mathit{f}}_1(t,{\mathit{x}}_1\left(t\right),{\mathit{x}}_2\left(t\right)), \left(3\mathrm{a}\right)\hfill \\ 0\hfill & =& {\mathit{f}}_2(t,{\mathit{x}}_1\left(t\right)),\left(3\mathrm{b}\right)\hfill \end{array}\right.$$

where ${\mathit{f}}_i$ are sufficiently differentiable, the variables ${\mathit{x}}_2$ of the algebraic part is absent from the constraint and the product of Jacobians $\frac{\partial {\mathit{f}}_2}{\partial {\mathit{x}}_1}\frac{\partial {\mathit{f}}_1}{\partial {\mathit{x}}_2}$ is nonsingular. Thus, all algebraic variables play the role of index-2 therefore this is a pure index-2 DAE.

Definition 2. 

Hessenberg Index-3

$$\left\{\begin{array}{ccc}{{\mathit{x}}_1}^{\prime }\left(t\right)\hfill & =& {\mathit{f}}_1(t,{\mathit{x}}_1\left(t\right),{\mathit{x}}_2\left(t\right),{\mathit{x}}_3\left(t\right)),\left(4\mathrm{a}\right)\hfill \\ {{\mathit{x}}_2}^{\prime }\left(t\right)\hfill & =& {\mathit{f}}_2(t,{\mathit{x}}_1\left(t\right),{\mathit{x}}_2\left(t\right)),\left(4\mathrm{b}\right)\hfill \\ 0\hfill & =& {\mathit{f}}_3(t,{\mathit{x}}_2\left(t\right)).\left(4\mathrm{c}\right)\hfill \end{array}\right.$$

In this case the product of three matrix functions $\mathbf{F}=\frac{\partial {\mathbf{f}}_3}{\partial {\mathbf{x}}_2}\frac{\partial {\mathbf{f}}_2}{\partial {\mathbf{x}}_1}\frac{\partial {\mathbf{f}}_1}{\partial {\mathbf{x}}_3}$ is nonsingular, where ${\mathbf{x}}_1$, ${\mathbf{x}}_2$ and ${\mathbf{x}}_3$ are $n_1$, $n_2$ and $n_3$ dimensions respectively. The index of a Hessenberg DAE is found, as in the general case, by differentiation. However, only the algebraic constraints must be differentiated.

Definition 3. 

Hessenberg Index-r

$$\left\{\begin{array}{ccc}{{\mathit{x}}_1}^{\prime }\hfill & =& {\mathit{f}}_1(t,{\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_r),\left(5\mathrm{a}\right)\hfill \\ {{\mathit{x}}_2}^{\prime }\hfill & =& {\mathit{f}}_2(t,{\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_{r-1}),\left(5\mathrm{b}\right)\hfill \\ ⋮\hfill & & \left(5\mathrm{c}\right)\hfill \\ {{\mathit{x}}_{\mathrm{i}}}^{\prime }\hfill & =& {\mathit{f}}_i(t,{\mathit{x}}_{i-1},{\mathit{x}}_i,\cdots ,{\mathit{x}}_{r-1}),3\le i\le r-1,\left(5\mathrm{d}\right)\hfill \\ ⋮\hfill & & \left(5\mathrm{e}\right)\hfill \\ 0\hfill & =& {\mathit{f}}_r(t,{\mathit{x}}_{r-1}).\left(5\mathrm{f}\right)\hfill \end{array}\right.$$

In this case the product of three matrix functions $\frac{\partial {\mathit{f}}_r}{\partial {\mathit{x}}_{r-1}}\frac{\partial {\mathit{f}}_{r-1}}{\partial {\mathit{x}}_{r-2}}\cdots \frac{\partial {\mathit{f}}_2}{\partial {\mathit{x}}_1}\frac{\partial {\mathit{f}}_1}{\partial {\mathit{x}}_r}$ is nonsingular.

2.2. Lie group structure of ODEs

A Lie group is a differentiable manifold equipped with a group structure that is consistent with the underlying topology of the manifold. The general linear group is a Lie group, whose manifold is an open subset $GL(n,\mathbb{R}):=\{\mathbf{G}\in {\mathbb{R}}^{n\times n}|det\mathbf{G}\ne 0\}$. Thus, $GL(n,\mathbb{R})$ is also an $n\times n$-dimensional manifold. The group composition is given by the matrix multiplication.

The general linear group $GL(n,\mathbb{R})$ gives uniquely a real Lie-algebra $gl(n,\mathbb{R})$. Consider a one-parameter subgroup $\mathbf{G}\left(t\right)$, $t\in \mathbb{R}$, of the general linear group $GL(n,\mathbb{R})$, which is a curve passing through the group identity at $t=0$, $\mathbf{G}\left(0\right)={\mathbf{I}}_n$, and which left acts on the n-dimensional Euclidean space ${\mathbb{R}}^n$, resulting in a congruence of curves in ${\mathbb{R}}^n$,

$$\mathbf{X}\left(t\right)=\mathbf{G}\left(t\right)\mathbf{X}\left(0\right)=\left[\mathbf{G}\left(t\right){\mathbf{G}}^{-1}\left(t_0\right)\right]\mathbf{X}\left(t_0\right),t_0,t\in \mathbb{R}.$$

(6)

Owing to the closure property of the Lie group, $\mathbf{G}\left(t\right){\mathbf{G}}^{-1}\left(t_0\right)$ also belongs to $GL(n,\mathbb{R})$. When $t_0$ is put very close to t, $\mathbf{G}\left(t\right){\mathbf{G}}^{-1}\left(t_0\right)$ is very close to the identity ${\mathbf{I}}_n$. Moreover,

$$\mathbf{A}\left(t\right):={\displaystyle \frac{\partial }{\partial t}}\left[\mathbf{G}\left(t\right){\mathbf{G}}^{-1}\left(t_0\right)\right]{|}_{t_0=t}=\dot{\mathbf{G}}\left(t\right){\mathbf{G}}^{-1}\left(t\right)$$

(7)

defines a string of tangent vectors on the tangent space at the group identity of the group manifold, more precisely, a continuously singly parametrized series of one-dimensional sub-algebra of the real Lie algebra $gl(n,\mathbb{R})$.

Differentiating Eq. (6), setting $t_0=t$ and then using Eq. (7) yields

$$\mathbf{X}{\left(t\right)}^{\prime }=\mathbf{A}\left(t\right)\mathbf{X}\left(t\right).$$

(8)

The flow generated by such a $gl(n,\mathbb{R})$ vector field is the congruence of curves resulting from solving the system (8). Essentially, it is very important that when one wants to establish a local coordinate on a smooth manifold on which a finite dimensional Lie-group of transformations is acting, the solution of the differential equations system is necessary.

The Lie-group method possesses a greater advantage than other numerical methods due to its Lie-group structure, and it is a powerful technique to solve ordinary differential equations (ODEs). Here we give a new form of the dynamics in Eq. (4a) or (4b) from the $GL(n,\mathbb{R})$ Lie-group structure. In order to fit the format in Eq. (8), the vector field ${\mathbf{f}}_i$ on the right-hand side of Eq. (4a) or (4b) can be written as

$$\left\{\begin{array}{ccc}{\mathbf{x}}_1^{\prime }\hfill & =& {\mathbf{A}}_1(t,{\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3){\mathbf{x}}_1 \left(9\mathrm{a}\right)\hfill \\ {\mathbf{x}}_2^{\prime }\hfill & =& {\mathbf{A}}_2(t,{\mathbf{x}}_1,{\mathbf{x}}_2){\mathbf{x}}_2\left(9\mathrm{b}\right)\hfill \end{array}\right.$$

where ${\mathbf{A}}_i={\displaystyle \frac{{\mathbf{f}}_i}{\left|\right|{\mathbf{x}}_i\left|\right|}}\otimes {\displaystyle \frac{{\mathbf{x}}_i}{\left|\right|{\mathbf{x}}_i\left|\right|}}$ is the coefficient matrix and $\mathbf{a}\otimes \mathbf{b}$ denotes the dyadic operation of $\mathbf{a}$ and b, i.e., $(\mathbf{a}\otimes \mathbf{b})\mathbf{c}=(\mathbf{b}·\mathbf{c})\mathbf{a}$. Here we suppose that $\left|\right|{\mathbf{x}}_i\left|\right|>0$.

2.3. An implicit $GL(n,\mathbb{R})$ Lie-group method

The equation (9) provides a new starting point for constructing the Lie-group $GL(n,\mathbb{R})$ approach. In order to develop a numerical scheme, we suppose that the coefficient matrix ${\mathbf{A}}_i$ is constant with

$${\overline{\mathbf{a}}}_i={\displaystyle \frac{{\overline{\mathbf{f}}}_i}{\left|\right|{\overline{\mathbf{x}}}_i\left|\right|}},{\overline{\mathbf{b}}}_i={\displaystyle \frac{{\overline{\mathbf{x}}}_i}{\left|\right|{\overline{\mathbf{x}}}_i\left|\right|}},$$

(10)

which can be obtained by taking the values of ${\mathbf{f}}_i$ and ${\mathbf{x}}_i$ at a suitable mid-point of $\overline{t}\in [t_0=0,t]$, where $t\le t_0+h$ and h is a small time stepsize.

By introducing new variables $w_i\left(t\right)={\overline{\mathbf{b}}}_i^T·{\mathbf{x}}_i\left(t\right)$, we obtain an augmented dynamical system

$${\displaystyle \frac{d}{dt}}\left(\begin{array}{c}{\mathbf{x}}_i\left(t\right)\\ w_i\left(t\right)\end{array}\right)=\left(\begin{array}{cc}\mathbf{0}& {\overline{\mathbf{a}}}_i\\ \mathbf{0}& {\overline{\mathbf{a}}}_i^T{\overline{\mathbf{b}}}_i\end{array}\right)\left(\begin{array}{c}{\mathbf{x}}_i\left(t\right)\\ w_i\left(t\right)\end{array}\right).$$

(11)

Because the state matrix in Eq. (11) can be decomposed as

$$\left(\begin{array}{cc}\mathbf{0}& {\overline{\mathbf{a}}}_i\\ \mathbf{0}& {\overline{\mathbf{a}}}_i^T{\overline{\mathbf{b}}}_i\end{array}\right)=\left(\begin{array}{cc}{\mathbf{I}}_n& {\displaystyle \frac{{\overline{\mathbf{a}}}_i}{{\overline{\mathbf{a}}}_i^T{\overline{\mathbf{b}}}_i}}\\ \mathbf{0}& 1\end{array}\right)\left(\begin{array}{cc}{\mathbf{I}}_n& \mathbf{0}\\ \mathbf{0}& {\overline{\mathbf{a}}}_i^T{\overline{\mathbf{b}}}_i\end{array}\right){\left(\begin{array}{cc}{\mathbf{I}}_n& {\displaystyle \frac{{\overline{\mathbf{a}}}_i}{{\overline{\mathbf{a}}}_i^T{\overline{\mathbf{b}}}_i}}\\ \mathbf{0}& 1\end{array}\right)}^{-1},$$

(12)

the solution of the linear dynamical system can be written analytically as follows

$$\left(\begin{array}{c}{\mathbf{x}}_i\left(t\right)\\ w_i\left(t\right)\end{array}\right)=exp\left[\left(\begin{array}{cc}\mathbf{0}& {\overline{\mathbf{a}}}_i\\ \mathbf{0}& {\overline{\mathbf{a}}}_i^T{\overline{\mathbf{b}}}_i\end{array}\right)(t-t_k)\right]\left(\begin{array}{c}{\mathbf{x}}_i\left(t_k\right)\\ w_i\left(t_k\right)\end{array}\right)=\left(\begin{array}{cc}{\mathbf{I}}_n& \rho (c_i,t-t_k){\overline{\mathbf{a}}}_i\\ \mathbf{0}& exp\left[{\overline{c}}_i(t-t_k)\right]\end{array}\right)\left(\begin{array}{c}{\mathbf{x}}_i\left(t_k\right)\\ w_i\left(t_k\right)\end{array}\right),$$

(13)

where $\rho (u,t)={\displaystyle \frac{e^{ut}-1}{u}}$,$c_i={\overline{\mathbf{a}}}_i^T{\overline{\mathbf{b}}}_i$. We reduce the augmented ${\mathbf{x}}_i-w_i$ system to the ${\mathbf{x}}_i$ system, obtaining

$${\mathbf{x}}_i\left(t\right)=[{\mathbf{I}}_n+\rho (c_i,t-t_k){\overline{\mathbf{a}}}_i{\overline{\mathbf{b}}}_i^T]{\mathbf{x}}_i\left(t_k\right)={\mathbf{G}}_i(t-t_k){\mathbf{x}}_i\left(t_k\right),$$

(14)

and $det{\mathbf{G}}_i(t-t_k)>0$. Note that ${\mathbf{G}}_i$ is indeed an element in the Lie group $GL(n,\mathbb{R})$. Therefore, by introducing the auxiliary variable $w_i$ and the augmented system, we are able to derive an explicit expression for the components of a one-parameter group $GL(n,\mathbb{R})$.

3. MELGDAE method for index-3 Hessenberg-DAEs

In this section, we present a novel modified extended Lie group DAE (MELGDAE) method for Hessenberg-DAEs with an index of 3,and provide a theoretical justification for the proposed method.

We first calculate the solution ${\mathbf{x}}_2\left(t_{k+1}\right)$ by the given variables ${\mathbf{x}}_1\left(t_k\right)$, ${\mathbf{x}}_2\left(t_k\right)$, then obtain the solution ${\mathbf{x}}_1\left(t_{k+1}\right)$ via given variables ${\mathbf{x}}_1\left(t_k\right)$, ${\mathbf{x}}_3\left(t_k\right)$ and ${\mathbf{x}}_2\left(t_{k+1}\right)$. Moreover, we compute the solution ${\mathbf{x}}_3\left(t_{k+1}\right)$ by the resulting value of ${\mathbf{x}}_1\left(t_{k+1}\right)$ and ${\mathbf{x}}_2\left(t_{k+1}\right)$. Hereafter, we denote the loop that solves for the variable ${\mathbf{x}}_i$ as the ${\mathbf{x}}_i$-loop. The MELGDAE method, which is presented in detail below, utilizes the notation ${\mathbf{x}}_{i\left[k\right]}={\mathbf{x}}_i\left(t_k\right)$.

 

For the ${\mathbf{x}}_2$-Loop:

S1:

Give $0\le {\theta }_2\le 1$, ${\mathbf{x}}_{1\left[k\right]}$, ${\mathbf{x}}_{2\left[k\right]}$.

S2:

Estimate ${\overline{\mathbf{x}}}_{2[k+1]}$ via Euler method as follow,

$${\overline{\mathbf{x}}}_{2[k+1]}={\mathbf{x}}_{2\left[k\right]}+h{\mathbf{f}}_2(t_k,{\mathbf{x}}_{1\left[k\right]},{\mathbf{x}}_{2\left[k\right]}).$$

(15)

S3:

Let ${\tau }_2=t_k+{\theta }_{2}h$.

S4:

Compute $\widehat{{\mathbf{x}}_{2[k+1]}}$:

S5:

If ${\widehat{\mathbf{x}}}_{2[k+1]}$ convergers according to a given stopping criterion

$$\left|\right|{\widehat{\mathbf{x}}}_{2[k+1]}-{\overline{\mathbf{x}}}_{2[k+1]}\left|\right|<{\epsilon }_2,$$

(17)

then let ${\mathbf{x}}_{2[k+1]}={\overline{\mathbf{x}}}_{2[k+1]}$; otherwise, let ${\overline{\mathbf{x}}}_{2[k+1]}={\widehat{\mathbf{x}}}_{2[k+1]}$ and goto S4.

 

For the ${\mathbf{x}}_1$-Loop:

S1:

Give $0\le {\theta }_1\le 1$, ${\mathbf{x}}_{1\left[k\right]}$, ${\mathbf{x}}_{3\left[k\right]}$, ${\mathbf{x}}_{2[k+1]}$ resulting from the ${\mathbf{x}}_2$-Loop.

S2:

Estimate ${\overline{\mathbf{x}}}_{1[k+1]}$ via Euler method as follow,

$${\overline{\mathbf{x}}}_{1[k+1]}={\mathbf{x}}_{1\left[k\right]}+h{\mathbf{f}}_1(t_k,{\mathbf{x}}_{1\left[k\right]},{\mathbf{x}}_{2[k+1]},{\mathbf{x}}_{3\left[k\right]}).$$

(18)

S3:

Let ${\tau }_1=t_k+{\theta }_{1}h$.

S4:

Compute $\widehat{{\mathbf{x}}_{1[k+1]}}$:

S5:

If ${\widehat{\mathbf{x}}}_{1[k+1]}$ convergers according to a given stopping criterion

$$\left|\right|{\widehat{\mathbf{x}}}_{1[k+1]}-{\overline{\mathbf{x}}}_{1[k+1]}\left|\right|<{\epsilon }_1,$$

(20)

then let ${\mathbf{x}}_{1[k+1]}={\overline{\mathbf{x}}}_{1[k+1]}$; otherwise, let ${\overline{\mathbf{x}}}_{1[k+1]}={\widehat{\mathbf{x}}}_{1[k+1]}$ and goto S4.

 

For the ${\mathbf{x}}_3$-Loop:

S1:

Let ${\tau }_3=t_k+{\theta }_{3}h$,${\overline{\mathbf{x}}}_{3[k+1]}={\mathbf{x}}_{3\left[k\right]}$.

S2:

Compute ${\widehat{\mathbf{x}}}_{3[k+1]}$ via Newton iterative method:

$${\widehat{\mathbf{x}}}_{3[k+1]}={\overline{\mathbf{x}}}_{3[k+1]}-{\mathbf{J}}^{-1}{\mathbf{f}}_3(t_{k+1},{\mathbf{x}}_2\left(\left[{\mathbf{x}}_1\left({\overline{\mathbf{x}}}_{3[k+1]}\right)\right]\right)),$$

(21)

where ${\tilde{\mathbf{x}}}_{1\left[k\right]}=(1-{\theta }_3){\mathbf{x}}_{1\left[k\right]}+{\theta }_3{\mathbf{x}}_{1[k+1]}$, ${\tilde{\mathbf{x}}}_{2\left[k\right]}=(1-{\theta }_3){\mathbf{x}}_{2\left[k\right]}+{\theta }_3{\mathbf{x}}_{2[k+1]}$,$\mathbf{J}=\frac{\partial {\mathbf{f}}_3}{\partial {\mathbf{x}}_3}=\frac{\partial {\mathbf{f}}_3}{\partial {\mathbf{x}}_2}\frac{\partial {\mathbf{x}}_2}{\partial {\mathbf{x}}_1}\frac{\partial {\mathbf{x}}_1}{\partial {\mathbf{x}}_3}{|}_{{\mathbf{x}}_1={\tilde{\mathbf{x}}}_{1\left[k\right]},{\mathbf{x}}_2={\tilde{\mathbf{x}}}_{2\left[k\right]},{\mathbf{x}}_3={\overline{\mathbf{x}}}_{3[k+1]}}$. We could obtain $\frac{\partial {\mathbf{x}}_1}{\partial {\mathbf{x}}_3}{|}_{{\mathbf{x}}_1={\tilde{\mathbf{x}}}_{1\left[k\right]},{\mathbf{x}}_2={\tilde{\mathbf{x}}}_{2\left[k\right]},{\mathbf{x}}_3={\overline{\mathbf{x}}}_{3[k+1]}}$, $\frac{\partial {\mathbf{x}}_2}{\partial {\mathbf{x}}_1}{|}_{{\mathbf{x}}_1={\tilde{\mathbf{x}}}_{1\left[k\right]},{\mathbf{x}}_2={\tilde{\mathbf{x}}}_{2\left[k\right]}}$ from Eq. (14).

S3:

If ${\widehat{\mathbf{x}}}_{3[k+1]}$ convergers according to a given stopping criterion

$$\left|\right|{\widehat{\mathbf{x}}}_{3[k+1]}-{\overline{\mathbf{x}}}_{3[k+1]}\left|\right|<{\epsilon }_3,$$

(24)

then let ${\mathbf{x}}_{3[k+1]}={\overline{\mathbf{x}}}_{3[k+1]}$; otherwise, let

$${\overline{\mathbf{x}}}_{3[k+1]}={\widehat{\mathbf{x}}}_{3[k+1]},$$

$${\mathbf{x}}_{2[k+1]}={\tilde{\mathbf{x}}}_{2[k+1]},$$

$${\mathbf{x}}_{1[k+1]}={\tilde{\mathbf{x}}}_{1[k+1]},$$

and goto S2.

Lemma 1. 

([23]) Let $\mathit{U}$ be an $n\times m$ matrix and $\mathit{V}$ be an $m\times n$ matrix. Then

$${\lambda }^m|\lambda {\mathit{I}}_n+\mathit{U}\mathit{V}|={\lambda }^n|\lambda {\mathit{I}}_m+\mathit{V}\mathit{U}|,$$

(25)

where λ is any scalar and ${\mathit{I}}_k$ denotes the $k\times k$ identity matrix.

Lemma 2. 

Let $\mathit{U}$ be an $n\times m$ matrix, $\mathit{V}$ be an $m\times n$ matrix, and $\mathit{a}$ and $\mathit{b}$ be non-zero $m\times 1$ vectors. Define the $m\times m$ matrix $\mathit{E}={\mathit{I}}_m+h\mathit{a}{\mathit{b}}^T$ for $h>0$. If $\mathit{U}\mathit{V}$ is nonsingular and $h<{\displaystyle \frac{1}{{\left|\right|\mathit{Q}\left|\right|}_2·{\left|\right|\mathit{a}\left|\right|}_2·{\left|\right|\mathit{b}\left|\right|}_2}}$, where $\mathit{Q}=\mathit{V}{\left(\mathit{U}\mathit{V}\right)}^{-1}\mathit{U}$, then $\mathit{U}$E$\mathit{V}$ is nonsingular.

Proof. 

To prove that $\mathbf{U}\mathbf{E}\mathbf{V}$ is nonsingular, we need to show that its determinant is nonzero.

To begin with, $\mathbf{U}\mathbf{V}$ is nonsingluar, that is, $n<=m$. And it is noted that ${\mathbf{Q}}^2=\mathbf{Q}$, then ${\left|\right|\mathbf{Q}\left|\right|}_2>=1$.

For case $n=m$, we obtain

$$\begin{array}{cc}\hfill det\left(\mathbf{U}\mathbf{E}\mathbf{V}\right)& =det\left(\mathbf{U}\right)det\left(\mathbf{E}\right)det\left(\mathbf{V}\right)\hfill \\ \hfill & =det\left(\mathbf{U}\right)det\left(\mathbf{V}\right)det\left(\mathbf{E}\right)\hfill \\ \hfill & =det\left(\mathbf{U}\mathbf{V}\right)det\left(\mathbf{E}\right)\hfill \end{array}$$

and using Lemma 1, we have

$$det\left(\mathbf{E}\right)=1+h{\mathbf{b}}^T\mathbf{a}>0,$$

where $h<{\displaystyle \frac{1}{{\left|\right|\mathbf{Q}\left|\right|}_2·{\left|\right|\mathbf{a}\left|\right|}_2·{\left|\right|\mathbf{b}\left|\right|}_2}}$.

Since both $det\left(\mathbf{U}\mathbf{V}\right)$ and $det\left(\mathbf{E}\right)$ are nonzero, then $det\left(\mathbf{U}\mathbf{E}\mathbf{V}\right)$ is also nonzero. Therefore $\mathbf{U}\mathbf{E}\mathbf{V}$ is nonsingular.

For case $n<m$, using Lemma 1, we obtain

$$\begin{array}{cc}\hfill det\left(\mathbf{U}\mathbf{E}\mathbf{V}\right)& =det(\mathbf{U}\mathbf{V}+h\mathbf{U}\mathbf{a}{\mathbf{b}}^T\mathbf{V})\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =det\left(\mathbf{U}\mathbf{V}\right)det({\mathbf{I}}_n+h{\left(\mathbf{U}\mathbf{V}\right)}^{-1}\mathbf{U}\mathbf{a}{\mathbf{b}}^T\mathbf{V})\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & =det\left(\mathbf{U}\mathbf{V}\right)(1+h{\mathbf{b}}^T\mathbf{V}{\left(\mathbf{U}\mathbf{V}\right)}^{-1}\mathbf{U}\mathbf{a})\hfill \end{array}$$

(28)

Moreover, $|{\mathbf{b}}^T{\mathbf{Q}\mathbf{a}|\le |\left|\mathbf{b}\right||}_2·{\left|\right|\mathbf{Q}\mathbf{a}\left|\right|}_2\le {\left|\right|\mathbf{Q}\left|\right|}_2·{\left|\right|\mathbf{a}\left|\right|}_2·{\left|\right|\mathbf{b}\left|\right|}_2$, then $1+h{\mathbf{b}}^T\mathbf{Q}\mathbf{a}>0$, which means that $\mathbf{U}\mathbf{E}\mathbf{V}$ is nonsingular from equation (28). □

Lemma 3. 

Let ${\mathit{U}}_1$ be an $n_1\times n_3$ matrix, ${\mathit{U}}_2$ be an $n_2\times n_1$ matrix, and ${\mathit{U}}_3$ be an $n_3\times n_2$ matrix. Let ${\mathit{a}}_1$ and ${\mathit{b}}_1$ be nonzero $n_1\times 1$ vectors, and let ${\mathit{a}}_2$ and ${\mathit{b}}_2$ be nonzero $n_2\times 1$ vectors. Define the matrices ${\mathit{E}}_1={\mathit{I}}_{n_1}+h{\mathit{a}}_1{\mathit{b}}_1^T$ and ${\mathit{E}}_2={\mathit{I}}_{n_2}+h{\mathit{a}}_2{\mathit{b}}_2^T$, where $h>0$. Suppose that ${\mathit{U}}_3{\mathit{U}}_2{\mathit{U}}_1$ is nonsingular and that h satisfies

$$h<min\{{\displaystyle \frac{1}{3·\left|\right|{\mathit{Q}}_1{\left|\right|}_2·\left|\right|{\mathit{a}}_1{\left|\right|}_2·\left|\right|{\mathit{b}}_1{\left|\right|}_2}},{\displaystyle \frac{1}{3·\left|\right|{\mathit{Q}}_2{\left|\right|}_2·\left|\right|{\mathit{a}}_2{\left|\right|}_2·\left|\right|{\mathit{b}}_2{\left|\right|}_2}},{\displaystyle \frac{1}{{2\left|\right|\mathit{M}\left|\right|}_2}}\},$$

where ${\mathit{Q}}_1={\mathit{U}}_1{\left({\mathit{U}}_3{\mathit{U}}_2{\mathit{U}}_1\right)}^{-1}{\mathit{U}}_3{\mathit{U}}_2$, ${\mathit{Q}}_2={\mathit{U}}_2{\mathit{U}}_1{\left({\mathit{U}}_3{\mathit{U}}_2{\mathit{U}}_1\right)}^{-1}{\mathit{U}}_3$, and $\mathit{M}={\mathit{U}}_1{\left({\mathit{U}}_3{\mathit{U}}_2{\mathit{U}}_1\right)}^{-1}{\mathit{U}}_3$ ${\mathit{a}}_2{\mathbf{b}}_2^T{\mathit{U}}_2$. Then, the matrix ${\mathit{U}}_3{\mathit{E}}_2{\mathit{U}}_2{\mathit{E}}_1{\mathit{U}}_1$ is nonsingular.

Proof. 

Firstly, ${\mathbf{U}}_3{\mathbf{U}}_2{\mathbf{U}}_1$ is nonsingular, that is $n_3\le n_1$ and $n_3\le n_2$.

Let

$${\mathbf{V}}_2={\mathbf{U}}_2{\mathbf{U}}_1,$$

that is ${\mathbf{U}}_3{\mathbf{V}}_2$ is nonsingular. And ${\mathbf{E}}_2={\mathbf{I}}_{n_2}+h{\mathbf{a}}_2{\mathbf{b}}_2^T$ and $h<{\displaystyle \frac{1}{3\left|\right|{\mathbf{Q}}_2{\left|\right|}_2·\left|\right|{\mathbf{a}}_2{\left|\right|}_2·\left|\right|{\mathbf{b}}_2{\left|\right|}_2}}$, where ${\mathbf{Q}}_2={\mathbf{U}}_2{\mathbf{U}}_1{\left({\mathbf{U}}_3{\mathbf{U}}_2{\mathbf{U}}_1\right)}^{-1}{\mathbf{U}}_3$. By Lemma 2, ${\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{V}}_2$ is nonsingular, which means ${\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{U}}_2{\mathbf{U}}_1$ is nonsingular.

Then, let

$${\mathbf{U}}^{\prime }={\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{U}}_2,$$

and ${\mathbf{U}}^{\prime }{\mathbf{U}}_1$ is nonsingular. In addtion, ${\mathbf{E}}_1={\mathbf{I}}_{n_1}+h{\mathbf{a}}_1{\mathbf{b}}_1^T$, if $h<{\displaystyle \frac{1}{\left|\right|{\mathbf{Q}}_1^{\prime }{\left|\right|}_2·\left|\right|{\mathbf{a}}_1{\left|\right|}_2·\left|\right|{\mathbf{b}}_1{\left|\right|}_2}}$, where ${\mathbf{Q}}_1^{\prime }={\mathbf{U}}_1{\left({\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{U}}_2{\mathbf{U}}_1\right)}^{-1}{\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{U}}_2$. By Lemma 2, ${\mathbf{U}}^{\prime }{\mathbf{E}}_1{\mathbf{U}}_1$ is nonsingular, which means ${\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{U}}_2{\mathbf{E}}_1{\mathbf{U}}_1$ is nonsingular.

Now, we need proof

$$h\left|\right|{\mathbf{Q}}_1^{\prime }{\left|\right|}_2·\left|\right|{\mathbf{a}}_1{\left|\right|}_2·\left|\right|{\mathbf{b}}_1{\left|\right|}_2<1.$$

(29)

It is noted that

$${\mathbf{Q}}_1^{\prime }={({\mathbf{I}}_{n_1}+h\mathbf{M})}^{-1}{\mathbf{Q}}_1({\mathbf{I}}_{n_1}+h\mathbf{M}),$$

(30)

where $h<{\displaystyle \frac{1}{{2\left|\right|\mathbf{M}\left|\right|}_2}}$. Moreover, we are able to obtain $\left|\right|{({\mathbf{I}}_{n_1}+h\mathbf{M})}^{-1}{\left|\right|}_2\le {\displaystyle \frac{1}{1-{\left|\right|h\mathbf{M}\left|\right|}_2}}<2$ and $\left|\right|{\mathbf{I}}_{n_1}+{h\mathbf{M}\left|\right|}_2\le \left|\right|{\mathbf{I}}_{n_1}{\left|\right|}_2+{\left|\right|h\mathbf{M}\left|\right|}_2<{\displaystyle \frac{3}{2}}$.

Then, we obtain

$$\begin{array}{cc}\hfill \left|\right|{\mathbf{Q}}_1^{\prime }{\left|\right|}_2& =\left|\right|{({\mathbf{I}}_{n_1}+h\mathbf{M})}^{-1}{\mathbf{Q}}_1({\mathbf{I}}_{n_1}+h\mathbf{M}){\left|\right|}_2\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & \le \left|\right|{({\mathbf{I}}_{n_1}+h\mathbf{M})}^{-1}{\left|\right|}_2·\left|\right|{\mathbf{Q}}_1{\left|\right|}_2·\left|\right|{\mathbf{I}}_{n_1}+h\mathbf{M}{\left|\right|}_2\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & <3\left|\right|{\mathbf{Q}}_1{\left|\right|}_2\hfill \end{array}$$

(33)

For $h<{\displaystyle \frac{1}{3\left|\right|{\mathbf{Q}}_1{\left|\right|}_2·\left|\right|{\mathbf{a}}_1{\left|\right|}_2·\left|\right|{\mathbf{b}}_1{\left|\right|}_2}}$, then we obtain the inequality (29). □

Theorem 1. 

For the ${\mathit{x}}_3$-Loop of Hessenberg DAEs with index 3, if h is sufficiently small, then the Jacobian matrix $\mathit{J}$ is nonsingular.

Proof. 

We first compute the $\frac{\partial {\mathbf{x}}_2}{\partial {\mathbf{x}}_1}$ by Eq.(22).

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{x}}_2}{\partial {\mathbf{x}}_1}& ={\tilde{d}}_{12\left[k\right]}\{{\tilde{\mathbf{a}}}_{12\left[k\right]}\frac{\partial }{\partial {\mathbf{x}}_1}{\tilde{\rho }}_{2\left[k\right]}+{\tilde{\rho }}_{2\left[k\right]}\frac{\partial }{\partial {\mathbf{x}}_1}{\tilde{\mathbf{a}}}_{12\left[k\right]}\}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & ={\tilde{d}}_{12\left[k\right]}\{{\tilde{\mathbf{a}}}_{12\left[k\right]}\left[{\displaystyle \frac{(ut-1)exp\left(ut\right)+1}{u^2}}{|}_{u={\tilde{c}}_{12\left[k\right]},t=h}\right]\frac{\partial }{\partial {\mathbf{x}}_1}{\tilde{c}}_{12\left[k\right]}+{\tilde{\rho }}_{2\left[k\right]}\frac{\partial }{\partial {\mathbf{x}}_1}{\tilde{\mathbf{a}}}_{12\left[k\right]}\}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & ={\tilde{d}}_{12\left[k\right]}\{{\tilde{\mathbf{a}}}_{12\left[k\right]}{\displaystyle \frac{({\tilde{c}}_{12\left[k\right]}h-1)exp\left({\tilde{c}}_{12\left[k\right]}h\right)+1}{{\tilde{c}}_{12\left[k\right]}^2}}\left[{\tilde{\mathbf{b}}}_{12\left[k\right]}^T\frac{\partial }{\partial {\mathbf{x}}_1}{\tilde{\mathbf{a}}}_{12\left[k\right]}\right]+{\tilde{\rho }}_{2\left[k\right]}\frac{\partial }{\partial {\mathbf{x}}_1}{\tilde{\mathbf{a}}}_{12\left[k\right]}\}\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & ={\tilde{d}}_{12\left[k\right]}\{{\tilde{\rho }}_{2\left[k\right]}{\mathbf{I}}_{n_2}+{\displaystyle \frac{({\tilde{c}}_{12\left[k\right]}h-1)exp\left({\tilde{c}}_{12\left[k\right]}h\right)+1}{{\tilde{c}}_{12\left[k\right]}^2}}{\tilde{\mathbf{a}}}_{12\left[k\right]}{\tilde{\mathbf{b}}}_{12\left[k\right]}^T\}\frac{\partial }{\partial {\mathbf{x}}_1}{\tilde{\mathbf{a}}}_{12\left[k\right]}\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & =\frac{{\tilde{d}}_{12\left[k\right]}}{\left|\right|{\tilde{\mathbf{x}}}_{2\left[k\right]}\left|\right|}\{{\tilde{\rho }}_{2\left[k\right]}{\mathbf{I}}_{n_2}+{\displaystyle \frac{({\tilde{c}}_{12\left[k\right]}h-1)exp\left({\tilde{c}}_{12\left[k\right]}h\right)+1}{{\tilde{c}}_{12\left[k\right]}^2}}{\tilde{\mathbf{a}}}_{12\left[k\right]}{\tilde{\mathbf{b}}}_{12\left[k\right]}^T\}\frac{\partial {\mathbf{f}}_2(t,{\mathbf{x}}_1,{\mathbf{x}}_2)}{\partial {\mathbf{x}}_1}.\hfill \end{array}$$

(38)

Similarly, we obtain the $\frac{\partial {\mathbf{x}}_1}{\partial {\mathbf{x}}_3}$ from Eq. (23).

$$\begin{array}{cc}\hfill \frac{\partial {\mathbf{x}}_1}{\partial {\mathbf{x}}_3}& ={\tilde{d}}_{11\left[k\right]}\{{\tilde{\mathbf{a}}}_{11\left[k\right]}\frac{\partial }{\partial {\mathbf{x}}_3}{\tilde{\rho }}_{1\left[k\right]}+{\tilde{\rho }}_{1\left[k\right]}\frac{\partial }{\partial {\mathbf{x}}_3}{\tilde{\mathbf{a}}}_{11\left[k\right]}\}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & ={\tilde{d}}_{11\left[k\right]}\{{\tilde{\rho }}_{1\left[k\right]}{\mathbf{I}}_{n_1}+{\displaystyle \frac{({\tilde{c}}_{11\left[k\right]}h-1)exp\left({\tilde{c}}_{11\left[k\right]}h\right)+1}{{\tilde{c}}_{11\left[k\right]}^2}}{\tilde{\mathbf{a}}}_{11\left[k\right]}{\tilde{\mathbf{b}}}_{11\left[k\right]}^T\}\frac{\partial }{\partial {\mathbf{x}}_3}{\tilde{\mathbf{a}}}_{11\left[k\right]}\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill & =\frac{{\tilde{d}}_{11\left[k\right]}}{\left|\right|{\tilde{\mathbf{x}}}_{1\left[k\right]}\left|\right|}\{{\tilde{\rho }}_{1\left[k\right]}{\mathbf{I}}_{n_1}+{\displaystyle \frac{({\tilde{c}}_{11\left[k\right]}h-1)exp\left({\tilde{c}}_{11\left[k\right]}h\right)+1}{{\tilde{c}}_{11\left[k\right]}^2}}{\tilde{\mathbf{a}}}_{11\left[k\right]}{\tilde{\mathbf{b}}}_{11\left[k\right]}^T\}\frac{\partial {\mathbf{f}}_1(t,{\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3)}{\partial {\mathbf{x}}_3}.\hfill \end{array}$$

(41)

Futhermore, we are able to get

$$\begin{array}{cc}\hfill \mathbf{J}& =\frac{\partial {\mathbf{f}}_3}{\partial {\mathbf{x}}_3}=\frac{\partial {\mathbf{f}}_3}{\partial {\mathbf{x}}_2}\frac{\partial {\mathbf{x}}_2}{\partial {\mathbf{x}}_1}\frac{\partial {\mathbf{x}}_1}{\partial {\mathbf{x}}_3}\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & =\frac{{\tilde{d}}_{12\left[k\right]}}{\left|\right|{\tilde{\mathbf{x}}}_{2\left[k\right]}\left|\right|}\frac{{\tilde{d}}_{11\left[k\right]}}{\left|\right|{\tilde{\mathbf{x}}}_{1\left[k\right]}\left|\right|}\frac{\partial {\mathbf{f}}_3}{\partial {\mathbf{x}}_2}{\mathbf{E}}_2\frac{\partial {\mathbf{f}}_2}{\partial {\mathbf{x}}_1}{\mathbf{E}}_1\frac{\partial {\mathbf{f}}_1}{\partial {\mathbf{x}}_3},\hfill \end{array}$$

(43)

where ${\mathbf{E}}_{i\left[k\right]}={\tilde{\rho }}_{i\left[k\right]}{\mathbf{I}}_{n_i}+{\displaystyle \frac{({\tilde{c}}_{1i\left[k\right]}h-1)exp\left({\tilde{c}}_{1i\left[k\right]}h\right)+1}{{\tilde{c}}_{1i\left[k\right]}^2}}{\tilde{\mathbf{a}}}_{1i\left[k\right]}{\tilde{\mathbf{b}}}_{1i\left[k\right]}^T$ ($i=1,2$), and $\mathbf{F}={\mathbf{U}}_3{\mathbf{U}}_2{\mathbf{U}}_1$( ${\mathbf{U}}_3=\frac{\partial {\mathbf{f}}_3}{\partial {\mathbf{x}}_2},{\mathbf{U}}_2=\frac{\partial {\mathbf{f}}_2}{\partial {\mathbf{x}}_1},{\mathbf{U}}_1=\frac{\partial {\mathbf{f}}_1}{\partial {\mathbf{x}}_3}$) is nonsingular.

From Lemma 1, it is noted that

$$|{\mathbf{E}}_{i\left[k\right]}|={\tilde{\rho }}_{i\left[k\right]}^{n_i-1}\{{\tilde{\rho }}_{i\left[k\right]}+{\displaystyle \frac{({\tilde{c}}_{1i\left[k\right]}h-1)exp\left({\tilde{c}}_{1i\left[k\right]}h\right)+1}{{\tilde{c}}_{1i\left[k\right]}^2}}{\tilde{\mathbf{b}}}_{1i\left[k\right]}^T{\tilde{\mathbf{a}}}_{1i\left[k\right]}\}=h{\tilde{\rho }}_{i\left[k\right]}^{n_i-1}exp\left({\tilde{c}}_{1i\left[k\right]}h\right)\ne 0.$$

For small $h<H$, ${\mathbf{E}}_{i\left[k\right]}=h{\mathbf{I}}_{n_i}+h^2{\tilde{\mathbf{a}}}_{1i\left[k\right]}{\tilde{\mathbf{b}}}_{1i\left[k\right]}^T+O\left(h^3\right)\approx h{\mathbf{I}}_{n_i}+h^2{\tilde{\mathbf{a}}}_{1i\left[k\right]}{\tilde{\mathbf{b}}}_{1i\left[k\right]}^T$. Moreover, by using Lemma 3, if

$$h<\underset{k}{min}\{{\displaystyle \frac{1}{3\left|\right|{\mathbf{Q}}_{1\left[k\right]}{\left|\right|}_2·\left|\right|{\tilde{\mathbf{a}}}_{11\left[k\right]}{\left|\right|}_2·\left|\right|{\tilde{\mathbf{b}}}_{11\left[k\right]}{\left|\right|}_2}},{\displaystyle \frac{1}{3\left|\right|{\mathbf{Q}}_{2\left[k\right]}{\left|\right|}_2·\left|\right|{\tilde{\mathbf{a}}}_{12\left[k\right]}{\left|\right|}_2·\left|\right|{\tilde{\mathbf{b}}}_{12\left[k\right]}{\left|\right|}_2}},{\displaystyle \frac{1}{2\left|\right|{\mathbf{M}}_k{\left|\right|}_2}}\},$$

where ${\mathbf{Q}}_{1\left[k\right]}={\mathbf{U}}_1{\left({\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{U}}_2{\mathbf{U}}_1\right)}^{-1}{\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{U}}_2$, ${\mathbf{Q}}_{2\left[k\right]}={\mathbf{U}}_2{\mathbf{U}}_1{\left({\mathbf{U}}_3{\mathbf{U}}_2{\mathbf{U}}_1\right)}^{-1}{\mathbf{U}}_3$ and

${\mathbf{M}}_k={\mathbf{U}}_1{\left({\mathbf{U}}_3{\mathbf{U}}_2{\mathbf{U}}_1\right)}^{-1}{\mathbf{U}}_3{\tilde{\mathbf{a}}}_{12\left[k\right]}{\tilde{\mathbf{b}}}_{12\left[k\right]}^T{\mathbf{U}}_2$, then ${\mathbf{U}}_3{\mathbf{E}}_2{\mathbf{U}}_2{\mathbf{E}}_1{\mathbf{U}}_1$ is nonsingular. Therefore, satisfying for small h with

$$h<\underset{k}{min}\{H,{\displaystyle \frac{1}{3\left|\right|{\mathbf{Q}}_{1\left[k\right]}{\left|\right|}_2·\left|\right|{\tilde{\mathbf{a}}}_{11\left[k\right]}{\left|\right|}_2·\left|\right|{\tilde{\mathbf{b}}}_{11\left[k\right]}{\left|\right|}_2}},{\displaystyle \frac{1}{3\left|\right|{\mathbf{Q}}_{2\left[k\right]}{\left|\right|}_2·\left|\right|{\tilde{\mathbf{a}}}_{12\left[k\right]}{\left|\right|}_2·\left|\right|{\tilde{\mathbf{b}}}_{12\left[k\right]}{\left|\right|}_2}},{\displaystyle \frac{1}{2\left|\right|{\mathbf{M}}_k{\left|\right|}_2}}\},$$

we are able to obtain $\mathbf{J}$ is nonsingular. □

4. Numerical Examples

In this section, we do numerical simulation experiments to illustrate the performance of the MELGDAE. All codes are written under Windows 10 system and run on a personal computer with Intel(R) Core(TM) i5-3570 CPU @ $3.40$ GHz, $4.00$ GB RAM and 64-bit operating system.

We consider a Hessenberg DAEs system with index 3 as follows, which can be found in [24]. where $t\in [0,1]$, $z_1\left(0\right)=z_2\left(0\right)=z_3\left(0\right)=z_4\left(0\right)=z_5\left(0\right)=1$ are the consistent initial values. And the exact solution is given by

$$z_1\left(t\right)=z_3\left(t\right)=e^{2t},z_2\left(t\right)=z_4\left(t\right)=e^{-t},z_5\left(t\right)=e^t.$$

Let ${\mathbf{x}}_1\left(t\right)={(z_1\left(t\right),z_2\left(t\right))}^T,{\mathbf{x}}_2\left(t\right)={(z_3\left(t\right),z_4\left(t\right))}^T,{\mathbf{x}}_3\left(t\right)=z_5\left(t\right)$, and ${\mathbf{f}}_1={(g_1,g_2)}^T,{\mathbf{f}}_2={(g_3,g_4)}^T,{\mathbf{f}}_3=g_5$, then $\mathbf{F}=\frac{\partial {\mathbf{f}}_3}{\partial {\mathbf{x}}_2}\frac{\partial {\mathbf{f}}_2}{\partial {\mathbf{x}}_1}\frac{\partial {\mathbf{f}}_1}{\partial {\mathbf{x}}_3}=2z_2\left(t\right)(2z_2{\left(t\right)}^2{z}_3\left(t\right)+1)$.

To obtain DAEs with index 2, we differentiate the constraint $g_5$ to obtain the following equation.

$$0=g_6=z_1\left(t\right)z_4\left(t\right)-z_2\left(t\right)z_3\left(t\right).$$

(45)

Then, we are able to obtain index 2 DAEs below.

We compared the numerical results of the MELGDAE method with those of RADAU, a code based on implicit Runge-Kutta methods, and MEBDFI, a code based on the MEBDF method that has better theoretical properties than BDF methods, to demonstrate the capability of the MELGDAE method. The comparison was made using the Hessenberg DAEs with index 2 or 3 mentioned above.

Without loss of generality, we assume the absolute error (${\epsilon }_1$,${\epsilon }_2$,${\epsilon }_3$) is ${10}^{-8}$, the relative error is ${10}^{-5}$, the (initial) step size is $h={10}^{-3}$, ${\theta }_1={\theta }_2={\theta }_3=0.5$. Based on Figure 1 and Figure 2, the MELGDAE method showed comparable accuracy with the RADAU method for the differential variables $z_1\left(t\right)$ and $z_3\left(t\right)$ for both index-2 and index-3 Hessenberg DAEs. In contrast, the MEBDFI method was less accurate for these variables. For the algebraic variable $z_5\left(t\right)$, both RADAU and MEBDFI methods provided more accurate results than the MELGDAE method. However, for the algebraic constraints $g_5\left(t\right)$ and $g_6\left(t\right)$, the MELGDAE method exhibited significantly higher accuracy compared to the RADAU and MEBDFI methods.

Upon observing Figure 1 and Figure 2, it is clear that all three methods exhibit a reduction in the absolute errors of $z_1\left(t\right)$, $z_3\left(t\right)$, and $z_5\left(t\right)$ as the index increases from 2 to 3. This reduction is not significant for the MELGDAE and RADAU methods. Therefore, we can conclude that the MELGDAE method has the potential to be extended for computing DAEs systems with an index higher than 3.

Furthermore, we conducted numerical experiments to verify the convergence accuracy of the MELGDAE method for index 3 DAEs. We assumed $h=\{2^{-10},2^{-9},2^{-8},2^{-7},2^{-6},2^{-5},2^{-4}\}$ and calculated their constants in $-lo{g}_2({max}_k|z_{i\left[k\right]}-z_i\left(t_k\right)\left|\right)=\nu (-lo{g}_{2}h)+\mu$ ($i=1,2,3,4,5$) and $-lo{g}_2({max}_k|g_5(t_k,z_{3\left[k\right]},z_{4\left[k\right]})\left|\right)=\nu (-lo{g}_{2}h)+\mu$, which represents ${max}_k|z_{i\left[k\right]}-z_i\left(t_k\right)|=O\left(h^{\nu }\right)$ and ${max}_k\left|g_5(t_k,z_{3\left[k\right]},z_{4\left[k\right]})\right|=O\left(h^{\nu }\right)$, using the standard least-squares method. The results are presented in Figure 3, which indicate that the MELGDAE method exhibits second-order convergence for $z_1\left(t\right)$, $z_2\left(t\right)$, $z_3\left(t\right)$, and $z_4\left(t\right)$ (all differential variables) and $g_5$ (the algebraic constraints), and first-order convergence for $z_5\left(t\right)$ (algebraic variable). Thus, this method may be extended to compute systems with an index of 4 or higher and is suitable for DAE systems that cannot be directly solved using other classical methods.

5. Conclusions

The Modified Extended Lie Group Differential Algebraic Equations (MELGDAE) method was developed in this paper to solve Hessenberg-DAEs with index 3. While the LGDAE method presented in [18] is only capable of handling index 2 systems, the MELGDAE method can solve higher index systems by taking into account the differential equation components with respect to the index in Hessenberg form. It is worth mentioning that the MELGDAE method demonstrates superior theoretical properties compared to the LGDAEs, as illustrated in Theorem 1.

We conducted a comparison of the MELGDAE method with the standard methods RADAU and MEBDFI on index 2 and 3 DAE systems. The MELGDAE method demonstrated promising results in terms of solution errors when compared to RADAU and MEBDFI, which are known to be powerful on different DAE systems. However, the MELGDAE method excelled in preserving the algebraic constraints. Moreover, all differential variables in index 3 Hessenberg-DAEs exhibited order 2 convergence with the MELGDAE method, whereas MEBDFI in [24] and RADAU in [7,25] demonstrated order reduction. Although the MELGDAE method has high stability and accuracy when solving strong nonlinear Hessenberg DAEs with index 3, there are still some theoretical gaps that need to be addressed in the future. Additionally, the MELGDAE method can be further extended to solve DAEs with an index of 4 or higher, which cannot be directly solved using standard methods.