A least-squares collocation method for nonlinear higher-index differential-algebraic equations

TL;DR

This paper presents a novel least-squares collocation method for solving nonlinear higher-index differential-algebraic equations, demonstrating promising numerical results and establishing initial convergence theory for nonlinear cases.

Contribution

It introduces a direct numerical approach using overdetermined polynomial least-squares collocation for nonlinear higher-index DAEs, with the first convergence proof for nonlinear problems.

Findings

Numerical experiments show impressive results.

The method is not much more expensive than standard collocation.

First convergence proof for nonlinear problems.

Abstract

We introduce a direct numerical treatment of nonlinear higher-index differential-algebraic equations by means of overdetermined polynomial least-squares collocation. The procedure is not much more computationally expensive than standard collocation methods for regular ordinary differential equations. The numerical experiments show impressive results. In contrast, the theoretical basic concept turns out to be considerably challenging. So far, quite recently convergence proofs for linear problems have been published. In the present paper we come up to a first convergence result for nonlinear problems.

Tables

Table 1. Table 1: Errror in L 2 ​ ( 0 , 1 ) superscript 𝐿 2 0 1 L^{2}(0,1) for N = 3 𝑁 3 N=3 for the pendulum example. n 𝑛 n equidistant grid points and M = N + 1 𝑀 𝑁 1 M=N+1 uniformly distributed collocation points have been used

$n$	$x$	$x^{\prime }$	$y$	$y^{\prime }$	$\lambda$
10	$4.42e-02$	$1.17e-01$	$1.83e-02$	$1.01e-01$	$6.25e-01$
20	$6.01e-03$	$1.76e-02$	$2.48e-03$	$1.98e-02$	$3.33e-01$
40	$8.28e-04$	$3.07e-03$	$3.41e-04$	$4.47e-03$	$1.72e-01$
80	$1.11e-04$	$6.26e-04$	$4.59e-05$	$1.07e-03$	$8.67e-02$
160	$1.42e-05$	$1.44e-04$	$5.87e-06$	$2.64e-04$	$4.34e-02$
320	$1.86e-06$	$3.50e-05$	$7.65e-07$	$6.58e-05$	$2.17e-02$
640	$2.32e-07$	$8.68e-06$	$9.57e-08$	$1.64e-05$	$1.08e-02$

Table 2. Table 2: Order estimate for N = 3 𝑁 3 N=3 for the pendulum example. n 𝑛 n equidistant grid points and M = N + 1 𝑀 𝑁 1 M=N+1 uniformly distributed collocation points have been used

$n$	$x$	$x^{\prime }$	$y$	$y^{\prime }$	$\lambda$
10	$2.9$	$2.7$	$2.9$	$2.3$	$0.9$
20	$2.9$	$2.5$	$2.9$	$2.1$	$1.0$
40	$2.9$	$2.3$	$2.9$	$2.1$	$1.0$
80	$3.0$	$2.1$	$3.0$	$2.0$	$1.0$
160	$2.9$	$2.0$	$2.9$	$2.0$	$1.0$
320	$3.0$	$2.0$	$3.0$	$2.0$	$1.0$

Table 3. Table 3: Errror in L 2 ​ ( 0 , 1 ) superscript 𝐿 2 0 1 L^{2}(0,1) for N = 5 𝑁 5 N=5 for the pendulum example. n 𝑛 n equidistant grid points and M = N + 1 𝑀 𝑁 1 M=N+1 uniformly distributed collocation points have been used

$n$	$x$	$x^{\prime }$	$y$	$y^{\prime }$	$\lambda$
10	$4.13e-04$	$1.28e-03$	$1.75e-04$	$1.99e-03$	$3.61e-02$
20	$4.59e-05$	$1.22e-04$	$1.88e-05$	$1.38e-04$	$4.90e-03$
40	$1.45e-06$	$4.44e-06$	$5.94e-07$	$6.94e-06$	$6.18e-04$
80	$3.43e-08$	$1.82e-07$	$1.41e-08$	$3.97e-07$	$7.74e-05$
160	$1.02e-09$	$1.04e-08$	$4.17e-10$	$2.45e-08$	$9.68e-06$
320	$5.57e-11$	$6.41e-10$	$2.28e-11$	$1.52e-09$	$1.21e-06$

Table 4. Table 4: Order estimate for N = 5 𝑁 5 N=5 for the pendulum example. n 𝑛 n equidistant grid points and M = N + 1 𝑀 𝑁 1 M=N+1 uniformly distributed collocation points

$n$	$x$	$x^{\prime }$	$y$	$y^{\prime }$	$\lambda$
10	$3.2$	$3.4$	$3.2$	$3.8$	$2.9$
20	$5.0$	$4.8$	$5.0$	$4.3$	$3.0$
40	$5.4$	$4.6$	$5.4$	$4.1$	$3.0$
80	$5.1$	$4.1$	$5.1$	$4.0$	$3.0$
160	$4.2$	$4.0$	$4.2$	$4.0$	$3.0$

Table 5. Table 5: Errors in H D 1 ​ ( 0 , 5 ) subscript superscript 𝐻 1 𝐷 0 5 H^{1}_{D}(0,5) for ( 5.2 ) using M = N + 1 𝑀 𝑁 1 M=N+1 for the Campbell-Moore example. n 𝑛 n equidistant grid points and M 𝑀 M Gaussian collocation points have been used

$n$	$N=1$	$N=2$	$N=3$	$N=4$	$N=5$
10	3.32e+1	4.53e+0	3.82e-1	7.02e-2	1.47e-3
20	3.32e+1	7.51e-1	1.02e-1	1.26e-2	1.24e-4
40	3.32e+1	3.03e-1	3.14e-2	2.52e-3	1.30e-5
80	3.32e+1	1.80e-1	1.22e-2	5.45e-4	1.54e-6
160	3.32e+1	1.17e-1	5.67e-3	1.25e-4	1.20e-6
320	3.32e+1	7.95e-2	2.73e-3	1.25e-4	1.20e-6

