Questions concerning differential-algebraic operators: Toward a reliable direct numerical treatment of differential-algebraic equations

TL;DR

This paper investigates the properties of differential-algebraic operators and their approximations, focusing on polynomial collocation methods to improve the direct numerical solution of higher-index differential-algebraic equations.

Contribution

It provides a detailed analysis of first-order differential-algebraic operators and justifies the use of overdetermined polynomial collocation for higher-index equations.

Findings

Justification of polynomial collocation for higher-index DAE

Analysis of differential-algebraic operator properties

Discussion of practical aspects for higher-order operators

Abstract

The nature of so-called differential-algebraic operators and their approximations is constitutive for the direct treatment of higher-index differential-algebraic equations. We treat first-order differential-algebraic operators in detail and contribute to justify the overdetermined polynomial collocation applied to higher-index differential-algebraic equations. Besides, we discuss several practical aspects concerning higher-order differential-algebraic operators and the associated equations.

Tables

Table 1. Table 1: Componentwise maximal error, polynomial degree N = 3 𝑁 3 N=3

$n$	Standard collocation			Least-squares collocation
	$i=1$	$i=2$	$i=3$	$i=1$	$i=2$	$i=3$
20	5.56e+006	3.03e+004	5.99e+004	2.09e-4	1.10e-06	2.18e-06
40	1.55e+017	4.23e+014	8.41e+014	5.03e-5	1.31e-07	2.65e-07
80	5.70e+038	7.76e+035	1.55e+036	1.23e-5	1.60e-08	3.20e-08
160	3.36e+082	2.29e+079	4.57e+079	3.06e-6	1.98e-09	4.00e-09
320	4.93e+170	1.68e+167	3.35e+167	7.68e-7	2.50e-10	5.00e-10

Table 2. Table 2: Error ( ∥ e 1 ∥ L 2 2 + ∥ e 2 ∥ H 1 2 + ∥ e 3 ∥ H 1 2 ) 1 2 (\rVert e_{1}\rVert^{2}_{L^{2}}+\rVert e_{2}\rVert^{2}_{H^{1}}+\rVert e_{3}\rVert^{2}_{H^{1}})^{\frac{1}{2}} of the collocation solution for N = 3 𝑁 3 N=3

	$M=2N+1$ uniform points		$M=N+1$ Gaussian points
$n$	error	order	error	order
10	6.31e-4		6.46e-4
20	1.44e-4	2.1	1.45e-4	2.2
40	3.47e-5	2.1	3.47e-5	2.1
80	8.53e-6	2.0	8.53e-6	2.0
160	2.12e-6	2.0	2.12e-6	2.0
320	5.27e-7	2.0	5.27e-7	2.0

Table 3. Table 3: Error ( ∥ e 1 ∥ L 2 2 + ∥ e 2 ∥ H 1 2 + ∥ e 3 ∥ H 1 2 ) 1 2 (\rVert e_{1}\rVert^{2}_{L^{2}}+\rVert e_{2}\rVert^{2}_{H^{1}}+\rVert e_{3}\rVert^{2}_{H^{1}})^{\frac{1}{2}} of the collocation solution for N = 1 𝑁 1 N=1

	$M=3$ uniform points		$M=2$ Gaussian points
$n$	error	order	error	order
10	5.65e-1		5.65e-1
20	3.93e-1	0.5	3.93e-1	0.5
40	2.49e-1	0.6	2.49e-1	0.7
80	1.85e-1	0.4	1.85e-1	0.4
160	1.42e-1	0.4	1.42e-1	0.4
320	1.12e-1	0.3	1.12e-1	0.3

Equations

$$x^{\prime }(t)+B(t)x(t)=g(t),t\in [a,b],x(a)=d\in {\mathbb{R}}^m,$$

$$x(t)=X(t,a)d+(Vg)(t),t\in [a,b],d\in {\mathbb{R}}^m,(Vg)(t)={\int }_a^{t}X(t,s)g(s)\mathrm{ds}.$$

$${\mathcal{T}}^{-1}(g,d)(t)=X(t,a)d+(Vg)(t),t\in [a,b].$$

$$x_2^{\prime }(t)+x_1(t)$$

$$t\eta x_2^{\prime }(t)+x_3^{\prime }(t)+(\eta +1)x_2(t)$$

$$t\eta x_2(t)+x_3(t)$$

$$x_1(t)=e^{-t}\sin t,x_2(t)=e^{-2t}\sin t,x_3(t)=e^{-t}\cos t.$$

$$Tx$$

$$\ker T=\{0\},\mathrm{im}T=\{g\in \mathcal{C}([0,1],{\mathbb{R}}^3):g_3\in {\mathcal{C}}^1([0,1],\mathbb{R}),g_2-g_3^{\prime }\in {\mathcal{C}}^1([0,1],\mathbb{R})\},$$

$$T^{-1}g=\left[\begin{array}{c}{g}_1-(g_2-g_3^{\prime }{)}^{\prime }\\ g_2-g_3^{\prime }\\ g_3-t\eta (g_2-g_3^{\prime })\end{array}\right],g\in \mathrm{im}T.$$

$$T^{-1}g=$$

$$g\in \{v\in \mathcal{C}([0,1],{\mathbb{R}}^3):v_2\in {\mathcal{C}}^1([0,1],\mathbb{R}),v_3\in {\mathcal{C}}^2([0,1],\mathbb{R})\}\subset \mathrm{im}T.$$

$$T$$

$$\mathring{T}\in \mathcal{L}(X,Y),\mathring{T}x:=E{x}^{\prime }+Fx,x\in \mathrm{dom}\mathring{T}:={\mathcal{C}}^1([a,b],{\mathbb{R}}^m)\subseteq X,$$

$$\Vert x{\Vert }_T=\Vert x\Vert +\Vert Tx\Vert ,x\in X_T:=\mathrm{dom}T,$$

$$\ker A(t)\oplus \mathrm{im}D(t)={\mathbb{R}}^k,t\in [a,b].$$

$$\mathring{T}x=E{x}^{\prime }+Fx$$

