Zeitschrift für Analysis und ihre Anwendungen

Journal for Analysis and its Applications Volume 18 (1999), No. 2, 393–406

The Discretization of

Neutral Functional Integro-Differential Equations by Collocation Methods

H. Brunner

Dedicated to Prof. L. von Wolfersdorf

Abstract. We study the approximate solution of certain neutral functional integro-differential equations describing the aeroelastic motions of certain airfoils, and of related Volterra equations with a delay argument. The focus of the paper is on questions concerning the (uniform) convergence and local superconvergence properties of discretization methods based on collocation techniques for equations with smooth kernels.

Keywords: Neutral functional integro-differential equations, Volterra integral equations with

delay, collocation methods

AMS subject classification: 45 E 10, 34 K 40, 65 R 20

1. Introduction

In mathematical modelling processes describing physical phenomena with memory effects, Volterra integral equations of the first kind often arise in non-standard form. Typical examples may be found in [13, 14, 18, 20] (inverse problems for the identification of memory kernels in heat conduction and viscoelasticity), [19] (qualitative behavior of semiconductor devices with abrupt pn-junctions), and [8, 9] (aeroelastic motion of airfoils with flap). In the following I shall focus on this last problem: it includes a neutral functional integro-differential equation which is equivalent to a Volterra integral equation of the first kind with a delay argument. We shall look at this problem from the viewpoint of a numerical analyst: why and how does the discretization of the integral equations differ from that used for ”standard” (first-kind) Volterra integral equations, and what are feasible methods (yielding high-order approximations) for their numerical solution? As we shall see in more detail below, the construction and analysis of high-order numerical methods for such equations is not trivial (in contrast to second-kind Volterra integral equations and integro-differential equations [5, 6]). A thorough understanding of this problem is crucial when solving, for example, ”mixed” systems consisting of Volterra integral equations of the second and first kind arising, for example, in the work of von Wolfersdorf (and his collaborators) cited above.

H. Brunner: Memorial University of Newfoundland, Dept. Math. & Stat., St. John’s, Newfoundland, Canada A1C 5S7; e-mail: hermann@math.mun.ca ISSN 0232-2064 / $ 2.50 \begin{matrix}c \\ ◯ \end{matrix} Heldermann Verlag Berlin

For ease of exposition I will focus on some appropriate generic form of the problem under consideration; the reader interested in detailed background information on the sources and the theory of the underlying general problem will find a list of appropriate references at the end of the paper.

The elastic motions of a 3-degree-of-freedom airfoil section with flap in a 2-dimen-sional incompressible flow can be described by a system of neutral functional integrodifferential equations of the form

\begin{matrix}\frac{d}{dt} & (A_{0}x(t) + ∫_{−r}^{0} A_{1}(s)x(t + s) ds) \\ & = B_{0}x(t) + B_{1}x(t − r) + ∫_{−r}^{0} K(s)x(t + s) ds + F (t) \end{matrix}  (t > 0)  (1.1)

with

x(s) = ϕ(s) (−r ≤ s ≤ 0; 0 < r = const).

Here, the matrices A_{0}, A_{1}(s) = [a_{ij}^{(1)}(s)], B_{0}, B_{1}, K(s) ∈ ℝ^{d×d} (d = 8) are given. The element a_{88}^{(1)} of A_{1}(s) has the form (see [9, 12] and their references)

a_{88}^{(1)}(s) = c(s)(−s)^{−α} + p(s)  (0 < α < 1);  (1.2)

c and p are smooth functions and, typically, α = \frac{1}{2}. The particular structure of these matrices, especially of A_{0} (whose last row consists of zeros), allows the decoupling of the weakly singular integral component; see [8, 9, 12] for details and further references. This integral component may be written generically as

\frac{d}{dt} (∫_{−r}^{0}(−s)^{−α}x(t + s) ds)  (1.3)

= b_{0}x(t) + b_{1}x(t − r) + ∫_{−r}^{0} k(s)x(t + s) ds + f(t) (t > 0)  (1.4)

where the kernel k is smooth. In this paper we shall consider an important special case of (1.11), namely

\frac{d}{dt} (∫_{−r}^{0}(−s)^{−α}x(t + s) ds) = f(t)  (t > 0; 0 < α < 1)  (1.5)

subject to the initial condition

x(s) = ϕ(s) (−r ≤ s ≤ 0).

Setting

D_{α}ϕ := ∫_{−r}^{0}(−s)^{−α}ϕ(s) ds  and  g(t) := D_{α}ϕ + ∫_{0}^{t} f(s) ds,

this initial-value problem may be written as

∫_{−r}^{0}(−s)^{−α}x(t + s) ds = g(t) (t > 0)

or, equivalently, as

\begin{matrix}∫_{t−r}^{t}(t − s)^{−α}x(s) ds & = g(t)  (t > 0) \\ x(s) & = ϕ(s)  (−r ≤ s ≤ 0). \end{matrix} \}  (1.6)