Table 6. Table 6: Order estimation for ( 5.2 ) using M = N + 1 𝑀 𝑁 1 M=N+1 for the Campbell-Moore example. n 𝑛 n equidistant grid points and M 𝑀 M Gaussian collocation points have been used. The row “theory” contains the expected orders. Note that Theorem  4.2 is only valid for N ≥ 5 𝑁 5 N\geq 5

$n$	$N=1$	$N=2$	$N=3$	$N=4$	$N=5$
20	0.0	2.6	1.9	2.5	3.6
40	0.0	1.3	1.7	2.3	3.3
80	0.0	0.7	1.4	2.2	3.1
160	0.0	0.6	1.1	2.1	0.6
320	0.0	0.6	1.1	0.0	0.0
theory		(0)	(1)	(2)	3

Equations

$$f((Dx{)}^{\prime }(t),x(t),t)$$

$$g(x(a),x(b))$$

$$g(u,v)=g(D^{+}Du,D^{+}Dv),u,v\in R^m,$$

$$A_{\ast }(t)(Dz{)}^{\prime }(t)+B_{\ast }(t)z(t)$$

$$G_{\ast a}z(a)+G_{\ast b}z(b)$$

$$A_{\ast }(t):=f_y((D{x}_{\ast }{)}^{\prime }(t),x_{\ast }(t),t),B_{\ast }(t):=f_x((D{x}_{\ast }{)}^{\prime }(t),x_{\ast }(t),t),t\in [a,b],$$

$$G_{\ast a}:=g_u(x_{\ast }(a),x_{\ast }(b)),G_{\ast b}:=g_v(x_{\ast }(a),x_{\ast }(b)).$$

$$\ker D\subseteq \ker G_{\ast a},\ker D\subseteq \ker G_{\ast b}.$$

$$A_{\ast }(t)(Dz{)}^{\prime }(t)+B_{\ast }(t)z(t)=0,t\in [a,b],G_{\ast a}z(a)+G_{\ast b}z(b)=d,$$

$$\underset{t\in [a,b]}{\max }|z(t)|\le {\kappa }_{BC}|d|,$$

$$\pi :a=t_0<t_1<\cdots <t_n=b,$$

$$X_{\pi }=\{$$

$$x_{\kappa }{|}_{[t_{j-1},t_j)}\in {\mathfrak{P}}_N,\kappa =1,\dots ,k,x_{\kappa }{|}_{[t_{j-1},t_j)}\in {\mathfrak{P}}_{N-1},\kappa =k+1,\dots ,m,j=1,\dots ,n\}.$$

$$0<{\tau }_1<\cdots <{\tau }_M<1$$

$$t_{ji}=t_{j-1}+{\tau }_i{h}_j,i=1,\dots ,M,j=1,\dots ,n,$$

$$f((Dx{)}^{\prime }(t_{ji}),x(t_{ji}),t_{ji})$$

$$g(x(t_0),x(t_n))$$

$$(R_{\pi ,M}w)(t_{ji})=w(t_{ji}),i=1,\cdots ,M,j=1,\cdots ,n,$$

$$W=\left[\begin{array}{c}{W}_1\\ ⋮\\ W_n\end{array}\right]\in {\mathbb{R}}^{mMn},W_j={\left(\frac{h_j}{M}\right)}^{1/2}\left[\begin{array}{c}w(t_{j1})\\ ⋮\\ w(t_{jM})\end{array}\right]\in {\mathbb{R}}^{mM},$$

$$\Vert R_{\pi ,M}w{\Vert }_{L^2}^2=W^T\mathcal{L}W,w\in {\mathcal{C}}_{\pi }([a,b],{\mathbb{R}}^m),$$

$${\psi }_{\pi ,M}(x)$$

$$=\Vert R_{\pi ,M}{w}_f{\Vert }^2+|g(x(a),x(b)){|}^2,x\in X_{\pi }.$$

$$c_{\mathcal{L}}|W{|}^2=c_{\mathcal{L}}{W}^{T}W\le W^T\mathcal{L}W\le C_{\mathcal{L}}{W}^{T}W=C_{\mathcal{L}}|W{|}^2,W\in {\mathbb{R}}^{mMn},$$

$$\psi (x)=\Vert w_f{\Vert }^2+|g(x(a),x(b)){|}^2\text{subject to }x\in X_{\pi },$$

$$H_D^1=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},x,\bar{x}\in H_D^1.$$

$$\Vert Dx{\Vert }_{\infty }\le \kappa \Vert Dx{\Vert }_{H^1}\le \kappa \Vert x{\Vert }_{H_D^1},x\in H_D^1.$$

$$F$$

$$(Fx)(t)$$

$$F_{BC}x$$

$$\mathcal{F}x$$

$$T$$

$$(Tx)(t)$$

$$T_{BC}x$$

$$\mathcal{T}x$$

Full text

1 Introduction

For regular ordinary differential equations and index-1 differential-algebraic equations standard collocation methods which rely on closed discretized systems111The number of unknowns equals the number of equations. are known to work well. Moreover, Hessenberg form index-2 differential-algebraic equations can be treated successfully by so-called projected collocation methods that complement standard collocation with an additional updating of the differential solution component by a projection step. This goes along with the well-posedness of the related initial and boundary value problems in natural settings; we refer to [12] for a detailed survey. In contrast, higher-index differential-algebraic equations lead to ill-posed222More precisely: Essentially ill-posed in Tichonov’s sense, that is, the related operators feature nonclosed ranges. initial and boundary value problems, and standard collocation methods necessarily fail unless an elaborate index-reducing preprocessing is incorporated, which utilizes derivative array systems.