$$=A(Dx{)}^{\prime }+Bx,\text{for all}x\in \mathrm{dom}\mathring{T}.$$

$$H_D^1((a,b),{\mathbb{R}}^m)=\{x\in L^2((a,b),{\mathbb{R}}^m):Dx\in H^1((a,b),{\mathbb{R}}^k)\}$$

$$(x,\bar{x}{)}_{H_D^1}=(x,\bar{x}{)}_{L^2}+((Dx{)}^{\prime },(D\bar{x}{)}^{\prime }{)}_{L^2}$$

$$\mathrm{dom}T=H_D^1((a,b),{\mathbb{R}}^m),Tx=A(Dx{)}^{\prime }+Bx,x\in \mathrm{dom}T.$$

$$\Vert (Dx{)}^{\prime }{\Vert }_{L^2}$$

$$\Vert x{\Vert }_{H_D^1}^2$$

$${\Pi }_{can}\in \mathcal{C}([a,b],\mathcal{L}({\mathbb{R}}^m,{\mathbb{R}}^m)),{\Pi }_{can}(t{)}^2={\Pi }_{can}(t),\mathrm{rank}{\Pi }_{can}(t)=l,t\in [a,b],$$

$$A(DX{)}^{\prime }+BX=0,X(c)={\Pi }_{can}(c),c\in [a,b],$$

$$T:H_D^1(\mathcal{I},{\mathbb{R}}^m)\to L^2(\mathcal{I},{\mathbb{R}}^m),Tx:=A(Dx{)}^{\prime }+Bx,x\in H_D^1(\mathcal{I},{\mathbb{R}}^m),\mathcal{I}:=(a,b).$$

$$\dim \ker T=l.$$

$$T_{IC}:H_D^1(\mathcal{I},{\mathbb{R}}^m)\to {\mathbb{R}}^l,T_{IC}x=G_{a}x(a),x\in H_D^1(\mathcal{I},{\mathbb{R}}^m),$$

$$\mathcal{T}:H_D^1(\mathcal{I},{\mathbb{R}}^m)\to L^2(\mathcal{I},{\mathbb{R}}^m)\times {\mathbb{R}}^l,\mathcal{T}x=(Tx,T_{IC}x),x\in H_D^1(\mathcal{I},{\mathbb{R}}^m).$$

$${\mathcal{T}}^{-1}(g,d)(t)$$

$$(Vg)(t)={\int }_a^{t}X(t,s)G_{\mu }(s{)}^{-1}g(s)\mathrm{ds},$$

$$\mathfrak{D}g=v_0+v_1+\cdots +v_{\mu -1},$$

$$v_{\mu -1}$$

Full text

1 Introduction

To a large extend, in the framework of numerical analysis, differential-algebraic equations (DAEs), in particular higher-index ones, are recognized as special ordinary differential equations (ODEs), and, accordingly, they are treated by means of derivative-array systems and an involved or preceded expensive index-reduction.

In contrast, the experiments and theoretical contributions reportet in [3, 10, 9] give rise to the conjecture that next to the existing derivative-array based methods there is further potential toward a reliable direct numerical treatment of DAEs. The main aim of this note is to fill the gap between the theoretical convergence results for least-squares collocation methods [10, 9] and its practical realization. Moreover, we will explore the relevant scope concerning higher-order DAEs.

Recap well-known facts concerning first-order ordinary differential operators: Let $B(t)\in\mathbb{R}^{m,m}$ be continuous. The initial value problem (IVP)

[TABLE]

has the unique solution

[TABLE]

Introduce the operators $Tx:=x^{\prime}+Bx$, $T_{IC}x:=x(a)$,  $\mathcal{T}x:=(Tx,T_{IC}x)$, such that the IVP is represented by the operator equation $\mathcal{T}x=(g,d)$ and the inverse of $\mathcal{T}$ is given by

[TABLE]

As it is well known, these operators have useful properties making ODE problems easily accessible for the numerical treatment: $T:\mathcal{C}^{1}([a,b],\mathbb{R}^{m})\rightarrow\mathcal{C}([a,b],\mathbb{R}^{m})$ is bounded and surjective, $\dim{\rm ker}\,T=m$, and hence, $T$ is fredholm, and $\mathcal{T}:\mathcal{C}^{1}([a,b],\mathbb{R}^{m})\rightarrow\mathcal{C}([a,b],\mathbb{R}^{m})\times\mathbb{R}^{m}$ is bounded and bijective, thus a homeomorphism. The same properties persist in the Hilbert space setting $T:H^{1}((a,b),\mathbb{R}^{m})\rightarrow L^{2}((a,b),\mathbb{R}^{m})$, $\mathcal{T}:H^{1}((a,b),\mathbb{R}^{m})\rightarrow L^{2}((a,b),\mathbb{R}^{m})\times\mathbb{R}^{m}$. Such properties will survive for index one differential-algebraic equations in an appropriately modified version. The situation in the higher index case is much more involved.

The paper is organized as follows: Section 2 deals with a small first-order index-three example and foreshadows the potential of overdetermined polynomial collocation. Section 3 is devoted to features of regular first-order arbitrary-index differential-algebraic operators acting in reasonable Hilbert spaces. We turn back to the overdetermined polynomial collocation in Section 4, now considered by operators representing finite-dimensional approximations and provide with the main Theorem 4.3 sufficient conditions justifying overdetermined collocation. Section 5 surveys higher-order differential-algebraic operators in this context.

Below, though using different norms we mark the related norms by extra tags merely on those places where confusions are actually imminent.

List of some symbols and abbreviations

[TABLE]

2 A symptomatic example

The DAE

[TABLE]

has index 3 uniformly for every $\eta\in\mathbb{R}$. To each sufficiently smooth $y$ there exists ecactly one solution $x$. The DAE is somewhat snaky, so that step by step integration methods generate waste unless an a priori or a posteriori incorporated regularization via derivative array systems is incorporated, e.g., [12, Example 8.5].

