DISCRETIZATION BASED INDICES FOR SEMILINEAR PARTIAL DIFFERENTIAL ALGEBRAIC EQUATIONS

W. LUCHT, K. STREHMEL

Martin-Luther-Universität Halle–Wittenberg Fachbereich Mathematik und Informatik Institut für Numerische Mathematik Postfach, D–06099 Halle, Germany email: lucht@mail.mathematik.uni-halle.de

Abstract. Semilinear partial differential algebraic equations (PDAEs) of the form Au_{t}(t, x)+ Bu_{xx}(t, x) + F (u) = f(t, x) (where F (u) ∈ ℝ^{n} is a nonlinear function of u ∈ ℝ^{n}) are studied where at least one of the matrices A, B ∈ ℝ^{n×n} is singular. For these systems we propose definitions of a differential time index and of a differential space index which are based on vertical and horizontal semidiscretization schemes, respectively. In these definitions an essential assumption is that the corresponding semidiscretizations converge. For two often used methods we show the convergence for index one systems.

AMS subject classification: 35E99, 35M10

Keywords: Differential algebraic equations, semilinear partial differential algebraic equations, coupled systems, indexes

1. Introduction

In this paper partial differential algebraic equations (PDAEs) of the diffusion reaction type

Au_{t} + Bu_{xx} + F (u) = f(t, x), t ∈ (0, t_{e}), t_{e} > 0, x ∈ Ω := (−l, l), l > 0,  (1.1)

are considered. A and B are assumed to be real constant (n, n)−matrices (n > 1) where at least one of these matrices is singular. F : D ⊂ ℝ^{n} → ℝ^{n} is, in general, some nonlinear (smooth) vector function which is bounded in the following sense. Let ξ ∈ ℝ^{n} and denote by |ξ| the Euclidean norm |ξ| := [∑_{i=1}^{n} ξ_{i}^{2}]^{1/2}. If D is a simply connected bounded set with 0 ∈ D, we assume throughout this paper that |F (ξ)| ≤ K_{0} for all ξ ∈ D, K_{0} being a finite positive constant. f ∈ ℝ^{n} is supposed to be a sufficiently smooth function of its arguments. PDAEs were investigated by several authors from different point of views (e.g. index definitions, Fourier transformation, method of lines), see above all [3], [5], [6], [7], [15], [20]. The linear constant coefficient version of the PDAE (1.1) has been considered also in [13], [14].

Systems of the type (1.1) occur frequently in applications, e.g. in reaction chemistry [9, Chapter 6] and in population dynamics [12, Chapter II], [18]. In typical examples, the components of F are polynomials in the variables u_{1}, ..., u_{n}, e.g. for n = 2

F_{1} = c_{1}u_{1}u_{2} + c_{2}u_{1}^{2}, F_{2} = c_{3}u_{2} − c_{4}u_{1}(1 − u_{2})u_{2}

with constant parameters c_{i}. The system (1.1) must be supplemented by suitable initial and boundary conditions (ICs and BCs) given below. As in the case of linear PDAEs it is a matter of fact that, in general, arbitrary initial and boundary values cannot be prescribed for all components of u being a solution of eq. (1.1). Here, we define a solution by the following

Definition 1.1. For any sufficiently smooth f, a function u = u(t, x), t > 0, x ∈ Ω, is a solution of the PDAE (in the classical sense) if it

(1) is sufficiently smooth,

(2) is uniquely determined by its initial and boundary values,

(3) solves (pointwise) the PDAE.

Remark 1.1. ”Sufficiently smooth” means that u and its derivatives entering the PDAE exist and are continuous for t > 0 and x ∈ Ω. Because at least one of the matrices A, B is singular it may be that the components of u need not have all the same smoothness. On the other hand the claim of smoothness of u requires also that the right hand vector f must be smooth. □ Guided by the linear case we introduce two index sets 𝔐_{IC} and 𝔐_{BC}, see [13].

Definition 1.2. If i ∈ 𝔐_{IC} ⊆ {1, ..., n}, then any IC can be given for u_{i}. If i ∉ 𝔐_{IC}, then an arbitrary IC cannot be prescribed for u_{i}.

Analogous is the

Definition 1.3. If j ∈ 𝔐_{BC} ⊆ {1, ..., n}, then any BC can be prescribed for u_{j}. j ∉ 𝔐_{BC} means the BC for u_{j} can only be calculated (from the PDAE).

After these definitions the BCs and ICs which can be chosen arbitrarily for an initial boundary value problem (IBVP) associated with eq. (1.1) can be specified. For simplicity we assume homogeneous Dirichlet BCs

R_{B} u_{j}(t, x) := u_{j}(t, ±l) = 0, j ∈ 𝔐_{BC}, t ∈ [0, t_{e}].  (1.2)

The given ICs are

u_{i}(0, x) = g_{i}(x) for x ∈ [−l, l], i ∈ 𝔐_{IC}.  (1.3)

The IC for u_{i}, i ∉ 𝔐_{IC}, and the BC for u_{k}, k ∉ 𝔐_{BC}, must be determined by means of the PDAE and the ICs for u_{j}, j ∈ 𝔐_{IC}, and/or the BCs for u_{l}, l ∈ 𝔐_{BC}. This yields the so called consistent ICs and BCs (for linear PDAEs this is outlined in [13]). Furthermore, we also require the compatibility conditions

R_{B} g(x) = R_{B} u(0, x)  (1.4)

between the BCs and the ICs.

Remark 1.2. We mention that there are PDAEs for which the two index sets are not unique. That means there are different possibilities to prescribe ICs and/or BCs. To give an example, we consider the case n = 2 and suppose the PDAE yields, e.g., a BC of the type γ_{1}u_{1}(t, ±l) + γ_{2}u_{2}(±l) = γ_{3} where γ_{1}, γ_{2} are coefficients different from zero. γ_{3} may be chosen arbitrarily. Then we can prescribe an arbitrary BC either for u_{1} or for u_{2}. An illustration is the simple linear PDAE

(\begin{matrix}1 & 0 \\ 0 & 1 \end{matrix}) u_{t} + (\begin{matrix}−1 & −1 \\ 0 & 0 \end{matrix}) u_{xx} + (\begin{matrix}0 & 0 \\ β & 0 \end{matrix}) u = f

(β is a constant different from zero). Here, it is 𝔐_{IC} = {1, 2}, but there are two possibilities for 𝔐_{BC}: either 𝔐_{BC} = {1} or 𝔐_{BC} = {2}. (For more details see [13, section 3].) □

2. Differential time and space index

