Numer. Math.

52, 129-145(1988) Numerische

Mathematik
9 Springer-Verlag 1988

Convolution Quadrature
and Discretized Operational Calculus. I.
C. Lubich
Institut ffirMathematik und Geometrie,Universit/itInnsbruck, Technikerstrasse 13,
A-6020 Innsbruck, Austria

Summary. Numerical methods are derived for problems in integral equations

(Volterra, Wiener-Hopf equations) and numerical integration (singular inte-
grands, multiple time-scale convolution). The basic tool of this theory is
the numerical approximation of convolution integrals
x

f*g(x)= S f(x--t)g(t)dt (x>O)

0

by convolution quadrature rules. Here approximations to f * g ( x ) on the

grid x = 0 , h, 2h ..... Nh are obtained from a discrete convolution with the
values of g on the same grid. The quadrature weights are determined with
the help of the Laplace transform of f and a linear multistep method. It
is proved that the convolution quadrature method is convergent of the order
of the underlying multistep method.
Subject Classifications: AMS(MOS): 65D30, 65R20, 65L05, 44A55; CR:
G1.9.

1. Introduction

1.1. Derivation of the Methods

We begin our investigations with the numerical approximation of convolution

integrals
x

f(t)g(x-t)dt (x>__O). (1.1)

0

We shall obtain methods of discrete convolution form

~, eo;(h)g(x-.jh) (1.2)
O<=jh<=x
130 C. L u b i c h

where h > 0 is the stepsize, and the quadrature weights coi(h) are the coefficients
of the power series

F(b(~)/h)= ~ ooj(h)~J. (1.3)

j=o

Here F is the Laplace transform of f and 6 ( ~ ) = ~ b j ~ J is the quotient of the

0
generating polynomials of a linear multistep method.
As will be pointed out in Part II of this work, quadrature methods of this
type are particularly suited to the approximation of convolution integrals whose
kernel f(t) is singular or has components at different time-scales; to the numeri-
cal evaluation of the expressions arising in classical operational calculus where
the Laplace transform F of the convolution kernel (and not the kernel f itself)
is known a priori; and to obtain stable discretizations of integral equations.
From the control engineer's viewpoint our results can be seen as results on
the discrete-time simulation of a continuous-time linear system given via its
transfer function F(s).
To explain this approach we start from the Laplace inversion formula

f ( t ) = ~ n / r~ F(2)e~td2 (t>0) (1.4)

which holds if, for example, we assume