Here we set $\eta=-2$ and determine $g_{1}$, $g_{2}$, $g_{3}$ sucht that the solution becomes

[TABLE]

We set $N\geq 1$ and approximate the solution components $x_{2}$ and $x_{3}$ by continuously connected piecewise polynomials of degree $N$ and the component $x_{1}$ by possibly discontinuous piecewise polynomials of degree $N-1$ on uniform partitions of the interval $[0,1]$ with stepsize $h=1/n$. On each subinterval we choose $M$ collocation points.

The classical or standard collocation procedures use $M=N$ collocation points per subinterval (cf. [13]). This results in $3nN$ equations to determine the $3nN+2$ parameters of the unknown collocation solution. In order to obtain a unique collocation solution, we pose additionally two consistent initial condition.

In contrast, choosing $M>N$ leads to an overdetermined collocation system which can be treated by a least-squares solver. Corresponding first experiments are reported in [10] with $M=2N+1$. Table 1 shows the componentwise maximal error in both versions: The standard collocation generates waste as expected, however, the overdetermined least-squares collocation provides surprisingly nice results.

Further experiments (cf. [10, 9]) give rise to the conjecture that a much smaller number $M\geq N+1$ will do in general and that the special position of the collocation points does not matter. Tables 2 and 3 are quoted from [9] to this effect. At this place it should be noted that the overdetermined collocation method is treated in [10, 9] against the background of a least-squares problem in Hilbert spaces. For this reason, the errors are now measured in $L^{2}$ and $H^{1}$ norms. The theoretically confirmed convergence order is $N-\mu+1=1$, but we observe order $2$. To date there is no theoretical recognition of this impressive, nice behavior!

Table 3 indicates the further interesting observation that even though no convergence is proved for $N=1$ the numerical approximations remain bounded and seem to converge with order 0.4.

The bounded operator $T:\{x\in\mathcal{C}([0,1],\mathbb{R}^{3}):x_{2},x_{3}\in\mathcal{C}^{1}([0,1],\mathbb{R})\}\rightarrow\mathcal{C}([0,1],\mathbb{R}^{3})$ associated with our test DAE reads in detail

[TABLE]

and one immediately checks that

[TABLE]

further $\mathcal{T}=T$, since $T$ itself is injective, and hence, no initial or boundary conditions are allowed. The inverse operator

[TABLE]

is unbounded in this setting since ${\rm im}\,T$ is a nonclosed subset in $\mathcal{C}([0,1],\mathbb{R}^{3})$. This is a fundamental contrast to the case of regular ODEs. Now we do not have closed range and fredholm properties, and the main ingredient of the inverse of $\mathcal{T}$ is not a nice Volterra operator but a higher-order differential operator. We have

[TABLE]

Analogous situations arise also, e.g., in the settings

[TABLE]

with their natural norms, [14].

3 Regular first-order DA operators in a Hilbert space setting

3.1 Setting

We begin this part by considering the operator, henceforth called a DA operator,

[TABLE]

so that the operator equation $\mathring{T}x=g$ represents the so-called standard form DAE $Ex^{\prime}+Fx=g$. The coefficient functions $E,F:[a,b]\rightarrow\mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{m})$ are at least continuous. The function spaces $X$ and $Y$ will be specified later.

We suppose a nontrivial leading coefficient $E$ the nullspace of which is a $\mathcal{C}^{1}$-subspace varying in $\mathbb{R}^{m}$. We are looking for a Hilbert-space setting with a bounded DA operator.

First we try $X=Y=L^{2}((a,b),\mathbb{R}^{m})$, with the usual norms. In this setting, $\mathring{T}$ is unbounded, but densely defined and closable, see [14]. The closure of $\mathring{T}$, $T:{\rm dom}\,T\subseteq X\rightarrow Y$, is densely defined and closed, but also unbounded. We apply the usual graph-norm approach: The space $X_{T}:={\rm dom}\,T$ equipped with graph-norm,

[TABLE]

is complete and $T:X_{T}\rightarrow Y$ is a bounded operator. How can we specify ${\rm dom}\,T$ and $T$? To answer this question we need so-called proper factorizations of $E$.

Definition 3.1**.**

The factorization $E=AD$ is called proper, if $0<k\leq m$,  $A:[a,b]\rightarrow\mathcal{L}(\mathbb{R}^{k},\mathbb{R}^{m})$ is continuous, $D:[a,b]\rightarrow\mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{k})$ is continuously differentiable, and

[TABLE]

There are many possible proper factorizations. We fix an arbitrary one $E=AD$, put $B=F-AD^{\prime}$, and observe that

[TABLE]

Observe that $A(Dx)^{\prime}+Bx=g$ is a DAE with so-called properly stated leading term. This indicates how the closure $T$ as well as ${\rm dom}\,T$ look like.

Theorem 3.2**.**

1. (i)

If $D:[a,b]\rightarrow\mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{k})$ is continuously differentiable and has constant rank, then the function space

[TABLE]

equipped with the inner product

[TABLE]

is a Hilbert space. 2. (ii)

For each proper factorization $E=AD$, it results that

[TABLE] 3. (iii)

The norm $\lVert\cdot\lVert_{H^{1}_{D}}$ is equivalent to the graph-norm $\rVert\cdot\lVert_{T}$.

Proof.

(i) is an immediate consequence of [14, Lemma 6.9].

(ii) The closure $T=\mathring{T}^{**}$ is provided by means of the biadjoint of $\mathring{T}$ in [14, Theorem 3.1 (3)].

(iii) Let $x\in{\rm im}\,T$ be given. From $A(Dx)^{\prime}=Tx-Bx$ we obtain $A^{-}A(Dx)^{\prime}=A^{-}Tx-A^{-}Bx$. Owing to the proper factorization we may deduce the representation $(Dx)^{\prime}=A^{-}Tx+((A^{-}A)^{\prime}D-A^{-}B)x$. This yields

[TABLE]

On the other hand we have simply $\lVert x\rVert_{T}=\lVert x\rVert_{L^{2}}+\lVert Tx\rVert_{L^{2}}\leq(1+\lVert T\rVert)\lVert x\rVert_{H^{1}_{D}}$. ∎