In the following we suppose that the IBVP (1.1)–(1.4) has an unique solution which is sufficiently smooth (this assumption is used extensively in sections 3 and 4). The notion of an index (or indexes) of a system of differential algebraic equations (DAEs) or of a system of linear PDAEs plays a fundamental role in the theory and numerical solution. For example, the convergence of a numerical solution u_{τ,h} of a linear IBVP (where τ and h are time and space step sizes of a full discretization scheme, respectively) to the exact solution u for τ → 0 and h → 0 depends strongly on its indexes [14].

There are several variants of indexes of general DAEs and linear PDAEs, see e.g. [1], [4], [6], [8], [10] and references therein. Of course, it is expected that semilinear PDAEs can also be characterized by indexes. In order to introduce these we go back to the index definition of linear constant coefficient PDAEs (considered also in [6], [7]). In [13] we defined an uniform differential time index ν_{dt} by a Fourier and a differential spatial index ν_{dx} by a Laplace transformation of the corresponding PDAE. Although these transformations do not work properly for semilinear PDAEs, the idea is to use for index definitions of semilinear PDAEs certain analogies between the Fourier and Laplace transformation of linear PDAEs on the one hand and of the vertical (method of lines, MOL) and horizontal (Rothe method) semidiscretization of semilinear PDAEs on the other hand, respectively. First we consider in the following subsection the MOL as the analogy of the Fourier transformation. This relationship between PDAEs and MOL is widely used in the literature, see e.g. [7], [1], [20]. The analogy of the Laplace transformation is the Rothe method which is the subject of subsection 2.2.

2.1. Differential time index. In this subsection a time index of the semilinear PDAE (1.1) is defined and for this we consider the space discretization of the IBVP (1.1) – (1.4) by means of finite differences. The PDAE (1.1) is space discretized on an equidistant grid

Ω_{h} = \{x_{k} : x_{k} = −l + kh, k = 1, …, N , h = \frac{2l}{N + 1}\}

for a given integer N ∈ ℕ^{+}. Let u(t, x_{k}) be the exact solution vector at time t and at the grid point x_{k}. Replacing u_{xx}(t, x_{k}) by

u_{xx}(t, x_{k}) ≈ \frac{1}{h^{2}} (u_{k+1}(t) − 2u_{k}(t) + u_{k−1}(t)) , k = 1, …, N

we obtain the semi-discretized equation (u_{k}(t) ≈ u(t, x_{k}))

A u′_{k}(t) + \frac{1}{h^{2}} B (u_{k+1}(t) − 2u_{k}(t) + u_{k−1}(t)) + F (u_{k}(t)) = f_{k}(t)  (2.1)

or in matrix representation (using the Kronecker product)

(I_{N} ⊗ A) U′(t) + (\frac{1}{h^{2}} P ⊗ B) U(t) + G(U(t)) = F̃(t) − r(t).  (2.2)

where the component u_{k}(t) of the vector U(t) ∈ ℝ^{nN} is the continuous time approximation to u(t, x_{k}) and the nN dimensional vectors G(U(t)), r(t), F̃(t) are given by

G(U(t)) = (F^{⊤}(u_{1}(t)), ..., F^{⊤}(u_{N}(t)))^{⊤} ,
r(t) = (\frac{1}{h^{2}} I_{N} ⊗ B) (u^{⊤}(t, −l), 0 …, 0, u^{⊤}(t, l))^{⊤},
F̃(t) = (f_{1}^{⊤}(t), …, f_{N}^{⊤}(t))^{⊤}, f_{k}(t) = f(t, x_{k}).

I_{N} is the (N × N)-unit matrix, and the matrix P is defined by

P = (\begin{matrix}−2 & 1 \\ 1 & −2 & & 1 \\ & & . & . & . \\ & & & 1 & & −2 \end{matrix}) ∈ ℝ^{N×N}.  (2.3)

We note for later use that the eigenvalues of \frac{1}{h^{2}} P are given by λ_{k} = −\frac{4}{h^{2}} sin^{2} (\frac{kπ}{2(N+1)}) for k = 1, ..., N.

Eq. (2.2) is (according to eq. (1.3)) completed by the consistent initial vector

U(0) = (g̃^{⊤}(x_{1}), ..., g̃^{⊤}(x_{N}))^{⊤} ∈ ℝ^{nN}  (2.4)

where the difference g̃ − g tends (component wise) to zero for h → 0. If A is singular, then eq. (2.2) is for fixed h a DAE (also called MOL–DAE) which may be handled by standard methods ([2], [10]). Especially, its usual differential index (which is a differential time index in our context) can be determined. To do this, the transformation to a semi–explicit form, as given in eqs. (3.14), (3.15), may be advantageous.

Let U_{h}(t) be the restriction of the exact solution u(t, x) to the spatial grid Ω_{h} and define the global space error of the semidiscretization by

η = η_{h} := U_{h}(t) − U(t).  (2.5)

Let furthermore for two vectors v, w ∈ ℝ^{nN} the discrete scalar product and the corresponding induced norm be defined by

(v, w)_{h} := ∑_{ν=1}^{nN} h v_{ν}w_{ν}, ∥ v ∥_{h}:= \sqrt{(v, v)_{h}}.

With this notation, the MOL-DAE (2.2) suggests the following

Definition 2.1. Suppose (1) lim_{h→0} ∥ η ∥_{h} = 0,

(2) there exists a step size h_{0} > 0 such that the differential index ν_{dt} of the DAE (2.2) is

the same for any h ∈ (0, h_{0}]. Then the differential time index of the semilinear PDAE (1.1) is (with respect to the space discretized scheme) defined to be ν_{dt}.

Remark 2.1. For the linear case F (u) = Cu (where C is a constant matrix) it is shown in [14] that under some reasonable assumptions ∥ η ∥_{h}= O(h^{2}) for h → 0 (i.e. the first condition in the index definition 2.1 is satisfied). The second condition is based on the idea that the index of the MOL–DAE should not change if h_{0} is sufficiently small (see also [3]). There are examples in which the MOL–DAE has for a specified step size h^{∗} an index which is different of the index for h < h^{∗} (see [14, section 4.4]). □

In section 3 we investigate the first condition of the definition for general semilinear PDAEs (1.1) of differential time index ν_{dt} = 1.

2.2. Differential space index. In order to define a differential space index of the semilinear PDAE (1.1) we consider the time discretization of the IBVP by means of finite differences (see [19]). We define for simplicity with a constant time step size τ = t_{e}/M, t_{e} < ∞, M ∈ ℕ^{+}, time grid points t_{m} = mτ, m = 0, ..., M. Let u^{m}(x) be an approximation for the exact solution u(t_{m}, x). Replacing in the PDAE (1.1) the time derivative by the backward Euler difference approximation we get on Ω the second order equation

