This is a reprint of the paper

The numerical solution of linear time-varying

DAEs with index 2 by IRK methods

by Ebroul Izquierdo

published in Rev. Colombiana Mat.

29 (1995), pp. 49–88

The numerical solution of

linear time-varying DAEs

with index 2 by IRK methods

Ebroul Izquierdo

Heinrich Hertz Institut, Berlin

Abstract. Differential-algebraic equations (DAEs) with a higher index can be approximated by implicit Runge-Kutta methods (IRK). Until now, a number of initial value problems have been approximated by Runge-Kutta methods, but all these problems have a special semi-explicit or Hessenberg form. In the present paper we consider IRK methods applied to general linear time-varying (nonautonom) DAEs tractable with index 2.

For some stiffly accurate IRK formulas we show that the order of accuracy in the differential component is the same nonstiff order, if the DAE has constant nullspace. We prove that IRK methods cannot be feasible or become exponentially unstable when applied to linear DAEs with variable nullspace. In order to overcome these difficulties we propose a new approach for this case. Feasibility, weak instability and convergence are proved. Order results are given in terms of the Butcher identities. Keywords and phrases. Ordinary differential equations, differential-algebraic equations, initial value problems, implicit Runge-Kutta methods. 1991 Mathematics Subject Classification. 65L06, 34D20

1. Introduction

In this paper we consider linear differential algebraic equations (DAEs) of the form

A(t)x′(t) + B(t)x(t) = q(t).  (1.1 a)

We assume that the coefficients A, B in (1.1 a) are continuous matrix functions A, B : [t_{0}, T] ⟶ L(ℝ⅁⋖∼^{⋗}) and the matrix A(t) has a smooth nullspace N(t) for t ∈ [t_{0}, T], i.e., there are continuously differentiable projection functions Q, P : [t_{0}, T] ⟶ L(ℝ^{⋗}) so that Q(t) is a projection onto the nullspace N(t) and P (t) ≡ I − Q(t). Note that a smooth nullspace N(t) has always constant dimension on [t_{0}, T] and A(t) has constant rank. Such equations arise in a variety of applications, e.g., electrical network theory, dynamical systems subjects to constraints, optimal control of lumped-parameter systems and reduced equations in singulary perturbed systems.

In the present paper, we are interested in index 2 equations (more precisely, DAEs being tractable with index 2). Besides the BDF implicit Runge-Kutta methods are popular for approximating initial value problems in DAEs. Since higher index DAEs are known to be ill-posed, some kind of instability should be expected. Nevertheless, some IRK methods are reported to be adequate for some special types of higher index equations. Especially, the so-called Hessenberg systems are well investigated ([1], [4], [9]). We will generalize these results to the general linear DAE (1.1 a). In order to overcome the inherent difficulties when the nullspace of A varies we introduce a new approach for IRK methods using the projection Q. Doing so we can show that these methods share all properties which have IRK methods when applied to Hessenberg systems even for general DAEs (1.1 a).

Following the idea of März [10] we use the matrix functions

G_{1} := A + BQ,
A_{1} := G_{1} − AP′Q,
G_{2} := A_{1} + BP Q_{1},

as well as the subspaces

N(t) := ker(A(t)),
S(t) := {x ∈ ℝ^{⋗} : 𝔹(≈)↶ ∈ im (𝔸(≈))},
N_{1}(t) := ker(A_{1}(t)),
S_{1}(t) := {x ∈ ℝ^{⋗} : 𝔹(≈)ℙ(≈)↶ ∈ im (𝔸_{⊮}(≈))}.

Here, t ∈ [t_{0}, T] and Q_{1}(t) denotes a projection onto N_{1}(t); further P_{1}(t) := I−Q_{1}(t). The subspaces N, S, N_{1} and S_{1} are called the canonical the subspaces of the DAE (1.1 a) (see Lemma A.1).

Definition 1.1 (März [11]) The DAE (1.1 a) is called transferable if G_{1}(t) is nonsingular for every t ∈ [t_{0}, T]. The DAE (1.1 a) is called tractable with index 2 if G_{1}(t) is singular, but G_{2}(t) is nonsingular for every t ∈ [t_{0}, T].

Note that G_{1}(t) is nonsingular iff so is A_{1}(t). Transferable DAEs are well understood and suitably modified numerical ODE methods work well for them. We refer to [1], [3] for this case. Nontransferable DAEs are essentially more complex. Tractability with index 2 characterizes an important class of nontransferable DAEs. Maybe the most important qualitative difference between transferable and nontransferable DAEs is that the first class of problems remain well-posed in the Hadamard sense, whereas initial and boundary value nontransferable problems are ill-posed in the sense of Tikhonov. This fact has to be taken into account for the numerical treatment.

The formulation of appropriate initial conditions for index 2 tractable DAEs depends on the DAE itself. The index 2 tractable DAE (1.1) subject to the condition

P (t_{0})P_{1}(t_{0})x(t_{0}) = b  (1.1 b)

represent an initial value problem (IVP) with consistent initial values.

According to the originally conceived method for the numerical solution of ordinary differential equations (Kutta, 1901), an IRK method can be realized for the DAE (1.1 a) in the following way: from an approximation x_{n−1} of the solution of the IVP (1.1) at t_{n−1}, these one-step methods construct an approximation x_{n} at t_{n} = t_{n−1} + h via the formulas

x_{n} = x_{n−1} + h ∑_{j=1}^{s} b_{j}X′_{j},  (1.2 a)

where X′_{j} is defined by

A(t_{nj})X′_{j} + B(t_{nj})X_{j} = q(t_{nj})  (1.2 b)

with t_{nj} := t_{n−1} + c_{j}h and internal stages X_{j} given by

X_{j} = x_{n−1} + h ∑_{l=1}^{s} a_{jl}X′_{l}, j = 1, 2, …, s.  (1.2 c)

Here a_{jl}, b_{j}, c_{j} are the coefficients determined by the method, and s is the number of stages. Usually, the matrix A = (a_{jl})_{j,l=1}^{s} and the vectors b = (b_{1}, b_{2}, …, b_{s})^{T}, c = (c_{1}, c_{2}, …, c_{s})^{T} are combined in the so-called Butcher Diagram [2] or Runge-Kutta scheme

c A b^{T}  .

The paper is organized as follows. In Section 2 we give conditions which ensure the feasibility of IRK methods. These conditions depend not only on the method, but essentially on the behaviour of the nullspace N(t). This is similar to what is known for the BDF [13]. In Section 3 we consider the case of a constant nullspace N(t) ≡ N. We show the weak instability and give order results in terms of the Butcher identities. In Section 4 we show that similar results are not true if the nullspace is variable. Therefore, for this case, we introduce a new approximation which makes use of the projection Q. In Section 5 we show that these new approximations have the same properties as the original IRK methods in the constant nullspace case. Section 6 contains some numerical examples.

Assumption: In the following, Q_{1}(1) denotes the canonical projection onto N_{1}(t) along S_{1}(t) and for any matrix function M(t) is