We emphasize that ${\rm dom}\,T$ and $T$ remain invariant under proper refactorizations $AD=\bar{A}\bar{D}$ since the closure of a densely defined closable operator is unique.111Note that this is at the same time consistent with [12, Theorem 2.21].

3.2 Regular DA operators

We adopt the regularity notion in [12, Sections 2.4.2 and 2.7]. Regularity of a linear DAE is solely a matter of its coefficients $E,F$ and $A,D,B$, respectively. The associated coefficient pair $(E,F)$ is regular on $[a,b]$ exactly if one (equivalently each) associated proper triple $(A,D,B)$ resulting from a proper factorization is regular on $[a,b]$. This is in full accordance with the fact that the closure $T$ does not depend on the special factorization.

Definition 3.3**.**

The DA operators $\mathring{T}$ and $T$ are regular with tractability index $\mu\in\mathbb{N}$ and characteristic values $0<r_{0}\leq\cdots\leq r_{\mu-1}<r_{\mu}=m$, if the associated DAE is regular on $[a,b]$ with these characteristics.

A regular DA operator is said to be fine, if the coefficients $A,D,B$ are smooth enough for the existence of completely decoupling projectors.

Definition 3.3 is consistent with [9, Definition 2.1]. The operator $\mathring{T}$ and its closure $T$ share their index and characteristic values. Recall that $l:=m-\sum^{\mu-1}_{j=0}(m-r_{j})\geq 0$ is the dynamical degree of freedom of the DAE. $l=0$ may happen, see, e.g., the example DAE in Section 2.

A regular on $[a,b]$ DAE is associated with the so-called canonical projector function

[TABLE]

which can be seen as generalization of the spectral projection for regular matrix pencils. Each IVP for the matrix-valued function $X$,

[TABLE]

is uniquely solvable. We refer to [12, Section 2.4], [13, Section 2.2],[14, Section 4.2] for details, but note that, owing to the continuity of the coefficients, the solution of the above IVP, $X(\cdot,c)$ is well-defined and continuous with continuously differentiable part $DX(\cdot,c)$. $X(t,c)$ is called maximal-size fundamental solution matrix normalized at $c$. It has constant rank $l<m$, and it may happen that $l=0$, see Section 2. This is in contrast to the case of regular ODEs.

In the following we focus on bounded DA operators given in their natural Hilbert spaces,

[TABLE]

If $T$ is regular, then

[TABLE]

Introduce the operator $T_{IC}$ to capture initial conditions as well as the composed operator $\mathcal{T}$ by

[TABLE]

Thereby we suppose that $G_{a}\in\mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{l})$ and ${\rm ker}\,D(a)\subseteq{\rm ker}\,G_{a}$. The latter condition ensures the relation $G_{a}=G_{a}D(a)^{+}D(a)$. Then owing to the continuous embedding $H^{1}\hookrightarrow\mathcal{C}$ the operator $T_{IC}$ is well-defined and bounded.

Next we adopt the notion of accurately stated initial condition [13, Definition 2.3] accordingly.

Definition 3.4**.**

The operator $T_{IC}$ is accurately stated if ${\rm im}\,\mathcal{T}={\rm im}\,T\times\mathbb{R}^{l}$ and the composed operator $\mathcal{T}$ is injective.

Owing to [13, Corollary 2.2], if $T$ is regular, then $T_{IC}$ is accurately stated exactly if ${\rm ker}\,G_{a}\cap{\rm im}\,\mathnormal{\Pi}_{can}(a)=\{0\}$.

Supposing $T$ to be fine and $T_{IV}$ to be accurately stated, the inverse of the composed operator $\mathcal{T}$ can be represented as (e.g.,[13, Section 2.2], [12, Section 2.6]).

[TABLE]

in which the functions $v_{i}$ are successively explicitly determined by simple multiplications with certain matrix coefficients, by differentiation of $Dv_{j}$, $j=i+1,\ldots,\mu-1$, and subsequent linear combinations,

[TABLE]

The coefficients $\mathcal{N}_{i,k},\mathcal{M}_{i,k},\mathcal{L}_{i}$, and $G_{\mu}$ are at least continuous. They are fully determined by $A,D,B$ via a sequence of admissible matrix functions corresponding to a complete decoupling.

With the representation (2) of the inverse $\mathcal{T}^{-1}$ we intend to emphasize, on the one hand, the partial resemblance to (1). On the other hand, the second term $\mathfrak{D}g$ emerge for DA operators only. It is a differential operator, with may be higher order, see Section 2. If $T$ has index $\mu>1$ then $\mathcal{T}^{-1}$ includes derivatives up to order $\mu-1$.

The representation (2) shows that the inverse of the composed operator actually decomposes into two parts. The first “good” part is close to (1), i.e., the case of regular ODEs. This part may disappear. The second part is always present and representative for DA operators. Unfortunately, the canonical projector function which separates the parts is practically available in a few special cases only.

We close this subsection by quoting further relevant results from [9, Section 2]

Theorem 3.5**.**

Let $T$ be fine with index $\mu\in\mathbb{N}$ and $T_{IC}$ be accurately stated.

1. (i)

If $\mu=1$, then ${\rm im}\,T=L^{2}(\mathcal{I},\mathbb{R}^{m})$, and the composed operator $\mathcal{T}$ is a homeomorphism. 2. (ii)

If $\mu>1$, then ${\rm im}\,T\subset L^{2}(\mathcal{I},\mathbb{R}^{m})$ is nonclosed and the inverse $\mathcal{T}^{-1}$ as well as the Moore-Penrose inverse $T^{+}$ are no longer continuous. Then the equation $\mathcal{T}x=(g,d)$ is essentially ill-posed in Tikhonov’s sense.

4 Justification of the overdetermind polynomial collocation