Bu_{xx}^{m}(x) + \frac{1}{τ}Au^{m}(x) + F (u^{m}(x)) = f(t_{m}, x) + \frac{1}{τ}Au^{m−1}(x), m = 1, ..., M,  (2.6)

where according to (1.3) the initial values are written for all components of u as u^{0}(x) = g(x), x ∈ \overline{Ω}. The BCs are

R_{B} u_{j}^{m}(x) = u_{j}^{m}(±l) = u_{j}(t_{m}, ±l) = 0, j ∈ 𝔐_{BC}.  (2.7)

If B is singular, we suppose that the consistent BCs for u_{j}^{m}, j ∉ 𝔐_{BC}, are known. Thus, for each m a boundary value problem (BVP) for a DAE must be solved.

Let U_{τ}(x) ∈ ℝ^{nM} and U(x) ∈ ℝ^{nM} be the vectors

U_{τ}(x) = (u^{⊤}(t_{1}, x), ..., u^{⊤}(t_{M}, x))^{⊤} ,  U(x) = (u^{1 ⊤}(x), ..., u^{M ⊤}(x))^{⊤}  (2.8)

such that

ξ(x) = ξ_{τ}(x) := U_{τ}(x) − U(x)  (2.9)

is the global time error of the semidiscretization. Let again a scalar product and a norm on ℝ^{nM} be given by

(v, w)_{τ} := ∑_{ν=1}^{nM} τv_{ν}w_{ν} and ∥ v ∥_{τ}:= \sqrt{(v, v)_{τ}},  (2.10)

respectively.

Definition 2.2. Assume (1) lim_{τ→0} ∥ ξ ∥_{τ}= 0,

(2) there is a time step size τ_{0} > 0 such that the differential index ν_{dx} of the DAE (2.6) is

independent of m and the same for any τ ∈ (0, τ_{0}]. Then the differential space index of the semilinear PDAE (1.1) is (with respect to the time discretized scheme) defined to be ν_{dx}.

As is obvious, this definition is analogous to the definition 2.1 of the differential time index. The question of convergence of ∥ξ∥_{τ} for τ → 0 for systems with ν_{dx} = 1 is considered in section 4.

3. Local convergence of the space discretized solution

According to the definition of the differential time index it is for a given problem necessary to verify the condition (1) of the definition 2.1. In order to do this we start from eq. (2.2) into which we insert the exact solution vector U_{h}. The result is the equation

(I_{N} ⊗ A)U′_{h}(t) + (\frac{1}{h^{2}} P ⊗ B) U_{h}(t) + G(U_{h}(t)) = F̃(t) − r(t) + α_{h}(t)  (3.1)

where α_{h}(t) is the space truncation error. It represents the defect of the grid function U_{h}(t) with respect to the DAE (2.2). Roughly speaking, α_{h} shows how good U_{h}(t) satisfies the DAE (2.2).

Definition 3.1. The semidiscretization is said to be consistent if

max_{t∈(0,t_{e})} ∥α_{h}(t)∥_{h} = o(1) as h → 0,

where consistency of order p^{∗} is defined again in an obvious way.

A componentwise Taylor expansion (remember that U_{h}(t) is smooth by assumption) yields

(\frac{1}{h^{2}} P ⊗ B) U_{h}(t) + r(t) = (I_{N} ⊗ B)U_{xx,h}(t) + h^{2}(I_{N} ⊗ B)γ_{h}(t),

γ_{h}(t) = \frac{1}{12} \frac{∂^{4}}{∂x^{4}} (u_{11}(t, ζ_{11}), ..., u_{1n}(t, ζ_{1n}), u_{21}(t, ζ_{21}), ..., u_{2n}(t, ζ_{2n}), ...,
u_{N1}(t, ζ_{N1}), ..., u_{Nn}(t, ζ_{Nn}))^{⊤},

ζ_{ij} ∈ (x_{i−1}, x_{i}), i = 1, .., N, j = 1, ..., n. Here, u_{ij} is the component u_{j} in the space grid point i. Therefore, the space truncation error can be written

α_{h}(t) = h^{2}(I_{N} ⊗ B)γ_{h}(t).  (3.2)

Now we assume that max_{−l≤x≤l} |\frac{∂^{4}}{∂x^{4}} u_{i}(t, x)| , i = 1, ..., n, is bounded for t ∈ (0, t_{e}) implying that

the estimate

∥α_{h}(t)∥_{h} = O(h^{2}) for h → 0, t ∈ (0, t_{e})  (3.3)

holds. Thus, the discretization in space is consistent of second order.

Subtracting eq. (2.2) from eq. (3.1) and taking into account the definition of the error function η, eq. (2.5), we obtain for η(t), t ∈ (0, t_{e}), the DAE

(I_{N} ⊗ A)η′(t) + (\frac{1}{h^{2}} P ⊗ B) η(t) + G(U_{h}(t)) − G(U(t)) = α_{h}(t)  (3.4)
η(0) = U_{h}(0) − U(0).  (3.5)

This equation shows that η(t) depends essentially on the truncation error α_{h}(t).

Our main interest is the case where A (and hence I_{N} ⊗ A) is singular. For applications it is advantageous to transform I_{N} ⊗ A with constant regular (nN × nN)–matrices S_{0} and S_{1} to diagonal form (e.g. by Gaussian elimination):

Ā := (\begin{matrix}I_{1} & 0 \\ 0 & 0 \end{matrix}) = S_{0}(I_{N} ⊗ A)S_{1}.  (3.6)

The dimension N_{1} of the unit matrix I_{1} is equal to the rank of I_{N} ⊗ A. If we define B̄ := S_{0}(P ⊗ B)S_{1}, Ū_{h} := S_{1}^{−1}U_{h}, Ū := S_{1}^{−1}U, η̄ := Ū_{h} − Ū = S_{1}^{−1}η, Ḡ(Ū_{h}) := S_{0}G(S_{1}Ū_{h}), Ḡ(Ū) := S_{0}G(S_{1}Ū) and ᾱ_{h} := S_{0}α_{h}, then eqs. (3.4), (3.5) are transformed into

Āη̄′(t) + \frac{1}{h^{2}} B̄η̄(t) + Φ̄ =ᾱ_{h}(t)  (3.7)
η̄(0) = Ū_{h}(0) − Ū(0).  (3.8)