7~
F(s) is analytic in a sector [ a r g ( s - c)[ < n -- ~o with ~o< ~, c elR and satisfies there
]F(s)] < m - I s ] - " for some M < oo, # > 0 . (1.5)
The contour F can then be chosen as running from oo-e -i("-~) to oo-e i~-~)
within the sector of analyticity of F(s). A condition equivalent to (1.5) is that
f(t) is analytic and of (at most) exponential growth in a sector containing the
positive real half-line and satisfies f ( t ) = O(t ~- 1) as t--, 0 within the sector. The
assumption/~ > 0 guarantees in particular that f(t) is locally integrable. Examples
include fractional powers, logarithms and most special functions.
The approach (1.2), (1.3) now comes about in the following way: Inserting
(1.4) into (1.1) and reversing the order of integration gives
x 1 x
S f(t)g(x-t)dt=~n i ~rF(2)~o eZtg(x--t)dtd2" (1.6)
0
We now approximate the inner integral by applying a linear multistep method
to the differential equation y' = 2y + g, y(0) = 0,
k k
Z aJYn+J-k-=h Z flJ(2Yn+J k+g((n+j-k)h)) (n>O), (1.7)
j=O j=O
C o n v o l u t i o n Q u a d r a t u r e a n d Discretized O p e r a t i o n a l Calculus. I. 131

with starting values Y-k . . . . . y _ ~ = 0 , and with g~C[0, oo) extended by 0 to

the negative real axis. Multiplying (1.7) by ~" and summing over n from 0 to
oe we obtain

(% (k + . . . + ~k)"Y(O : (flO (k +... + ilk)" (h 2' y(() + h-g(O),

co co
with the generating (formal) power series Y(0 = ~Y, (", g ( 0 = ~ g ( n h) (". Solving
0 0
this equation for Y(0 we find that y, is the n-th coefficient of the power series
( b ( O / h - 2) -1 g((), where

6 (0 = (~o ~k +... + C~k _ i ~ + C~k)/(flO~k + fig - I ~ + fig)" (t .8)

Hence the resulting approximation of (1.6) at x = nh is the n-th coefficient of

1 /6(0 2)-lg(Od2=F(6(~h) ) (1.9)

where the equality holds by Cauchy's integral formula. The coefficients of the
right-hand side of (1.9) are the Cauchy product of the two sequences {~o;(h)}
of (1.3) and {g(jh)}. This gives finally (1.2).
The above derivation indicates also how to obtain error estimates: One
studies first the error introduced by approximating the inner integral in (1.6)
by the multistep method (1.7), then multiplies the obtained error bound by
the bound (1.5) of F(2) and integrates along the contour F. This will actually
be carried out in Sects. 2 and 3.
Of the linear multistep method we shall assume that it is A(a)-stable with
a>(p of (1.5), stable in a neighbourhood of infinity, strongly zero-stable and
consistent of order p. In terms of 6(0, these conditions can be expressed as:

3(0 is analytic and without zeros in a neighbourhood of the

closed unit disc ]~1__<1, with the exception of a zero at ~ = 1. (1.10a)

[arg6(0l<n-c~ for [ ~ l < l , forsome c~>q~. (1.10b)

1
7 6 ( e - h ) = l + O ( h p) forsome p > l . (1.10c)
n

Well-known examples are the backward differentiation formulas of order p < 6,

given by 3(~)= ~-(1--0 i, with e = 90 ~ 90 ~ 88 ~ 73 ~ 51 ~ 18~ for p--1, ..., 6,
respectively, i= 1
Condition (1.5) on the transfer function (or symbol) F(s) of the convolution
(1.1), and condition (I.10) on the discretization method will serve as the main
assumption for all results of this paper.
Remark. (1.3) is well-defined if 5o/h is in the domain of analyticity of F(s). Since
30 > 0 by (1.10), this is satisfied at least for sufficiently small h >0. The integrals
in (1.6) and (1.9) are absolutely convergent, because ~ > 0 in (1.5).
132 C. Lubich

1.2. Symbolic Notation and Historical Remarks

It will be convenient to use transfer function notation for the convolution (1.1),
x

F(s)g(x)=f*g(x)= ~ f ( t ) g ( x - t ) d t (s: complex variable). (1.i1)

0

We use the abbreviation Sh=(~(~)/h and employ discrete transfer function nota-
tion for (1.2), (1.3),

F(Sh)g(x)= ~, coi(h)g(x--jh). (1.12)

O<jh<x

Note that F(Sh)=~'c%(h)( ~ is the discrete Laplace transform of the sequence

0

{coj(h)}.)

With this notation we obtain from (1.9) the "Cauchy integral formula"

1
F (s,) g (x) = T ~ ~ F (2~). (sh - .~) -1 g (x). d ,~, (t. 13 )

and the same relation, with s formally in place of sh, holds by (1.6), since ( s - 2)-
is the Laplace transform of e ~'.
We shall use the notation (1.11) also for distributional convolution, in partic-
ular for differentiation

sg=~'.g=g'+g(0).~ (Dirac's 6). (1.14)

As usual, we identify distributions with functions outside their singular support

(which wilt be at most {0} in our applications), so that (sg)(x)=g'(x) for x > 0 .
The derivative (1.14) is discretized by the backward difference approximation

1
Shg(X)=~ Z 6~g(x--jh). (1.15)
O<=jh<=x
With this in mind, the discretization process leading from (t.11) to (1.12) has
a simple interpretation: The symbol of differentiation, s, is replaced by the symbol
of a backward difference quotient, Sh. This idea of a discretized operational calcu-
lus has actually very old origins and seems to have been rediscovered several
times. Of particular historical interest in this context is the work of Post [12]
who has shown pointwise convergence of F(Sh) g(x) as h ~ 0 for the simple
difference quotient sh g(x)=(g(x)-g(x--h))/h, thereby generalizing earlier work
of Liouville [10] and Gr/inwald [8] on fractional integrals, i.e. the case F(s)= s -u.
Post's work appears to have gone unnoticed in numerical analysis (in which
he himself showed no interest), and the author of the present paper is not
aware of any results in the literature going beyond those of Post.
Convolution Quadrature and Discretized Operational Calculus. I. 133

1.3. Outline of the Paper

In Sect. 2 we study the convergence properties of the linear multistep method
(1.7) (i.e., of (Sh--2)- 1g(x)) for 2 varying in a sector of the complex left half-plane.
Such problems have previously been studied in the numerical analysis of para-
bolic differential equations by Ztfimal [15] and Crouzeix and Raviart [2], see
also Gekeler [7]. Here we give (and require) stronger estimates in terms of
the input g, using different techniques.
Section 3 contains the first main results. The convergence properties of con-
volution quadrature rules F(Sh)g(x) are now immediately obtained by integra-
ting the error bounds of Sect. 2 along the contour F. The convergence order
of the underlying multistep method, p, is obtained for geCP[0, oo) by adding
a few correction terms to F(Sh) g(x),
p-2
F(Sh)~g(x)=F(sh)g(x)+ ~ w,j(h)g(jh) at x=nh. (1.16)
j=0

These are required because of the special choice of starting values in (1.7), unless
g near 0 or f near x=nh is "small".
In Sect. 4 we turn our attention to the coefficients {o,(h) of (1.3) themselves,
or rather to f,(h)=oo,(h)/h which is shown to approximate the inverse Laplace
transform f(t) of F(s) at t = nh (away from 0) with the full order of the multistep
method (1.7). This result (which comes rather unexpected) extends the classical
Post-Widder inversion formula. Its numerical importance is perhaps not so much
in the numerical inversion of Laplace transforms, but in the fact that it justifies
replacing the weights co,(h) by hf(nh) for nh bounded away from 0, and in
a very simple choice of the starting weights w,j(h) of (1.16).
In Sect. 5 we give some extensions of the results in Sect. 3. These concern
the approximation of derivatives
(d~k x
dx] Jof ( t ) g ( x - t ) d t (x>O)

("finite part integrals") and the convergence of (1.2) under weakened smoothness
assumptions on g, in particular for the important case g (t)= t ~- 1~(t) with fl > 0,
smooth, for which p-th order convergence can again be preserved by a modifi-
cation of the form (1.16).
Topics which have been omitted are an extension to the non-sectorial,
"hyperbolic" case, where F(s) is analytic and suitably bounded only in some
half-plane Re s>c. For A-stable linear multistep methods (1.7) one can still
obtain convergence results (of order at most 2, according to Dahlquist's [4]
well-known order barrier.) Another omission concerns convolution quadratures
2
in which the underlying multistep method is the trapezoidal rule, Sh=~(1
- ~ ) / ( i + ~ ) , which does not satisfy (1.10a). It appears, however, in the control
literature where it is known as Tustin's method (after Tustin [13]; see, e.g.,
Franklin and Powell [6]).
134 c. Lubich

2. L i n e a r M u l t i s t e p M e t h o d s A p p l i e d to y ' = 2y+g
with ~ Varying in a Sector

In this section we derive some technical estimates which will be required for
the subsequent development. We begin with some preparation.
Since F(Sh)g is the convolution of g with the sequence {co,(h)}~, the associati-
vity of convolution yields

F(Sh)(g~ * g2) = (F(Sh)ga)* g2, (2.1)

and similarly

F(s)(gl * g2)=(F (s) gl)* g2 . (2.2)

We shall use these relations in place of the usual Peano kernel technique. They
permit to reduce the study of the error F(Sh) g(x) F(s) g(x) to the polynomial
case g (x)= x q, since

p- 1 g(q)(o) 1 1
g (x) = Y' ~ x q+ ~ ( t ' - * g('))(x), (2.3)
q=O

the Taylor expansion of g at 0.

Here we have denoted by t p-I the function t~--~tp-~ on 1-0, oe), and the
same notation will also be used for other exponents in the following.
We now study the error of (s,- 2)- 1 g(x), i.e., of the linear multistep method
(1.7).
[,emma 2.1. Let Sh= 6(()/h satisfy (1.10).
a) For l a r g ( - - 2 ) l = e <c~ and IhRl_>r>0 we have

l(Sh--R)-ltq(x)--(s--R)-ltq(x)l< C~.p x/h.hq ( q = 0 , 1 . . . . p)

=lh21 I,ll
with C=C(r) and p=p(r)< 1 independent of h~(O, 1], x > h and 2.
TC
b) To arbitrary 0 < ~ : < 1 there exists t o > 0 such that for l a r g ( - 2 ) l < e ' < ~
and Ih21<ro
IC[Cile'ax['hq+l ( q = 0 , 1 ..... p - l )
I(Sh -- 2 ) - 1 t q (X) - - ( S - - 2 ) - 1 t q (X) [ < [ e ~ ax 1" h;
121 (q=P)
with C independent of he(O, 1], x > 0 and 2.
Remark. Except for q = 0 the estimate in a) holds also for 0__<x<h. (In this
interval F(sh)g(x)=F(6o/h).g(x), and error estimates are easily derived.)
It is easily verified that

(Sh--2)- 1g(X)--(S-- 2)- I g(x)= --(Sh-- 2)- 1(Shy--y')(X) (2.4)

Convolution Quadrature and Discretized Operational Calculus. I. 135

where y = ( s - 2 ) - 1 g , the solution of y ' = 2y + g, y ( 0 ) = 0. W e study Sn y-Y' in Lem-

ma 2.2 (consistency), and ( s h - 2 ) -1 in L e m m a 2.3 (stability). L e m m a 2.1 is de-
rived from these estimates without new difficulty.

L e m m a 2.2. a) F o r Iarg ( - 2) 1< c~'< 2 and I h 21 > r > 0 we have

C
[(S h -- S) (S -- 2) - 1 t q (X) I ~ ]h~T" px/h. h q (q = O, 1 . . . . , p)

with C = C (r) and p-=p(r) < t independent o f h ~(O, I], x > h and 2.
b) There exists r > O such that f o r ] h 2 1 < r

I(sh--s)(s--2) -1 t~(x)l< C1.12P-qeZXl.hP+C2.px/h.h q (q=O, 1. . . . . p)

with C t , C 2 and p < 1 independent o f he(O, 1], x >O and 2.

Proof. First we show that

](Sh--S)tq(x)l<=C.pX/h.h q - t ( q = 0 , 1. . . . . p) (2.5)

with C and p < 1 independent of h and x > 0.

Taylor expansion ia (I.10c) shows (cf. L e m m a 5.3 of Henrici [9]) that
oO

1 ~=oSi(x_jh)q=q.xq_ 1 for q = 0 , 1,.. p, for all h and x.

hj

By (1.10a), 5 i = O(p~) for some p < 1. Let n o w x - ~ n h + O h with 0 < 0 < 1. Then

[(Sh--s) tq(x)]= h ~ ( S j ( x - - j h ) q - - q ' x q - 1

j=o

i co oo

This proves (2.5).

a) By (2.1) and (2.2),

(s h - s)(s - 2) -1 t q = (s h _ s)(e a'* t q) = e ~' * (Sh -- S) t q.

For ~ = e -~ with ~ > p the estimate (2.5) yields thus

x

I(Sh - - S)(S-- 2) - 1 tq(x) t ~ C. h q- 1. e - ~x/h S e'tRe (h~) + Wh d t.

0

(For q = 0 there appears an additional term ]e~l on the right-hand side of

this estimate, caused by s 1 = ~ (Dirac's t~). F o r x > h this term does, however,
not affect the final estimate, and is therefore omitted in the following.) Fix
136 c. Lubich

now r > 0 . We can choose ~7='7(r)>0 such that

7~
Re(h2)+<--c<0 for l a r g ( - s [hR[>=r.

Then the upper limit of the above integral can be extended to infinity, and
we obtain the estimate

C.h .e-~X/h.hq- I<

C .~x/h.hq"
[(sh -- s) (s -- 2)- 1 t q(x)] <
- [Re(h2)+~ --]h2I
This gives a).
b) By partial integration,
1 H tq + 1 tq + 2 tp tp
~. (s-2)-ltq q!
:__,cAt __
( q + I~.V+ 2 ~ + ... _~,~p--q-- 1 _p!
_ _~_,~p-q ~ . , egt.

Using (2.5) with q + 1, ..., p instead of q we obtain thus for [h2[ <const.

[(Sh--S)(S--A)-ltq(x)[<Cz.pX/h.hq+C1.]2P-q].hV-l.(fh*]ea'[)(x). (2.6)

Let now p = e -~ and choose r > 0 so small that

Re(h2)+7>c>0 for [ h 2 [ < r .

Then
x h
(p,/n, l ea, l)(x) = [eXX[. ~ e - t [Re(h2)+ 7]/h d t <=- IeX~ I.
e
0

Inserting this estimate in (2.6) yields the result. []

We have ( S n - 2 ) - ~ = h ( 6 ( O - h 2)-1. F o r the coefficients of this power series
we have the following stability estimate.
L e m m a 2.3. a) For [arg ( -- z) I _-<e' < e and [z [> r > 0 the power series (8 (0 - z)- 1
MOO
is rnajorized by ~ ~o p" (n with M = M (r) and p = p(r) < 1 independent of z.

7~
b) Let 0 < < 1. There exists r = r(tz) > 0 such that for ] a r g ( - z)[ __<e' < ~ and
oo
]z[<r the series ( 6 ( O - z ) - 1 is majorized by M.~le"KZ[ ~n with M independent
ofz. o
Proof a) By (1.10) the series 6 ( Q - z is for [zl>r, [ a r g ( - z ) l < ~ ' < a analytic
and without zeros in a disc ]([__<1/p with p < 1, where p is independent of
z (but depends on r). Also

M
sup j ( b ( O - z ) - l [ < j ~ T
I~l_-<l/p
Convolution Quadrature and DiscretizedOperationalCalculus. I. 137

with M < oo independent of z. Hence Cauchy's estimate (e.g. Ahlfors [-1, p. 98])
yields that the n-th coefficient of (6(~)--z)-1 is bounded by Mp"/Iz[. This gives
a).
b) Since 3(1)=0, the estimates in a) deteriorate as r ~ 0 . For small [z[ we
use instead (1.10c). We define R(z) implicitly by
6(1/R(z))= z, R(0)= 1. (2.7)
We split
1 1 --R(z)~ 1
6([)--z 6([)--z 1--R(z)[

From (1.10c) we obtain

R(z)=e~+O(zV+l).
Hence there exists r--r(~c)> 0 such that

~) = [ e (1-K)~+O(z p+I)I=]I + ( 1 - g ) z + O ( z 2 ) l < = l

r 7/:
for [z] __<r and [ a r g ( - z)[ __<a < ~. Therefore

1 1
is majorized by
1-R(z)~ 1-[e~Z[~"

On the other hand, using (1.10a) Cauchy's estimate yields

1--R(z)( M
6(()--z is majorized by l-p('

with M and p < ! independent of z. Hence

t M 1
~(~)--z is majorized by 1-p~ 1-[e~Z[~"

We can choose r so small that e~'p<l. For the n-th coefficient of the above
series we then obtain

pJle~Zl"-J<:leK~["~ pJleKZ]-J<:le~Zl"--
j=O j=O
1 - p e ~r"

This proves b). []

3. Convergence of Operational Quadratures

We are now in a position to show convergence of F(Sh)g(x).

Theorem 3.1. Under the assumptions (1.5), (1.10) we have
138 C. Lubich

IF (sh) g (x) - F (s) g (x)[< C - x ~- 1. {h ]g (0) l + ... + h p- 1[ g(p- 2)(0)1

+ hP.([g ~p- 1)(0)1 + x - m a x [gr
O__<t__<x

where the constant C does not depend on he(O, [;f], x~[h, 2] with f i x e d ,2< oo,
and g ~ C p [0, ~].
Remark. If additionally F(s) is exponentially stable (i.e., F(s) is analytic in some
half-plane Re s > Y with negative 7), T h e o r e m 3.1 extends to a result of uniform
convergence on the half-line x => 1, with exponential decay in the low-order error
terms caused by g(0) . . . . , g(p-2)(0). (The p r o o f is even slightly simpler, since
no contour-shifting is necessary on intervals b o u n d e d away from 0.)
Proof. By (2.1)-(2.3),
1
E g~q)(O) F (sh) tq(x) + T~ pS -- l-)~ ' ( F (sh) t p- l * g(P))(x),
F (Sh) g (X) = p-1 (3.1)
q=O

and the same holds with s formally in place of Sh. By (1.13),

1
F (Sh) t q(X) -- F (s) t q(x) = ~ : F (,l)- [(Sh -- 2)-1 t q(X)-- (S-- 2)-1 t q(X)]- d 2.

We now substitute w = 2 x in the integral. Since the c o n t o u r x F comes arbitrarily

close to the origin for small x, we replace it by an equivalent c o n t o u r / ' 1 which
is independent of x~[h, 2], has a positive distance to the origin and, a p a r t
from a c o m p a c t subset, is again contained in a sector ]arg(-w)[=< a ' < a.
For[w]=[2x[<const. we obtain from the classical convergence theory of
linear multistep m e t h o d s [3, 9]

] ( s h - 2 ) - l t q ( x ) - ( s - 2 ) - l t q ( x ) l < = C . h q+i for x e [ h , 2].

F o r t arg ( - w) 1= [a r g ( - - 2) ] __<c( < ~ L e m m a 2.1 gives

/ px/h
I(Sh-- 2)- I t q (X)-- (S-- 2) - 1 t q(x)[ _--<C. (I eZX/2 [" h q + 1 + zha) 9
We observe p~/h<=C ' h / x and insert the a b o v e estimates and the b o u n d in (1.5)
into the integral representation o v e r / ' 1 . This yields the estimate

If(sh) ta(x)-- F(s) ta(x) I _--<C - x " - l - h a + 1

for x e [ h , 2], q = 0 , 1. . . . . p - 1 . (3.2)

(The factor x " - 1 comes from [F(w/x)l < M - x " . I w l - " and d2 = d w / x . The integral
remains absolutely convergent, because # > 0.)
Further, for 0 < x < h we have

F(sh ) t p - 1 (x)=Coo(h)x p- 1 = F(fo/h) x p- 1 =O(hU+ p- 1)

where we have used (1.5) in the last estimate. Since f ( t ) = O ( t "-~) as t ~ O we

have also
F ( s ) t P - l ( x ) = O ( h "+p-i) for O < x < h .

Using (3.1) we obtain finally the estimate of the theorem. []

Convolution Quadrature and Discretized Operational Calculus. I. 139

We can get rid of the low order error terms in Theorem 3.1 by a simple
modification:
p--2
F(sh)- g(x)=F(Sh)g(x)+ ~, w.~(h)g(jh) at x=nh, (3.3)
j=O

where the correction quadrature weights w,~(h) are determined such that the
quadrature formula becomes exact for polynomials up to degree p - 2 :

f(Sh) ~ t q ( x ) = F ( s ) tq(x) for q = 0 , 1, ..., p--2. (3.4)

This gives a Vandermonde system of linear equations for w,i(h ). An alternative

to (3.4) (for x bounded away from 0) will be given in Corollary 4.2 below.
Corollary 3.2. Under the assumptions (1.5), (1.10) the method (3.3), (3.4) satisfies
for geCP[O, ,23
1F(sn)- g(x) - F(s) g(x) l <=C. x u- 1 h p

with C independent of hE(0, h-] and x~[h, ~].

1
Proof With r(x)= ( p - 2 ) ~ (tp-2*g(P-1))(x)' the remainder in the Taylor expan-
sion of g at 0, we have by (3.4) at x = n h

F (Sh)~ g(x)-- F (s) g(x) = F (Sh)- r(x) -- F (s) r(x)

p-2
= F(Sh) r(x)-- F(s) r(x) + ~, w,2(h) r(jh). (3.5)
j=o

The weights w.j(h) are determined from the Vandermonde system

p--2
w,j(h).jqhq=V(s)tq(x)-F(Sh)tq(x) (q=0, 1..... p - 2 )
j=0

where the right-hand side is bounded by C. x u- 1. hq+ 1 because of Theorem 3.1.

Cancelling the factor hq we see that w.~(h) are bounded by

]w,,~(h)l<C.xU-l-h (x=nh).
Moreover,

I r ( j h ) l ~ C . h p-1 for j = 0 , 1, . . . , p - - 2 .

Inserting these estimates in (3.5) and applying Theorem 3.1 (with r instead of
g) gives the result. []

4. Approximation of the Inverse Laplace Transform

We define f,(h) by GO

F(b(O/h) = h ~ f.(h) (" (4.1)

n=O
140 C. Lubich

so that
f,(h)=co,(h)/h (4.2)

with e),(h) given by (1.3). We show that f,(h) approximate the inverse Laplace
transform of F(s). For the special case 6(~)= 1 --~ (backward Euler) this is known
as the Post-Widder inversion formula [12, 14].
Theorem 4.1. Assume (1.5), (1.10), and let f (t) denote the inverse Laplace transform
of F(s). Then
[L(h)-f(nh)l<=C.xU-l-V.h p (x=nh)
where the constant C does not depend on he(O, hi and xe [h, ~] with fixed ,2 < oo.
Remarks. a) For exponentially stable F(s) (i.e., exponentially decaying f ) Theo-
rem 4.1 can be extended to a uniform p-th order estimate on the half-line x > 1,
with exponential decay of the (absolute) error.
b) For F (s) = (s-- 2)- 1 the coefficients f, (h 2) are the solution of a linear mult-
istep method applied to y'=2y with the special choice of starting values Y0
= 1 / ( 8 o - h 2), y_ k . . . . . Y 1 = 0. It is rather unexpected that f , (h 2) are p-th order
approximations to e nh;~ for n h bounded away from 0.
Proof a) By Cauchy's integral formula,
I F(,~)
F(6(~)/h)=~ni ~r (6(~)/h--2) d2

where F is again the path of (1.4). We denote the n-th coefficient of the power
series ( 6 ( ~ ) - z ) - 1 by e.(z) and thus have

f.( h ) = ~ /Ir ~ F(2) e.(h2) d2 (4.3)

which closely resembles the inversion formula

f(nh)=2~t ~rF(2)e"h~d2. (4.4)

We are thus led to study the difference en(h2)-e "ha. As in the proof of Theo-
rem 3.1, we substitute w = 2x ( x = nh) and replace the contour x F by the contour
F 1 (independent of x) given there.
b) For z=h2 with ]arg(-z)] =<~'<a and ]z] =>r>0 Lemma 2.3a) shows

e,(z)=O =e"z+O for some p < l . (4.5)

c) Let now Iz] < r, r sufficiently small. As in the proof of Lemma 2.3 b) let
R (z) be defined by

6 (I/R (z)) = z (4.6)

C o n v o l u t i o n Q u a d r a t u r e a n d Discretized O p e r a t i o n a l Calculus. I. 141

so that
R(z)=e~+O(zP+l). (4.7)
7"C
Let first Ia r g ( - z)[ < a' < ~. Using IR (z) l < Ie~l for some 0 < K< 1, the relation

n-1
R (z)" -- e "~ = ~ R (z)" - ' - j e jz [R (z) -- e =]
j=O
yields
, 7Z
IR(z)"-e"ZI<=C.z~le~"Z I for Iarg(-z)l_<~ < ~ , Izl<=r. (4.8)

Using [R (z) l< eelzl for some ff > 1, one has the well-known estimate

[R(z)"-en=l<C.z p for Inz]<const., lzl<=r. (4.9)

d) We split again

1 1-R(z)( 1
6(O--z- 6(O--z 1--R(z)~

and recall that by Cauchy's estimate the coefficients of(1 - R(z)O/(a(O--z), which
we denote by ?, (z), satisfy

?,(z)=O(p") forsome p < l (uniformlyin[zl<r). (4.10)

We write

e,(z)= ~ ?j(z)R(z)"-J=R(z) ". ~, ?j(z)R(z) - ~ - ~ ?j(z)R(z) "-j. (4.11)

j=0 j=0 j=n+l

The last term is O(p") for p of (4.10), provided that pea'r< 1 (~7 as above). By
definition of ~/,(z) and by de l'Hospital's rule

Tj(z)R(z) - j 1-R(z)_( _ -R(z)

6 ( O - - z ~=gm-' 6'(1/R(z))"
j=O

Differentiating (4.6) we obtain that this expression equals

R'(z) e =+ 0 (z p)
~ YJ(z)R(z)-J- R(z) e z + O ( z p+I) = 1 + O(zP). (4.12)
j=0

Hence we obtain from (4.8) and (4.9)

7~
le,(z)-e"=l ~ C . ( z p le~"Z I + p") for [ a r g ( - z ) l < e ' < ~ , Izl~r, (4.13)

le.(z)-e"zl < C.(zP + p ") for Inzl_<const., Izl~r. (4.14)

142 c. Lubich

e) Inserting the estimates (4.5), (4.13), (4.14) and the bound of (1.5) yields

~C.xU-l.n-P=C.xU-a-P.hP. []
Theorem 4.1 is also of practical interest in the approximation of F(s)g(x).
Let
p-2
F(Sh)~g(x)=F(Sh)g(x)+ ~ co,_i(h)cjg(jh ) (x=nh) (4.15)
j=O

where cj do not depend on n and h (cf. (3.3).)

Corollary 4.2. Assume (1.5), (1.10), and choose ci in (4.15) as the correction weights
of the p-th order Newton-Gregory formula (end-point correction of the trapezoidal
rule). Then the method (4.15) satisfies for g~CP[0, oo)
[ F(sh) ~ g(x)-- F (s) g(x) t <=C. hp

where the constant C does not depend on h ~ (0, h-] and x = n h ~ [Xo, x] with fixed
Xo>0.
Examples. p = 2: Co= -- 1/2 (trapezoidal rule)
p = 3 : Co= - 7 / 1 2 , c~=1/12
p = 4 : C o = - 5 / 8 , c1=1/6, c 2 = - 1 / 2 4 ,

Proof Let u(t) a smooth function with

0 for t <-Xo/3
u(t) =
1 for t>2Xo/3,

and let v(t)= 1 - u ( t ) . By Theorem 3.1,

F (sh)~ (u g)(x) = F (s) (u g)(x) + 0 (hV).

Further, at x = n h,
p--2
F(Sh) ~ (vg)(x) = h ~ f,_j(h)(vg)(jh) + h ~ f,_j(h) c~(vg)(jh)
j=o j=o
p-2
= h ~ f ( ( n - j ) h ) ( v g ) ( j h ) + h ~, f ( ( n - j ) h) c~(vg)(jh) + 0 (h;)
j=O j=0

= V(s)(vg)(x) + O(hP).

Here we have used Theorem 4.1 in the second line, and p-th order convergence
of the Newton-Gregory rule in the last line. Since g - - u g + v g , this yields the
Newton-Gregory result. []
Convolution Quadrature and DiscretizedOperationalCalculus. I. 143

5. Extensions: Derivatives of Convolution Integrals, Non-Smooth Input

In this section we extend Theorem 3.1 in two directions. First we consider the
approximation of

d k~of(t)g(x_t)dt=skF(s)g(x ) (X > 0).

Theorem 5.1. Under the assumptions (1.5), (1.10) we have for integer k >=0
]sk F(sh) g(X)-- sk F(s) g(X)] ~ C" Xu- 1 -k. {h [g(0)] + h p- 1 [g(P- z)(0)]
+ h p.([gr 1)(0)] + x [g<P)(0)]+ . . . + x k ]gtp+k- 1)(0)] +Xk+ 1 max IgtP+*)(t)])}
O<_t<_x

where the constant C does not depend on he(0, h-I, x~[h, Y:] with ~ < o c , and
g ~ C p+k [0, 2].
Remarks. a) Theorem 5.1 can be interpreted as an extension of Theorem 3.1
to the case p < 0 in (1.5). (Take F(s)=skF(s), fi =/~--k.)
b) For exponentially stable F(s) the result can again be extended to the
whole half-line, with exponential decay of the error terms involving
g(O). . . . . g(P+k-1)(0).
c) Corollary 3.2 remains valid with # - k instead of #.
d) Theorem 4.1 can be extended similarly.
The proof of Theorem 5.1 is documented in LuNch [11]. For the sake of
brevity we shall here only sketch the main arguments.
Outline of the proof We write

s k F (Sh) g - - sk F (s) g = (sk - - s k) F (s) g + s k (F (Sh)g -- F (s) g).

The first term on the right-hand side can be shown to satisfy an estimate of
the desired type by using analyticity and growth properties (as t ~ 0) off(t).
The estimate for the second term follows from (1.13) and the error bound
I Skh [(Sh __ ,~,) - 1 t q __ (S -- 2) - I t o] (X)]

<C.(12k e~ZXl hq+14 l + [ Ih 2 l pX/h.hq+l_k)

= (q=O, 1..... p - - l )

(
<=C. ]2k+P-q-le~XlhP+ l+lh2~PX/h.h q+l-k ' ) ( q = p ..... p + k - 1 )

where ~c>0, p < l , and C are independent of he(0, 1], x > h and 2 with ]arg(
-,t)t<~'<~.
This estimate is obtained using techniques similar to those of Sect. 2 and
the following asymptotic expansion for R(z) defined in (2.7):

R (z)" = e"Z [ 1 + Pl (n z) z~ + ... + Pu + 1 - p(n z) zu] + Rem (n, z)

144 C. Lubich

where Pj are polynomials (of degree j), Pj(0)=0, and the remainder satisfies for
7~
[ a r g ( - z)] < a' < ~- and [z] < r (r sufficiently small)

[Rem(n,z)]<C[e~"Z.z N+I} forsome x > 0 . []

Next we turn to the situation where g(t) is smooth on t > 0 , but has an asymptotic
expansion in fractional powers of t at t = 0. Here we have the following extension
of Theorem 3.1 (or equivalently, of (3.2)). Convergence of order p for F(Sh)g
can be restored by eliminating low-order error terms as in Corollary 3.2.
Theorem 5.2. Under the assumptions (1.5), (1.10) we have

1 f C'Xu-I'h# for O<fl<=p

[F(Sh)ta-l(x)--F(s)ta- ( x ) l < ~ C . x U - l + a - P - h p for fl>p (fl real)

where the constant C is independent of he(0, h-] and x ~ [h, if] with ~ < oo.
Remark. For 0 < fl < 1 the above estimate is understood to hold for

F(Sh) g(x)= ~ eoj(h)g ( x - j h )

O<_jh<x-h

which differs from (1.12) in the omission of the last term of the sum, in order
that g is not sampled arbitrarily close to the singularity at 0. All the previous
results remain valid also for this modified definition.
The proof uses Theorem 5.1 and the following lemma.
Lemma 5.3. Let #, v ~ O , - 1,--2 .... real numbers. The convolution of two
sequences u, = O(n u- 1) and vn = O(n ~- 1) satisfies

~ ujv._j=O(n~-l), where ? = m a x { ~ , v,#+v}.

j=O

This result is easily derived using the relation [5, p. 47]

n --I ~ n u- 1
(-1) ( n )=~(~[l+O(n-1)].

Remark. This formula is the special case F(s)=s -u, 6 ( 0 = 1--(, nh= 1 of Theo-
rem 4.1.
Proof of Theorem 5.2. a) If fl > p, we use
t ~- 1 t p- 1 t #-p- 1

r(13) - r ( p ) * ~

By (2.1) and (2.2),

[F(Sh) -- F(s)] (tp-1 , t/~-p- 1) = ([f(sh) -- t ( s ) ] tP- 1), t p- p- 1,

and the result follows from (3.2), since t#-p- 1 is locally integrable.
Convolution Quadrature and Discretized Operational Calculus. I. 145

b) For fl < p we use

t~ - 1 tp - 1
:S p a __
r(fl) r(p)
and write
[f(sh)-- F(s)] s "-a tp- i

:[S~(shaF(sh))--sP(S-aF(S))]tP-I--F(Sh)[S~Sh#--SPS-#]tP-1. (5.1)

By Theorem 5.1 we have for x~[h, ~]

[[s~ ( s ; a F (Sh)) -- S p ( S - a F (s))] t p - 1 (x)[ < C . x ~ + a - 1 - p. h p < C . x u ~. h a (5.2)

and
] [ Sp Sh fl _ _ S p S - fl] t p - 1 (X)[ ~ C " X fl - 1 - p . h p. (5.3)

Further, by (l.5) and Theorem 4.1,

[o%(h)]<C.h ~, [e),(h)l<C.h".n "-1 (n>l, nh<~). (5.4)

By Lemma 5.3 the estimates (5.3), (5.4) yield

IF(Sh)[Sg-a--sP-a]tP-l(x)[=l ~ ogj(h)'[s~-fl--sP-a]tP-l(x--jh)[
O<jh<x-h
< C x " - 1. h a. (5.5)

Now (5.1), (5.2) and (5.5) give the result. []

References
I. Ahlfors, L.: Complex analysis. New York: McGraw-Hill (1953)
2. Crouzeix, M., Raviart, P.A.: Approximation des probl~mes d'~volution, Premiere Partie. Report,
Universit6 de Rennes (1980)
3. Dahlquist, G.: Convergence and stability in the numerical integration of ordinary differential
equations. Math. Scand. 4, 33-53 (1956)
4. Dahlquist, G. : A special stability problem for linear multistep methods. BIT 3, 27~43 (1963)
5. Erd61yi, A. (ed.): Higher transcendental functions I. New York: McGraw-Hill 1953
6. Franklin, G.F., Powell, J.D.: Digital control of dynamic systems. Reading: Addison-Wesley 1980
7. Gekeler, E.: Discretization methods for stable initial value problems. SLNM 1044. Berlin, Heidel-
berg, New York: Springer 1984
8. Gri.inwald, A.K.: Ober ,,begrenzte" Derivationen und deren Anwendung. Z. Math. Phys. 12, 441
480 (1867)
9. Henrici, P.: Discrete variable methods in ordinary differential equations. New York: Wiley 1962
t0. Liouville, J. : M6moire sur le calcul des diff6rentielles ~ indices quelconques. J. l'Ecole Polytechnique
13, 71-162 (1832)
11. Lubich, Ch.: Discretized operational calculus I. Technical Report, Inst. f. Math. u. Geom., Univer-
sit/it Innsbruck, 1984
12. Post, E.L.: Generalized differentiation. Trans. Am. Math. Soc. 32, 723-78t (1930)
13. Tustin, A.: A method of analysing the behaviour of linear systems in terms of time series. J.
lEE 94, 130-142 (1947)
14. Widder, D.V.: The Laplace transform. Princeton: Princeton University Press 1941
15. Zlfimal, M.: Finite element multistep discretizations of parabolic boundary value problems. Math.
Comput. 29, 350-359 (1975)

Received September 5, 1985/September 20, 1987