The above equation is a first-kind Volterra integral equation with weakly singular kernel and constant delay r > 0.

It is well known (see, for example, [4, 5, 6, 16]) that the discretization of (non-delay) Volterra integral equations of the first kind with regular or weakly singular kernel by quadrature or collocation methods does in general not yield a (uniformly) convergent numerical method. We shall discuss the analogous convergence problems encountered in the discretization of (1.5) and (1.6) in Sections 3 and 4.

2. Discretization based on collocation

We begin by describing the framework for the piecewise polynomial collocation method which underlies the discretization for equation (1.6) and its generalization (2.4) below. Let I := [a, b] be the (compact) interval on which a given Volterra integral equation is to be solved, assume that Π_{N} denotes a mesh for I,

Π_{N} : a = t_{0} < t_{1} < … < t_{N} = b,

and set h_{n} := t_{n+1} − t_{n} (n = 0, 1, …, N) and h := max_{(n)}{h_{n}}. For our purposes the mesh points {t_{n}} will be given by

t_{n} = a + (\frac{n}{N} )^{q} · (b − a)  (n = 0, 1, …, N),  (2.1)

where the grading exponent q satisfies q ≥ 1; it will depend on the degree of regularity of the solution at t = a.

An approximation u to the solution of the given Volterra integral equation will be sought in the linear space S_{m−1}^{(−1)}(Π_{N}) of discontinuous (real) piecewise polynomials of degree not exceeding m − 1 ≥ 0,

S_{m−1}^{(−1)}(Π_{N}) := {u : u|_{(t_{n},t_{n+1}]} ∈ π_{m−1} (n = 0, 1, …, N − 1)},

whose dimension is Nm. This collocation solution u ∈ S_{m−1}^{(−1)}(Π_{N}) is to satisfy the integral equation on the set

X_{N} := {t_{n} + c_{i}h_{n} : 0 < c_{1} < … < c_{m} ≤ 1 (n = 0, 1, …, N − 1)} ⊂ I  (2.2)

of collocation points. This set is completely determined by the given mesh and the collocation parameters {c_{i}}.

A suitable, and commonly used, local representation of the collocation solution u ∈ S_{m−1}^{(−1)}(Π_{N}) is given by

u(t_{n} + sh_{n}) = ∑_{j=1}^{m} L_{j}(s)U_{n,j}  (s ∈ (0, 1]),  (2.3)

where U_{n,j} := u(t_{n} + c_{j}h_{n}); the L_{j} denote the Lagrange fundamental polynomials with respect to the set {c_{j}}. Detailed background information on the discretization of various types of differential and integral equations by spline collocation methods may be found for example in [5, 6].

We now turn to the first-kind Volterra integral equation with constant delay r > 0,

\begin{matrix}∫_{t−r}^{t} k(t − s)x(s) ds & = g(t)  (t ∈ (0, T]) \\ x(s) & = ϕ(s)  (−r ≤ s ≤ 0) \end{matrix} \}  (2.4)

where the kernel k satisfies either k ∈ C^{ν}(I) for some integer ν ≥ 1, with k(0) > 0, or is of the form k(t) = t^{−α} (0 < α < 1). It is easily seen that, due to the delay argument t−r, solutions to (2.4) will in general be discontinuous at the points ξ_{μ} := μr (μ ∈ ℕ_{0}), in analogy to those of delay differential equations with constant delay ([10, 21]; see also [2]). We shall refer to these points as primary discontinuities for (2.4) (in contrast to secondary discontinuities generated by discontinuities in the given initial function ϕ); we assume without loss of generality that T = ξ_{M+1} for some positive integer M.

Thus, the role of the interval [a, b] in the above description of the collocation method will be assumed successively by the intervals [ξ_{μ}, ξ_{μ+1}] (μ = 0, 1, …, M), while that of S_{m−1}^{(−1)}(Π_{N}) is taken by the polynomial spline spaces S_{m−1}^{(−1)}(Π_{N}^{(μ)}) corresponding to the meshes

Π_{N}^{(μ)} : 0 = t_{0}^{(μ)} < t_{1}^{(μ)} < … < t_{N}^{(μ)} = ξ_{μ+1} (0 ≤ μ ≤ M)

the mesh for the subintervals [ξ_{μ}, ξ_{μ+1}], with h_{n}^{(μ)} := t_{n+1}^{(μ)} − t_{n}^{(μ)}. The choice of these meshes obviously corresponds to the use of so-called constrained meshes in the numerical solution of delay differential equations with constant delay r > 0 (see, for example, [21]).

On (ξ_{μ}, ξ_{μ+1}] the collocation solution u = u_{μ} ∈ S_{m−1}^{(−1)}(Π_{N}^{(μ)}) is determined from the collocation equation

∫_{t−r}^{t} k(t − s)u(s) ds = g(t)  (t ∈ X_{N}^{(μ)}; μ = 0, 1, …, M),  (2.5)