D_{M} = diag (M(t_{n1}), …, M(t_{ns})).

Moreover, we assume that all matrix functions and projections are sufficiently smooth.

2. Existence and uniqueness of the Runge-Kutta solution

First we study the existence and uniqueness of a Runge-Kutta solution, i.e., we try to answer the question: when does the system (1.2 b) have a unique solution X_{j} (or X′_{j}) for j = 1, 2, …, s?

Suppose that the Runge-Kutta matrix A is nonsingular. Then, the IRK method (1.2) is equivalent to

x_{n} = ρx_{n−1} + ∑_{j=1}^{s} ∑_{l=1}^{s} b_{j}â_{jl}X_{l},  (2.1 a)

A(t_{nj}) ∑_{l=1}^{s} â_{jl}(X_{l} − x_{n−1}) + hB(t_{nj})X_{j} = hq(t_{nj}), j = 1, 2, …, s,  (2.1 b)

where A^{−1} = (â_{jl})_{j,l=1}^{s} and ρ := 1 − ∑_{i=1}^{s} ∑_{j=1}^{s} b_{j}â_{ij} ([6], [8]).

We now consider the IRK method defined by (2.1). Let us introduce some notation.

D_{A} = diag (A(t_{n1}), …, A(t_{ns})), D_{B} = diag (B(t_{n1}), …, A(t_{ns})),
Dq = diag (q(t_{n1}), …, q(t_{ns})), 11_{s} = (1, …, 1)^{T} ∈ ℝ^{∼},
\overline{X} = (X_{1}^{T} , X_{2}^{T} , …, X_{s}^{T})^{T} ∈ ℝ^{⋗∼}.

Now, system (2.1 b) can be written in the following form

[D_{A}(A^{−1} ⊗ I_{m}) + hD_{B}]\overline{X} = hDq11_{s} + D_{A}(A^{−1}11_{s} ⊗ x_{n−1}).  (2.2)

The IRK method (2.1) is feasible if the matrix M := D_{A}(A^{−1} ⊗ I_{m}) + hD_{B} is nonsingular. Unfortunately, simple examples show that M is not always so.

Example 2.1. Consider

A(t) = ( \begin{matrix}0 & 0 & 0 \\ t & e^{−t} & 0 \\ 0 & 0 & t \end{matrix} ) , B(t) = ( \begin{matrix}t & e^{−t} & 0 \\ 0 & 0 & 3 \\ 0 & 0 & t + 1 \end{matrix} ) .

The related DAE is transferable with index 2 for t ∈ (0, ∞) and

Q = ( \begin{matrix}0 & −e^{−t}/t & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{matrix} )

is a projector onto N(t) = ker (A(t)). A generic 1-stage IRK method yields

M = A(t)(a^{−1} ⊗ I_{3}) + hB(t) = ( \begin{matrix}ht & he^{−t} & 0 \\ a^{−1}t & a^{−1}e^{−t} & 3h \\ 0 & 0 & a^{−1}t + ht + h \end{matrix} ) .

Obviously, this matrix is singular for all t, h ∈ (0, ∞) and all a ∈ ℝ \ {⊬}. □

The following theorem provides a necessary and sufficient condition for the existence of M^{−1}. We present here a detailed proof which makes frequent use of the statements of Lemma A.2. For a more detailed proof we refer the reader to [7], [8].

Theorem 2.1 Suppose that the IVP (1.1) is tractable with index 2 and moreover, that A and I_{ms} − ED_{QQ_{1}P′} are nonsingular for all t ∈ [t_{0}, T], where Q_{1}(t) denotes the canonical projection onto N_{1}(t) along S_{1}(t). Then the system (2.2) is uniquely solvable for sufficiently small h. Here

E := I_{ms} + [(C ⊙ A^{−1}) ⊗ I_{m}], C = (c_{jl})_{j,l=1}^{s} ∈ L(ℝ^{∼}), c_{jl} = c_{j} − c_{l}, and

C ⊙ A^{−1} = (c_{jl}â_{jl})_{j,l=1}^{s} ∈ L(ℝ^{∼}) is the Hadamard product of matrices.

Proof. First we decouple (2.2) in the P P_{1}–, P Q_{1}– and Q-components, using the projection matrices defined in Section 1, and Lemma A.2. Finally, we calculate M^{−1}.

System (2.2) is equivalent to

[D_{AP} (A^{−1} ⊗ I_{m}) + hD_{AQ} + D_{BQP} (A^{−1} ⊗ I_{m}) + hD_{BQ} + hD_{BP} ]\overline{X}

= hDq11_{s} + D_{A}(A^{−1}11_{s} ⊗ x_{n−1}).

Using G_{1}, A_{1} and G_{2} (see Section 1) we can write this equation as:

[(D_{A_{1}} + D_{BP Q_{1}})(D_{P_{1}P} (A^{−1} ⊗ I_{m})+hD_{P_{1}Q}+hD_{Q_{1}})+hD_{AP′Q} + hD_{BP P_{1}}]\overline{X}

= hDq11_{s} + D_{A}(A^{−1}11_{s} ⊗ x_{n−1}).

Since (1.1) is tractable with index 2, there exists G_{2}^{−1} = (A_{1} + BP Q_{1})^{−1}. Hence, after multiplying the latter equation by G_{2}^{−1} we obtain:

[D_{P_{1}P} (A^{−1} ⊗ I_{m}) + hD_{P_{1}P P′Q}+hD_{P_{1}Q}+hD_{Q_{1}} +hD_{G_{2}^{−1}AP′Q}+hD_{G_{2}^{−1}BP P_{1}}]\overline{X}

= hD_{G_{2}^{−1}q}11_{s} + D_{G_{2}^{−1}A}(A^{−1} ⊗ I_{m})(11_{s} ⊗ x_{n−1}).

We multiply the latter equation by D_{P P_{1}}, \frac{1}{h}D_{QP_{1}} and \frac{1}{h}D_{Q_{1}} and use Lemma A.2 together with the following identities:

D_{Q_{1}Q} = 0, D_{Q_{1}} = D_{Q_{1}G_{2}^{−1}BP} , D_{P P_{1}P} = D_{P P_{1}}, D_{P P_{1}Q} = 0,
D_{QP_{1}Q} = D_{Q}, D_{QP_{1}P} = −D_{QQ_{1}}, D_{Q_{1}G_{2}^{−1}A} = 0,

Note that most of these relations are evident. For a detailed proof we refer the reader to [8]. See also Lemma A.1.

Thus, we obtain the system

((A^{−1} ⊗ I_{m}) + hD_{P P_{1}G_{2}^{−1}B} + R_{P P_{1}})D_{P P_{1}} \overline{X}
+ R_{P P_{1}}D_{P Q_{1}} \overline{X} + (hD_{P P_{1}P′} + R_{P P_{1}})D_{Q}\overline{X}  (2.3 a)
= hD_{P P_{1}G_{2}^{−1}q}11_{s} + D_{P P_{1}}(A^{−1} ⊗ I_{m})(11_{s} ⊗ x_{n−1}),
D_{Q_{1}} \overline{X} = D_{Q_{1}G_{2}^{−1}q}11_{s},  (2.3 b)