In this section we deal with regular higher-index DA operators $T$, the related composed operator $\mathcal{T}$ and their approximations $R_{\pi,M}TU_{\pi}$ and $\mathcal{R}_{\pi,M}\mathcal{T}U_{\pi}$. On the background of the corresponding properties we provide new sufficient convergence conditions for the overdetermined polynomial collocation.

4.1 Basic technicalities

We consider the linear IVP or BVP,

[TABLE]

with the constant matrix $D=[I\;0]\in\mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{k})$, ${\rm rank}\,D=k$, and at least continuous matrix coefficients $A:[a,b]\rightarrow\mathcal{L}(\mathbb{R}^{k},\mathbb{R}^{m}),\;B:[a,b]\rightarrow\mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{m})$. The DAE (5) is supposed to be fine in the sense of [12, Section 2.6] , with tractability index $\mu\in\mathbb{N}$ and dynamical degree of freedom $l\leq k$. Recall that $\mu>1$ necessarily implies $l<k$. The matrices $G_{a},G_{b}\in\mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{l})$ are supposed to satisfy the conditions

[TABLE]

Condition (6) is further supposed to be accurately stated in the sense of [13, Definition 2.3]), so that the problems

[TABLE]

are uniquely solvable for each $d\in\mathbb{R}^{l}$. In particular, the homogeneous linear BVP has the trivial solution only.

The function $g$ is assumed to be admissible, so that the DAE (5) is solvable. Then the BVP (5), (6) has exactly one solution $x_{*}$ to be approximated later on.

Following the ideas of [10, 9] we represent the BVP (5), (6) as operator equation $\mathcal{T}x=y:=(g,d)$ in Hilbert spaces by introducing the spaces

[TABLE]

equipped with the inner products

[TABLE]

and operators

[TABLE]

The operator $T_{BC}$ is well defined and bounded owing to condition (7) and the continuous embedding $H^{1}((a,b),\mathbb{R}^{k})\hookrightarrow\mathcal{C}([a,b],\mathbb{R}^{k})$. Then, the DA operator $T$ as well as the composed operator $\mathcal{T}$ are obviously bounded. Moreover, $\mathcal{T}$ is injective and ${\rm im}\,\mathcal{T}={\rm im}\,T\times\mathbb{R}^{l}$. At this place let us emphasize again that we focus our interest on higher-index DAEs, $\mu\geq 2$, but then ${\rm im}\,T$ is a nonclosed subset of $L^{2}$ and $\mathcal{T}^{-1}$ is an unbounded operator, cf. Section 3, also [10, 9, 14].

Given the partition

[TABLE]

with stepsizes $h_{j}=t_{j}-t_{j-1}$, minimal stepsize $h_{\pi,\,min}$, and maximal stepsize $h_{\pi}$, we denote by $\mathcal{C}_{\pi}=\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})$ the space of piecewise continuous functions having breakpoints merely at the mesh points. Note that the supremum-norm $\lVert\cdot\rVert_{\infty}$ is well-defined for the elements of $\mathcal{C}_{\pi}$.

Next we fix a number $N\geq 1$ and introduce the space of ansatz functions to approximate $x_{*}$ by piecewise polynomial functions,

[TABLE]

The finite-dimensional space $X_{\pi}$ is a closed subspace of $H^{1}_{D}$, and the latter decomposes into the topological sum $X_{\pi}\oplus X_{\pi}^{\bot}=H^{1}_{D}$. We agree upon that

[TABLE]

For later reference, the following norm in $X_{\pi}$ will be needed:

[TABLE]

The ansatz space $X_{\pi}$ has dimension $nNm+k$. Choosing values

[TABLE]

we specify $M\geq N+1$ collocation points per subinterval, i.e.,

[TABLE]

and are then confronted with the overdetermined collocation system of $nMm+l>nNm+k$ equations for providing an approximation $x\in X_{\pi}$, namely,

[TABLE]

As a matter of course, the choice $M>N$ goes along with an overdetermined system (11),(12) comprising more equations than unknowns. This is different from standard collocation methods for ODEs and index-1 DAEs, e.g., [13]. Here we treat the overdetermined collocation system in a least-squares sense.

Let $R_{\pi,M}:\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})\rightarrow\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})$ denote the restriction operator which assigns to $w\in\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})$ the piecewise polynomial $R_{\pi,M}w\in\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})$ of degree less than or equal to $M-1$ such that the interpolation conditions,

[TABLE]

are satisfied. We also assign to $w\in\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})$ the vector $W\in\mathbb{R}^{mMn}$,

[TABLE]

which yields (cf. [10, Subsection 2.3])

[TABLE]

with the matrix $\mathcal{L}$ being positive definite, symmetric and independent222The entries of $\mathcal{L}$ are fully determined by the corresponding $M$ Lagrangian basis polynomials, thus, by $M$ and $\tau_{1},\ldots,\tau_{M}$. of $h_{\pi}$. There are further constants $\kappa_{l},\kappa_{u}>0$ such that

[TABLE]

If $w\in\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})$ is of class $\mathcal{C}^{M}$ on each subinterval of the partition $\pi$, then

[TABLE]

Additionally, we introduce the restriction operator $\mathcal{R}_{\pi,M}:\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})\times\mathbb{R}^{l}\rightarrow\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})\times\mathbb{R}^{l}$ by

[TABLE]

The overdetermined least-squares collocation means now that we seek an element $\tilde{x}_{\pi}\in X_{\pi}$ minimizing the functional

[TABLE]

With $w=Tx-g$ we may represent

[TABLE]

which reveals that by minimizing $\psi_{\pi,M}(x)$ subject to $x\in X_{\pi}$ we actually provide a least-squares solution of the collocation system (11),(12). The mathematics behind is closely related to special properties of the restriction operator $R_{\pi,M}$ on the one hand, but on the other hand, to the problem to minimize the functional

[TABLE]

for which (16) serves as approximation.

4.2 The operators $R_{\pi,M}TU_{\pi}$ and $\mathcal{R}_{\pi,M}\mathcal{T}U_{\pi}$