Recently ([7, 8]) first promising experiments concerning an least-squares overdetermined polynomial collocation directly applied to the DAE without any preprocessing have been reported. The theoretical justification appears to be quite challenging. So far, only sufficient convergence conditions are obtained for linear problems [7, 6, 8]. In the present paper we provide a first proof for nonlinear problems.

The paper is organized as follows: In Section 2 we state the problem in detail. Then we provide a Hilbert space setting in Section 3. This setting is more comfortable for the treatment of the given ill-posed problems. In Section 4, we introduce and investigate a kind of Newton-iteration related to a single partition, which uses bounded outer inverses as discussed in [15] and which serves in the end as background for the Gauss-Newton iteration applied to an overdetermined collocation system. Then, we consider nested multiple partitions to ensure convergence of the iteration-projection method in Section 6. The examples in Section 5 confirm the capability of the approach, but, having said that, they also indicate that our sufficient convergence conditions seem to be too unsubtle still. Finally, we provide some remarks and conclusions.

We use the symbol $\rVert\cdot\lVert$ for different function and operator norms. In general, in the given context things will be unambiguous. Only on certain places, to prevent maybe imminent confusions we indicate the special norms by the corresponding subscripts, e.g., $\rVert\cdot\lVert_{L^{2}}$.

Some notations and abbreviations

[TABLE]

2 The issue and basic technicalities

We deal with IVPs and BVPs given in the form

[TABLE]

with $[a,b]$ being a compact interval, $D=[I\;0]\in\mathcal{L}(\mathbb{R}^{m},\mathbb{R}^{k})$, ${\rm rank}\,D=k$, and data $f(y,x,t)\in\mathbb{R}^{m},y\in\mathbb{R}^{k},x\in\mathcal{D}_{x}\subseteq\mathbb{R}^{m},t\in\mathcal{D}_{t}\subset\mathbb{R}$, $\mathcal{D}_{t}\supset[a,b]$, $g(u,v)\in\mathbb{R}^{l},u,v\in\mathbb{R}^{m}$. The functions $f$ and $g$ are supposed to be at least continuous together with their partial derivatives $f_{y},f_{x},g_{u},g_{v}$.

We assume that the BVP (1), (2) has the solution $x_{*}:[a,b]\rightarrow\mathbb{R}^{m}$ to be approximated. $x_{*}$ is supposed to be continuous with continuously differentiable part $Dx_{*}$. Later on, among others for obtaining convergence orders, additional smoothness will be required.

Moreover, the DAE (1) is supposed to be regular with (tractability) index $\mu\in\mathbb{N}$ and characteristics $0<r_{0}\leq\cdots\leq r_{\mu-1}<r_{m}=m$ around $x_{*}$, that means, the graph $\{(x_{*}(t),t):t\in[a,b]\}$ resides within a regularity region having these characteristics (e.g., [11, Definition 3.28]). Note that then the derivative $(Dx)^{\prime}$ is properly involved in the DAE (1) so that $f_{y}(y,x,t)$ has full column-rank $k$.

Furthermore, in condition (2), we apply $l=m-\sum_{i=0}^{\mu-1}(m-r_{i})\geq 0$ which is the dynamical degree of freedom of the DAE. Recall that regular ODEs are indicated by $l=k=m$, regular index-1 DAEs by $l=k<m$, but higher-index DAEs by $l<k<m$. We are mainly interested in the last case. We further suppose the function $g$ to satisfy the relation

[TABLE]

so that the initial or boundary condition (2) actually applies to the differentiable component $Dx$ only.

Together with the BVP (1),(2) we consider the linear BVP,

[TABLE]

with

[TABLE]

We assume the solution $x_{*}$ and possibly the data $f$ to be sufficiently smooth so that the linearized DAE (4) is fine in the sense of [11, Section 2.6]. Since the solution $x_{*}$ resides in a regularity region of the DAE (1), the linear DAE (4) inherits the characteristic values and the index $\mu$ of the nonlinear DAE, see [14, Page 279]. Furthermore, owing to condition (3) it holds that

[TABLE]

Condition (2) is supposed to be stated in such a way that the linear BVP (4),(5) features accurately stated boundary condition in the sense of [12, Definition 2.3]), meaning that the problems

[TABLE]

are uniquely solvable for each $d\in\mathbb{R}^{l}$, and the solutions satisfy the inequality

[TABLE]

with a constant $\kappa_{BC}$. In particular, the homogeneous linear BVP, that is, the so-called variational problem, has then the trivial solution only.

Given the partition

[TABLE]

with stepsizes $h_{j}=t_{j}-t_{j-1}$, maximal stepsize $h_{\pi}$, and minimal stepsize $h_{\pi,\min}$. Denote by $\mathcal{M}_{[r]}$ the set of all partitions $\pi$ the ratio of the maximal stepsize by the minimal stepsize of which is uniformly bounded by the constant $r<\infty$.

Let $\mathcal{C}_{\pi}([a,b],\mathbb{R}^{m})$ denote the space of piecewise continuous functions having breakpoints merely at the mesh points.

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

[TABLE]

This ansatz space has dimension $nNm+k$. Choosing values

[TABLE]

we specify $M$ collocation points per subinterval, i.e.,

[TABLE]

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

[TABLE]

The choice $M=N$ corresponds to the standard polynomial collocation yielding $nNm+l$ equations, which works well for regular ODEs and index-1 DAEs, with dynamical degree $l=k=m$ and $l=k<m$, respectively (cf. [12]). In contrast, higher-index DAEs feature always a dynamical degree $0\leq l<k<m$. As it is well-known, completing the collocation system by additional $k-l$ consistent boundary conditions does not result in a suitable method owing to the ill-posedness of the higher-index problem, e.g., [7, Example 1.1]. As a matter of course, the choice $M>N$ goes along with an overdetermined system (10),(11) comprising more equations than unknowns.