D_{QP_{1}G_{2}^{−1}BP P_{1}} \overline{X} − \frac{1}{h}(A^{−1} ⊗ I_{m})D_{QQ_{1}} \overline{X} + (I_{ms} − D_{QQ_{1}P′})D_{Q}\overline{X}

= D_{QP_{1}G_{2}^{−1}q}11_{s} − \frac{1}{h}D_{QQ_{1}}(A^{−1} ⊗ I_{m})(11_{s} ⊗ x_{n−1}) + \frac{1}{h}R_{QQ_{1}} \overline{X},  (2.3 c)

where

R_{P P_{1}} = h[((C ⊙ A^{−1}) ⊗ I_{m})D_{(P P_{1})′} + O(h)],
\frac{1}{h}R_{QQ_{1}} \overline{X} = ((C ⊙ A^{−1}) ⊗ I_{m})D_{(QQ_{1})′} + O(h))\overline{X}
= ((C ⊙ A^{−1}) ⊗ I_{m})D_{(Q′Q_{1}+QQ′_{1})}((D_{P P_{1}} + D_{P Q_{1}} + D_{Q} + O(h))\overline{X}.

Inserting (2.3 a) and (2.3 b) into the last equation and using the relations Q′Q_{1}P P_{1} = 0 and Q_{1}Q = 0 we obtain

\frac{1}{h}R_{QQ_{1}} \overline{X} = V + ((C ⊙ A^{−1}) ⊗ I_{m})D_{QQ′_{1}Q}\overline{X} + O(h∥\overline{X}∥),  (2.4)

where

V = ((C ⊙ A^{−1}) ⊗ I_{m})D_{QQ′_{1}}[h(A ⊗ I_{m})D_{P P_{1}G_{2}^{−1}q}11_{s}
+ (A ⊗ I_{m})D_{P P_{1}}(A^{−1} ⊗ I_{m})(11_{s} ⊗ x_{n−1})]
+ ((C ⊙ A^{−1}) ⊗ I_{m})D_{(QQ_{1})′P Q_{1}G_{2}^{−1}q}11_{s} .

We insert (2.4) into (2.3 c). The equation obtained together with (2.3 a-b) can be written

F [ \begin{matrix}D_{P P_{1}} \overline{X} \\ D_{Q_{1}} \overline{X} \\ D_{Q}\overline{X} \end{matrix} ] = J + [ \begin{matrix}O(h) & O(h) & −hS + O(h^{2}) \\ 0 & 0 & 0 \\ O(h) & O(h) & O(h) \end{matrix} ] [ \begin{matrix}D_{P P_{1}} \overline{X} \\ D_{P Q_{1}} \overline{X} \\ D_{Q}\overline{X} \end{matrix} ] ,  (2.5)

where

F = [ \begin{matrix}I_{m}s & 0 & 0 \\ 0 & I_{ms} & 0 \\ D_{QP_{1}G_{2}^{−1}B} & −\frac{1}{h}(A^{−1} ⊗ I_{m})D_{Q} & D \end{matrix} ] ,

J = ( \begin{matrix}h(A ⊗ I_{m})D_{P P_{1}G_{2}^{−1}q}11_{s} + (A ⊗ I_{m})D_{P P_{1}}(A^{−1} ⊗ I_{m})(11_{s} ⊗ x_{n−1}) \\ D_{Q_{1}G_{2}^{−1}q}11_{s} \\ D_{QP_{1}G_{2}^{−1}q}11_{s} − \frac{1}{h}D_{QQ_{1}}(A^{−1} ⊗ I_{m})(11_{s} ⊗ x_{n−1}) + V \end{matrix} ) ,

S = [(A ⊗ I_{m}) + (A(C ⊙ A^{−1}) ⊗ I_{m}]D_{(P P_{1})′Q}.

Note that in (2.5) the identity QQ′_{1}Q = QQ′_{1}QQ = QQ_{1}P′Q has been used.

The matrix F is nonsingular iff D is nonsingular and in this case we can compute F^{−1}:

F^{−1} = [ \begin{matrix}I_{ms} & 0 & 0 \\ 0 & I_{ms} & 0 \\ −D^{−1}D_{QP_{1}G_{2}^{−1}B} & \frac{1}{h}D^{−1}(A^{−1} ⊗ I_{m})D_{Q} & D^{−1} \end{matrix} ] .

Hence, from (2.5) we have

[ \begin{matrix}D_{P P_{1}} \overline{X} \\ D_{P Q_{1}} \overline{X} \\ D_{Q}\overline{X} \end{matrix} ] = [ \begin{matrix}I_{ms} & 0 & 0 \\ 0 & D_{P} & 0 \\ 0 & 0 & I_{ms} \end{matrix} ] F^{−1}J

+ [ \begin{matrix}O(h) & O(h) & O(h) \\ 0 & 0 & 0 \\ O(h) & O(h) & O(h) \end{matrix} ] [ \begin{matrix}D_{P P_{1}} \overline{X} \\ D_{P Q_{1}} \overline{X} \\ D_{Q}\overline{X} \end{matrix} ] .  (2.6)

Using Banach’s lemma we can solve system (2.6) for sufficiently small h. Then we use the identity \overline{X} = D_{P P_{1}} \overline{X} + D_{P Q_{1}} \overline{X} + D_{Q}\overline{X} to construct the unique solution of (2.2). □

Remarks

• If the nullspace N(t) of A(t) is constant, the IRK method is always feasible, for h sufficiently small, because in this case P′ = 0, i.e., Theorem 2.1 ensures the feasibility of IRK methods applied to DAEs with index 2 and constant nullspace. This makes this theorem very valuable because many important classes of DAEs have constant nullspace, e.g., systems in Hessenberg or semi-explicit form. Note that for all these general DAEs, A(t) has the form

A(t) = ( \begin{matrix}I & 0 \\ 0 & 0 \end{matrix} ) ,

so that P = A and trivially P′ = 0.

• It is clear that before we try to solve numerically a DAE with index 2 we must check the nullspace N(t). If N(t) is variable, we have little hope that the method remain convergent (see next Section). In this case we advise to approximate the DAE using additional strategies: e.g., regularization methods [5], or the modified method proposed in Section 5.

• For 1-stage IRK methods it holds that I_{m} − ED_{QQ_{1}P′} = I_{m} − QQ_{1}P′, because in this case C = 0. Hence, to verify the solvability of the system (2.2) it is sufficient to check the nonsingularity of the matrix I_{m} − QQ_{1}P′. In example 2.1,

I_{3} − QQ_{1}P′ = ( \begin{matrix}1 & e^{−t}/t & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{matrix} )

is always singular.

• In [11] it has been proved that BDF methods for the numerical solution of (1.1) are feasible if the matrix I_{m} − QQ_{1}P′ is nonsingular and h is sufficiently small. For IRK methods with s-stages, s ≥ 2, the feasibility of the method depends on the method and the IVP.

• Summarizing, the statements of Theorem 2.1 present an important tool for practical numerical applications as well as for investigations about stability and convergence.

Example 2.2. Let

A(t) = ( \begin{matrix}0 & 0 \\ 1 & ηt \end{matrix} ) , B(t) = ( \begin{matrix}1 & ηt \\ 0 & 1 + η \end{matrix} )

with η ∈ ℝ. The related DAE is transferable with index 2

Q(t) = ( \begin{matrix}0 & −ηt \\ 0 & 1 \end{matrix} ) , A_{1}(t) = ( \begin{matrix}0 & 0 \\ 1 & 1 + ηt \end{matrix} )

and

Q_{1}(t) = ( \begin{matrix}1 + ηt & ηt(1 + ηt) \\ −1 & −ηt \end{matrix} ) .

It is known that the backward Euler method is not feasible if η = 1 (see [11]). If we instead apply the 2-stage Lobatto IIIC-formula [2] to this DAE, we have

A^{−1} = ( \begin{matrix}1 & 1 \\ −1 & 1 \end{matrix} ) , C = ( \begin{matrix}0 & −1 \\ 1 & 0 \end{matrix} ),

E = ( \begin{matrix}1 & 0 & −1 & 0 \\ 0 & 1 & 0 & −1 \\ −1 & 0 & 1 & 0 \\ 0 & −1 & 0 & 1 \end{matrix} )

I_{ms} − ED_{QQ_{1}P′} = ( \begin{matrix}1 & −η^{2}t_{n1} & 0 & η^{2}t_{n2} \\ 0 & 1 + η & 0 & −η \\ 0 & η^{2}t_{n1} & 1 & −η^{2}t_{n2} \\ 0 & −η & 0 & 1 + η \end{matrix} ) .

This latter matrix is singular iff η = −\frac{1}{2}. Moreover,

M = ( \begin{matrix}h & hηt_{n1} & 0 & 0 \\ 1 & ηt_{n1} + h(1 + η) & 1 & ηt_{n1} \\ 0 & 0 & h & hηt_{n2} \\ −1 & −ηt_{n2} & 1 & ηt_{n2} + h(1 + η) \end{matrix} ) .

and M is singular iff η = −\frac{1}{2}. □

Now we will study the stability and convergence of the IRK method (1.2). To that end we consider two cases: constant nullspace and variable nullspace.

3. Stability and convergence for the constant nullspace case

IRK methods become unstable when they are applied to DAEs with index 2. This fact is well-known in the simple case of linear DAEs with constant coefficients and constant stepsize. We refer to [1],

[8] for this case. For linear DAEs with variable coefficients but constant nullspace N(t) some terms of the system (2.3) are lost, because P′ = 0. Solving this system with respect to D_{P P_{1}} \overline{X}, D_{P Q_{1}} \overline{X} and D_{Q}\overline{X}, we obtain

D_{P P_{1}} \overline{X} = h((A ⊗ I_{m}) + O(h))D_{P P_{1}G_{2}^{−1}q}11_{s} − h[(A(C ⊙ A^{−1}) ⊗ I_{m})D_{P P′_{1}}
+ O(h)]D_{P Q_{1}G_{2}^{−1}q}11_{s} + (D_{P P_{1}} + O(h))D_{P} (11_{s} ⊗ x_{n−1}),  (3.1 a)
D_{P Q_{1}} \overline{X} = D_{P Q_{1}G_{2}^{−1}q}11_{s}  (3.1 b)
D_{Q}\overline{X} = [\frac{1}{h}(A^{−1} ⊗ I_{m})D_{QQ_{1}} + D_{QP_{1}} + ((C ⊙ A^{−1}) ⊗ I_{m})D_{QQ′_{1}P Q_{1}}
+ O(h)]D_{G_{2}^{−1}q}11_{s} + [\frac{1}{h}(A^{−1} ⊗ I_{m})D_{QP_{1}P} − D_{QP_{1}G_{2}^{−1}BP P_{1}}  (3.1 c)
+ ((C ⊙ A^{−1}) ⊗ I_{m})D_{QQ′_{1}}(D_{P P_{1}} − I_{ms}) + O(h)]D_{P} (11_{s} ⊗ x_{n−1}).

Clearly, the nontrivial term 1/h(A^{−1} ⊗ I_{m})D_{QP_{1}P} in (3.1 c) grows unboundedly if h tends to zero, i.e., the method remains unstable. Nevertheless we can prove convergence. To make this fact more transparent, we show that in the constant nullspace case the instability is weak and it refers only to the Q-component. To that purpose we consider, in addition to the vectors x_{n} obtained from the formulas (2.1), the x̃_{n} resulting from the perturbed systems

x̃_{n} = ρx̃_{n−1} + ∑_{j=1}^{s} ∑_{l=1}^{s} b_{j}â_{jl}X̃_{l} + hδ_{n,s+1},  (3.2 a)

\frac{1}{h}A(t_{nj}) ∑_{l=1}^{s} â_{jl}(X̃_{l} − x̃_{n−1}) + B(t_{nj})X̃_{j} = q(t_{nj}) + δ_{nj}; j = 1, 2, …, s.  (3.2 b)

Assume ∥δ_{ij}∥ ≤ δ for all i = 0, 1, …, n, j = 1, 2, …, s + 1.

Analogously to (2.2), we can write the system (3.2 b) in the following form

M\overline{X̃} = hD_{(q+δ_{n})}11_{s} + D_{A}(A^{−1}11_{s} ⊗ x_{n−1}),

where \overline{X̃} = (X̃_{1}^{T} , X̃_{2}^{T} , …, X̃_{s}^{T})^{T}.

According to theorem 2.1, M is nonsingular for sufficiently small h, i.e., the system (3.2 b) is uniquely solvable.

Let us introduce some new notations:

△x_{k} := x̃_{k} − x_{k}, △P x_{k} := P (x̃_{k} − x_{k}), △P P_{1}x_{k} := P P_{1}(t_{k})(x̃_{k} − x_{k}),

△P Q_{1}x_{k} := P Q_{1}(t_{k})(x̃_{k} − x_{k}), △Qx_{k} := Q(x̃_{k} − x_{k}),

△P P_{1}X_{j} := P P_{1}(t_{nj})(X̃_{j} − X_{j}), △P Q_{1}X_{j} := P Q_{1}(t_{nj})(X̃_{j} − X_{j}),

△QX_{j} := Q(X̃_{j} − X_{j}) ∈ ℝ^{⋗},

△\overline{X} := \overline{X̃} − \overline{X} ∈ ℝ^{⋗∼} for k = n − 1, n, j = 1, 2, …, s.

For the difference △\overline{X} we obtain (see 3.1):

D_{P P_{1}}△\overline{X}=h((A ⊗ I_{m}) + O(h))D_{P P_{1}G_{2}^{−1}δ_{n}}11_{s} − h[(A(C ⊙ A^{−1}) ⊗ I_{m})D_{P P′_{1}}
+O(h)]D_{P Q_{1}G_{2}^{−1}δ_{n}}11_{s} + (D_{P P_{1}} + O(h))D_{P} (11_{s} ⊗ △x_{n−1}),  (3.3 a)
D_{P Q_{1}}△\overline{X}=D_{P Q_{1}G_{2}^{−1}δ_{n}}11_{s} ,  (3.3 b)
D_{Q}△\overline{X} = [\frac{1}{h}(A^{−1} ⊗ I_{m})D_{QQ_{1}} + D_{QP_{1}} + ((C ⊙ A^{−1}) ⊗ I_{m})D_{QQ′_{1}P Q_{1}}
+O(h)]D_{G_{2}^{−1}δ_{n}}11_{s} + [\frac{1}{h}(A^{−1} ⊗ I_{m})D_{QP_{1}P} − D_{QP_{1}G_{2}^{−1}BP P_{1}}  (3.3 c)
+((C ⊙ A^{−1}) ⊗ I_{m})D_{QQ′_{1}}(D_{P P_{1}} − I_{ms}) + O(h)]D_{P} (11_{s} ⊗ △x_{n−1}).

Using Taylor’s expansion about t_{n−1}, and equating like powers of h, it is easy to see that the differences △X_{j}, j = 1, 2, …, s, satisfy:

∥△P P_{1}X_{j}∥ ≤ ∥△P P_{1}x_{n−1}∥ + hC(∥△P x_{n−1}∥ + δ),  (3.4 a)
∥△P Q_{1}X_{j}∥ ≤ Cδ,  (3.4 b)
∥△QX_{j}∥ ≤ \frac{C}{h} (∥△P Q_{1}x_{n−1}∥ + h∥△P x_{n−1}∥ + δ),  (3.4 c)

with C a constant.

These estimates reflect the influence of pertubations on the solutions of (2.1 b). To prove the main results of this seccion we will need the following lemma. For the proof we refer to [4], [8].

Lemma 3.1. Suppose that the two sequences of positive numbers (y_{n}), (z_{n}) satisfy

y_{n} ≤ y_{n−1} + hC(y_{n−1} + z_{n−1} + D_{1}),
z_{n} ≤ | ρ |z_{n−1} + hC(y_{n−1} + z_{n−1} + \frac{D_{2}}{h} ).

Then we have, for hn ≤ C_{0},

y_{n} ≤ C_{1}(y_{0} + hy_{0} + hz_{0} + D_{1} + D_{2}), if | ρ | ≤ 1,
z_{n} ≤ C_{1}(hy_{0} + (| ρ |^{n} + h)z_{0} + hD_{1} + D_{2}), if | ρ | < 1,
z_{n} ≤ C_{1}(hy_{0} + (1 + h)z_{0} + hD_{1} + \frac{D_{2}}{h} ), if | ρ | = 1.

Here C, C_{0}, C_{1}, D_{1}, D_{2} are positive constants.

Theorem 3.1 Suppose the assumptions of Theorem 2.1 are satisfied and furthermore, P′ = 0. If the coefficients of the IRK method (2.1) satisfy the relation |ρ| = |b^{T} A^{−1}11_{s} − 1| < 1, then the IRK method (2.1) applied to the DAE (1.1) is weakly unstable. This weak instability refers to the Q-component of the numerical solution. The P -component is stable. Proof. Subtracting (2.1a) from (3.2 a), we obtain

△x_{n} = ρ △x_{n−1} + (b^{T} A^{−1} ⊗ I_{m})△\overline{X} + hδ_{n,s+1}.  (3.5)

Multiplying this equation by P P_{1}P (t_{n}) = P P_{1}(t_{n}), P Q_{1}P (t_{n}) = P Q_{1}(t_{n}) and Q and using Taylor’s expansion about t_{n−1} on the right-hand side of the so obtained equations, we have

△P P_{1}x_{n} =ρ△P P_{1}x_{n−1} + hρ(P P′_{1}(t_{n−1})O(h))△P x_{n−1}
+ (b^{T} A^{−1} ⊗ I_{m})D_{P P_{1}}△\overline{X} + h((b^{T} A^{−1} ⊗ I_{m})D_{γP P′_{1}}
+ O(h))D_{P} △\overline{X} + hP P_{1}(t_{n})δ_{n,s+1},  (3.6 a)
△P Q_{1}x_{n} =ρ△P Q_{1}x_{n−1} + hρ(P Q′_{1}(t_{n−1}) + O(h))△P x_{n−1}
+ (b^{T} A^{−1} ⊗ I_{m})D_{P Q_{1}}△\overline{X}
+ h((b^{T} A^{−1} ⊗ I_{m})D_{γP Q′_{1}} + O(h))D_{P} △\overline{X}
+ hP Q_{1}(t_{n})δ_{n,s+1},  (3.6 b)
△Qx_{n} =ρ△Qx_{n−1} + (b^{T} A^{−1} ⊗ I_{m})D_{Q}△\overline{X} + hQδ_{n,s+1},  (3.6 c)

where γN(t_{nj}) = (1 − c_{j})N(t_{nj}) for any matrix function N.

Inserting (3.3 a) and (3.3 b) into (3.6 a) and (3.6 b) and using Taylor’s expansion about t_{n−1} on the right-hand side of the so obtained equations, we get

∥△P P_{1}x_{n}∥ ≤ ∥ρ△P P_{1}x_{n−1} + (b^{T} A^{−1} ⊗ I_{m})(11_{s} ⊗ △P P_{1}x_{n−1})∥
+ hC_{1}^{(n)}( ∥△P x_{n−1}∥ + δ)
= ∥(1 − b^{T} A^{−1}11_{s})△P P_{1}x_{n−1} + b^{T} A^{−1}11_{s}△P P_{1}x_{n−1}∥
+ hC_{1}^{(n)}( ∥△P x_{n−1}∥ + δ)
= ∥△P P_{1}x_{n−1}∥ + hC^{(n)}(∥△P x_{n−1}∥ + δ).  (3.7 a)
∥△P Q_{1}x_{n}∥ ≤ | ρ | ∥△P Q_{1}x_{n−1}∥ + ∥(b^{T} A^{−1} ⊗ I_{m})D_{P Q_{1}G_{2}^{−1}δ_{n}}11_{s}∥
+ hC_{1}^{(n)}(∥△P x_{n−1}∥ + δ)
≤ | ρ | ∥△P Q_{1}x_{n−1}∥ + hC^{(n)}(∥△P x_{n−1}∥ + \frac{δ}{h}),  (3.7 b)

where C^{(n)} and C_{1}^{(n)} are constants.

Applying lemma 3.1 to the recursions (3.7 a-b) with ∥△P P_{1}x_{k}∥, ∥△P Q_{1}x_{k}∥, δ, max_{i=1}^{n} C^{(i)} instead of y_{k}, z_{k}, D_{j}, C, respectively, k = n − 1, n, j = 1, 2, we obtain: If |ρ| ≤ 1, then

∥△P P_{1}x_{n}∥ ≤ C_{1}( ∥△P P_{1}x_{0}∥ + h∥△P P_{1}x_{0}∥ + h∥△P Q_{1}x_{0}∥ + δ).  (3.8 a)

If |ρ| < 1, then

∥△P Q_{1}x_{n}∥ ≤ C_{1}(h∥△P P_{1}x_{0}∥ + (| ρ |^{n} + h)∥△P Q_{1}x_{0}∥ + δ).  (3.8 b)

If |ρ| = 1, then

∥△P Q_{1}x_{n}∥ ≤ C_{1}(h∥△P P_{1}x_{0}∥ + (1 + h)∥△P Q_{1}x_{0}∥ + \frac{δ}{h}).  (3.8 c)

With x̃_{0} − x_{0} := δ_{0,s+1}, we arrive at

∥△P x_{n}∥ ≤ ∥△P P_{1}x_{n}∥ + ∥△P Q_{1}x_{n}∥ ≤ Cδ,  (3.9)

where C is constant and |ρ| < 1.

Now we consider the Q-component. Inserting (3.3 c) into (3.6 c) and using Taylor’s expansion about t_{n−1} on the right-hand side of the so obtained equation, we have

∥△Qx_{n}∥ ≤ | ρ | ∥△Qx_{n−1}∥ + \frac{C_{1}}{h} (∥△P Q_{1}x_{n−1}∥ + h∥△P x_{n−1}∥ + δ).

Inserting (3.9) into this equation, we get

∥△Qx_{n}∥ ≤ | ρ |∥△Qx_{n−1}∥ + \frac{C_{2}}{h} δ
≤ | ρ |^{n} ∥△Qx_{0}∥ + \frac{C_{2}}{h} δ ∑_{i=0}^{n−1} |ρ|^{i} ≤ \frac{C}{h} δ  (3.10)

with C constant. I.e., the IRK method is weak unstable on the Q-component.

□

Remarks

• If |ρ| = 1, then the P -component is also unstable because in this case we can only obtain an estimate of the form

∥△P x_{n}∥ ≤ \frac{C}{h} δ.

Furthermore, for the Q-component we obtain the estimate

∥△Qx_{n}∥ ≤ ∥△Qx_{0}∥ + \frac{C}{h^{2}} (h∥△P P_{1}x_{0}∥ + ∥△P Q_{1}x_{0}∥ + \frac{δ}{h}).

• A very important class of IRK methods, the superconvergent symmetric methods, may exhibit severe order reduction when applied to DAEs with index 2. For example, the socalled Gauss-Legendre methods, which have order 2s for nonstiff ODEs, exhibit very bad results in numerical experiments. These methods yield ρ = (−1)^{s}.

Now we study the convergence properties of the IRK method (2.1). We derive necessary and sufficient conditions for the local error of an IRK method when applied to the DAE (1.1). Then we study error propagation and derive estimates for the global error. Before we can state the main result of this subsection, we need some definitions.

Definition 3.1.

• The local error d_{n} of the method (2.1) is given by

d_{n} = ρx(t_{n−1}) + ∑_{j=1}^{s} b_{j} ∑_{l=1}^{s} â_{jl}X_{l} − x(t_{n}).

Here X_{l} is obtained from (2.1 b) with x_{n−1} = x(t_{n−1}).

• The j-th internal local truncation error δ_{j}^{(n)} at t_{n} of an s-stage IRK method (2.1) is given by

δ_{j}^{(n)} = x(t_{n−1}) + h ∑_{l=1}^{s} a_{jl}x′(t_{nl}) − x(t_{nj}),  (3.11 a)

δ_{s+1}^{(n)} = x(t_{n−1}) + h ∑_{l=1}^{s} b_{l}x′(t_{nl}) − x(t_{n}).  (3.11 b)

• The algebraic conditions B(p) and C(q) (Butcher identities) of an s-stage IRK method are given by

B(p) : ∑_{i=1}^{s} b_{i}c_{i}^{k−1} = \frac{1}{k} , for k = 1, 2, …, p,

C(q) : ∑_{j=1}^{s} a_{ij}c_{j}^{k−1} = \frac{c_{i}^{k}}{k} , for k = 1, 2, …, q and all i.

• The internal stage order K_{I} of an s-stage IRK method is given by

K_{I} = max{k : B(k) and C(k)}.

The numerical solution by the IRK method (2.1) satisfies (3.2 b) where δ_{nj}, j = 1, 2, …, s, represent errors in the solution of the system (2.1 b). The true solution x(t) of (1.1) satisfies

x(t_{n}) = ρx(t_{n−1}) + ∑_{j=1}^{s} ∑_{l=1}^{s} b_{j}â_{jl}(x(t_{nl}) + δ_{l}^{(n)}) − δ_{s+1}^{(n)}.  (3.12 a)

\frac{1}{h}A(t_{nj}) ∑_{l=1}^{s} â_{jl}(x(t_{nl}) − x(t_{n−1}) + δ_{j}^{(n)}) + B(t_{nj})x(t_{nj}) = q(t_{nj}).  (3.12 b)

where δ_{j}^{(n)}, j = 1, 2, …, s + 1, are the internal local truncation errors at t_{n}.

Subtracting (3.12) from (3.2), we obtain

e_{n} = ρe_{n−1} + (b^{T} A^{−1} ⊗ I_{m})(\overline{E} − δ̄^{(n)}) + δ_{s+1}^{(n)},  (3.13 a)
\frac{1}{h}A(t_{nj}) ∑_{l=1}^{s} â_{jl}E_{l} + B(t_{nj})E_{j} = \frac{1}{h}A(t_{nj}) ∑_{l=1}^{s} â_{jl}(δ_{l}^{(n)} + e_{n−1}), j = 1, 2, …, s, (3.13 b)

where e_{k} := x_{k} − x(t_{k}), k = n − 1, n, E_{j} := X_{j} − x(t_{nj}), j = 1, 2, …, s, \overline{E} =

(E_{1}^{T} , E_{2}^{T} , …, E_{s}^{T})^{T} ∈ ℝ^{⋗∼}, δ̄^{(n)} = ((δ_{1}^{(n)})^{T} , (δ_{2}^{(n)})^{T} , …, (δ_{s}^{(n)})^{T} )^{T} ∈ ℝ^{⋗∼},

ē_{n−1} = (e_{n−1}^{T}, e_{n−1}^{T}, …, e_{n−1}^{T})^{T} ∈ ℝ^{⋗∼} and δ̄_{n} = (δ_{n1}^{T}, δ_{n2}^{T}, …, δ_{ns}^{T})^{T} ∈ ℝ^{⋗∼}.

According to theorem 2.1, for sufficiently small h, and P′ = 0 the system (3.13 b) is uniquely solvable and the solution \overline{E} satisfies (see 3.1)

D_{P P_{1}} \overline{E} =(D_{P P_{1}} + O(h))D_{P} (δ̄^{(n)} + ē_{n−1}) + O(h∥δ̄_{n}∥),  (3.14 a)
D_{P Q_{1}} \overline{E} =D_{P Q_{1}G_{2}^{−1}} δ̄_{n},  (3.14 b)
D_{Q}\overline{E} =[\frac{1}{h}(A^{−1} ⊗ I_{m})D_{QP_{1}P} + ((C ⊙ A^{−1}) ⊗ I_{m})D_{QQ′_{1}}(D_{P P_{1}} − I_{ms})
− D_{QP_{1}G_{2}^{−1}BP P_{1}} + O(h)]D_{P} (δ̄^{(n)} + ē_{n−1}) + O(\frac{1}{h}∥δ̄_{n}∥).  (3.14 c)

Multiplying (3.13 a) by P P_{1}(t_{n}), P Q_{1}(t_{n}) and Q and using Taylor’s expansion about t_{n−1} on the right-hand side of the equations, we obtain

P P_{1}e_{n} =ρP P_{1}e_{n−1} + hρ(P P′_{1}(t_{n−1}) + O(h))P e_{n−1}
(b^{T} A^{−1} ⊗ I_{m})D_{P P_{1}}(\overline{E} − δ̄^{(n)})
+ h[((b^{T} A^{−1} ⊗ I_{m})D_{γP P′_{1}} + O(h))]D_{P} (\overline{E} − δ̄^{(n)})
+ P P_{1}(t_{n})δ_{s+1}^{(n)},  (3.15 a)
P Q_{1}e_{n} =ρP Q_{1}e_{n−1} + hρ(P Q′_{1}(t_{n−1}) + O(h))P e_{n−1}
+ (b^{T} A^{−1} ⊗ I_{m})D_{P Q_{1}} \overline{E} + h[((b^{T} A^{−1} ⊗ I_{m})D_{γP Q′_{1}} + O(h))]D_{P} \overline{E}
− P Q_{1}(t_{n})[(b^{T} A^{−1} ⊗ I_{m})δ^{(n)} − δ_{s+1}^{(n)}],  (3.15 b)
Qe_{n} =ρQe_{n−1} + (b^{T} A^{−1} ⊗ I_{m})D_{Q}(\overline{E} − δ̄^{(n)}) + Qδ_{s+1}^{(n)} ,  (3.15 c)

where P P_{1}e_{k} := P P_{1}(t_{k})e_{k}, P Q_{1}e_{k} := P Q_{1}(t_{k})e_{k}, k = n − 1, n. Inserting (3.14) into (3.15) we thus arrive at

P P_{1}e_{n} = P P_{1}e_{n−1} + O(h(∥P e_{n−1}∥ + ∥δ̄^{(n)}∥ + ∥δ̄_{n}∥)) + P P_{1}(t_{n})δ_{s+1}^{(n)}  (3.16 a)

P Q_{1}e_{n} =ρP Q_{1}e_{n−1} + O(h∥P e_{n−1}∥ + h∥δ̄^{(n)}∥ + ∥δ̄_{n}∥)
− P Q_{1}(t_{n})[(b^{T} A^{−1} ⊗ I_{m})δ̄^{(n)} − δ_{s+1}^{(n)}],  (3.16 b)

Qe_{n} =ρQe_{n−1} + (b^{T} A^{−1} ⊗ I_{m})[\frac{1}{h}(A^{−1} ⊗ I_{m})D_{QP_{1}P} − D_{QP_{1}G_{2}^{−1}BP P_{1}}
+ ((C ⊙ A^{−1}) ⊗ I_{m})D_{QQ′_{1}}(D_{P P_{1}} − I_{m}) + O(h)]D_{P} (δ̄^{(n)} + ē_{n−1})
− (b^{T} A^{−1} ⊗ I_{m})δ̄^{(n)} + Qδ_{s+1}^{(n)} + O(\frac{1}{h}∥δ̄_{n}∥).  (3.16 c)

The local error of the IRK method (2.1) is given by (3.16) for x_{n−1} = x(t_{n−1}), i.e., e_{n−1} = 0.

Theorem 3.2. Let the assumptions of Theorem 2.1 be satisfied and furthermore assume that, P′ = 0. Suppose the IRK method satisfies the Butcher identities B(p) and C(q) for p ≥ q ≥ 1 and the numerical errors in the solution of the system (2.1 b) are δ_{nj} = O(h^{q+α}). Then the local error satisfies:

(a) P d_{n} = O(h^{q+1}), Qd_{n} = O(h^{q}), if α = 1.

(b) P P_{1}(t_{n})d_{n} = O(h^{q+2}), for p ≥ q + 1 and α = 1.

(c) P Q_{1}(t_{n})d_{n} = O(h^{q+2}), if α = 2 and the coefficients of the IRK method satisfy

b^{T} A^{−1}c^{q+1} = 1

where c^{q+1} = (c_{1}^{q+1}, c_{2}^{q+1}, …, c_{s}^{q+1})^{T}.

(d) P d_{n} = O(h^{q+2}), if α = 2, b^{T} A^{−1}c^{q+1} = 1 and p ≥ q + 1.

(e) Qd_{n} = O(h^{q+1}), if α = 2 and the coefficients of the IRK method satisfy

b^{T} A^{−2}c^{q+1} = q + 1.

Proof.

(a) It follows directly from (3.16) with e_{n−1} = 0, δ̄^{(n)} = O(h^{q+1}), δ_{s+1}^{(n)} = O(h^{q+1}), δ̄_{n} = O(h^{q+1}) and P e_{n} = P P_{1}e_{n} + P Q_{1}e_{n}.

(b) It follows from (3.16 a) with e_{n−1} = 0, δ̄^{(n)} = O(h^{q+1}), δ̄_{n} = O(h^{q+1}) and δ_{s+1}^{(n)} = O(h^{q+2}) (because B(q + 1)).

(c) Inserting e_{n−1} = 0, δ̄^{(n)} = O(h^{q+1}) and δ̄_{n} = O(h^{q+2}) into (3.16 b), we obtain

P Q_{1}(t_{n})d_{n} = −P Q_{1}(t_{n})[(b^{T} A^{−1} ⊗ I_{m})δ̄^{(n)} − δ_{s+1}^{(n)}] + O(h^{q+2}).  (3.17)

For q ≤ p, the following relations

δ̄^{(n)} = h^{q+1}[(Ac^{q} ⊗ x^{(q+1)}(t_{n−1})) − \frac{1}{q + 1}(c^{q+1} ⊗ x^{(q+1)}(t_{n−1}))]
+ O(h^{q+2})  (3.18 a)

and

δ_{s+1}^{n} = h^{q+1}(b^{T} c^{q} − \frac{1}{q + 1})x^{(q+1)}(t_{n−1}) + O(h^{q+2}).  (3.18 b)

hold. This yields

(b^{T} A^{−1} ⊗ I_{m})δ̄^{(n)} − δ_{s+1}^{(n)} =

h^{q+1}(b^{T} A^{−1} ⊗ I_{m})[(Ac^{q}) ⊗ x^{(q+1)}(t_{n−1}) − \frac{1}{q + 1}(c^{q+1} ⊗ x^{(q+1)}(t_{n−1}))]

− h^{q+1}(b^{T} c^{q} − \frac{1}{q + 1})x^{(q+1)}(t_{n−1}) + O(h^{q+2})

= \frac{h^{q+1}}{q + 1}(1 − b^{T} A^{−1}c^{q+1})x^{(q+1)}(t_{n−1}) + O(h^{q+2}).

Inserting this equation in (3.17) and using b^{T} A^{−1}c^{q+1} = 1 we obtain

P Q_{1}(t_{n})d_{n} = O(h^{q+2}).

(d) Follows directly from (b) and (c).

(e) Inserting (3.18 a) into (3.16 c) and using δ̄^{(n)} = O(h^{q+1}), δ_{s+1}^{(n)} = O(h^{q+1}), δ̄_{n} = O(h^{q+2}) and e_{n−1} = 0 we obtain

Qd_{n} = h^{q}QP_{1}P (t_{n})((b^{T} A^{−1}c^{q} − \frac{1}{q + 1}b^{T} A^{−2}c^{q+1}) ⊗ x^{(q+1)}(t_{n−1})) + O(h^{q+1}).

Inserting the identity b^{T} A^{−2}c^{q+1} = q + 1 into this equation and using B(q) and C(q) we arrive at

Qd_{n} = O(h^{q+1}) . □

Theorem 3.2 gives conditions on the coefficients of an IRK method so that it attains a given order for the local error when applied to the DAE (1.1) with constant nullspace. The conditions b^{T} A^{−1}c^{q+1} = 1 and b^{T} A^{−2}c^{q+1} = q + 1 have been studied in [1] when IRK methods are applied to linear constant coefficient systems of index 1 and 2. Furthermore, in [1] Brennan, Campbell and Petzold derive conditions for the local error of an IRK method applied to linear constant coefficient systems of arbitrary index k.

Now we use the recursions (3.16) to derive estimates for the global error.

Theorem 3.3. Let the assumptions of theorem 2.1 be satisfied and furthermore assume that P′ = 0. Suppose the IRK method satisfies |ρ| < 1. For q ≥ 1, let e_{0} = O(h^{q+α}) be the errors in the initial conditions and δ_{ij} = O(h^{q+α}) the numerical errors in the solution of the system (2.1 b) at the i-th integration step. Then the global error e_{n} satisfies:

• ∥P e_{n}∥ ≤ Ch^{q}, if B(q), C(q) and α = 0,

• ∥P e_{n}∥ ≤ Ch^{q+1}, if B(q + 1), C(q) and α = 1,

• ∥Qe_{n}∥ ≤ Ch^{q}, if B(q), C(q) and α = 1,
• ∥Qe_{n}∥ ≤ Ch^{q+1}, if B(q + 1), C(q), b^{T} A^{−1}c^{q+1} = 1, b^{T} A^{−2}c^{q+1} = q + 1 and α = 2.

Proof. From (3.16 a-b) we immediatly obtain the estimates

∥P P_{1}e_{n}∥ ≤ ∥P P_{1}e_{n−1}∥ + hC^{(n)}(∥P P_{1}e_{n−1}∥ + ∥P Q_{1}e_{n−1}∥ + D_{1}), (3.19 a)

∥P Q_{1}e_{n}∥ ≤ | ρ | ∥P Q_{1}e_{n−1}∥ + hC^{(n)}(∥P P_{1}e_{n−1}∥ + ∥P Q_{1}e_{n−1}∥ + \frac{D_{2}}{h} ) (3.19 b)

with C^{(n)} constants and

D_{1} = \{ \begin{matrix}O(h^{q}), if α = 0, B(q) and C(q), \\ O(h^{q+1}), if α = 1, B(q + 1) and C(q), \end{matrix}  (3.20 a)

D_{2} = \{ \begin{matrix}O(h^{q}), if α = 0, B(q) and C(q), \\ O(h^{q+1}), if α = 1, B(q) and C(q), \\ O(h^{q+2}), if α = 2, b^{T} A^{−1}c^{q+1} = 1, B(q) and C(q). \end{matrix}  (3.20 b)

Applying Lemma 3.1 to the recursions (3.19) we obtain

∥P P_{1}e_{n}∥ ≤ C_{1}(∥P P_{1}e_{0}∥ + h∥P P_{1}e_{0}∥ + h∥P Q_{1}e_{0}∥ + D_{1} + D_{2}),  (3.21 a)
∥P Q_{1}e_{n}∥ ≤ C_{1}(h∥P P_{1}e_{0}∥ + (|ρ|^{n} + h)∥P Q_{1}e_{0}∥ + hD_{1} + D_{2}).  (3.21 b)

Now, the estimates for P e_{n} follow from (3.20), (3.21) and the inequality

∥P e_{n}∥ ≤ ∥P P_{1}e_{n}∥ + ∥P Q_{1}e_{n}∥.

In order to prove the first estimate for Qe_{n}, we consider (3.16 c) and (3.21 b). This last inequality together with (3.20) yield

P Q_{1}e_{n} = O(h^{q+1}), if B(q), C(q) and α = 1.

Inserting this last equation together with P e_{n} = O(h^{q}), δ̄_{n} = O(h^{q+1}), δ̄^{(n)} = O(h^{q+1}) and δ_{s+1}^{(n)} = O(h^{q+1}) into (3.16 c), we obtain

Qe_{n} = ρQe_{n−1} + O(h^{q}).

The first estimate for Qe_{n} follows directly from this recursion. In order to prove the second estimate for Qe_{n}, we insert (3.20) with α = 2, B(q + 1), C(q) and b^{T} A^{−1}c^{q+1} = 1 into (3.21 b). This implies

P Q_{1}e_{n} = O(h^{q+2}).

Inserting this equation together with P e_{n} = O(h^{q+1}), δ̄_{n} = O(h^{q+2}), δ_{s+1}^{(n)} = O(h^{q+2}) and δ̄^{(n)} = O(h^{q+1}) into (3.16 c) we obtain

Qe_{n} = ρQe_{n−1} + \frac{1}{h}(b^{T} A^{−2} ⊗ I_{m})D_{QP_{1}P} δ̄^{(n)} + O(h^{q+1}).

Using the identity b^{T} A^{−2}c^{q+1} = q + 1 we thus arrive at

Qe_{n} = ρe_{n−1} + O(h^{q+1}).

The second estimate for Qe_{n} follows directly from this recursion. □

For many important methods such as the Gauss-Legendre method, |ρ| = 1. This means that the statements of theorem 3.3 are not valid for this class of IRK. The following theorem shows how strong the order reduction is when IRK methods with |ρ| = 1 are applied to linear DAEs with index 2 and constant nullspace.

Theorem 3.4. Assume that the hypothesis of Theorem 2.1 are satisfied. Suppose that the IRK method satisfies |ρ| = 1, and further assume that, for q ≥ 1, e_{0} = O(h^{q+β}) is the error in the initial

16 EBROUL IZQUIERDO

condition and δ_{ij} = O(h^{q+α}) is the numerical error in the solution of the system (2.1 b) at the ith integration step. Then the global error e_{n} satisfies:

• ∥P e_{n}∥ ≤ Ch^{q}, ∥Qe