We begin this section by providing useful norm inequalities. Regarding convergence properties for $h_{\pi}$ tending to zero we have in mind sequences of partitions. In favor for an easier reading we drop an extra labeling, but we thoroughly assure that the indicated constants do not depend on the partitions and stepsizes $h_{\pi}$ in fact. We allow partitions $\pi$ having quotients $h_{\pi}/h_{\pi,min}\leq r$, with a global bound $1\leq r<\infty$. A consequence of [7, Theorem 3.2.6] is that there exists a constant $c_{K}=c_{K}(r)$ such that

[TABLE]

for all functions $z\in\mathcal{C}_{\pi}([a,b],\mathbb{R}^{s})$ being a polynomial of degree less than or equal to $K$ on each subinterval of the partition $\pi$.

Lemma 4.1**.**

There is a constant $\kappa>0$ such that

[TABLE]

Proof.

The inequality $\lVert x\rVert_{H^{1}_{D}}^{2}\leq(b-a)\lVert x\rVert_{\mathcal{C}^{1}_{D}}^{2}$ is evident for all $x\in X_{\pi}$. On the other hand, owing to (19) each arbitrary $x\in X_{\pi}$ satisfies

[TABLE]

and finally

[TABLE]

∎

If $A$ and $B$ are constant matrices, then $Tx$ is piecewise polynomial with degree less than or equal to $N$ for $x\in X_{\pi}$. Owing to $M\geq N+1$ this leads to

[TABLE]

In general, the operators $R_{\pi,M}TU_{\pi}$ and thus $\mathcal{R}_{\pi,M}\mathcal{T}U_{\pi}$ are well-defined on $H^{1}_{D}$ since $TU_{\pi}x$ belongs to $\mathcal{C}_{\pi}$ for all $x\in H^{1}_{D}$.

Proposition 4.2**.**

Let the DA operator $T$ be fine with index $\mu\in\mathbb{N}$ and $T_{BC}$ be accurately stated. Let $M\geq N+1$ and let the entries of $A$ and $B$ be of class $\mathcal{C}^{M}$. Then the following assertions are valid for all sufficiently fine partitions $\pi$:

1. (i)

There is a constant $C_{AB1}$ such that

[TABLE] 2. (ii)

There is a constant $C_{AB2}$ such that

[TABLE] 3. (iii)

If additionally $M\geq N+\mu$, then there is a constant $C$ such that

[TABLE] 4. (iv)

If the entries of $A$ and $B$ are polynomials of degree less than or equal to $N_{AB}$ and $M\geq N+1+N_{AB}$, then (20) and (22) remain valid, too.

Proof.

(i) We choose a number $K\geq 1$, $K+N\leq M$, $K$ interpolations nodes $0<\tilde{\tau}_{1}<\cdots<\tilde{\tau}_{K}<1$ and provide piecewise entrywise polynomial approximations $\tilde{A}$ and $\tilde{B}$ of degree $K-1$ such that

[TABLE]

Then there are constants $c_{1},c_{2}$ such that

[TABLE]

The auxiliary operator $\tilde{T}:H^{1}_{D}\rightarrow L^{2}$

[TABLE]

is bounded, $\lVert\tilde{T}\rVert\leq\sqrt{2}c_{2}$ and owing to $M\geq N+K$ it holds that $R_{\pi,M}\tilde{T}U_{\pi}x=\tilde{T}U_{\pi}x$ for all $x\in H^{1}_{D}$.

Next, to each arbitrary $x\in X_{\pi}$ we form $v:=(\tilde{T}-T)x\in\mathcal{C}_{\pi}$ and the vector $V\in\mathbb{R}^{mMn}$,

[TABLE]

which yields (cf. (13))

[TABLE]

On the other hand, regarding Lemma 4.1 and $x=U_{\pi}x$ we estimate

[TABLE]

For each arbitrary $x\in H^{1}_{D}$ it follows that

[TABLE]

which verifies the inequality (21) with a suitable bound $C_{AB1}$. Then it also results that

[TABLE]

and the assertion is verified.

(ii) To each arbitrary $x\in X_{\pi}$ we set $w=A(Dx)^{\prime}+Bx$ and derive on each subinterval of the partition $\pi$ that

[TABLE]

Since $(Dx)^{(N+1)}$ and $x^{(N+1)}$ vanish identically, we obtain the inequality

[TABLE]

with a constant $c_{3}$ being determined by the coefficients $A$ and $B$, and their involved derivatives. By Lemma 4.1 this yields

[TABLE]

Owing to [9, Lemma 4.2] it follows that

[TABLE]

We emphasize that the $C^{*}_{i}$ are also constants independent of the partition and stepsize. Then there is a constant $c_{4}$ such that

[TABLE]

Finally we arrive at

[TABLE]

and further

[TABLE]

(iii) The Moore-Penrose inverse of $\mathcal{T}U_{\pi}$ satisfies, by [9, Theorem 4.1], the inequality

[TABLE]

with a positive constant $c_{\gamma}$.

Denote by $\mathcal{V}_{\pi}$ and $\mathcal{V}_{\pi,M}$ the orthoprojectors of $L^{2}\times\mathbb{R}^{l}$ onto ${\rm im}\,\mathcal{T}U_{\pi}$ and ${\rm im}\,\mathcal{R}_{\pi,M}\mathcal{T}U_{\pi}$, respectively. Assertion (ii) implies now

[TABLE]

and further, for sufficiently fine partitions,

[TABLE]

Next we represent

[TABLE]

The subspace ${\rm im}\,\mathfrak{A}={\rm im}\,\mathcal{T}U_{\pi}$ has finite dimension, ${\rm im}\,\mathfrak{E}\subseteq{\rm im}\,\mathfrak{A}$. Furthermore, it holds that ${\rm ker}\,\mathfrak{A}={\rm ker}\,\mathcal{T}U_{\pi}={\rm ker}\,U_{\pi}$ and ${\rm ker}\,\mathfrak{E}\supseteq{\rm ker}\,U_{\pi}$. By [8, Lemma A.2] it follows that

[TABLE]

Now the inclusions

[TABLE]