Here we always set $M>N$ and treat the overdetermined collocation system in a least-squares sense. More precisely, 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 $M$ 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. [8, Subsection 3])

[TABLE]

with a positive definite, symmetric matrix $\mathcal{L}$. The entries of $\mathcal{L}$ do not at all depend on the partition $\pi$. They are fully determined by the corresponding $M$ Lagrange basis polynomials.

Letting $w_{f}(t)=f((Dx)^{\prime}(t),x(t),t),\;t\in[a,b]$, we introduce the functional

[TABLE]

The overdetermined least-squares collocation means now that we seek an element $\tilde{x}_{\pi}$ making the value $\psi_{\pi,M}(\tilde{x}_{\pi})$ as small as possible. Note that there are positive constants $c_{\mathcal{L}}$, $C_{\mathcal{L}}$ such that

[TABLE]

which justifies the labeling least squares collocation. We refer to [8, 7] for a number of promising numerical experiments, see also Section 5. Expression (12) serves to indicate the basic numerical procedure, whereas formula (13) suggests that 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 (13) serves as approximation. We refer to [6] for properties of the restriction operator in this context. The objective of the present paper is to contribute to the background problem (14).

3 Hilbert space setting

Following the ideas of [8, 7] concerning linear problems, we investigate also the nonlinear problem (1),(2) described in Section 1 as operator equation $\mathcal{F}x=0$ in a Hilbert space setting, which is most comfortable for treating ill-posed problems. Besides standard function spaces such as $L^{2}$, $H^{1}$, $\mathcal{C}$, etc., equipped with usual inner products and norms, we use the space

[TABLE]

equipped with the inner product

[TABLE]

$H^{1}_{D}$ is a Hilbert space, [14, Lemma 6.9]. Owing to the continuous embedding $H^{1}((a,b),\mathbb{R}^{k})\hookrightarrow\mathcal{C}([a,b],\mathbb{R}^{k})$, e.g., [1, Theorem 0.4], $x\in H^{1}_{D}$ implies $Dx\in C([a,b],\mathbb{R}^{k})$, and it holds

[TABLE]

We introduce the nonlinear operators $F,F_{BC},$ and $\mathcal{F}$,

[TABLE]

[TABLE]

as well as the linear operators $T,T_{BC},$ and $\mathcal{T}$,

[TABLE]

[TABLE]

We are merely interested in the local behavior of $F$ and $\mathcal{F}$ and suppose

[TABLE]

Regarding condition (3) as well as (15), we find the operators $F_{BC}$ and $T_{BC}$ well defined. $F_{BC}$ is Fréchet-differentiable, which can be checked by straightforward computation. In particular, $F_{BC}(x_{*})=T_{BC}$. Moreover, supposing the partial derivatives $g_{u},g_{v}$ to be Lipschitz continuous, there is a constant $L_{BC}$ such that

[TABLE]

The linear operators $T$ and $\mathcal{T}$ are obviously bounded. The operator $F$ is closely related to a certain Nemyckij operator as Proposition 3.1 below indicates. In the convergence proofs we will need that $F$ and thus $\mathcal{F}$ are Gâteaux-differentiable on their domain with uniformly bounded Gâteaux-derivatives,

[TABLE]

Proposition 3.1 provides sufficient conditions to justify these assumptions.

Moreover, we will need the inequality

[TABLE]

to be valid with a constant $L_{F}$ for the Gâteaux-derivative $F^{\prime}$ where $X_{\pi}$ is given by (2). Proposition 3.1 provides conditions also for this property to hold. Having (23), we are provided with a constant $L$ such that

[TABLE]

Note that $L_{F}$ and $L$ depend on the stepsize ratio $r$.

Now the BVP (1),(2) is represented by the operator equation $\mathcal{F}x=0$ and the least-squares functional (14) we are mainly interested in reads now

[TABLE]

By construction, one has $T=F^{\prime}(x_{*})$ and $\mathcal{T}=\mathcal{F}^{\prime}(x_{*})$. The equation $\mathcal{F}^{\prime}(x_{*})z=0$ represents the homogeneous variational BVP (7), with $d=0$, which has the trivial solution only. Therefore, the operator $\mathcal{F}^{\prime}(x_{*})$ is injective. At this place we emphasize again, that higher-index DAEs lead to ill-posed problems. In the context here this means that ${\rm im}\,F^{\prime}(x_{*})$ and ${\rm im}\,\mathcal{F}^{\prime}(x_{*})$ are nonclosed subsets in $L^{2}$ and $L^{2}\times\mathbb{R}^{l}$, respectively, see [7, Theorem 2.4], and the inverse $\mathcal{F}^{\prime}(x_{*})^{-1}$ is unbounded.

Proposition 3.1**.**

Let $f$ and $D$ be as described in Section 2, with $\mathcal{D}_{x}=\mathbb{R}^{m}$, $\mathcal{D}_{y}=\mathbb{R}^{k}$ and bounded partial derivatives $f_{y}$ and $f_{x}$.

1. (i)

Then, $x\in H^{1}_{D}$ implies $Fx\in L^{2}$, and $F$ is Gâteaux-differentiable, with the Gâteaux-derivative $F^{\prime}(x)$,

[TABLE]

Moreover, $F^{\prime}(x)$ is uniformly bounded. 2. (ii)

If, additionally, the partial derivatives $f_{x}$ and $f_{y}$ satisfy the inequalities

[TABLE]

for all $x_{1},x_{2}\in\mathbb{R}^{m}$ and $y_{1},y_{2}\in\mathbb{R}^{k}$, then there is a constant $L_{F}=L_{F}(r)$ such that

[TABLE]

that is (23).

Proof.