with values u(s) = ϕ(s) if s ∈ [−r, 0]. Here, in analogy to (2.2),

X_{N}^{(μ)} := {t_{n}^{(μ)} + c_{i}h_{n}^{(μ)} : 0 < c_{1} < … < c_{m} ≤ 1 (n = 0, 1, …, N − 1)}.

Equation (2.5) may be written as

∫_{ξ_{μ}}^{t} k(t − s)u(s) ds = G_{μ}(t)  (t ∈ X_{N}^{(μ)})  (2.6)

with

G_{μ}(t) := g(t) − ∫_{t−r}^{ξ_{μ}} k(t − s)u(s) ds.  (2.7)

The collocation method for (2.4) is then described by equations (2.8) and (2.9) below: for given μ = 0, 1, …, M, let

u_{μ}(t_{n}^{(μ)} + sh_{n}^{(μ)}) = ∑_{j=1}^{m} L_{j}(s)U_{n,j}^{(μ)}  (s ∈ (0, 1]),  (2.8)

with U_{n,j}^{(μ)} := u_{μ}(t_{n}^{(μ)} + c_{j}h_{n}^{(μ)}), be the local representation of u = u_{μ} ∈ S_{m−1}^{(−1)}(Π_{N}^{(μ)}) on the subinterval (t_{n}^{μ}, t_{n+1}^{(μ)} ∈ (ξ_{μ}, ξ_{μ+1}] (cf. (2.3)). The values U_{n}^{(μ)} := (U_{n,1}^{(μ)}, …, U_{n,m}^{(μ)})^{T} ∈ ℝ^{m} are given by the solution of the linear algebraic system

\begin{aligned}h_{n}^{(μ)} ∑_{j=1}^{m} ∫_{0}^{c_{i}} k(h_{n}^{(μ)}(c_{i} − s))l_{j}(s) ds & = G_{μ}(t_{n}^{(μ)} + c_{i}h_{n}^{(μ)}) − Φ_{n,i}^{(μ)} \\ (i = 1, …, m; n & = 0, 1, …, N − 1) \end{aligned}  (2.9)

where

Φ_{n,i}^{(μ)} := ∑_{ℓ=0}^{n−1} h_{ℓ}^{(μ)} ∫_{0}^{1} k(t_{n}^{(μ)} + c_{i}h_{n}^{(μ)} − t_{ℓ}^{(μ)} − sh_{ℓ}^{(μ)})u(t_{ℓ}^{(μ)} + sh_{ℓ}^{(μ)})ds  (2.10)

represents the values of the ”history term” for [ξ_{μ}, ξ_{μ+1}],

Φ_{n}^{(}μ)(t) := ∫_{ξ_{μ}}^{t_{n}^{(μ)}} k(t − s)u(s) ds

at the collocation points t = t_{n}^{(μ)} + c_{i}h_{n}^{(μ)}. On the initial interval (ξ_{0}, ξ_{1}] we have

G_{0}(t) = g(t) − ∫_{t−r}^{0} k(t − s)ϕ(s) ds.  (2.11)

It is readily verified that for kernels k = k(t) satisfying either k ∈ C^{1}(I), k(0) > 0, or k(t) = t^{−α} (0 < α < 1), the linear system (2.9) has a unique solution U_{n}^{(μ)} ∈ ℝ^{m} for each n = 0, 1, …, N − 1 and each μ = 0, 1, …, M.

3. Uniform convergence of collocation solutions

In this section we shall analyze the uniform convergence of collocation solutions u ∈ S_{m−1}^{(−1)}(Π_{N}^{(μ)}) for delay Volterra integral equations (2.4) with smooth k: k ∈ C^{m+1}(I), with k(0) > 0. Analogous results corresponding to equations with weakly singular kernels will be given elsewhere (but compare also [4] for related open problems).

3.1 Regularity of solutions. The first-kind Volterra integral equation (2.6) has a solution that is continuous at the primary discontinuity points t = ξ_{μ} = μr only if

G_{μ}(ξ_{μ}) = g(ξ_{μ}) − ∫_{ξ_{μ−1}}^{ξ_{μ}} k(ξ_{μ} − s)x(s) ds = 0

(compare (2.7)). In particular, for μ = 0 this necessary condition reduces to

g(0) − ∫_{−r}^{0} k(−s)ϕ(s) ds =: η − D_{0}ϕ = 0.  (3.1)

The proof of the regularity result in Lemma 3.1 is straightforward: it is based on classical Volterra theory (see also [4, 6]). A similar continuity result for the solution of (2.4) with weakly singular kernel k(t) = t^{−α} (0 < α < 1) is given in [15]; its regularity properties will be discussed elsewhere.

Lemma 3.1. Assume that k ∈ C^{d+1}(I), with k(0) > 0, g ∈ C^{d+1}(I), and ϕ ∈ C^{d}[−r, 0], where d ≥ 0. Then:

(i) The (unique) solution x = x(t) of (2.6) is continuous on I = [0, T] if and only if (3.1) holds. Moreover, for d ≥ 1 we have

x ∈ C^{d}[ξ_{μ}, ξ_{μ+1}] (μ = 0, 1, …, M; ξ_{M+1} = T)

but in general x ∉ C^{1}(I).

(ii) Whenever η−D_{0}ϕ ≠ 0, x is discontinuous at t = ξ_{μ} (μ = 0, 1, …, M); however, we still have x ∈ C^{m}(ξ_{μ}, ξ_{μ+1}].

3.2 Convergence results. In the following analysis we will assume that the meshes Π_{N}^{(μ)} are uniform, that is, h^{(μ)} = h = \frac{r}{N} for all μ = 0, 1, …, M.

Theorem 3.2. Assume that k ∈ C^{d+1}(I), with k(0) > 0, g ∈ C^{d+1}(I), and ϕ ∈ C^{d}[−r, 0] for some d ≥ 0. Let x be the (unique) solution to (2.4) and denote by u ∈ S_{m−1}^{(−1)}(Π_{N}^{(μ)}) (μ = 0, 1, …, M) the collocation solution for the first-kind Volterra integral equation with delay r > 0, (2.4). Suppose that (3.1) holds. Then

||x − u||_{∞} := sup_{t∈I} |x(t) − u(t)| → 0  as N → ∞ (Nh^{(μ)} = Nh = r)  (3.2)

(with Nh^{(μ)} = const) if and only if the collocation parameters {c_{i}} characterizing the collocation points X_{N}^{(μ)} satisfy

−1 ≤ ρ_{m} := (−1)^{m+1} ∏_{i=1}^{m} \frac{1 − c_{i}}{c_{i}} ≤ 1.  (3.3)

If d = m and (3.4) holds, the order of convergence is given by

||x − u||_{∞} ≤ \{ \begin{matrix}Ch^{m} & if ρ_{m} ∈ [−1, 1] \\ Ch^{m−1} & if ρ_{m} = 1. \end{matrix}

Here (and in Section 4) C denotes a generic (finite) constant independent of h.

Remark. An analogous result holds if the collocation space is S_{m}^{(0)}(Π_{N}^{(μ)}), that is, the space of continuous piecewise polynomials of degree not exceeding m ≥ 1, provided we have 0 < c_{1} < … < c_{m} = 1. However, in this case condition (3.2) has to be replaced by

∏_{i=1}^{m−1} \frac{1 − c_{i}}{c_{i}} ≤ 1.

If the collocation parameters are the m Gauss points in (0, 1) (i.e. the zeros of the shifted Legendre polynomial P_{m}(2s−1)) where c_{m} < 1 but for which (3.2) is true (with = replacing ≤), then the corresponding collocation solution in S_{m}^{(0)}(Π_{N}^{(μ)}) is divergent. This can easily be verified by adapting the results in [16] for classical (non-delay) Volterra integral equations of the first kind to the present situation.

Proof of Theorem 3.2. We shall sketch the proof for the case when d = m; the corresponding analysis is readily modified to encompass data with less regularity (i.e. 1 ≤ d ≤ m − 1). Consider (2.6) on the interval I_{μ} := [ξ_{μ}, ξ_{μ+1}] and observe that, by Lemma 3.1/(i), its (unique) solution has the regularity x ∈ C^{m}(I_{μ}). The collocation error e := x − u solves

∫_{ξ_{μ}}^{t} k(t − s)e(s) ds = − ∫_{t−r}^{ξ_{μ}} k(t − s)e(s) ds for t ∈ X_{N}^{(μ)} (μ = 0, 1, …, M)

where e(s) = 0 for s ∈ [−r, 0]. Thus, on I_{0} = [0, r] the classical convergence theory for collocation in S_{m−1}^{(−1)}(Π_{N}^{(0)}) for non-delay Volterra integral equations of the first kind (cf. [1, 6]) shows immediately that (3.3) holds for any m ≥ 0 if and only if (3.4) is satisfied. (An analogous result for collocation in S_{m}^{(0)}(Π_{N}^{(0)} can be derived by using the techniques of [16].)

Let now 1 ≤ μ ≤ M. For ease of exposition, and without any loss of essential features of the problem, we will assume that k(t) ≡ 1 and h_{n}^{(μ)} = h = \frac{r}{N} for all n and μ. The above error equation can then be written as

∫_{ξ_{μ}}^{t_{n}^{(μ)}+c_{i}h} e(s) ds = − ∫_{t_{n}^{(μ−1)}+c_{i}h}^{ξ_{μ}} e(s) ds (i = 1, …, m).

If n ≥ 1, replace n by n − 1 and c_{i} by c_{m}, and subtract the resulting error equation from the one above, to obtain (after an obvious change of variables)

∫_{0}^{c_{i}} e_{μ,n}(s) ds = − ∫_{c_{m}}^{1} e_{μ,n−1}(s) ds + ∫_{c_{m}}^{1} e_{μ−1,n−1}(s) ds + ∫_{0}^{c_{i}} e_{μ−1,n−1}(s) ds (3.4)

where we have set e_{μ,n} := e|_{(t_{n}^{(μ)},t_{n+1}^{(μ)}]}. On (t_{n}^{(μ)}, t_{n+1}^{(μ)}], the regularity of x and Taylor’s formula allow us to write

e_{μ,n}(s) := e(t_{n}^{(μ)} + sh) = ∑_{j=1}^{m} β_{n,j}^{(μ)}s^{j−1} + h^{m}R_{μ,n}(s)  (s ∈ (0, 1])  (3.5)

where R_{μ,n}(s) := \frac{s^{m}}{m!} \frac{dx^{m}}{dt^{m}} z_{μ,n} t_{n}^{(μ)} < z_{μ,n} < t_{n+1}^{(μ)}.

Let the m × m matrices P and Q be defined by

P := (∫_{0}^{c_{i}} s^{j−1}ds)  and  Q := (∫_{c_{m}}^{1} s^{j−1}ds)

and set β_{n}^{(μ)} := (β_{n,1}^{(μ)}, …, β_{n,m}^{(μ)})^{T} ∈ ℝ^{m}. Substitution of (3.6) (and the expressions corresponding to other indices) into (3.5) leads to

P β_{n}^{(μ)} = −Qβ_{n−1}^{(μ)} + [P β_{n}^{(μ−1)} + Qβ_{n−1}^{(μ−1)}] + O(h^{m})  (3.6)

(n = 1, …, M; μ = 1, …, M). The matrix P is non-singular (it is essentially a Vandermonde matrix), and Q has rank one. The non-trivial eigenvalue λ_{m} of P^{−1}Q is found to be

λ_{m} = (−1)^{m+1} ∏_{i=1}^{m} \frac{1 − c_{i}}{c_{i}} .

It thus follows from the elementary theory of difference equations (see, for example, [11: Section 3.3]) or [17: Section 4.2]) that (3.6) has bounded solutions (uniformly in n, as N → ∞ and Nh = r for all μ = 0, 1, …, M if and only if ρ(P^{−1}Q) = |λ_{m}| ≤ 1.

The order reduction from m to m − 1 when λ_{m} = 1 (implying that m must be odd, and hence m ≥ 1) is due to the fact that the expression for the general solution of (3.6) now contains a term of the form N · O(h^{m}), where Nh = r

Remark. If the continuity condition (3.1) does not hold for the given data k, g and ϕ in (2.4), then (3.3) and the order result will be valid only on compact subintervals of the form [ξ_{μ} + ε, ξ_{μ+1}] (μ = 0, 1, …, M) where ε > 0.

4. Local superconvergence of collocation solutions

The main purpose of this section consists in

(i) studying the advantages of applying the collocation method directly to a given neutral functional integro-differential equation or to some equivalent form, and

(ii) pointing to an important question whose answer is not yet known.

As for (i) we shall show that in certain cases (but not in others) it is possible to achieve so-called local superconvergence on some finite subset of the interval on which the equation is to be solved. This survey of results that either have been established before or can readily be proved by the adaptation of classical arguments will allow us to see the above question in a wider perspective.

In order to do this, we shall consider a neutral functional integro-differential equation which is more general than (1.5), namely the scalar analogue of the system (1.1) (with B_{0} = B_{1} = K = 0):

\frac{d}{dt} (a_{0}x(t) + ∫_{−r}^{0} a_{1}(s)x(t + s) ds) = f(t)  (0 < t ≤ T)  (4.1)

with x(s) = ϕ(s) (−r ≤ s ≤ 0). Suppose that u ∈ S_{m+d}^{(d)}(Π_{N}^{(μ)}) (μ = 0, 1, …, M; ξ_{M+1} = T) denotes a collocation approximation to the solution x of (4.1); depending on whether we use (4.1) directly or employ its differentiated form,

a_{0} x′(t) + \frac{d}{dt} (∫_{−r}^{0} a_{1}(s)x(t + s) ds) = f(t)  (0 < t ≤ T)  (4.2)

or its integrated form,

a_{0}x(t) + ∫_{−r}^{0} a_{1}(s)x(t + s) ds = g(t)  (0 < t ≤ T)  (4.3)

with

g(t) := ∫_{0}^{t} f(s) ds + a_{0}x(0) + Dϕ

where Dϕ := ∫_{−r}^{0} a_{1}(s)ϕ(s) ds (see also (1.5) and (1.6)) as the basis for the numerical solution, the value of d will either be d = −1 or d = 0.

Assume that for sufficiently smooth solutions this collocation solution satisfies

||x − u||_{∞} := sup_{t∈I} |x(t) − u(t)| = O(h^{p})

for some p ≥ 1 (usually – but not always: compare Theorem 3.2 – we have p = m). If there exists a finite subset T_{N}^{∗} of I so that

max_{t∈T_{N}^{∗}} |x(t) − u(t)| = O(h^{p^{∗}})  for some p^{∗} > p,  (4.4)

whenever x has sufficient regularity, then u is said to be (locally) superconvergent on T_{N}^{∗} of order p^{∗}. Obvious candidates for T_{N}^{∗} are given by the mesh points in (ξ_{μ}, ξ_{μ+1}],

T_{N}^{∗} = I_{N} := {t_{n}^{(μ)} : 1 ≤ n ≤ N and 0 ≤ μ ≤ M}  (4.5)

and the collocation points,

T_{N}^{∗} = X_{N} := ⋃_{μ=1}^{M} X_{N}^{(μ)}

with X_{N}^{(μ)} as introduced in (2.5).

Consider first the set I_{N}: for which equivalent forms of (4.1) is local superconvergence possible?

Theorem 4.1. Suppose that a_{1} ∈ C^{m+κ}[−r, 0], with a_{1}(0) > 0, f ∈ C^{m+κ−1}(I), and ϕ ∈ C^{m+κ}[−r, 0] for some κ > 0. Assume that the integrated form (4.3) of (4.1) with a_{0} = 0 is solved by collocation in S_{m−1}^{(−1)}(Π_{N}^{(μ)}) (μ = 0, 1, …, M). Then:

(i) There do not exist points {c_{i}} with 0 < c_{1} < … < c_{m} ≤ 1 for which

max_{t∈I_{N}} |x(t) − u(t)| = O(h^{p^{∗}}) with p^{∗} > m.

(ii) If g(0) = Dϕ = 0, and if the collocation parameters {c_{i}} are given by the zeros of the polynomial (s−1)P′_{m}(2s−1) (the m Lobatto points in the left-open interval (0, 1]; P_{m}(2s − 1) denotes the shifted Legendre polynomial of degree m), then

max_{t∈P_{N}} |x(t) − u(t)| ≤ Ch^{m+1}  (4.6)

whenever m is odd and κ ≥ 1. Here, the set P_{N} is given by

P_{N} := \{t = t_{n}^{(μ)} + \frac{h^{(μ)}}{2} : 0 ≤ n ≤ N − 1 and 0 ≤ μ ≤ M\}.

The order p^{∗} = m + 1 in (4.6) is best possible.

The proof of these assertions is based on an essentially straightforward modification of the arguments in the proof of Theorem 3.2 and those given in [6: Section 5.5]. Details are left to the reader.

Example 4.1. If u ∈ S_{0}^{(−1)}(Π_{N}^{(μ)}) (m = 1) and c_{1} = 1 (hence t_{n}^{(μ)} + c_{1}h = t_{n+1}^{(μ)}), and if we denote the restriction of u to the subinterval (t_{n}^{(μ)}, t_{n+1}^{(μ)}] by u_{n+1}^{(μ)}, then the collocation equation corresponding to (4.3) with a_{0} = 0 is given by

h ∫_{0}^{1} a_{1}(h(s − 1)) ds · u_{n+1}^{(μ)} + h ∑_{ℓ=1}^{n−1} ∫_{0}^{1} a_{1}(t_{ℓ}^{(μ)} − t_{n+1}^{(μ)} + sh) ds · u_{ℓ+1}^{(μ)} = G_{μ}(t_{n+1}^{(μ)}) (4.7)

(n = 0, 1, …, N − 1 and μ = 0, 1, …, M), with

G_{μ}(t_{n+1}^{(μ)}) := g(t_{n+1}^{(μ)}) − h ∑_{ℓ=n+1}^{N−1} ∫_{0}^{1} a_{1}(t_{ℓ}^{(μ−1)} − t_{n+1}^{(μ)} + sh) ds · u_{ℓ+1}^{(μ−1)}.

Since a_{1} is continuous and satisfies a_{1}(0) > 0, u_{n+1}^{(μ)} is uniquely defined for all sufficiently small h = h^{(μ)} > 0. According to Theorem 4.1(ii) we then have

max_{(n,μ)} |x(t_{n+ \frac{1}{2}}^{(μ)} ) − u_{μ,n+1}| ≤ Ch^{2}

whenever the solution x is sufficiently regular: at the midpoints of each subinterval [t_{n}^{(μ)}, t_{n+1}^{(μ)}] the collocation solution has an O(h^{2})-error, while at the mesh points the collocation error is in general only O(h).

Is local superconvergence on I_{N} possible if a_{0} ≠ 0 in (4.3)? If collocation is at the Gauss points (i.e. if the {c_{i}} are given by the zeros of the shifted Legendre polynomial P_{m}(2s − 1)), then the answer is again in the negative: a local order p^{∗} > p = m (e.g. p = 2m) can be attained only by applying a certain postprocessing method to u.

Theorem 4.2. In (4.1) assume that a_{0} ≠ 0, a_{1} ∈ C^{m+κ}(I), f ∈ C^{m+κ−1}(I), and ϕ ∈ C^{m+κ}[−r, 0] for some integer κ > 0. Let u ∈ S_{m−1}^{(−1)}(Π_{N}^{(μ)}) (μ = 0, 1, …, M) be the collocation solution to the integrated form (4.3) of (4.1).

(i) If collocation is at the Gauss points, then

||x − u||_{∞} ≤ Ch^{m}  and  max_{t∈I_{N}} |x(t) − u(t)| ≤ Ch^{m}.

In other words, collocation at the Gauss points does not lead to local superconvergence at the mesh points for (4.3).

(ii) The iterated collocation solution corresponding to the above collocation solution u,

u_{it}(t) := \frac{1}{a_{0}} (g(t) − ∫_{−r}^{0} a_{1}(s)u(t + s) ds) (t ∈ I)

satisfies

||x − u_{it}||_{∞} ≤ Ch^{m+1}  and  max_{t∈I_{N}} |x(t) − u_{it}(t)| ≤ Ch^{2m}

whenever κ ≥ m.

Proof. Part (i) and the local superconvergence result in (ii) (p^{∗} = 2m) were established in [2]. In [7] the global superconvergence result (p^{∗} = m + 1) was proved for second-kind Volterra integral equations without delay. The technique used in that proof is readily extended to encompass linear Volterra integral equations with constant delay r > 0

Consider now (4.2), the differentiated form of (4.1), with a_{0} ≠ 0, and assume that a_{1} is smooth on I. We may then write the equation as

x′(t) + a_{1}(0)x(t) − a_{1}(−r)x(t − r) + ∫_{t−r}^{t} a′_{1}(s − t)x(s) ds = f(t) (t ∈ I)

where we have assumed that a_{0} = 1.

Theorem 4.3. Let the assumptions of Theorem 4.2 with κ = m hold, and assume that the solution of (4.7) is approximated by collocation in S_{m}^{(0)}(Π_{N}^{(μ)}) (μ = 0, 1, …, M). If the collocation parameters {c_{i}} are the Gauss points, then

max_{t∈I_{N}} |x(t) − u(t)| ≤ Ch^{2m},

that is, the collocation solution exhibits local superconvergence of (optimal) order p^{∗} = 2m > p = m at the mesh points.

Proof. For the special case where a′_{1}(t) ≡ 0 equation (4.7) reduces to a delay differential equation for which the above assertion is well known (see, for example, [21] and its references, especially to the work of Bellen (1984)).

If a′_{1}(t) ≢ 0 but a_{1}(−r) = 0, the superconvergence results of [2] apply. The local superconvergence property (4.8) for the general case follows by combining the ideas behind the proofs for the two particular cases

5. Concluding remarks

Rewriting the given neutral functional integro-differential equation (4.1) in one of the forms (4.2) or (4.3) requires a integration or differentiation step which in certain concrete practical situations may not be carried out analytically. However, the ”direct” discretization of equation (4.1) (and its more general counterpart (1.1)) remains to be studied. The problem is even more difficult when a_{0} = 0, since necessary and sufficient conditions (analogous to (3.3) in Theorem 3.2 for regular kernels) for the uniform convergence of collocation solutions to (1.5) or (1.6) are not yet known (see also [4]).

The numerical solution of the weakly singular neutral functional integro-differential equation (1.5) has been studied in [12], using the semigroup framework ([8, 9]; also [15]) underlying this functional equation: its underlying idea is the rewriting of (1.5) as a linear first-order hyperbolic partial differential equation with non-local boundary conditions. The low regularity of the analytical solution (which, in general, has unbounded derivatives at the points ξ_{μ} = μr) leads to low-order numerical approximations if the discretization is based on uniform meshes.

This is also true for collocation solutions u ∈ S_{m−1}^{(−1)}(Π_{N}^{(μ)}) (μ = 0, 1, …, M): for any set {c_{i}} of collocation parameters for which u converges uniformly to the solution x of (1.6) (subject to the continuity condition (3.1)), the result on the optimal order of convergence in Theorem 3.2 reduces to

||x − u||_{∞} ≤ Ch^{α}

regardless of the value of m ≥ 1 (see also [4]). This is due to the fact that at t = ξ_{μ}^{+} the solution x of (1.6) behaves like (t − ξ_{μ})^{α} and thus has an unbounded derivative at ξ_{μ}. If m = 1 and c_{1} = 1 (see Example 4.1), we have uniform convergence of u to x, and hence the orders of global and local convergence (at the midpoints P_{N} = \{t_{n}^{(μ)} + \frac{h}{2}\}) are given by

||x − u||_{∞} ≤ Ch^{α}  and  max_{t∈P_{M}} |x(t) − u(t)| ≤ h^{2−α},

respectively.

For collocation solutions in the space S_{m−1}^{(−1)}(Π_{N}^{(μ)}) these low orders of convergence can only be raised to that given by the approximation power of the space if graded meshes of the form (2.1), with grading exponent q = \frac{m}{α}, are employed. The use of such meshes in collocation methods for the equation (1.5) and the corresponding first-kind Volterra integral equation (1.6) is currently being studied by T. L. Herdman and the author.

Acknowledgements. This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC Research Grant No. OGP0009406). Part of the research was carried out in May 1998 at the Mittag-Leffler Institute in Djursholm/Sweden (during the special year on Computational Methods for Differential Equations); the author gratefully acknowledges the hospitality extended to him by the Institute and the director of the program, Professor Vidar Thomeé, during his visit.

Discretization of Integro-Differential Equations 405

References

[1] Brunner, H.: Discretization of Volterra integral equations of the first kind (II). Numer. Math. 30 (1978), 117 – 136.

[2] Brunner, H.: Iterated collocation methods for Volterra integral equations with delay arguments. Math. Comp. 62 (1994), 581 – 599.

[3] Brunner, H.: The numerical solution of neutral Volterra integro-differential equations with delay arguments. Ann. Numer. Math. 1 (1994), 309 – 322.

[4] Brunner, H.: On systems of Volterra difference equations associated with collocation methods for weakly singular Volterra integral equations. In: Proc. Third Intern. Conf. Difference Equ. Appl., Taipei 1997 (to appear).

[5] Brunner, H.: The use of splines in the numerical solution of differential and Volterra integral equations. In: Spline Functions and the Theory of Wavelets, Vol. I and II, Montréal 1996 (ed.: S. Dubuc; CRM Proceedings and Lecture Notes: Vol. 17). Providence (RI): Amer. Math. Soc. 1999 (to appear).

[6] Brunner, H. and P. J. van der Houwen: The Numerical Solution of Volterra Equations (CWI Monographs: Vol. 3). Amsterdam: North-Holland 1986.

[7] Brunner, H. and N. Yan: On global superconvergence of iterated collocation solutions to linear second-kind Volterra integral equations. J. Comput. Appl. Math. 67 (1996), 185 – 189.

[8] Burns, J. A., Herdman, T. L. and H. W. Stech: Linear functional differential equations as semigroups on product spaces. SIAM J. Math. Anal. 14 (1983), 98 – 116.

[9] Burns, J. A., Herdman, T. L. and J. Turi: Neutral functional integro-differential equations with weakly singular kernels. J. Math. Anal. Appl. 145 (1990), 371 – 401.

[10] Hale, J. K. and S. M. Verduyn Lunel: Introduction to Functional Differential Equations (Appl. Math. Sci.: Vol. 99). New York: Springer-Verlag 1993.

[11] Henrici, P.: Elements of Numerical Analysis. New York: John Wiley & Sons 1964.

[12] Ito, K. and J. Turi: Numerical methods for a class of singular integro-differential equations based on semigroup approximation. SIAM J. Numer. Anal. 28 (1991), 1698 – 1722.

[13] Janno, J. and L. v. Wolfersdorf: Inverse problems for identification of memory kernels in viscoelasticity. Math. Methods Appl. Sci. 20 (1997), 291 – 314.

[14] Janno, J. and L. v. Wolfersdorf: Identification of weakly singular memory kernels in heat conduction., Z Angew. Math. Mech. (ZAMM) 77 (1997), 243 – 257.

[15] Kappel, F. and K. P. Zhang: On neutral functional differential equations with nonatomic difference operator. J. Math. Anal. Appl. 113 (1986), 311 – 343.

[16] Kauthen, J.-P. and H. Brunner: Continuous collocation approximations to solutions of first kind Volterra equations. Math. Comp. 66 (1997), 1441 – 1459.

[17] Ortega, J. M.: Numerical Analysis: A Second Course. New York: Academic Press 1972.

[18] Unger, F.and L. v. Wolfersdorf: Inverse Probleme zur Identifikation von Memory-Kernen. Freiberger Forschungshefte, C 458 Mathematik, Technische Universität Bergakademie Freiberg 1995.

[19] Unterreiter, A.: Volterra integral equation models for semiconductor devices. Math. Methods Appl. Sci. 19 (1996), 425 – 450.

[20] Wolfersdorf, L. von: On identification of memory kernels in linear theory of heat conduction. Math. Methods Appl. Sci. 17 (1994), 919 – 932.

[21] Zennaro, M.: Delay differential equations: theory and numerics. In: Theory and Numerics of Ordinary and Partial Differential Equations (ed.: M. Ainsworth; Advances in Numerical Analysis: Vol. 4). Oxford: Clarendon Press 1995, pp. 291 – 333.

Received 03.09.1998