hold and, hence, ${\rm ker}\,\mathcal{R}_{\pi,M}\mathcal{T}U_{\pi}={\rm ker}\,U_{\pi}$. In the end we compute

[TABLE]

and hence $\lVert(\mathcal{R}_{\pi,M}\mathcal{T}U_{\pi})^{+}\rVert\leq 2\lVert\mathfrak{A}^{+}\rVert$.

(iv) This assertion is a direct consequence of the fact that $TU_{\pi}x$ is piecewise polynomial of degree less than or equal to $N+N_{AB}$. ∎

4.3 Error estimations

Recall that $x_{*}$ denotes the sought solution, i.e., $\mathcal{T}x_{*}=y$ for given $y=(g,d)\in{\rm im}\,T\times\mathbb{R}^{l}$. As descibed in Subsection 4.1 the overdetermined polynomial collocation actually means that we generate the minimizer of the functional $\phi_{\pi,M}$, cf. (17), that is,

[TABLE]

to approximate $x_{*}$. Now we provide a corresponding error estimation.

Suppose that the solution is smooth, $x_{*}\in\mathcal{C}^{N}([a,b],\mathbb{R}^{m})$, $Dx_{*}\in\mathcal{C}^{N+1}([a,b],\mathbb{R}^{m})$. With the $N$ nodes $0<\tau_{*\,1}<\cdots<\tau_{*\,N}<1$, the interpolating function $p_{*}\in X_{\pi}$ uniquely defined by

[TABLE]

satisfies the inequalities

[TABLE]

in which the constant $c_{*}$ is determined by $x_{*}$. Next, owing to Proposition 4.2(i), we may estimate

[TABLE]

Denoting $w_{*}=T(x_{*}-p_{*})\in\mathcal{C}_{\pi}$ and using the Lagrange basis polynomials we further derive

[TABLE]

Because of

[TABLE]

we obtain

[TABLE]

such that

[TABLE]

Proposition 4.2(iii) implies, for $M\geq N+\mu$, that $U_{\pi}=(\mathcal{R}_{\pi,M}\mathcal{T}U_{\pi})^{+}\mathcal{R}_{\pi,M}\mathcal{T}U_{\pi}$, which gives rise to the error representation

[TABLE]

so that

[TABLE]

At the end we arrive at

[TABLE]

We summarize the result in the following theorem.

Theorem 4.3**.**

Let the DA operator $T$ be fine with index $\mu\in\mathbb{N}$ and $T_{BC}$ be accurately stated. Let $M\geq N+\mu$ and let the entries of $A$ and $B$ be of class $\mathcal{C}^{M}$. Let $g\in{\rm im}\,T$, and the solution $x_{*}=\mathcal{T}^{-1}(g,d)$ be of class $\mathcal{C}^{N}$ with $Dx_{*}$ of class $\mathcal{C}^{N+1}$. Then there is a constant $\bar{C}$ such that the error estimation

[TABLE]

is valid for all sufficiently fine partitions $\pi$.

In contrast to Theorem 4.3, the earlier error estimation from [9, Theorem 3.1 (a)] is given for the least-squares approximation $x_{\pi}$ (cf. (18)),

[TABLE]

Then, supposing the entries of $A$ and $B$ to be polynomials of degree less than or equal to $N_{AB}$ and letting $M\geq N+1+N_{AB}$ the estimation (23) is derived by a different technique ([9, Theorem 5.1(a)]). In the present context of operator properties this means that then the operators $R_{\pi,M}TU_{\pi}$ and $TU_{\pi}$ coincide, see Proposition 4.2.

So far we know only sufficient convergence and order conditions. The question concerning an appropriate or even optimal choice of $N\geq 1$ and $M\geq N+1$ remains open. Of course, smaller $N$ and $M$ are associated with less computational effort. So far, in experiments $M=N+1$ works quite well. In general, the practical performance is much better than we can substantiate till now. Much further analysis is needed.

5 Higher-order DA operators

General linear order-$s$ DA operators have the form (e.g., [6, 4, 15])

[TABLE]

with at least continuous matrix-coefficients $E_{i}$ and a singular leading coefficient $E_{s}$. So far the consolidated knowledge of higher-order DAEs and the related operators is rather poor. Naturally the class of higher-order DA operators is much more complex than the class of first-order ones, nevertheless the numerical treatment of the corresponding DAEs might be often easier than expected. In particular, systems of ODEs of mixed order, which can be handled by traditional approved software packages such as COLNEW, COLDAE and BVPSUITE ([1, 2, 11]), actually apply to higher-order DAEs directly. For instance, the simple system

[TABLE]

corresponds to the DA operator $\mathring{T}:{\rm dom}\,\mathring{T}=\mathcal{C}^{2}([a,b],\mathbb{R}^{2})\subset\mathcal{C}([a,b],\mathbb{R}^{2})\rightarrow\mathcal{C}([a,b],\mathbb{R}^{2})$

[TABLE]

and its extension $T$, with ${\rm dom}\,T=\{z\in\mathcal{C}([a,b],\mathbb{R}^{2}):z_{1}\in\mathcal{C}^{2}([a,b],\mathbb{R}),z_{2}\in\mathcal{C}^{1}([a,b],\mathbb{R})\}$,

[TABLE]

One has ${\rm im}\,T=\mathcal{C}([a,b],\mathbb{R}^{2})$ and $\dim{\rm ker}\,T=3$, and hence $T$ is fredholm, so that IVPs and BVPs can be stated in a well-posed way. It arises the question which further DA operators allow a reliable direct treatment of the corresponding DAE. For the time being we are not able to present a general answer. Below, we survey certain related aspects.

In the early paper [3] higher-order Hessenberg DAEs arising from higher-order ODEs subject to constraints are introduced and analyzed with respect to their perturbation index. Written in the form (24) the associated operators are

[TABLE]

and, with properly involved derivatives by the additional matrix $D=\begin{bmatrix}I&0\end{bmatrix}$, ${\rm rank}\,D=m_{1}$,

[TABLE]