(i) Consider the operators $J:H^{1}_{D}\rightarrow L^{2}((a,b),\mathbb{R}^{k})\times L^{2}$ given by $Jx=((Dx)^{\prime},x)$ and $\tilde{F}:L^{2}((a,b),\mathbb{R}^{k})\times L^{2}\rightarrow L^{2}$ defined as the Nemytskij operator

[TABLE]

Under the stated conditions on $f$, $\tilde{F}$ is well-defined [1, Theorem 1.2.2]. Moreover, it is Gâteaux-differentiable and its Gâteaux-differential is given by [1, Theorem 1.2.7]

[TABLE]

Now, $F=\tilde{F}\circ J$. Hence,

[TABLE]

where we used $u=J(x)$ and $v=J(z)$.

The norm of the derivative can be estimated by

[TABLE]

where $C$ denotes a bound on the partial derivatives $f_{y}$ and $f_{x}$. Hence, the Gâteaux-derivative is uniformly bounded.

(ii) We will need an inverse inequality for functions from $X_{\pi}$. A consequence of [3, Theorem 3.2.6] is the estimate

[TABLE]

for a constant $c$ independent of $\pi\in\mathcal{M}_{[r]}$.

Let $\pi\in\mathcal{M}_{[r]}$ and $x,\bar{x}\in X_{\pi}$. Then ist holds

[TABLE]

Applying (28), we arrive at

[TABLE]

which proves the assertion. ∎

Remark 3.2**.**

According to Propsition 3.1, the Gâteaux-derivative $F^{\prime}$ is contiuous on each $X_{\pi}$. Hence, it is Fréchet-differentiable there. However, $F$ is in general not Fréchet-differentiable on $H^{1}_{D}$ unless it has a very special structure. A discussion of related question can be found in [1, Section 1.2].

Corollary 3.3**.**

Let the partial derivatives $f_{y}$ and $f_{x}$ satisfy the Lipschitz condition in Proposition 3.1(ii) locally. Let $\tilde{x}\in{\rm dom}\,F$ be a sufficiently smooth function, possibly not belonging to $X_{\pi}$, $\lVert\tilde{x}-x_{*}\rVert\leq\rho/2$. Then it holds, for all $\tau\in[0,1]$,

[TABLE]

with a constants $\hat{L}$. In particular, for $\tau=0$, we obtain

[TABLE]

and

[TABLE]

Proof.

Let $I_{\pi}:H_{D}^{1}\cap C_{\pi}([a,b],\mathbb{R}^{m})\rightarrow X_{\pi}$ be a piecewise polynomial interpolation operator. In order to be specific, consider node sequences

[TABLE]

and define $I_{\pi}$ componentwise. For a component $x_{\kappa}\in C[a,b]$, $1\leq\kappa\leq k$, $I_{\pi,\kappa}x_{\kappa}$ is the piecewise polynomial interpolation using the nodes $\bar{t}_{ji}=t_{j-1}+\sigma_{i}^{d}h_{j}$, $i=0,\ldots,N$, $j=1,\ldots,n$. Analogously, for $x_{\kappa}\in C_{\pi}[a,b]$, $I_{\pi,\kappa}$, $k<\kappa\leq m$ is the piecewise polynomial iterpolation using the nodes $\bar{t}_{ji}=t_{j-1}+\sigma_{i}^{a}h_{j}$, $i=1,\ldots,N$, $j=1,\ldots,n$. Then we set $I_{\pi}=[I_{\pi,1},\ldots,I_{\pi,m}]^{T}$.

Let $R_{\pi}=\tilde{x}-I_{\pi}\tilde{x}$ be the remainder. Standard interpolation results provide the estimate

[TABLE]

For all sufficiently fine partitions $\pi\in\mathcal{M}_{[r]}$, $I_{\pi}\tilde{x}$ belongs also to ${\rm dom}\,F$.

Since $I_{\pi}$ is the identity on $X_{\pi}$, we have, for each $x\in X_{\pi}\cap{\rm dom}\,F$,

[TABLE]

Following the lines of the proof of Proposition 3.1(ii) we arrive at the estimate

[TABLE]

hence

[TABLE]

Then we obtain

[TABLE]

This proofs the assertion. ∎

4 Properties related to individual sufficiently fine partitions $\pi$

This section is to provide an approximation of the solution $x_{*}$ by means of an iteration residing in $X_{\pi}$ for an arbitrary sufficiently fine individual partition $\pi\in\mathcal{M}_{[r]}$.

The space of ansatz functions $X_{\pi}$ is defined by (2) as before. Below we frequently apply the topological decompositions

[TABLE]

and the associated orthoprojectors

[TABLE]

in which $x\in{\rm dom}\,F$. $F^{\prime}(x_{*})$ is a fine DAO with index $\mu$ and $\mathcal{F}^{\prime}(x_{*})$ is injective, but its inverse is unbounded if $\mu>1$.

Lemma 4.1**.**

Let $x_{*}$ be sufficiently smooth. Let $N>\mu-1$. Choose $s\in\mathbb{R}$ with $s>\mu-1/2>0$ and $\rho_{\pi}:=c_{\rho}h^{s}_{\pi}$, with a constant $c_{\rho}>0$. Then there is a constant $c_{\gamma}>0$, such that the following relations become valid:

[TABLE]

for each arbitrary mesh $\pi\in\mathcal{M}_{[r]}$ with sufficiently small $h_{\pi}$.

Proof.

The existence of $c_{\gamma}>0$ as well as the inequality (33) are ensured by [7, Theorem 4.1] concerning the instability threshold. $c_{\gamma}$ may depend on the ratio $r$. The injectivity of $\mathcal{F}^{\prime}(x_{*})$ immediately implies ${\rm ker}\,\mathcal{F}^{\prime}(x_{*})U_{\pi}={\rm ker}\,U_{\pi}$.

For $x\in\bar{\mathfrak{B}}(x_{*},\rho_{\pi})\cap X_{\pi}$, $\rho_{\pi}<\rho$, we have