where for short Φ̄ := Ḡ(Ū_{h}(t)) − Ḡ(Ū(t)).

The aim is now to estimate the error norm ∥η∥_{h} or equivalently ∥η̄∥_{h} for h → 0 (note that the transformation from η to η̄ is essentially only a change of rows and columns in the system (3.4)). In order to do this we must handle the nonlinearity of the IBVP (1.1). We take a point of view which considers the IBVP ”local in time”. The basic idea of the method can be found also at other places in the literature, e.g. for nonlinear parabolic problems in [11].

The nonlinear term Φ̄ can be written (according to (2.5))

Φ̄ = Φ̄(η̄) = Ḡ(Ū_{h}) − Ḡ(Ū_{h} − η̄)  (3.9)

Note that ∥U_{h}∥_{h} is bounded by assumption (see definition 1.1 and the statement at the beginning of section 2). Furthermore, we also assume that the solution of the MOL–DAE (2.2) is smooth such that the solutions U_{h}(t) and U(t) will stay for small t > 0 close to the initial states U_{h}(0) and U(0). Therefore, we consider the vector function η̄(t) being for t ∈ (0, t_{ϵ}) a member of the set

Y_{ϵ} := {y ∈ ℝ^{nN}, ∥y − η̄(0)∥ ≤ ϵ, ϵ > 0.}  (3.10)

Defining a scalar characteristic function ψ by