We set $X_{T}:={\rm dom}\,T$ and introduce the norm $\lVert x\rVert=\lVert x\rVert_{\infty}+\lVert Dx\rVert_{H^{s}}$, $x\in X_{T}$, so that $T:X_{T}\rightarrow\mathcal{C}([a,b],\mathbb{R}^{m})$ is bounded.

If $E_{0,22}$ remains nonsingular, then it is evident that

[TABLE]

i.e., $T$ is also fredholm. Therefore the operator $T$ is marked as index-1 operator independently of $s\geq 1$. The corresponding IVPs and BVPs are well-posed and can be treated even by standard polynomial collocation. By introducing new variables $v_{i}=(Dx)^{i},\;i=1,\ldots,s-1$, one obtains an associated first-order formulation of the DAE in the variables $v_{1},\ldots,v_{s-1},x$, which has the same index 1.

If $E_{0,22}$ vanishes identically, but $E_{s-1,21}E_{0,12}$ remains nonsingular, then

[TABLE]

and $T$ is marked as index-2 operator independently of $s\geq 1$. Although ${\rm im}\,T$ is a nonclosed subset in $\mathcal{C}([a,b],\mathbb{R}^{m})$, utilizing the special problem structure and incorporating special projections, corresponding BVPs with accurately stated boundary conditions can be directly treated by COLDAE, see [3]. Again, substituting the derivatives $(Dx)^{(i)}$ by new variables in the DAE leads to a first-order index-2 system which can be solved by overdetermined least-squares collocation, too.

If $E_{0,22}$ as well as $E_{1,21},\ldots,E_{s-1,21}$ vanish identically, but $E_{0,21}E_{0,12}$ remains nonsingular, then the operator $T$ has index $s+1$, and

[TABLE]

Also here, substituting the derivatives $(Dx)^{(i)}$ by new variables leads to an index-$(s+1)$ first-order DAE, and the latter can be treated by overdetermined least-squares collocation.

The index notion used in [3] is the perturbation index $\mu\in\mathbb{N}$ of a suitable first-order formulation333Not surprisingly, the classical procedure of turning a higher-order ODE into a first-order system applied to a DAE increases the differentiation index and leads to different solvability results and smoothness requirements. For details and examples we refer to [15].. In the operator context we consider the extension $T$ of $\mathring{T}:{\rm dom}\,\mathring{T}\subset\mathcal{C}([a,b],\mathbb{R}^{m})\rightarrow\mathcal{C}([a,b],\mathbb{R}^{m})$, turn then to the bounded operator $T:X_{T}\rightarrow\mathcal{C}([a,b],\mathbb{R}^{m})$ and show that the elements $g$ of ${\rm im}\,T$ are involved therein together with parts of derivatives up to order $\mu-1$, and $\mu$ is the smallest such number.

In contrast to the standard form (24), the version with properly involved derivatives can be seen as source of a reasonable first-order formulation. We conjecture that this idea applies also to further classes of DAEs. To emphasize the capabilities of this idea we use [15, Example] for a demonstration.

Example 5.1**.**

We refactorize the leading term of the second order operator $\mathring{T}$ given by

[TABLE]

by

[TABLE]

The resulting operator T reads

[TABLE]

*Introducing the new variable $v=(Dx)^{\prime}=x_{1}^{\prime}+(t+1)x_{2}^{\prime}+x_{2}$ leads to a first-order formulation with index $\mu=3$.444In [15] merely an index-4 first-order formulation was obtained with $v=x_{1}^{\prime}+(t+1)x_{2}^{\prime}$. *

A completely different approach to operators (24) and the corresponding DAEs is proposed in [6] in the context of a ring of operators acting on $\mathcal{C}^{\infty}([a,b],\mathbb{R}^{m})$, among them arbitrary-order differential operators of the form (24) with real-analytic coefficients555Note that the usual semi-norm family defining the topology of $\mathcal{C}^{\infty}$ is much too strong to measure practically relevant approximation errors, e.g., [14, Subsection 2.4.2]. This serves as background of the following index notion by means of left regularizing operators

Definition 5.2**.**

Each operator $Ly=\sum_{i=0}^{k}F_{i}y^{(i)}$, $y\in\mathcal{C}^{k}([a,b],\mathbb{R}^{m})$, such that the superposition $L\circ\mathring{T}$ is a regular ODE of the same order as $\mathring{T}$, i.e.,

[TABLE]

is called left regularizing operator of $\mathring{T}$. The smallest possible $k$ is said to be the index of $\mathring{T}$.

Here we denote the index in the sense of Definition 5.2 by $\mu_{C}$. The construction of left regularizers is closely related to the evaluation of derivative arrays, see [6]. For $s=1$, constructing a left regularizer is equivalent to providing a so-called completion ODE ([5]) and $\mu_{C}$ equals the differentiation index. In turn, for regular first-order DAEs the differentiation index equals the perturbation index and also the tractability index as well.

However, for $s>1$ things are completely different, so that $\mu_{C}$ is no longer helpful in view of the practical treatment of higher-order DAEs. The following two simple examples allow a first insight. Both examples have the form (25) resp. (26). We compare the perturbation index $\mu$ applied in [3] and $\mu_{C}$.

Example 5.3** ($\mu=1$,  $\mu_{C}=s$).**

[TABLE]

Example 5.4** ($\mu=s+1$,  $\mu_{C}=2s$).**

[TABLE]

We finish by mentioning that in [15, 16] the given DAE is transformed via derivative arrays to a so-called strangeness-free mixed-order system which then can be handled by standard methods. Under additional quite special conditions, such a system is provided by evaluating involved matrix polynomials in [4].

6 Conclusions

We have explored properties of regular first-order DA operators and their finite-dimensional counterparts associated with the polynomial overdetermined least-squares collocation and provided an new convergence result.

We notice substantial progess in view of the consolidation of the polynomial overdetermined least-squares collocation. Nevertheless, there are essential open questions, e.g., concerning the choice of $N$ and $M$.

Furthermore we have surveyed corresponding results concerning higher-order operators and the direct treatment of higher-order DAEs.