[TABLE]

and

[TABLE]

Making the stepsize $h_{\pi}$ small enough and regarding Corollary 3.3, (3), and (33) yields

[TABLE]

Applying Lemma A.2 of the appendix it results that

[TABLE]

and further

[TABLE]

Taking into account (37) we have

[TABLE]

and, in particular, (34). It also follows that

[TABLE]

Multiplying the last identity from the right by ($\mathcal{F}(x)U_{\pi})^{+}$ yields

[TABLE]

that means (35), and (36) follows immediately. ∎

It should be noted that $s$ in the previous lemma is not restricted to be an integer.

As previously agreed upon, there exists $x_{*}$ such that $\mathcal{F}x_{*}=0$, thus $\psi(x_{*})=0$, $F^{\prime}(x_{*})$ is a fine DAO, the varionational problem $\mathcal{F}^{\prime}(x_{*})z=0$ features accurately stated boundary condition, and the composed operator $\mathcal{F}^{\prime}(x_{*})$ is injective. Assuming the solution $x_{*}$ to be smooth enough we apply the estimates (cf. [7])

[TABLE]

in which $N$ is again the polynomial degree used for the ansatz space $X_{\pi}$.

Since the inverse $\mathcal{F}^{\prime}(x_{*})^{-1}$ is unbounded, standard Newton-like iterations cannot be expected to work well here. Instead we apply a kind of projected Newton iteration using the bounded Moore-Penrose inverse333Note that $(\mathcal{F}^{\prime}(x)U_{\pi})^{+}$ is a bounded outer inverse of $\mathcal{F}^{\prime}(x_{*})$. $(\mathcal{F}^{\prime}(x)U_{\pi})^{+}$ against the background of Lemma 4.1.

More precisely, supposing that $h_{\pi}$ is small enough, we take an initial guess $x_{0}\in\bar{\mathfrak{B}}(x_{*},\rho_{\pi})\cap X_{\pi}$ and provide the correction $z_{1}$ by means of the least-squares problem

[TABLE]

and then put $x_{1}=x_{0}+z_{1}$, and so on. By construction, $z_{1}$ is well defined and belongs to $X_{\pi}$, and so does the new iteration $x_{1}$. Notice that $z_{1}=U_{\pi}z_{1}$ serves as descent direction of the functional $\psi$ at $x_{0}$, as long as $\mathcal{V}_{\pi}(x_{0})\mathcal{F}x_{0}\neq 0$, because of

[TABLE]

Next we ask if $x_{1}$ belongs to the ball $\bar{\mathfrak{B}}(x_{*},\rho_{\pi})$ . For this aim we derive

[TABLE]

Then

[TABLE]

hence, applying Corollary 3.3 for $\tilde{x}=\tau x_{0}+(1-\tau)x_{*}$, $\rVert\tilde{x}-x_{*}\lVert\leq\rho_{\pi}\leq\frac{1}{2}\rho$, and supposing $N>s$,

[TABLE]

for sufficiently small $h_{\pi}$, cf. (38). Next, for

[TABLE]

we obtain a constant $c_{*}$ such that

[TABLE]

Now, to ensure that $x_{1}$ belongs to the ball $\bar{\mathfrak{B}}(x_{*},\rho_{\pi})$, we are confronted with the requirement

[TABLE]

which becomes valid by choosing $N$ so that

[TABLE]

for all sufficiently fine meshes $\pi\in\mathcal{M}_{[r]}$. Then we continue the iterations by providing

[TABLE]

The sequence $\{x_{k}\}$ remains in $\bar{\mathfrak{B}}(x_{*},\rho_{\pi})$. Furthermore we have

[TABLE]

There is a number $k_{\pi}\in\mathbb{N}$ so that one has $\left(\frac{1}{2}\right)^{k+1}\leq\frac{c_{*}}{c_{\rho}}h_{\pi}^{N-\mu+1-s}$ for all $k\geq k_{\pi}$, and hence

[TABLE]

We summarize what we get:

Theorem 4.2**.**

*Let $\mathcal{F}x=0$ denote the operator formulation from Section 3 associated with the BVP (1),(2), $\mathcal{F}x_{*}=0$, ${\rm ker}\,\mathcal{F}^{\prime}(x_{*})=\{0\}$, and $x_{*}$ be sufficiently smooth for (39) to hold.

Let the radius $\rho_{\pi}$ and the bound $\Gamma_{\pi}$ be as introduced in Lemma 4.1, and*

[TABLE]

and the mesh $\pi\in\mathcal{M}_{[r]}$ be sufficiently fine. Then the iteration (42) starting from $x_{0}\in\bar{\mathfrak{B}}(x_{*},\rho_{\pi})\cap X_{\pi}$ remains therein and there is a number $k_{\pi}\in\mathbb{N}$ such that the estimate (44) is valid and

[TABLE]

Proof.

It only remains to verify (46) which is a simple consequence of (3) and (44):

[TABLE]

∎

Let us emphasize that the constants $c_{\gamma},c_{\rho},c_{\alpha}$, $c_{*}$, and $M_{*}$ are global bounds for all partitions $\pi\in\mathcal{M}_{[r]}$.

5 Numerical experiments

In this section, we present the results of some experiments in order to illustrate the properties of the proposed method.

The nonlinear least-squares method (25) has been implemented in Matlab. Instead of (25), its approximation $\psi_{\pi,M}$ of (13) has been used. The finite-dimensional problems have been solved using a Matlab implementation of a Gauss-Newton method following the lines of [4, Section 4.3]. The iteration has been stopped if no further improvement in $\psi_{\pi,M}(x_{k})$ could be observed. For the purposes of investigating the convergence of the method, an interpolation of the exact solution has been used as an initial guess.

5.1 The mathematical pendulum