ψ(y) = \{ \begin{matrix}1 for & y ∈ Y_{ϵ} \\ 0 for & y ∉ Y_{ϵ} \end{matrix}

and a new interaction term Φ̃ by Φ̃(η̄) := ψ(η̄)Φ̄(η̄) we consider beside eqs. (3.7), (3.8) also η̃ as solution of the IVP

Āη̃′(t) + \frac{1}{h^{2}} B̄η̃(t) + Φ̃(η̃(t)) =ᾱ_{h}(t)  (3.11)
η̃(0) = Ū_{h}(0) − Ū(0)  (3.12)

By continuity, there must be a time t_{ϵ} > 0 (t_{ϵ} depends on ϵ) such that for t ∈ (0, t_{ϵ}) η̄(t) can be identified with η̃(t). Remark 3.1. According to the foregoing construction we can obtain (with respect to t) only local results (e.g. estimates of ∥η̄(t)∥_{h}). The idea is, of course, to continue the procedure, i.e. to take the values of η̄(t_{ϵ}) as new initial values etc. □ The essential point is that the nonlinear term in eq. (3.11) can be estimated by

∥Φ̃(y)∥_{h} = \{ \begin{matrix}≤ K_{ϵ}∥y∥_{h} for & y ∈ Y_{ϵ} \\ 0  for & y ∉ Y_{ϵ} \end{matrix}  (3.13)

where K_{ϵ} is finite for ϵ > 0. Because

Φ̃(η̃) = ψ(η̃)Φ̄(η̃) = ψ(η̃)[Ḡ(Ū_{h}) − Ḡ(Ū_{h} − η̃)] = J̃(η̃)η̃

where J̃ is defined by

J̃(η̃) := ψ(η̃) ∫_{0}^{1} Ḡ_{U}(U_{h} + θη̃) dθ,

the DAE (3.11) takes the (semi explicit) form

Āη̃′ + \frac{1}{h^{2}} B̄η̃ + J̃(η̃)η̃ = ᾱ_{h}  (3.14)

whose second equation is explicitly

Q̃_{21}(η̃)η̃_{1} + Q̃_{22}(η̃)η̃_{2} = ᾱ_{2h}.  (3.15)

Here, we introduced for later use the abbreviation

Q̃_{2j}(η̃) := \frac{1}{h^{2}} B̄_{2j} + J̃_{2j}(η̃(t)), j = 1, 2.  (3.16)

Furthermore, we have split the vector η̃ according to eq. (3.6) η̃ = (η̃_{1}^{⊤}, η̃_{2}^{⊤})^{⊤}, η̃_{1} ∈ ℝ^{N_{1}}, η̃_{2} ∈ ℝ^{N_{2}} with N_{1} + N_{2} = nN, and the matrix B̄ has been decomposed corresponding to the dimension of the vectors η̃_{1} and η̃_{2}.

The representation (3.15) can be used to define a differential index of the DAE (3.14) which is equivalent to (3.11) and to the original error equation (3.4). We say the DAE (3.11) (and also the MOL–DAE (2.2)) is of differential index 1 for 0 < t < t_{ϵ} if the matrix Q̃_{22}(·) (which is bounded for fixed h) is regular on Y_{ϵ}. In the following we split a vector v ∈ ℝ^{nN} and its norm into two parts: v = (v_{1}^{⊤}, v_{2}^{⊤})^{⊤}, v_{1} ∈ ℝ^{N_{1}}, v_{2} ∈ ℝ^{N_{2}},

∥v∥_{h}^{2} = ∑_{ν=1}^{N_{1}} hv_{1ν}^{2} + ∑_{ν=1}^{N_{2}} hv_{2ν}^{2} =: ∥v_{1}∥_{1h}^{2} + ∥v_{2}∥_{2h}^{2}, N_{1} + N_{2} = nN.

This is the usual discrete L^{2}−norm which is independent of N. If M is a rectangular matrix of the form, e.g., (N_{2}, N_{1}), then we denote by ∥M∥_{2h} the natural matrix norm

∥M∥_{2h} = sup_{∥v∥_{1h}>0} \frac{∥Mv∥_{2h}}{∥v∥_{1h}} , v ∈ ℝ^{N_{1}}.

After these preliminaries we can prove

Theorem 3.1. Suppose (1) the unique solution of the PDAE (1.1) is sufficiently smooth,

(2) the matrix B̄^{⊤}Ā + Ā^{⊤}B̄ is positive semidefinit,

(3) ᾱ_{h} is for h → 0 of second order, see eq. (3.3),

(4) the error DAE (3.14) is of differential index 1, i.e. the matrix Q̃_{22}(y) is regular, and ∥Q̃_{22}^{−1}(y)∥_{2h} ≤ K_{22} for h > 0, y ∈ Y_{ϵ} (K_{22} independent of h),

(5)

∥Q̃_{22}^{−1}(y) Q̃_{21}(y)∥_{2h} ≤ K_{21}, y ∈ Y_{ϵ},

for h > 0, even in the limit h → 0 (K_{21} independent of h),

(6) ∥Ū_{h}(0) − Ū(0)∥ = O(h^{2}) for h → 0. Then there exists t_{ϵ} > 0 such that

∥η̃(t)∥_{h} + Z_{h}(t) = O(h^{2}) for h → 0 and t ∈ (0, t_{ϵ})  (3.17)

where Z_{h}(t) is defined by

Z_{h}(t) := \frac{1}{h}[ ∫_{0}^{t} e^{C(t−τ)} (η̃(τ), (B̄^{⊤}Ā + Ā^{⊤}B̄)η̃(τ))_{h} dτ ]^{1/2},

C := 1 + K_{ϵ}(2 + 2K_{21} + K_{22}).

Proof. With eq. (3.11) and the Cauchy–Schwarz inequality it is for η̃ ∈ Y_{ϵ}

\frac{d}{dt}(Āη̃, Āη̃)_{h} = (Āη̃′, Āη̃)_{h} + (Āη̃, Āη̃′)_{h} = −\frac{1}{h^{2}} [(B̄η̃, Āη̃)_{h} + (Āη̃, B̄η̃)_{h}] − 2(Φ̃, Āη̃)_{h}

+2(ᾱ_{h}, Āη̃)_{h} ≤ −\frac{1}{h^{2}} (η̃, (B̄^{⊤}Ā + Ā^{⊤}B̄)η̃)_{h} + 2∥Φ̃∥_{h} · ∥Āη̃∥_{h} + ∥ᾱ_{1h}∥_{1h}^{2} + ∥η̃_{1}∥_{1h}^{2}.

The estimate (3.13) yields ∥Φ̃∥_{h} ≤ K_{ϵ} (∥η̃_{1}∥_{1h} + ∥η̃_{2}∥_{2h}), and from eq. (3.15) follows by means of the fourth assumption

∥η̃_{2}∥_{2h} ≤ ∥Q̃_{22}^{−1}∥_{2h} · ∥ᾱ_{2h}∥_{2h} + ∥Q̃_{22}^{−1}Q̃_{21}∥_{2h} · ∥η̃_{1}∥_{1h}.  (3.18)

These two estimates are inserted into the first inequality which then can be written

\frac{d}{dt}∥η̃_{1}∥_{1h}^{2} ≤ −\frac{1}{h^{2}} (η̃, (B̄^{⊤}Ā + Ā^{⊤}B̄)η̃)_{h} + [2K_{ϵ} + 2K_{ϵ}∥Q̃_{22}^{−1}Q̃_{21}∥_{2h} + K_{ϵ}∥Q̃_{22}^{−1}∥_{2h} + 1] ∥η̃_{1}∥_{1h}^{2}
+∥ᾱ_{1h}∥_{1h}^{2} + K_{ϵ}∥Q̃_{22}^{−1}∥_{2h} · ∥ᾱ_{2h}∥_{2h}^{2}.

Using assumptions (4) and (5) and a simple version of the Gronwall lemma (see Appendix) ∥η̃_{1}∥_{1h}^{2} can be estimated:

∥η̃_{1}(t)∥_{1h}^{2} + Z_{h}^{2}(t) ≤ e^{Ct}[∥η̃_{1}(0)∥_{1h}^{2} + ∫_{0}^{t} e^{−Cτ} (∥ᾱ_{1h}(τ)∥_{1h}^{2} + K_{ϵ}K_{22}∥ᾱ_{2h}(τ)∥_{2h}^{2}) dτ ].

With this relation an estimate for ∥η̃_{2}∥_{2h} can now be obtained from inequality (3.18):

∥η̃_{2}(t)∥_{2h} + K_{21}Z_{h}(t) ≤ K_{22}∥ᾱ_{2h}∥_{2h} + \sqrt{2}K_{21}e^{Ct/2}[∥η̃_{1}(0)∥_{1h}^{2}
+ ∫_{0}^{t} e^{−Cτ} (∥ᾱ_{1h}(τ)∥_{1h}^{2} + K_{ϵ}K_{22}∥ᾱ_{2h}(τ)∥_{2h}^{2}) dτ ]^{1/2}.

Combining these two estimates with assumptions (3) and (6) and the continuity of η̃(t) we obtain the result stated in (3.17). ⊠

As explained above, it is η̃(t) = η̄(t) for t ∈ (0, t_{ϵ}). Since η = S_{1}η̄ we have also an estimate for η(t) needed in the definition of the differential time index (definition 2.1).

Remark 3.2. The order relation (3.17) makes sense since the second assumption guarantees that the function Z_{h}(t) is not negative. In the special case that the matrices A and B would be positive definite and negative definite, respectively (i.e. the system (1.1) is parabolic), assumption (2) is satisfied. □

4. Convergence of the time discretized solution

Because the definition 2.2 of the differential space index is based on the convergence of the error norm ∥ξ∥_{τ} → 0, we consider in this section the convergence problem of the system (2.6) for τ → 0 in the case that the differential space index ν_{dx} is 1. Let

D = (\begin{matrix}1 & 0 \\ −1 & 1 & 0 \\ . & . & . \\ & 0 & −1 & 1 \end{matrix} ) ∈ ℝ^{nM×nM},  (4.1)

E(U(x)) = (F^{⊤}(u^{1}(x)), ..., F^{⊤}(u^{M}(x)))^{⊤} ∈ ℝ^{nM},

a(x) = \frac{1}{τ} (I_{M} ⊗ A) (u^{0⊤}(x), 0, ..., 0)^{⊤} ∈ ℝ^{nM},

F_{R} = (f^{⊤}(t_{1}, x), ..., f^{⊤}(t_{M}, x))^{⊤} ∈ ℝ^{nM}.

Then the system (2.6) can be written compactly

(I_{M} ⊗ B) U_{xx} + \frac{1}{τ} (D ⊗ A) U + E(U) = F_{R} + a(x)  (4.2)

where the vector of unknowns U is defined in (2.8), and the BCs for u_{j}, j ∈ 𝔐_{BC}, are given in eq. (2.7). The consistent BCs for u_{j}, j ∉ 𝔐_{BC}, are assumed to be determined by means of the system (2.6) or (4.2). In order to obtain an equation for the error ξ (defined in eq. (2.9)) we apply the technique already used in section 3. The exact solution vector U_{τ} ∈ ℝ^{nM} is inserted into eq. (4.2) with the result

(I_{M} ⊗ B)U_{τ,xx} + \frac{1}{τ}(D ⊗ A)U_{τ} + E(U_{τ}) = F_{R} + a(x) + β_{τ}(x).  (4.3)

β_{τ}(x) is the time truncation error

β_{τ} = τ(I_{M} ⊗ A)γ_{τ}(x)  (4.4)

where

γ_{τ}(x) = −\frac{1}{2} \frac{∂^{2}}{∂t^{2}} (u_{11}(ζ_{11}, x), ..., u_{1n}(ζ_{1n}, x), u_{21}(ζ_{21}, x), ..., u_{2n}(ζ_{2n}, x), ...,
u_{M1}(ζ_{M1}, x), ..., u_{Mn}(ζ_{Mn}, x))^{⊤},

ζ_{ij} ∈ (t_{i−1}, t_{i}), i = 1, ..., M, j = 1, ..., n. u_{ij} is the component j at the time grid point i. Let 𝔐_{BC}^{R} be the index set (according to definition 1.3) for the DAE (4.2). Comparing eqs. (4.2) and (4.3) we see that we can choose the BCs

U_{j}(±l) = U_{τj}(±l) for j ∈ 𝔐_{BC}^{R}.  (4.5)

(Note that the exact solution U_{τ} is considered to be known here.) This can also be verified by transforming both eqs. (4.2) and (4.3) to a semi explicit representation by means of the same transformation matrices R_{0}, R_{1} defined below. Subtracting eq. (4.2) from eq. (4.3) we obtain the desired equation for ξ:

(I_{M} ⊗ B)ξ_{xx} + (\frac{1}{τ}D ⊗ A + J_{R}(ξ)) ξ = β_{τ}(x).  (4.6)

Here, we have written E(U_{τ}) − E(U) = E(U_{τ}) − E(U_{τ} − ξ) = J_{R}(ξ)ξ with the Jacobian

J_{R}(ξ) = ∫_{0}^{1} E_{U}(U_{τ} + θξ) dθ.  (4.7)

The BCs for ξ should satisfy

ξ_{j}(±l) = 0 for j ∈ 𝔐_{BC}^{R}.  (4.8)

To obtain an error estimate we first transform eq. (4.6) to a semi–explicit form. Then we recognize those components of the transformed error for which arbitrary BCs (taken as zero) can be prescribed. These we consider to belong to the Sobolew space H_{0}^{1}(Ω) such that the Poincaré inequality

∥v∥_{L^{2}(Ω)} ≤ α_{1} ∥\frac{dv}{dx}∥_{L^{2}(Ω)}, v ∈ H_{0}^{1}(Ω),  (4.9)

(α_{1} independent of v) can be applied to such component of the error vector. From this we obtain (under some further assumptions) an estimate of the global time discretization error.

Let R_{0}, R_{1} be regular (nM, nM)−matrices which transform I_{M} ⊗ B into

(\begin{matrix}I_{1} & 0 \\ 0 & 0 \end{matrix}) = R_{0}(I_{M} ⊗ B)R_{1}

where the unit matrix I_{1} is of order M_{1}. M_{2} is then used to denote M_{2} = nM − M_{1}. Define

ω := R_{1}^{−1}ξ,  D̄(ω) := R_{0} (\frac{1}{τ}D ⊗ A + J_{R}(R_{1}ω)) R_{1},  d̄:= R_{0}β_{τ}.  (4.10)

Furthermore, we decompose ω = (ω_{1}^{⊤}, ω_{2}^{⊤})^{⊤} , ω_{i} ∈ ℝ^{M_{i}}, i = 1, 2. Other vectors in ℝ^{nM} and some (nM, nM)−matrices needed in the sequel are decomposed in a similar way. Then the system (4.6) can be written in the semi–explicit form

ω_{1,xx} + D̄_{11}(ω)ω_{1} + D̄_{12}(ω)ω_{2} = d̄_{1}  (4.11)
D̄_{21}(ω)ω_{1} + D̄_{22}(ω)ω_{2} = d̄_{2}.  (4.12)

We see that the matrices R_{0}, R_{1} (which only depend on I_{M} ⊗ B and not of J_{R}) lead to the result that arbitrary BCs can be prescribed for ω_{1}. Thus we take the BCs for ω_{1} as

ω_{1}(±l) = 0.  (4.13)

The consistent BCs for ω_{2} must be determined from eq. (4.12). Nonlinear BVPs for DAEs of the type (4.11), (4.12) (of indexes 1 and 2) are considered in [16], [17]. At this point we assume now that the system (4.11), (4.12) is of differential index 1. Therefore, we suppose that eq. (4.12) has a solution for ω_{2} which we write ω_{2} = φ(ω_{1}). If the PDAE (1.1) is linear, then D̄ is independent of ω (J_{R}, eq. (4.7), does not depend on ξ), and ν_{dx} = 1 means that the square matrix D̄_{22} is regular. The function φ is in this case simply ω_{2} = φ(ω_{1}) = D̄_{22}^{−1} (d̄_{2} − D̄_{21}ω_{1}). With this information an estimate for ω can be obtained from the system (4.11), (4.12).

In the nonlinear case of index ν_{dx} = 1 we insert ω_{2} = φ(ω_{1}) into eq. (4.11) and write the resulting equation

ω_{1,xx} + H(ω_{1}) = d̄_{1},  (4.14)
H(ω_{1}) := D̃_{11}(ω_{1})ω_{1} + D̃_{12}(ω_{1})φ(ω_{1})  (4.15)

where D̃_{ij}(ω_{1}) := D̄_{ij}(ω_{1}, φ(ω_{1})), i, j ∈ {1, 2}. Now we suppose that each component ω_{1j}, j = 1, ..., M_{1}, belongs to the Sobolow space H_{0}^{1}(Ω). To exploit this we multiply eq. (4.14) from the left by ω_{1}^{⊤} and integrate with respect to x from −l to l. A partial integration and the BC (4.13) yield

∫_{−l}^{l} ω_{1,x}^{⊤}ω_{1,x} dx − ∫_{−l}^{l} ω_{1}^{⊤}H(ω_{1}) dx = − ∫_{−l}^{l} ω_{1}^{⊤}d̄_{1} dx.  (4.16)

Let scalar products and norms (according to the discrete L^{2}−scalar product (2.10)) be defined by

(v, w)_{iτ} := ∑_{ν=1}^{M_{i}} τv_{ν}w_{ν}, ∥v∥_{iτ} := \sqrt{(v, v)_{iτ}}, i = 1, 2,  (4.17)

(v, w)_{iτ,L^{2}} := ∫_{−l}^{l}(v(x), w(x))_{iτ} dx, ∥v∥_{iτ,L^{2}} := \sqrt{(v, v)_{iτ,L^{2}}}, i = 1, 2.  (4.18)

Note that these definitions imply that the natural matrix norm of the form ∥I_{M_{i}} ⊗ Q∥_{iτ,L^{2}} is independent of M_{i} (Q ∈ ℝ^{n×n}).

We multiply eq. (4.16) by τ, apply Poincarés inequality (4.9) for each component of ω_{1}, use the Cauchy–Schwarz inequality and the definitions (4.18) to get the basic inequality

α_{1}∥ω_{1}∥_{1τ,L^{2}} − \frac{(ω_{1}, H(ω_{1}))_{1τ,L^{2}}}{∥ω_{1}∥_{1τ,L^{2}}} ≤ ∥d̄_{1}∥_{1τ,L^{2}}.  (4.19)

As is obvious, with this relation we get under some assumptions a convergence result summarized in

Theorem 4.1. Suppose

(1) the solutions of the PDAE (1.1) and the time discretized equation (4.2) are sufficiently smooth,

(2) eq. (4.12) can be solved for ω_{2},

(3) (ω_{1}, H(ω_{1}))_{1τ,L^{2}} ≤ 0 for τ → 0. Then

lim_{τ→0} ∥ω_{1}∥_{1τ,L^{2}} = 0.  (4.20)

Note added in proof: The first assumption guarantees that the components ω_{1j}, j = 1, ..., M_{1}, can really be seen as elements of the space H_{0}^{1}(Ω). Assumption (2) is the condition for a system with differential space index 1. The third assumption (discussed further below) implies that the second term (including the minus sign) on the left hand side of ineq. (4.19) is not negative. The right hand side goes to zero for τ → 0 because d̄∼ β_{τ} = O(τ).

To investigate the meaning of the assumption 3 it is instructive to consider linear constant coefficient PDAEs with F (u) = Cu where C ∈ ℝ^{n×n} is a constant matrix. In this case D̃ = D̄ is independent of ω_{1}, and D̄_{22} must be regular. Inserting φ(ω_{1}) = D̄_{22}^{−1} (d̄_{2} − D̄_{21}ω_{1}) into eq. (4.15) we obtain

H(ω_{1}) = H_{lin}ω_{1} + D̄_{12}D̄_{22}^{−1}d̄_{2}, H_{lin} := D̄_{11} − D̄_{12}D̄_{22}^{−1}D̄_{21},  (4.21)

such that eq. (4.14) can be rewritten

ω_{1,xx} + H_{lin}ω_{1} = d̄_{1} − D̄_{12}D̄_{22}^{−1}d̄_{2}.  (4.22)

If now ∥D̄_{12}D̄_{22}^{−1}∥_{1τ,L^{2}} is bounded for τ → 0 (which is a realistic assumption in many applications), then it is also H(ω_{1}) → H_{lin}ω_{1} for τ → 0. Under which conditions is assumption (3) of the theorem satisfied? One part of the answer can be found by means of the simple

Lemma 4.1. If the eigenvalues of the square matrix H_{lin} are not positive, then the BVP (4.22) with ω_{1}(±l) = 0, has an unique solution.

The statement can be proved by transforming H_{lin} to Jordan normal form. Then it is easily seen that the homogeneous BVP (zero right hand side) has only the trivial solution for non positive eigenvalues of H_{lin}. If furthermore the matrix H_{lin} is diagonalizable by an orthogonal transformation matrix, then assumption (3) is fulfilled. Therefore, we have for linear constant coefficient PDAEs the following result: if ∥D̄_{12}D̄_{22}^{−1}∥_{1τ,L^{2}} is bounded for τ → 0 and if the BVP (4.14), (4.13) (or equivalently the BVP (4.22), (4.13)) has in any case an unique solution and if the matrix H_{lin} is diagonalizable (by an orthogonal transformation), then assumption (3) of the theorem is satisfied.

To discuss assumption (3) of the theorem for semilinear PDAEs we restrict ourselves to the case where again the matrix

D̃_{22}(ω_{1}) := D̄_{22}(ω_{1}, φ(ω_{1})) is regular.  (4.23)

Under this assumption, eq. (4.11) can be written

ω_{1,xx} + H_{nl}(ω_{1})ω_{1} = r_{1}(ω_{1})  (4.24)

with the abbreviations

H_{nl}(ω_{1}) :=D̃_{11}(ω_{1}) − D̃_{12}(ω_{1})D̃_{22}^{−1}(ω_{1})D̃_{21}(ω_{1})  (4.25)
r_{1}(ω_{1}) :=d̄_{1} − D̃_{12}(ω_{1})D̃_{22}^{−1}(ω_{1})d̄_{2}.  (4.26)

It is seen immediately that the BVP (4.24)–(4.26) (with BC (4.13)) is formally in its structure similar to the BVP (4.22), (4.21) (with the same BC). Therefore, a linearization, say about ω_{1} = 0, gives in lowest order the special case already considered. In general, under the condition (4.23) we can say that assumption (3) of the theorem is satisfied, e.g., if the eigenvalues of the matrix H_{nl}(ω_{1}) are negative for τ → 0 and if it can be diagonalized by an orthogonal transformation.

The third assumption of theorem 4.1 is fulfilled also when H can be seen as a strong monotone operator on a (suitably chosen) domain D_{H} (a possibility is D_{H} ⊆ (H_{0}^{1}(Ω))^{M_{1}} with 0 ∈ D_{H}), i.e.

(ω_{1} − ω̄_{1}, H(ω_{1}) − H(ω̄_{1}))_{1τ,L^{2}} ≤ C_{0} < 0 for ω_{1}, ω̄_{1} ∈ D_{H},  (4.27)

under the condition

(ω_{1}, H(0))_{1τ,L^{2}} = 0 for τ → 0.  (4.28)

C_{0} is a given constant independent of the elements of D_{H}. The last relation is fulfilled, e.g., in the linear case (4.21) under the condition that D̄_{12}D̄_{22}^{−1} is bounded for τ → 0 (there is a class of PDAEs where this holds).

Lemma 4.2. Suppose that the relations (4.27) and (4.28) are valid. Then

(1) the solution of the BVP (4.14), (4.13) is unique,

(2) the third assumption of theorem 4.1 is satisfied.

Proof. To prove the first statement we suppose that there are two different solutions ω_{1}, ω_{1}^{∗} ∈ D_{H} of the BVP (4.14), (4.13). Subtracting first the equation for ω_{1}^{∗} from the equation for ω_{1} and multiplying the resulting equation with ω_{1} − ω_{1}^{∗} we get after integration ∫_{−l}^{l} ... dx

− (ω_{1,x} − ω_{1,x}^{∗}, ω_{1,x} − ω_{1,x}^{∗})_{1τ,L^{2}} + (ω_{1} − ω_{1}^{∗}, H(ω_{1}) − H(ω_{1}^{∗}))_{1τ,L^{2}} = 0.

Eq. (4.27) implies that this relation is a contradiction unless ω_{1} = ω_{1}^{∗}. The second statement follows from eq. (4.27) when we set ω̄_{1} = 0 and take eq. (4.28) into account. This gives (ω_{1}, H (ω_{1}))_{1τ,L^{2}} ≤ C_{0} < 0 for τ → 0, i.e. the third assumption of theorem 4.1 is satisfied. ⊠

An error estimate for ω_{2} can be obtained from eq. (4.12) (or from ω_{2} = φ(ω_{1})) and theorem 4.1. In the constant coefficient case, ∥ω_{2}∥_{2τ,L^{2}} can be estimated by

∥ω_{2}∥_{2τ,L^{2}} ≤ ∥D̄_{22}^{−1}∥_{2τ,L^{2}} · ∥d̄_{2}∥_{2τ,L^{2}} + ∥D̄_{22}^{−1}D̄_{21}∥_{2τ,L^{2}} · ∥ω_{1}∥_{1τ,L^{2}}

such that convergence can be shown if D̄_{22}^{−1}D̄_{21} is bounded for τ → 0 (provided ∥ω_{1}∥_{1τ,L^{2}}

converges). On the whole we have an error estimate for ω = (ω_{1}^{⊤}, ω_{2}^{⊤})^{⊤} = R_{1}^{−1}ξ where the estimate for ξ is needed in the index definition 2.2.

5. Conclusion

It was the aim of this paper to define for semilinear PDAEs indexes which can be used in theory and numerical solution. The method proposed in this paper for the definition of differential time and space indexes is based on special (often applied) semidiscretization schemes and uses certain analogies between the Fourier and Laplace transformation on the one hand and of the vertical and the horizontal line method on the other hand. Although we have (up to now) no proof that the indexes are independent of the semidiscretization chosen it is expected that the indexes of a given PDAE are the same at least for certain classes of discretizations. An essential assumption in our index definitions is the convergence of the semidiscretized solution towards the exact solution (assumed to exist and to be unique). We have given two theorems containing conditions for convergence (for the semidiscretization scheme considered here). In the special case that the given PDAE is linear with constant coefficients, our index definitions for semilinear PDAEs are essentially identical with the traditional index definitions of linear PDAEs.

Appendix

The following version of the Gronwall lemma is used in section 3: Lemma: Suppose y ∈ C^{1}([0, T]), ψ ∈ C([0, T]), satisfies

y′(t) ≤ Cy(t) + ψ(t), C = const. ≥ 0.

Then

y(t) ≤ e^{Ct} [y(0) + ∫_{0}^{t} e^{−Cτ}ψ(τ) dτ ] .

Proof. Let z(t) := e^{−Ct}y(t). Then

z′(t) = −Ce^{−Ct}y(t) + e^{−Ct}y′(t) ≤ e^{−Ct}ψ(t).

Integration yields

z(t) − z(0) ≤ ∫_{0}^{t} e^{−Ct}ψ(τ) dτ.

Inserting here the definition of z(t) we get the lemma. ⊠

Acknowledgement

The authors are indebted to Prof. Dr. S. L. Campbell for his helpful comments on a first version of this paper.

References

1. M. Arnold. A note on the uniform perturbation index. interner Bericht, Universität Rostock, 1995.

2. K.E. Brenan, S.L. Campbell, and L.R. Petzold. Numerical solution of initial-value problems in differentialalgebraic equations. North-Holland Publ. Co., Amsterdam, 1989.

3. S. L. Campbell. DAE approximations of PDE modeled control problems. Proc. IEEE mediterranean symposium on new directions in control and automation, Crete, pages 407–414, 1994.

4. S.L. Campbell and C.W. Gear. The index of general nonlinear DAEs. Numer. Math., 72:173–196, 1995.

5. S.L. Campbell and W. Marszalek. DAEs arising from travelling wave solutions of PDEs. 1996.

6. S.L. Campbell and W. Marszalek. The index of an infinite dimensional implicit system. Math. Mod. of Syst., 1(1):1–25, 1996.

7. S.L. Campbell and W. Marszalek. ODE/DAE integrators and MOL problems. ZAMM, 76:S1, 251–254, 1996.

8. E. Griepentrog and R. März. Differential-algebraic equations and their numerical treatment. Teubner-Texte zur Mathematik, Band 88, Leipzig, 1986.

9. P. Grindrod. The theory and applications of reaction–diffusion equations. Clarendon Press, Oxford, 1996.

10. E. Hairer and G. Wanner. Solving Ordinary differential equations II. Springer-Verlag, 1996.

11. H.O. Kreiss. Initial-boundary value problems and the Navier-Stokes equations. Academic Press, Inc., Boston, 1989.

12. A.W. Leung. Systems of nonlinear partial differential equations. Kluwer Acad. Publ., Dortrecht, Boston, London, 1989.

13. W. Lucht, K. Strehmel, and C. Eichler-Liebenow. Linear partial differential algebraic equations, Part I: Indexes, consistent boundary/initial conditions. Report 17, Martin-Luther-Universität Halle, Fachbereich Mathematik und Informatik, submitted for publication, 1997.

14. W. Lucht, K. Strehmel, and C. Eichler-Liebenow. Linear partial differential algebraic equations, Part II: Numerical solution. Report 18, Martin-Luther-Universität Halle, Fachbereich Mathematik und Informatik, submitted for publication, 1997.

15. W. Marszalek. Analysis of partial differential algebraic equations, Ph. D. Thesis, North Carolina State University, Raleigh. 1997.

16. R. März. On linear differential-algebraic equations and linearizations. Appl.Num.Math., (18):267–292, 1995.

17. R. März. Extra-ordinary differential equations. Attempts to an analysis of differential-algebraic systems. Preprint Nr. 97-8, Humboldt-Universität zu Berlin, Mathematisch-Naturwissenschaftliche Fakultät II, Institut für Mathematik, 1997.

18. J.D. Murray. Mathematical biology. Springer–Verlag Berlin, 1990.

19. K. Rektorys. The method od discretization in time and partial differential equations. D. Reidel Publ. Comp. Dordrecht, 1982.

20. B. Simeon. Modelling a flexible slider crank mechanism by a mixed system of DAEs and PDEs. Math. Mod. of Syst., 2(1):1–18, 1996.