This problem has been used in many publications for demonstrating properties of algorithms for the solution of differentail algebraic systems. We use the formulation

[TABLE]

The underlying interval is $(0,1)$. The parameters are chosen to be $g=16$, $L=\sqrt{8}$. We consider the initial values $y(0)=2$ and $y^{\prime}(0)=0$. This problem has index 3. Therefore, the results of Theorem 4.2 are only valid if $N\geq 5$. For $N=5$, $s=\mu-1/4$ can be chosen. However, the expected orders are observed in all cases $N\geq 2$. The case $N=1$ is rather surprising since we observed bounded solutions instead of diverging ones.

In Tables 1 and 2 as well as Tables 3 and 4 results for $N=3$ and $N=5$, respectively, are presented. In both cases, uniform grids and $M=N+1$ uniformly distributed collocation points per subinterval have been used.

5.2 An example proposed by S.L. Campbell and E. Moore

In [2], the following system is used as an example:

[TABLE]

The solution considered in the reference is

[TABLE]

yielding

[TABLE]

In [2], the inequality $r>\rho$ is supposed and the numerical experiments are carried out for $\rho=5$ and $r=10$. We use the same parameters in the following experiment. Under these conditions, the problem has index 3.

A thorough discussion as well as numerical experiments of the version linearized in the solution $x_{*}$ is given in [7]. In order to stimulate discussions of the least-squares method for nonlinear problems, also results for the original nonlinear version have been provided in this reference. We cite the results in Tables 5 and 6. Theorem 4.2 is only valid for $N\geq 5$ in this example and, thus, the corresponding order is strictly proven. However, the expected orders are observed in allowed cases $N\geq 2$. The case $N=1$ is rather surprising besause we observe bounded solutions even if we expecteddiverging ones.

6 Multilevel approach

We use $N$ and $s$ as previously agreed, that is $N-\mu+1>s>\mu-1/2>0$. Given an additional constant $q$ with $0<q<1$ we now deal with a sequence of partitions $\pi_{i}\in\mathcal{M}_{[r]}$,

[TABLE]

such that the associated ansatz spaces are nested,

[TABLE]

and $h_{\pi_{i}}\rightarrow 0$ if $i\rightarrow\infty$. Let $\pi_{0}$ be fine enough for Lemma 4.1 and Theorem 4.2 to hold. This means that

[TABLE]

to ensure the applicability of Lemma 4.1 and to make the iterations on the level $\pi_{0}$ to stay in $\bar{\mathfrak{B}}(x_{*},\rho_{\pi_{0}})\cap X_{\pi_{0}}$. Both conditions are satisfied correspondingly a fortiori on the further levels due to the smaller stepsizes $h_{\pi_{i}}$. In the consequence, Theorem 4.2 applies on each level, i.e., for $x_{0}^{[\pi_{i}]}\in\bar{\mathfrak{B}}(x_{*},\rho_{\pi_{i}})\cap X_{\pi}$ the sequence

[TABLE]

remains in $\bar{\mathfrak{B}}(x_{*},\rho_{\pi_{i}})\cap X_{\pi_{i}}$ and there exists a number $k_{\pi_{i}}\in\mathbb{N}$ such that $\left(\frac{1}{2}\right)^{k+1}\leq\frac{c_{*}}{c_{\rho}}h_{\pi_{i}}^{N-\mu+1-s}$ for all $k\geq k_{\pi_{i}}$, and hence

[TABLE]

Since the ansatz spaces are nested, $x_{k_{\pi_{i}+1}}^{[\pi_{i}]}$ belongs to $X_{\pi_{i+1}}$. Replacing the condition $h_{\pi_{i}}^{N-\mu+1-s}\leq\frac{1}{2}\frac{c_{\rho}}{c_{*}}$ by the stronger one

[TABLE]

yields

[TABLE]

Then $x_{k_{\pi_{i}}+1}^{[\pi_{i}]}$ belongs to $\bar{\mathfrak{B}}(x_{*},\rho_{\pi_{i+1}})\cap X_{\pi_{i+1}}$ and we are allowed to choose at the next level

[TABLE]

We summarize our result:

Theorem 6.1**.**

*Let $\mathcal{F}x=0$ denote the operator formulation from Section 3 associated with the BVP (1),(2), $\mathcal{F}x_{*}=0$, ${\rm ker}\,\mathcal{F}^{\prime}(x_{*})=\{0\}$, and $x_{*}$ be sufficiently smooth for (39).

Let (45) be given and $0<q<1$.*

Let the sequence of partitions $\pi_{i}\in\mathcal{M}_{[r]}$, $i\geq 0$, be such that the ansatz spaces are nested and the maximal stepsizes are related by $qh_{\pi_{i}}=h_{\pi_{i+1}}$. Let the the mesh $\pi_{0}$ be sufficiently fine,

[TABLE]

Then the iteration (47),(48),(51), with the initial guess $x_{0}^{[\pi_{0}]}\in\bar{\mathfrak{B}}(x_{*},\rho_{\pi_{0}})\cap X_{\pi_{0}}$ is well defined and yields

[TABLE]

7 Remarks and conclusions

We have presented and investigated a nonlinear least-squares method for approximating higher index differential-algebraic equations. The idea consists of discretizing the preimage space $H^{1}_{D}$ by piecewise polynomials and to form an overdetermined collocation system to determine an approximating solution. The resulting overdetermined system is solved in a least-squares sense. In the numerical experiments, the method behaved very well despite its simplicity. In particular, the method is not much more expensive than the standard collocation method applied to explicit ordinary differential equations and index-1 differential-algebraic equations.

The main tool both for the convergence proof and for the numerical solution of the discretized problems is a variant of the Newton method. For a large class of nonlinear index-$\mu$ tractable equations, this method applied to the discretized system is shown to deliver appropriate approximations provided that the polynomial order is large enough. The numerical experiments indicate, however, that the strong condition on the polynomial order does not seem to be necessary. In particular, the order of convergence corresponds to that of linear index-$\mu$ tractable differential-algebraic equations. So the present result should be considered as a first step towards a theoretical foundation of the method.

Remark 7.1**.**

1. (i)

Under the conditions of Theorem 4.2 we could not show that the sequence $\{x_{k}\}$ converges. 2. (ii)

If there is a minimizer $x_{\pi,*}$ of (14) in ${\mathfrak{B}}(x_{*},\rho_{\pi})\cap X_{\pi}$, then it holds $(\mathcal{F}^{\prime}({x}_{\pi,*})\mathcal{U}_{\pi})^{+}\mathcal{F}({x}_{\pi,*})=0$.

Since $\bar{\mathfrak{B}}(x_{*},\rho_{\pi})\cap X_{\pi}$ is compact, the sequence $\{x_{k}\}$ has a convergent subsequence. However, we were not able to show that, for an accumulation point $\hat{x}_{\pi}$, it holds $(\mathcal{F}^{\prime}(\hat{x}_{\pi})\mathcal{U}_{\pi})^{+}\mathcal{F}(\hat{x}_{\pi})=0$. 3. (iii)

In the context of regularization methods for nonlinear illposed problems, the so-called Scherzer, or tangential cone, condition is often used **[16, 9, 10]**. However, in the context of differential-algebraic equations, this conditions requires very hard conditions on the structure of the system. Therefore, it is of minor use here.

Appendix A An auxillary result

The convergence proof for the Gauss-Newton method requires an estimation of the norm and distance of Moore-Penrose inverses of derivatives of a nonlinear operator. In the case of finite dimensional spaces, such results are well-known and can be found, for example, in [13]. However, we need similar statements in the case of infinite dimensional spaces. This appendix provides the necessary lemmas.

Let $X$ and $Y$ be Hilbert spaces (not necessarily finite dimensional) and $A:X\rightarrow Y$ a linear and compact operator. Both operators $A^{\ast}A$ and $AA^{\ast}$ are selfadjoint compact operators. Their spectra consist only of nonnegative eigenvalues with finite multiplicity (with the possible exception of $\lambda=0$). If the eigenvalues have an accumulation point, then it is 0. The nonzero eigenvalues are identical (even with respect to their multiplicity) for both $A^{\ast}A$ and $AA^{\ast}$. Let the nonzero eigenvalues be sorted according to

[TABLE]

Let then $\{u_{i}\}\subset X$ and $\{v_{i}\}\subset Y$ be a complete orthonormal system of eigenvalues444The systems are not necessarily complete in $X$ and $Y$, respectively! for the operators $A^{\ast}A$ and $AA^{\ast}$,

[TABLE]

We set $\sigma_{i}=\sqrt{\lambda_{i}}>0$. This provides us with

[TABLE]

The system $\{\sigma_{i},u_{i},v_{i}\}$ is called a singular system of $A$ with the singular values $\sigma_{i}$. In particular, we have the representations

[TABLE]

Here, $(\cdot,\cdot)$ denotes the scalar products in $X$ and $Y$, respectively. Note that these sums can be both finite and infinite.

We are interested in perturbation results for the singular values of an operator $A$. The following lemma is proven in [5, Corollary VI.1.6].

Lemma 1**.**

Let $X$ and $Y$ be Hilbert spaces. Let $A,B:X\rightarrow Y$ be compact linear operators. Let $\sigma_{i}(A)$, $i=1,\ldots,\nu(A)$ and $\sigma_{i}(B)$, $i=1,\ldots,\nu(B)$ be the singular values of $A$ and $B$, respectively.555Both $\nu(A)$ and $\nu(B)$ may be finite or infinite. Assume without loss of generality $\nu(A)\leq\nu(B)$. Then it holds

[TABLE]

The next step consists of the establishement of bounds for the Moore-Penrose pseudoinverse. For compact operators as above with the singular system $\{\sigma_{i},u_{i},v_{i}\}$ it has the representation

[TABLE]

for all $y\in{\rm dom}\,(A^{+})$. An immediate consequence is:

1. (i)

$\lVert A\rVert=\sigma_{1}$. 2. (ii)

$A^{+}$ is bounded if and only if the number of singular values is finite. In that case it holds $\lVert A^{+}\rVert=\sigma_{\nu(A)}^{-1}.$

The following lemma presents modifications of [13, Theorem (8.15)].

Lemma 2**.**

Let $A,B:X\rightarrow Y$ be compact linear operators acting in the Hilbert spaces $X,Y$.

1. (i)

Assume $\dim X<\infty$ and $\nu(B)\leq\nu(A)=r$, $r\geq 1$. Moreover, let

[TABLE]

to hold and set $\epsilon=\lVert B-A\rVert$. Then it holds $\nu(B)=\nu(A)$ and

[TABLE]

where $\sigma_{\nu(A)}$ is the smallest singular value of $A$. 2. (ii)

Assume that $X$ decomposes in $X=X_{f}\oplus X_{f}^{\bot}$, $\dim X_{f}<\infty$, and

[TABLE]

Then it follows that

[TABLE]

Proof.

It holds $\lVert A^{+}\rVert=\sigma_{\nu(A)}^{-1}$ such that, by assumption, $\sigma_{i}-\epsilon>0$, $i=1,\ldots,\nu(A)$. Hence, $B=A+(B-A)$ has at least $\nu(A)$ nonvanishing singular values because of Lemma A.1. Hence, $\nu(B)\geq\nu(A)$. Together with the assumption, this provides $\nu(B)=\nu(A)$. Consequently, $\lVert B^{+}\rVert\leq 1/(\sigma_{\nu(A)}-\epsilon)$. (ii) is a consequence of (i). ∎