M.......

»'-,‘—.

7-... —- ya u.-

_..—..g

huh-usxvnu

..ulvlnlrufhell

H

< ‘

NUMERICAL.

umromw ACCURATE

EQUATIONS A

USING EXTRAPOLATiON VAND

v

ERENTIAL

INTERzPOLATION

SOLUTIONS T0 DIFF

 

 

 

                  

t .
a. . A
, «emf;
( 11.91.35
$3.15: . Int}:
0 w"
.

 

 

 

 

 

 

.‘1f12l3uwf
5,113.. $50.! 1

$31.31 x .:
.5 £34. 0......

  

 

 

 
 

» ,vffiflvvuwr. m ..

 

   
     

  

LIBRARY

Michiga :1 State.
University

This is to certify that the

thesis entitled

"Uniformly Accurate NUmerical Solutions
to Differential Equations Using
Extrapolation and Interpolation"

presented by

Richard A. Rogers

has been accepted towards fulfillment
of the requirements for

_EhJ._Jegee in Mathematics

Gerald D. Taylor/
Mary J. Winter

Major professor

 

Date JUJJL 25, 1974 fluflJ-‘D. 29%., A

04639

E? y
BINDING a! .

HMS &; '
800K BINDE INC.

LIBRARY amotns
SPHIIBPOflI

 

 

UNIFORM
m DIFFERENT

 

In this v
for solving ordir
those methods the
all powers of hC
afixed integer.

with such methods

grid points be 101

we develop the u]

conb'm
es extr
apes

2c ‘ -
effluent func-

Prod -
uce a highly
”est mesh

lie a -

Zesh,

ABSTRACT

UNIFORMLY ACCURATE NUMERICAL SOLUTIONS
TO DIFFERENTIAL EQUATIONS USING EXTRAPOLATION
AND INTERPOLATION

BY

Richard Allan Rogers

In this work we are concerned with numerical methods
for solving ordinary differential equations. We consider
those methods that have asymptotic error expansions involving
all powers of hq, where h is the steplength and q is
a fixed integer. The process of extrapolation can be employed
with such methods to obtain highly accurate solutions at
grid points belonging to the coarsest mesh. In Chapter I
we develop the "pullback interpolation method". This method
combines extrapolation with Hermite interpolation of the
coefficient functions for the asymptotic error expansion to
produce a highly accurate solution at all grid points of the
finest mesh. When q is l or 2 pullback interpolation
yields uniform accuracy at all grid points of the finest

mesh.

In Chapter II the pullback interpolation method
is modified so as to be applicable to boundary value
prdblems. In addition, an elementary proof of the stability
Of V. Pereyra's finite difference scheme for solving two

Point boundary value prdblems is given.

 

.00
1/

00

(i In (Shop
eoaations with c
are shown to be
Because of the 1.
acy obtained thr

for these proble

 

 

 

A finite
second order dele

in Chapter III .

save an asymptot

StePlength.

0 Richard Allan Rogers

'l/
00

G In Chapter III we consider difference differential

equations with constant retardation. The methods of Chapter I

are shown to be applicable to first order delay equations.

Because of the presence of the delay term, the uniform accur-

acy obtained through pullback interpolation is indispensible

for these problems.

A finite difference scheme for directly solving

second order delay equations is constructed and analyzed

in Chapter III. The global discretization error is shown to

have an asymptotic error expansion in even powers of the

s teplength.

 

 

UNI FORM
TO DIFFEREN’

 

 

 

n Partial

 

Mic]

UNIFORMLY ACCURATE NUMERICAL SOLUTIONS
TO DIFFERENTIAL EQUATIONS USING EXTRAPOLATION
AND INTERPOLATION

BY

Richard Allan Rogers

A THESIS

Submitted to
Michigan State University
in partial fulfullment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mathematics

1974

TO MY WIFE
KAREL

ii

' an"!

 

I would 1
3:. Gerald D. Ta;
and encouragement

toDr. Lee M. Son

aivice and guidan

 

 

 

 

 

ACKNOWLEDGEMENTS

I would like to thank my major professors,
Dr. Gerald D. Taylor and Dr. Mary J. Winter for their help
and encouragement. Also, I want to extend my appreciation
to Dr. Lee M. Sonneborn and Dr. Gerald D. Taylor for their

advice and guidance over the course of my graduate career.

iii

 

INTROD'

CERER I: THE p
INITI

Section 1

Section 2

Section

LA)

Section 1i
5

SECtiOn

33% .ER II: mg ‘l

Section 1

Section

SeCtiOn *
CiimR 111‘ Tm
DI;
RE?
SectiOn .
section ,
SeetiOn ;
Section 1

Section ;

 

Section ;

sectiOn .

l

l
l

TABLE OF CONTENTS

INTRODUCTION . . . . . . . . . . . . . . . 1

CHAPTER I: THE PULLBACK INTERPOLATION METHOD FOR

INITIAL VALUE PROBLEMS . . . . . . . . . 10
Section 1: Statement of the Problem . . . . 10
Section 2: The Pullback Interpolation

Method . . . . . . . . . . . . 16
Section 3: The Pullback Method for q=1

and 2 o o o o o o o o o o o o o 27
Section 4: Determination of e$(a) . . . . . 31
Section 5: Numerical Results for Initial

Value Problems . . . . . . . . . 56

CHAPTER II: TWO POINT BOUNDARY VALUE PROBLEMS . . . . 77

Section 1: The Problem and its

Discretization . . . . . . . . . 77
Section 2: Stability . . . . . . . . . . . . 84
Section 3: The Numerical Method . . . . . . 90

CHAPTER III: THE NUMERICAL SOLUTION OF DIFFERENCE
DIFFERENTIAL EQUATIONS WITH CONSTANT

RETARDATION . . . . . . . . . . . . . . 95
Section 1: First Order Equations . . . . . . 95
Section 2: Second Order Equations . . . . . 108
Section 2.1: Consistency . . . . . . . . . . 112
Section 2. 2: Stability . . . . . . . . . . . 121
Section 2. 3: The Global Discretization

Error . . . . . . . . . . . . . 132
Section 2.4: Implementing the Second
Order Method . . . . . . . . . 141
Section 2.5: Numerical Results for Second
Order Equations . . . . . . . . 147
BIBLIOGRAPHY . . . . . . . . . . . . . . . 152

iv

Table
Table
Table
Table
Table
Table
Table
Table
Table
Table
Table
Table
Table
Table
Table

Table

10 .

ll
12
13

14
15
16

LIST OF TABLES

l4
14
45
52
63
65
67
68
7O
71
72
74
105

107
148

150

 

 

In this
:iqtes for obtai
initial value pr
solutions by men:
nllalso consid‘

computed solutio:

 

Consider

x'(t)

N

(1)

\¢

x(a)

H

 

INTRODUCTION

In this section we will discuss some standard tech-
niques for obtaining numerical solutions to first order
initial value problems. The process of refining computed
solutions by means of extrapolation will be explained. We
will also consider the question of uniform accuracy of the

computed solution.

Consider the first order initial value prdblem

X'(t)

f(t,x(t) ),
(1)

x(a) a. a<tib.

We will assume that (1) has a unique solution w(t) which
depends continuously on the initial condition x(a) = a.
Conditions on f(t,x(t)) which will guarantee this are well

known and are given in Chapter I.

The majority of procedures for solving (1) numerically
are based on discretization. As such, they yield an approx—
imate solution to (1) on a discrete point set contained in
[a,b]. we will only consider the case where the discrete
points are equally spaced in [a,b]. In this case the discrete

point set can be conveniently represented as a grid

G = [tn=a+nh:n=O,l, ...,N},

where the parameter h = b§a and is called the steplength.

 

 

 

The thee
tit) and the a?

for
b‘. X(tn'h)

A compUt
viich takes the
Xitmj'h) and
to be a linear

is referred to a.

Three of
rule (sometimes (

rule and Gragg ' s

Euler's 1

l2)

50‘: t‘
.anon (2) exhi

‘5 Obtained given

0 ution at

 

The theoretical solution to (1) will be denoted by
m(t) and the approximate numerical solution will be denoted

by X(tn,h) for each grid point tn.

A computational method for determining X(tn+k'h)
which takes the form of a linear relationship between
X(tn+j'h) and f(tn+j’x(tn+j'h))' j=O,l,...,k 18 said

to be a linear k-step method. The class of all such methods

is referred to as the class of linear multistep methods.

Three of the simpler multistep methods are Euler's
rule (sometimes called the Euler—Cauchy method), the trapezoid

rule and Gragg's modified midpoint rule.
Euler's rule is summarized by the equation

(2) X(t h) = X(tn,h)+hf(tn,X(tn,h)), n=O,l,...,N—l.

n+l'

Equation (2) exhibits how the approximate solution at tn+l

is obtained given that one has already obtained an approx-
imate solution at tn. Since information is required only

at the preceeding grid point in order to obtain an approx-
imate solution at the next grid point, Euler's rule is

a one-step method. To initialize or start the method requires
one piece of information which is given by the initial con-
dition in (1). That is, take X(to,h) = x(a) = a. Since
everything on the right hand side of (2) is known when we

are trying to compute X(t h), the relationship (2) is

n+l'

 

 

 

m

explicit and Euler's rule is an explicit multistep method.
The trapezoid rule is given by

(3) X(t h) =X(tn,h) +%[f(tn,x(tn,h)) + f(t (t h))],

n+1' n+l'X n+l’

n=O,l,...,N-l.

This is again a one-step method and is initialized by taking
X(t0,h) = x(a) = a. However, the right hand side of (3) is
not completely known. Indeed, when trying to compute the sol-

ution at tn+1 the term f(tn+l'x(t h)) is unknown

n+1'

because it involves the solution we are trying to compute.
Whenever f(t,x(t)) is linear in x(t), and in some other
special cases, equation (3) can be solved explicitly for

X(t h). However, in general this is not possible and the

n+l'

trapezoid rule is said to be an implicit multistep method
because of this behavior. When (3) cannot be solved explic-

itly for X(t h) a root finding procedure or functional

n+l'

iteration can be used to find the solution.

The midpoint rule is an explicit two-step method

defined by

(4) X(tn 2.h) =X(tn,h) +2hf(t (t h)), n=o,1,...,N-2.

+ n+l'X n+l'

To compute X(tn 2,h) we need to know both X(tn,h) and

+

X(t h). Thus, the midpoint rule requires two starting

n+l'
values X(t0,h) = a and X(tl,h). The second of these can

be determined in a variety of ways.

f.\ Flt.
m... .o.. a It a do tr
r . . .1. -l/ \l m... .l hu
~ 3 .3 13 Kl r 0 l: e
_ a: O z 0 at X.
N: c.~ Ce ts ”PW n l. .3 x
T. c .

 

 

 

 

Gragg's modification of the midpoint rule [11,12]
is Uwofold. First, to obtain the additional starting value

use Euler's rule. That is,
X(t1,h) = d-+hf(a,a).

Second, at the grid point tN==b use a smoothing procedure
to make the computed solution more stable. The smoothing
procedure consists of averaging three computed solutions

and is similar to a device originally employed by Milne and
Reynolds [20,21]. Gragg's modified midpoint rule is formally

given by
X(tlflfl =c1+5hf(a,a):

(5) X(t h) =X(tn,h) +2hf(t (t h)), n=o,1,...,N-1:

n+2' n+l’X n+l'

l l

X(b,h) =71IX (t h) +-2-x (tN,h) +Zx (t h).

N-l' N+l’

The three numerical methods given by equations (2)
(3) and (5) have an important similarity. Each has an
asymptotic error expansion that involves all powers of hq
for a fixed integer q. That is, for Euler's rule, the trap—
ezoid rule, and Gragg's modified midpoint rule, the computed
solution satisfies
(6) X(t ,h) = co(t ) + 2e (t )hqk.

n n k=1 k n

The functions ek(t) are independent of h and q is a

fixed integer. For Euler's rule q=l and for the trapezoid

yu‘
OI.

 

MHz-i
"sb'u
4
.
b :4
Q ‘~
!
z m .
, . vat.

.1_

 

 

 

nl‘m

Hang
\

rule and Gragg's modified midpoint rule q=2. The expansion
(6) for Gragg's modified midpoint rule is only valid at

tN==b while for the other two methods (6) is valid at all
grid points tn. The expansion (6) is valid only when

f ec°°[[a,b] x (-oo,co)].

If f has only a finite number of continuous deri-
vatives, a truncated version of (6) is valid, namely

M
. _ qk q(M+l)
(6 ) X(tn.h) — CNtn) +k§lek(tn)h + 0(h ).

The notation 0(hj) means that the function being suppressed

behaves like a constant multiplied by h3 as h-eo. Formally,

we say that g(t) = 0(h3) if -L1£%LL §;C as h-+O, where

h
C is a constant. The length, M, of the expansion (6')
depends on q and the number of continuous derivatives of f

that exist.

The existence of asymptotic expansions of the form
(6) or (6') for Euler's rule, the trapezoid rule and Gragg's
modified midpoint rule was originally proved by Gragg [11,12].
Stetter [31] and Pereyra [24] have also studied the existence
of such expansions. Gragg's results for Euler's rule and the

trapezoid rule are presented in more detail in Chapter I.

The existence of such expansions is important because
numerical methods which have error expansions of the form

(6) or (6') are amenable to extrapolation. Basically, the

 

 

prscl

rn~r~
sva
A

U
‘n
fi
3 u
,var
olsvt
..
..

 

 

c!
v: CE

process of extrapolation is a means of combining several
computed solutions, each of low accuracy, in such a manner

as to obtain a computed solution with high accuracy.

Extrapolation is not a recent development, dating
back to at least 1654 when Huygens [14] used it to improve
Archimedes polygonal approximation to v. Extrapolation was
first systematically studied by Richardson [29] early in this
century and has often been referred to as either Richardson
extrapolation or the deferred approach to the limit, the
latter being the title of Richardson's 1927 paper. An
excellent survey article on extrapolation and its applications

has been written by Joyce [l6].

Extrapolation,vflmn1applied to Obtaining an accurate
solution to a differential equation.i£5most easily explained
as follows: Let HZ>O be a fixed basic steplength and
suppose an accurate solution to (l) is desired at the point
a-FH. Define a sequence of steplengths hk==H/2k

k=O,l,...,K and grids

k _ . . ._ k
Gk — [ti — a+ihk. i—O,l,...,2 }.
All grids contain a and a-+H and the grids are nested
With Gk<:Gk+l Vk. On each grid Gk compute a numerical
solution to (1) using a method which has an asymptotic error
expansion of the form (6') with MT>K. At a-+H we have

I<+]. computed solutions X(tN,hk) for k=O,l,...,K.

 

 

 

“0-th1

F ‘
.A.‘ #hl
U'. 5..

rf
,1
(D
'11

:o'c:

us: 5

sm-

«gc

Extrapolation is the process of forming a linear combination
of these K-kl solutions in such a manner so as to eliminate
the first K error terms of the expansion (6'). That is,
constants ck are determined so that

K
(7) '2': c X(a+H. ) = CPO?) +o(Hq(K+1’).
k=0 k h‘k

Aitken [1.] and Neville [22] independently devised
an iterative scheme by means of which extrapolation can be

performed without explicit computation of the constants ck

in (7). The convergence of this iterative scheme under

suitable hypotheses was established by Gragg [11,12].

(K+1)

Note that in order to obtain 0(Hq ) accuracy

at a grid point you must have K+l computed solutions avail—

(K+1)

able to work with. Thus, extrapolation will yield 0(Hq )

accuracy only at the point a+H which is common to all grids

Gk' Extrapolation can be performed at other grid points.

Hq(K+l))

However, it will not yield 0( accuracy at these

points. For instance, using the fact that GkCIGk+1 Vk.
extrapolation will yield 0(HqK) accuracy at the midpoint
H

a-eiu At other grid points, extrapolation will yield even

less accuracy.

Lindberg [19] has developed a method based on
extrapolation and Lagrange interpolation which can be used to
Obtain 0(Hq(K+l)-1) accuracy all at grid points of the finest

grid G

K

 

 

T5 nfl'

“gin";
' F 's .
l.'. .e,
cyp-
nae

Utes.
Q

In Chapter I we present what we have termed "the
pullback interpolation method". It utilizes extrapolation

and Hermite interpolation to obtain 0(Hq(K+1))

accuracy

at all points of the finest grid when q=l or 2. The first
four sections of Chapter I are devoted to develOping the method
and establishing a theoretical basis for it. In Section 5
Lindberg's method is presented and compared with pullback'
interpolation both theoretically and numerically. Extensive

numerical tests were performed and these results are also

presented in Section 5 of Chapter I.

The focus in Chapter I is entirely on first order
ordinary differential equations. In Chapter II we consider

two point boundary value problems of the form

x"(t) = f(t,x(t).X'(t)).
(8)

x(a) = A, x(b) = B.

Pereyra [24.25.26.27,28] has developed a finite difference
scheme which yields an approximate solution to (8) that

has an asymptotic error expansion of the form (6') with q=2.

Pereyra's results are summarized in Chapter II and
pullback interpolation is modified so as to be applicable to
boundary value problems of the form (8). In addition, in

Section 2 of Chapter II, we present a new proof of the

D.

AU
t

 

stability of Pereyra's difference scheme. In comparison

to Pereyra's proof.cnnx;is considerably more elementary.

In Chapter III we consider the numerical solution
of difference differential equations with constant retard-

ation. First order equations

x(t) = f(t.X(t).X(t-r))

are analyzed in Section 1 and pullback interpolation is
shown to be a viable solution technique for these prdblems.

Numerical results are presented to support this contention.

The remainder of Chapter III is devoted to the
develOpment and analysis of a finite difference method for

directly solving second order equations of the form

(9) x(t)-+f(t,x(t),x(t-r),x(t),x(t-r)) = o.

The approximate solution to (9) is shown to have an asymptotic
expansion of the form (6') with q=2. A modification of

pullback interpolation is shown to be applicable.

"“L

Vi.
§
V
I

 

\}

I.\

Ch

Au..-

.3"

CHAPTER I

THE PULLBACK INTERPOLATION METHOD
FOR INITIAL VALUE PROBLEMS

Section 1. Statement of the Problem

 

In this chapter we consider the first order initial

value problem.

Y'(t) f(t.y(t))

(l)
y(a) = a aétgb.

We shall assume that f(t,y(t)) is a continuous function
of t and satisfies a uniform Lipschitz condition with
respect to its second argument. Under these assumptions
it is well known (see Keller [17]) that (1) has a unique
solution, m(t), which depends Lipschitz continuously on
the initial data y(a) = a. y(t) may be either a scalar-
valued or a vector-valued function. If y(t) is vector—
valued then f(t,y(t)) will be a vector—valued function
of the variable t and the vector y(t) and (1) will be
a system of first order differential equations. This case
also arises when we reduce a mth order differential
equation to a system of m first order differential
equations. The standard technique for accomplishing this
can be found in Lambert [18]- The numerical methods to

be considered for solving (1) will work for either the

10

 

‘ .I
h“ (”.H
.4 writ/.4

refine
CH a i:

COIiis

ll

scalar or vector-valued case. In the vector-valued case
extrapolation and the pullback method, to be explained in
this chapter, can be done independently for each component
of the solution vector to obtain a refined solution. Since
the same work is done for each component, one at a time,
without using any results involving other components. the
refinement process appears well suited to implementation

on a parallel processing computer such as ILLIAC IV (see

Corliss [ 5]).

Turning our attention to the numerical solution of
(1), let h>O be a fixed basic steplength and for each
k=O,l,...,K define steplengths hk = h/2k and grids

k k . . k

Gk — [t..ti — a+1hk, 1—O,l,...2 }.
2k+l points: all grids contain the two points t3 = a and
k . .
t — a+h, and the grids are nested, that is, Gk CLGk+1Vk.

2k

Each grid Gk contains

In what follows we shall be concerned with methods
for the numerical solution of (l) on the grids Gk which
give an approximation, Y(t§,hk), to the theoretical sol-

. k .
ution, m(ti), such that the error has an asymptotic expan-

sion

(2) y(tk h ) - (tk) + :3: jqe (tk) + o( (“*1”)
1' k “9 i jzlhk j i hk

for each i=O,l,...,2k and for each k=O,l,...,K. The

coefficient error functions, endt)' are independent of

hk' for each k, and q is a fixed integer (usually 1 or 2)

0"

 

12

peculiar to the method being employed. The exact length

of the expansion (2) will depend on both the method being
employed and the number of derivatives of f(-,-) which
exist. In general, (M+l)q continuous derivatives of
f(t,y(t)) with respect to t are required for an expansion
with M error terms. This is equivalent to the theoretical

solution m(t) having (M+l)q+1 continuous derivatives.

As mentioned before, Gragg [11,12], Pereyra [24],
and Stetter [31] have investigated the existence of such
expansions. Examples of methods which yield such error
expansions are Euler's method (q=l). the usual generalization
of the trapezoid rule (q=2) and Gragg's modified midpoint
rule (q=2). An important result obtained in each of the
above mentioned studies is that the coefficient error
functions em(t) satisfy an inhomogeneous linear variational
equation on agtgb, of the form

I

em(t) — J(t)em(t)

ll
U0
8

(3)

II
.0
3
II
H
3

em(a)

The arguments of the inhomogeneous terms, bm(-), involve

the theoretical solution m(t), previous error functions
e1,...,em_1, the function f(t,y(t)), and various derivatives
thereof. The differentiability of the various error

functions depends on the differentiability of f(t,y(t)).

The left hand side of (3) is the Frechet derivative of (l),

 

 

 

 

grid

of
1

13

considered as a differential operator, operating on em(t).
Alternately, J(t) is the Jacobian of f(t,y(t)) evaluated
at the theoretical solution, ¢(t), of (l). The left hand
side of (3) may be obtained formally by assuming y(t)
depends on a parameter 1, differentiating (l) with respect

I O
to this parameter and setting em = %%f and em =

If we compute the numerical solution of (l) on the
grids Gk k=O,l,...,M using a method which has an expansion
of the form (2), we may then employ extrapolation to obtain
a solution Y(a+h) = wia+h) + O(h(M+l)q). However we are
not able to obtain comparable accuracy at the intermediate
points. For instance, using extrapolation, the solution
at the midpoint satisfied Y(atg) = m(a+%) + 0(th). And,
as the following example illustrates, we cannot interpolate
several extrapolated values to obtain a solution with equi—

valent accuracy.

Example 1: y' = y2; y(0) = .2.- ogt_<_3.

The theoretical solution to this problem is
®(t) = §%?-. Using the trapezoid rule and extrapolation
with h=1, M:3 we compute the solution, Y(t), at the
points t=l,2, and 3 by resolving the same problem three
times with initial conditions determined by the computed
solution at t=1 and 2. The results are given in

Table 1 below.

p
or"-
.‘v‘

Per:
the
t
O
l
2

t. li\2 .«a\.-‘ .\J\../\

 

 
 

14

A variety of interpolation schemes are possible
for determining the solution at intermediate grid pohits.
Two such are summarized in Table 2 below. L(t) is the
Lagrange cubic interpolation polynomial for the data
(t,Y(t)) given in Table l and H(t) is the quartic

Hermite interpolation polynomial for the same data with

the added condition Y'(O) = (Y(0))2 = (.2)2 = .04.
TABLE 1
t cp (t) Y(t) Wt) - Y(t)
0 .200 000 000 .200 000 000 0
1 .250 000 000 .250 000 000 0
2 .333 333 333 .333 333 330 3 x 10"9
3 .5000 000 000 .499 999 762 2.38 x 10'7
TABLE 2
t. wlt) L(t) wlt)-L(t) H(t) wtt)-H(t)
%-.222222222 .223958320 -1.736x10’3 .222395831 —1.736x10-4
%-.285714285 .284375013 1.339x10'3 .285312506 4.018x10'4
§-.400000000 .403124923 --3.125x10"3 .401562434 —1.562x10"3

As an examination of Tables 1 and 2 quickly
reveals, neither of the interpolation schemes gives accuracy

that is comparable to that of the computed solution.

15

In the next section we will present and analyze
the pullback interpolation method. It is based on a .
Hermite interpolation scheme for approximating the error
functions em(t) successively, beginning with the last
term of the error expansion. The details are worked out
for a general expansion of the form (2). However in
practice q is usually either one or two and the reader

may find it helpful to bear this fact in mind.

A v
A»..
n\u

ta
\Au

,|x

 

[“0

16

Section 2. The Pullback Interpolation Method
Still assuming that our numerical method and the
problem at hand are such that the expansion (2) is valid,
we now take K = M. At the point t = a+h<EGO we have the
M+l computed solutions and error expansions
M

(4) Y(t.hk) = w(t) + Zihgqej(t) + 0th;
3:

M+l)q), k=O,l,...,M.

Equation (4) can be regarded as a system of M+l poly-
nomial equations in h. Since the hk's are distinct these
polynomials are linearly independent and (4) may be solved
as a matrix system for the unknown vector

(w(t).e1(t).e2(t)....,eM(t))T.

In matrix form (4) is given by

(4') AU = Y + E
where
C q 2g (M-1)q Mqu
1, ho. ho , . . ., ho , ho
q 2q (M-llq Mq
A = 1! hll hl I ' . ° ' hl I hl
q 2q (M-l)q Mq
Ll'lhr hM , .. . hM , hM _J

 

 

is (M+l) x (M+l); and U,Y and E are (M+1) X l

vectors given by

 

 

5P}
,_“\

l7

 

 

 

 

 

 

r 1 r '
rCp(t)‘i Y(t,h0) 0(héM+l)qyi
el(t)l Y(t,hl) O(h£M+l)q)
U = . : Y = . : and E =
respectively.

Define (A,I) to be the (M+l) x 2(M+l) matrix
obtained by adjoining the identity matrix to A. A
standard theorem from linear algebra (see Cullen [ 6])
tells us that if (A,I) is row equivalent to (1,3) ——

denoted by (A,I) ~ (1.8) —— then B = A-l.

Since A is nonsingular we can find a sequence
of elementary row Operations such that (A,I) ~'(D,C)
where D is the diagonal matrix whose (j,j) entry is
h‘j-1)q, j=l,...,M+l. Moreover, since V k hk = h/2k.
the row operations which accomplish this reduction are
independent of h, in that each involves multiplication
or division by a constant only. (D,C) can be further

1) by dividing the jth row by h(3-1)q.

reduced to (I,A-
Multiplying both sides of (4') by A.1 we have

1 l

1 — . .
Y + A E. The vector U is the prec1se

U = A5 AU = A’
solution to the system (4) and the vector A-lY is the
numerical solution we can obtain. The vector A-lE is

the vector of errors for the numerical solution.

18

Now vk hkgh, so we can write

 

 

 

 

)- ‘1 r' a
0(héM+l)q) 0(h(M+l)q)
0(h{“”’q’ 0(h(M+1)q)
E = =
0(héM+l)q 0(h(M+l)q)
L b _ A
Since each component of the jth row of A-1 has a factor

of h-(j-l)q and this is the only dependence of these
entries on h, the jth component of A—lE is
0(h(M+l)q-(j-l)q) = 0(h(M7j+2)q). Thus our computed
solution 25% = (S(t),'él(t),é'2(t)....,'é’M(t))T will be
equal to the precise solution U = (cp(t),el(t),e2(t),...,eM(t))T

with an error of magnitude
A-lE = (0(b‘”+l’q).o(h"q).0(h‘M‘1’q)....,0(bq))T.

It should be emphasized that what we actually
Obtain when computing a: solution to (4) is the vector
A‘lY which is only an approximation to the actual solution

U. This can be conveniently summarized by saying the

solution to (4), (m(t),el(t),...,eM(t))T, is known with
accuracy (0(h(M+l)q),O(th),...,O(hq))T.

Since eM(a) = O and m(a) is known exactly:
we can Obtain efi(a) exactly from (3). That is,

eM(a) = 0 implies that efi(a) = bM(a). To evaluate

‘ ‘2} 71“ I

 

Ma) requi
Pi

:lt). J(t) .
evaluated at
to us. Thee
iifferential
details wilj

this chapte]

Thus
know eMiaH
interpolatir
Pieces of da
Hermite inte
known to 9
37°“ vk=1...

any positiw

l9

bM(a) requires knowledge of higher order derivatives of
m(t), J(t), f(t,m(t)) and the errors em(t), m=1,...,M—l,
evaluated at the point t=a, all of which are available

to us. These derivatives may be obtained by successively
differentiating (l) and (3) with respect to t. The
details will be worked out for specific methods later in

this chapter.

Thus we know eM(a) and efi(a) exactly and we
know eM(a+h) to 0(hq). Construct PM(t) the Hermite
interpolating polynomial of degree 2 to the above three
pieces of data. It is well known that the error in this
Hermite interpolation is 0(h3). Thus, since eM(a+h) is
known to 003) we have PM(s) = eM(s) + 009) + oniq).
NOW"Vk=1,...,M, hk<h implies hk<(hy' where y is

any positive integer and therefore O(h£)::0(hy). Thus

Vk=l,...,M and Vs 6[a,a+h] we have
Mq Mq BM
(5) hk PM(S) = hk eM(S) + Oih )
where BM 5 min((M+l)q, Mq+3). Note that for qg3
BM = (M+1)q.

Proceeding to the grid G1. let t be any

point in G1\GO (since GO contains 20-kl = 2 points

and G contains 21-+l = 3 points, there is only one

1
such t, namely t = a-kh/Z). Since the grids are nested
this: t.€Gk 'Vk=l,...,M and we have the M computed

solutions and error expansions

 

(6) mm]

Substituting
0f M lineal

unknown vectc

(7) Y(t'hk)

20

M—l .
(6) Y(t.hk) = o(t) + . hflqej(t) + hfiqu(t) + OlhfiM+l)q)

J=1

for k=l,...,M.

Substituting (5) into (6), we obtain the following system

of M linearly independent equations in h for the

unknown vector (cp(t),el(t),...,eM_l(t))T
Mq “'1 'q (n+1) BM
(7) Y(t’hk) - hk PM(t) = co(t) + 2 hi ej(t)+ 0011‘ q) + 001 )
i=1
M-l . B
= co(t) + .231 hflqejfi) + 0(h M),
J:
k=l,...,Mr
because 5M S_(M+l)q and hk < h.
Solving (7) we obtain the solution vector
(cp(t),e1(t),...,eM_l(t))T with accuracy
6 B - 6 -2q 6 -(M-l)q
(oth M).O(hM ). om“ )..... om” ))T

We now have the following information about eM_1(s).

For t = a+h EGO we know eM-l(t) with accuracy 0(h2q)
and for t = a+h/2 6G1\GO we know eM_l(t) with accuracy

(6 -(M-1)q
0(h M )

As before eM_1(a) = O and we can deter-
mine eM_1(a) exactly. Thus we can construct PMF1(S)
the Hermite interpolation polynomial of degree three
interpolating these four data points. Since the error

of interpolation will be 0(h4) and hk<h we have

vs €[a,a+h],

 

lj

 

 

 

 

~')
(9) h}?! J qu.
The Humber El
(10) Emj 5]
and EM ls de
Let t
nJ+l Conta

21

(e) héM-l)qPM-1(S) = héM—l)qu-l(s) + 0(h(M-l)q+4)
B
+ o(h‘“+1’q) + o(h M)
(M-l)q BM-l
= hk eM-l(S) + 0(h )I

where BM—l min (BM,(M—l)p+4). Here we have used the

fact that BM«£ (M+1)q. We note that for q 3.2,

BMel = (M+1)q.
Proceeding by induction, suppose PM(s), PM_1(s).
...,PM_J(S) for J < M—l have all been constructed. where

Vj=0,l,...,J, PM_j(s) is a Hermite interpolation poly—

nomial of degree 23+l with
(M-j)q _ (M-j)q 6M—3
(9) hk PMej(S) — hk eM_j(s) + 0(h ) vs €[a,a+h].

The numbers BM-j are defined by

(M-j)q + 23 + 2), j=1,...,J,

and BM 15 defined by (5). Note that BM-j S-BM-j+l V3.
Let t be any point in GJ+l \ GJ. Since

GJ+1 contains 2J+1 + 1 points and GJ contains 2J + 1

points there are 2J+1 — 2J = 2J such points. Since

J+l \GJ 18 also an element of

Gk for k=J+2,J+3,...,M. Thus at each such t we have

GRCGk-i-l Vk each tEG

the following system of M-J linearly independent equations

in h:

 

(ll) Y(t.hk)

 

since 41 ,3 E
. '“M-J _-

the System (11;

we Obtain the 1

i l
on M a

M‘Jel it) at
"it
s.)

22

(11) Y(t ) %; (M'j)qp (t)
'hk . hk Maj
3:0
M-J-l . J B .
= om + Z 151%. (t) +o(b1£M+l’q) + 2001 M'3)
':l 3 '=0
3 J
M-J-l . B
= w(t)-+ Z) hflqe-(t)+-0(h M’J). k=J+l,J+2,...,M:
i=1 3
since BM—J g-BM-J+1 g_... g_BM_1 S-BM-S (M+1)q. Solving
the system (11) for the unknown (w(t),el(t),...,eM_J_1(t))T,
we obtain the solution with accuracy
6 _ B _ -q B _ -2q 6 -(M—J-l)q
meJLomMJ).MhMJ LnnomMJ H?
From solv1ng (11) at each p01nt 1n GJ+1 \ GJ we
obtain knowledge of eM_J_1(t) at 2J points with
BM_J—(M-J-l)q
accuracy 0(h ). From our work on G we

0

know eM_J_l(a+h) with accuracy 0(h(M+l)q-(MPJ-l)q)

J+2)q).

OKh< From our work on G1 \ GO we know

BM-(M-J—1)q
eM_J_l(a+h/2) with accuracy 0(h ). By our

inductive hypothesis Vj < J, we have obtained knowledge

j . . .
of eM—J—1(t) at 2 p01nts 1n GJ+1\Gj *w1th accuracy

BMFJ-(MeJ-l)q
0(h ). Thus we know
J j J+1 J+1 .
e (t) at l + 232 = l + 2 — l = 2 points
M—J—l ._
3—0
with varying accuracies. In addition we also have

eM_J_1(a) = O and we can determine eM_J_1(a) exactly.

 

 

and we know 6

Thus we have k

interpolation I
above data. T

J+1
5012 +2) an.

(12) héM-J-lh

where .52.
P M- J— l

the error in ir

 

5(..iM‘J‘1)q+2J*

mint-'nEM-J-l) qn

since 3
M~J i (
Hence b
PQl
l ',M‘l
G]
' and (1
O) a
r
On the
Exbres
. Sign
:59 the CO

23

1

Thus we have knowledge of (t) at 2J+ +1 points

eM-J-l

and we know eM_J_l(a). Construct PM_J_1(s) the Hermite
J+1

interpolation polynomial of degree 2 +1 based on the

above data. The error of interpolation will be

J+1
0(h2 +2) and Vs €[a,a+h],
(M-J—l)q _ (M—J-1)q 6were
(12) bk PM_J_l(s) — hk eM_J_1(S) + o(h ).
where aM—J—l is determined in the following manner:

J+1
the error in interpolation is héM-J-1)q0(h2 +2)

1+

2).

J+
OKh(M-J-l)q+2 2) and the error in the given data is
_ _ B _.-(M-J—l)q _ _
min({hfin J l)qo( M 3 )}g=0' héM J l)qo(h(J+2)q))
B .
= min ((oni M'Jngzo. 0(h(M+1)q))
B
= min (0(h M'J). 0(h(M+l)q))*
B
= 0(h M-J)
since B < (M+1)q Thus B = min (b (M—J—1)q+2J+l+
M-J — ' M—J—l M—J'

Hence by induction we can determine PM_j(s) for

j=O,l,...,M—l such that PMej(s) has degree 23—1 and

(9) and (10) are valid for j=O,l,...,M—l.

On the last grid G we have the following

M

expression for m(t) for each t.€GM\GM_l.

*see the comment following (11)

 

[13) Y(t,hM

24

_. (M- DC:
(13) Y(t.hM) jZ):hM PM_ j(t)
M-l B .
= cp(t) + o<hb§M+1’q) + 20m 1"‘3)
i=0

61
= «J(t) + OH)

since B1£P233"SPMS(M+1)Q°

To obtain the final solution, Y(t), at all grid

points in G the finest grid, we proceed as follows.

MI
For t = a+h, Y(t) is the first component, g(t), of

the computed solution to the system (4). For

t<EGM \ GM-l’ Y(t) is the solution of (13). For any
M
t'EGMpl' Since kLng = GM and the Gk s are nested,

there exists an index J, depending on t, such that
t EGJ+1 \GJ. The solution at this t, Y(t), will be
the first component, E(t), of the computed solution to

the system (11).

Since 61 is the smallest of the Bj's, the

solution Y(t) constructed above satisfies
31
(14) Y(t) = cp(t) + O(h ).

Before we examine the Bj's in detail we should
point out that the pullback method does not require use
of the iterative scheme of Aitken and Neville to perform

extrapolation. In the pullback method extrapolation is

a
. . “.1

 

actually ac

systems.

It
for the sys
matrix for
matrix of 1
the matrix

PIOgreSSiv.

 

 

25

actually accomplished when we solve the various matrix

systems.

It should also be pointed out that the matrices
for the systems (11) are "nested" in the sense that the
matrix for the Jth system is an easily obtained sUb-
matrix of the (J-l)st system. Thus we need only define
the matrix for (4) and the others can be obtained by

progressively deleting rows and columns of this matrix.

Specifically, recall that the matrix for the

system (4) - the case j=O - is given by

’ q 2q (M-l)q Mq‘
1, ho, hO , . . . . ., hO , hO
q 2q (M—l)q Mq
l, h1' h1 I o o o o o, hl , hl
(15)
. s 2q (M-1)q Mq
_1, 15f hM , . . . .., hM , hM_J

 

 

The matrix for the system (7) - the case j=l - is given by

" q 2q (M-1)qq
1' h1, hl ' o o o o o, hl i
q 2q (M—1)q
1' h2' h2 I 0 o o o o, h2
. q 2q (M—1)q
_1, 15V hM , . . . . H hM J

 

 

This matrix can be obtained from (15) by deleting the

first row and the last column of (15).

 

For 96

is given by

F1, 1
1. l
l, 1

 

MMS and the

COnSe<

itis Only he<

26

For general J the matrix of the system (11)

is given by

i
F q 2q (M—J-1)q
1' hJ+l' hJ+l' . . . .. hJ+1
q 2q (M-J-llq
1' hJ+2' hJ+2' ' ° ° " hJ+2
q 2q (M—J—l)q
L_l, hM, hM , , hM 3

 

 

which can easily be obtained by deleting the first J+1

rows and the last J+1 columns of the matrix (15).

Consequently, when using the pullback technique

it is only necessary to define one matrix.

For large M the Vandermonde matrix (15) is
known to be ill-conditioned. BjOrck and Pereyra [2 ]
have developed an efficient algorithm for solving such
Vandermonde systems, which can be utilized for solving our

system.

Turning our attention to the Bj's, we shall
examine the cases q=l and q=2 in detail in the next

section.

 

Section

lerra l .

q:
PLO;
Hence c

with

the v;

this

(16)

and

27

Section 3. The Pullback Method for q=1 and 2

we will need the following rather obvious lemma.

1

Lemma 1. a) Let f(x) -2x + 2X+ then Vx€[1,¢), f(x)22.

l

b) Let g(x) —x + 2X+ then Vxe[l,m), g(x)22.

Proof. a) f'(x) = -2 + 2X+l

1n2>-2+2"20 for x21.
Hence on [1,«) f(x) is a monotone increasing function

with f(1) = 2 and therefore f(x).2 2.

b) g(x) = f(x) + x 2_f(x) 2_2 by part a) since

x is positive. Cl

Our previous inductive argument has established
the validity of equation (10) for j=1,...,Me1. From

this and (5) we can conclude for q=2 that

BM = min((M+l)2,2M+3) = (M+l)2
(16)

= min(B (M-j)2 + 23 + 2), j=l,...,M—l:

BM-j M-j+1'

and for q=1 that

BM = min(M+1,M+3) = M~+ 1
(17)

= min(B Mrj + 2j + 2). j=l,...,M~l.

BM-l M-j+l'

Theorem 1. For q=1 and 2 we have

BM-j = (M+1)ql 3:031, OOOIM-lo

 

Proof.

and (17)

(16) am

from 011

fOr

=1 OI

28

Proof. The proof is by induction on j. From (16)

and (17) we have BM = (M+l)q and for j=1 we have

BM—l min(BM,(M-l)q + 2 + 2)

min((M+l)q,(M-l)q + 2 + 2)

(M+1)q.

Assume BM-j = (M+l)q where l<j<M-1. Then by

(16) and (17)
BM—(j+1) = BM-j-i = min(BM_jv(M-j-1)q + 2j+1 + 2)

min((M+l)q,Mq-qj + 2J+1 - q + 2)

from our inductive assumption. By Lemma 1 -qj + 23+1.2 2

for q=1 or 2 and Obviously -q + 2.2 0. Thus for

q=1 or 2
(M+1lq S.Mq + 2 S.Mq - qj + 2J+1 - q + 2
which establishes that B = (M+1)q. C]

M - (j+l)

Theorem 1 shows that for the cases q=1 and 2
the pullback interpolation method yields the same accuracy
at every point of the finest grid as that Obtained by

extrapolation at the endpoint, a+h.

If M is large then the polynomials P P etc.;

0! 1'
will be of high degree and "Runge's phenomenon“ may destroy

the accuracy of our computed solution. This can be

 

circumvented

 

 

 

nomials. Exa
limitation on
to which we k
need only int
accuracy Comp
affect the 0v)

"Runge ' S Phen

BefOr
PM‘J' we need
Ma) 13 Siv

29

circumvented to some degree by using lower degree poly-
nomials. Examining the proof of Theorem l'we see that the
limitation on the accuracy we can Obtain is the accuracy

to which we know Thus for large M and j 'we

eM—j'
need only interpolate on enough points to guarantee

accuracy comparable to that of our data. This will not
affect the overall accuracy of our scheme and may avoid

"Runge's phenomenon".

Before we can actually compute the polynomials

P we need to know e' (a). The method fer determining

M—j M-j
e$(a) is given in the next section but at this point a
few more comments on the actual implementation of the

pullback method are in order

In the case where y(t) = (y1(t)...,y‘(t)) is
an fi-dimensional vector, the functions Y(t,hk). ¢(t)
and ej(t) are themselves L—dimensional. Consequently
the applicable system - (4), (7) or (11) - must be
solved L times. once for each component of the vectors.
Also a complete set of interpolation polynomials PM-j
must be computed for each component. We again wish to

point out the potential for using parallel processing

computers in this case.

There are cases where y is L—dimensional and

not all the above work is required. Specifically, suppose

 

the L—dimens
order differe
order differe.
one is intere.

equation and :

If th.
need utilize .
of the soluti<
0“ [0.13] wh.
the problem (j
equation on
as the initia;

may be repeat‘

In th:

at

30

the L—dimensional prOblem resulted from reducing an 1th
order differential equation to a system of 1 first
order differential equations. Furthermore, suppose that
one is interested only in the solution to the original

equation and not in its various derivatives.

If the solution is desired on [O,h] then one
need utilize the pullback technique on the first component
of the solution vector only. If the solution is desired
on [O,b] where b>h the standard procedure is to solve
the problem (1) on [O,h], then solve the same differential
equation on [h,h+h] using the computed solution at h
as the initial condition for the new prOblem. This procedure

may be repeated as many times as necessary to reach b.

In this latter case. the system (4) must be
solved 2 times since we need to know all the derivatives
at h accurately. However the other systems and the
polynomials PM—j need to be computed only for the first

components of the vectors.

 

 

M D)
In or

P .. of the

M-3

determine e1;1

is given in t

carried out f

In or
results, we w
are for the c
Plifies the a

valid fOr' an

ensional Case
I

31

Section 4. Determination of e$(a)
In order to construct the Hermite polynomials,
PM-j' of the previous section it is necessary to first

determine efiyj(a). The method for computing e$(a)
is given in this section. The actual computations are

carried out for m=1,...,4 and examples are given.

In order to make effecient use of Gragg's [11.12]
results, we will follow his notation. The examples given
are for the case when y is l-dimensional. This sim—
plifies the actual computations. However,1flmatheory is
valid for, and will be presented for, the general t-dim-

. l L
enSional case, y=(y ,...,y ).

Denote by J the Jacobian matrix of f, evaluated
at the theoretical solution w,

Bftt.gn(t)) a<t<b

J(t) BY _ _

and define symmetric k-linear Operators f(k)(t,m(t)),
agtgb,by
é . 0 2% akflgcoun ‘1 ‘k

(k) _
f (t:CP(t))Yl---Yk — ]_ ° =1 1'1 5‘ y]. °°°Yk
1 4k by ...ay

L

for k=2,3,...

fck)(t,cp(t))yl...yk is the kth Fréchet deri-
vative of f(t,m(t)) operating on the vectors yl...yk.

In the case where y is l-dimensional, J(t) is the

 

 

w
W’.

32

first partial with respect to y of f and f(k)(t,¢(t))yl..

(k)
_ .____(_1_CL(_L . ' k
_ B fay: t )yl...yk, that is, f( )(t,Cp(t))ylo--Yk

is the kth partial of f with respect to y multiplied
by Yl""’yk' The yj's are functions of t on which

f(k)(t,m(t)) Operates.

Let
(18a) eO(t) E w(t)
and, for m=l,2,..., let em(t) satisfy
e$(t) = J(t)em(t) + am(t) + bm(t)
(18b)
emh)==o agtib,
where
m ( k+1)
(18c) am(t> = - Zakem‘jk (t)
k=1
and
(18d) Zlbm(t)zm a Z) fi%f(k)(t.m(t))(m23em(t)zm)k.
m=l k=2 ' =1

The integer q and the constants ak are determined

from a generating function

(18e) A(z) = Zlakzqk

k=O

Gragg [11,12] has shown that both Euler's rule
and the trapezoid rule have asymptotic error expansions of

the form (2) with the coefficient functions determined by

.Yk

33

 

(18). For Euler's rule the generating function is
a) Z
l k _ e -1

For the trapezoid rule

 

_ 2 .2
(20) A(z) — z tanh<2>

_ __1.§_2 _2__§4____17(_z_)6 ___62(28

' 1 3(2) + 15(2) 315 2 + 2835 2)

_ + 1:1)“+122“(22”—1)n (§)2n-l i ...
(2n): ”n 2
. th .
where Bn is the n Bernoulli number.
From (18a) eo(t) E m(t) and from (1)

m'(t) = f(t,m(t)). m(a) = a; therefore we can find higher
order derivatives of e0(t) at t=a by successively
computing total derivatives of f and evaluating them

at t=a:

dp’lf(t.m(t1)

l

(p) _ (p)
(21) e (a) = (a) =
O m dtp- t=a,

p=l,2,...

From (18b) ei(t) = J(t)el(t) + a1(t) + b1(t)

and from (18d) bl(t) E 0. Using (18c), we find

 

al(t) = -a1eo(qk+l)(t). Thus
dp‘1[J(t)e (t)]
(p) _ l (qk+p)

The second term on the right of (22) is known from (21)

and since J(t) is known, the first term on the right

 

34

can be computed because it involves only lower order

derivatives of e1(t).

Proceeding inductively we assume that all derivatives

of e0,e1. . - at the point t=a can be evaluated,

.'e.
j-l
then from (18b)

83(t) = J(t)ej(t) + aj(t) + bj(t)

so that

dp‘1[J(t)ej(§)]

dtp—1 t=a

 

(23) egp)(a) = +agp'l1t)

J —a

+bgp'l1t)‘
t=a t-

 

Let us now consider each term on the right of (23) in

turn.

The derivatives of ej(t) appearing in

dP'1[J(t)e.(t)]
11 y are all of order less than p and can be

 

dtp'
evaluated at t=a successively from the lowest to the
highest. Also J(t) and all its derivatives can be
computed. Thus, the first term on the right of (23)

is known.

From (18c) we see that aj(t) depends on various
derivatives of the error functions e0,e1,...,ej_1 all
of which are known at the point t=a by our inductive

hypothesis. Consequently agp-l)(t) t=a is known.

35

Finally, bj(t) can be determined by collecting
the proper coefficients of zj on the right hand side
of (18d). bj(t) will depend on various Frechet derivatives
of f Operating on various error functions. The Frechet
derivatives can be computed. The error functions appearing
as arguments of the Frechet derivatives are from the
collection eo,el,...ej_1. This is easily seen by Observing
that the outer sum on the right hand side of (18d) begins
with k=2, which precludes the possibility of ej(t)

I
appearing as an argument of a Frechet derivative in the

expression for bj(t).

Once bj(t) is determined it is necessary to
find its (p-l)st derivative. The Frechet derivatives
can be differentiated with respect to t and evaluated
at t=a and the derivatives of the error functions are
known at t=a by our inductive hypothesis. Thus,

bgp-1)(t) t=a can be computed and equation (23) is valid.

Actually, equation (23) can be expressed entirely
/
in terms of eO(t), J(t), the Frechet derivatives and
various derivatives of these functions evaluated at t=a.

This will be done later for the case q=1.

Theoretically, at least, we can Obtain e$(a)
for any m=O,1,... . Computationally, the larger m is,

the more complicated the expression to be evaluated. As we

36

shall see, the computations Of ei(a), eé(a) and eé(a)
are fairly easy to perform for q=1 or q=2. In addition,

e4(a) is reasonable.

We consider the cases q=1 and q=2 separately.

In what follows we will use a superscript in parentheses

to denote differentiation with respect to t. The sole
exception to this notation will be that f(k)

to denote the kth Frechet derivative of f.

will continue

Assuming the validity Of (18), with q=1, (18c)

and (18e) become

(18c)' amu) = —k2ake§j;1’(t)

and

(18e)' A(z) = {501ka
k=O

respectively.

From (18b), with t=a, we have e£1)(a) = am(a) + bm(a).
However, since em(a) = O, Vm, we can conclude from (18d)

that bm(a) = 0, Vin. Thus,
(24) ergl) (a) = am(a)

For m=1 a1(t) = -a1 e(2)(t), so that

(25) aip)(a) = —aleép+2)(a) for any p.

37

In particular, we have

ei1)(a) = -a

e(2)
1eO (a).

91511)

(26)

In order to determine (a) for m=2,3,4 we
will need the second through fourth derivatives of el(t)
evaluated at t=a. TO obtain these we proceed as follows:

differentiating (18b) with m=1, we Obtain

d[J(t)el(t)]

 

 

 

e{2)(t) = dt + a{1’(t) + bf1’(t).
Now b1(t) a o and al(t) = ale eéz)(t ). Therefore
(27a) (2) (t) = d[J(::el(t)]-oc ale eté”( )

J(l) (t)el(t) + J(t)e](_1) (t) - Al en) (t).
Since e1(a) = o and e{1)(a) = —a1eeéz)(a ) we have
(27b) e{2)(a) = —a1J(a)e(2)(a ) — aleé3)(a)
= -al[eé3)(a) + J(a)eéZ)(a)].

Differentiating (27a) we have

(28a) e{3)(t) = d2[J(t):1(t)]-«xleé4)(t )

dt

J‘Z)(r)el(r)+2J(1’(t)e{1’(r)+J(r)e{2’(t)

e(4)
aleo (t)

by Liebnitz's rule. Using (26), (27b) and e1(a) = O,

we can evaluate (28a) for t=a and Obtain

38

-2J(l)

(28b) ef3)(a) (a)ale(2)(a)-J(a)al[e (3)(a)+J(a)e(2)(a)]

e(4)
0180 (a)

= -a1{[2J(1)(a)+(J(a))2]e (2)(a)+J(a)e(3)(a)+e(4)(a)}.

Differentiating (28a) we have
d3[J(t)el(t)]

dt3

e(5)
le 0 (t t)

 

(29a) el(t) - a

J(3)(t)el(t)+3J(2)(t)e(l)(t )+3J(l)(t)e{2’(t)

+J(t)e{3)(t)-a1fa (5)(t).

0

Using (26), (27b), (28b) and el(a) = O we have

(29b) ei4)(a) = -3a J(2)(a)e(2)(a)-3a

l J(l)(a)[eé3)(a)+J(a)eé2)(a)l

l
-alJ(a){[2J(l)(a)+(J(a))2]eéZ)(a)+J(a)eéB)(a)

e(4) (5)(a )

(a)}- —ale0

(2) (2)

= -e1{[3a (a)+5J(a>J(1’(ar+qun)3 ]e (a) +

[3J(1)(a)+(J(a))2]eé3)(a)+J(a)eé4)(a)+eé5)(a)}.

Note that all expressions for the derivatives of el
evaluated at t=a involve only the derivatives of eO

and the derivatives of J with respect to t. The
derivatives of J with respect to t are a straight-

forward calculation and the derivatives of eO can be

Obtained by successively differentiating and evaluating

)

39

the original equation (1).

(2) e(3)

For m=2: a2(t) = -d1e1(t)-OL2 (t) so that
(30) asp)(a) = -a1e{p+2)(a) - aze eép+3)(t) for any p.

In particular using (24) and (27b) we have

(31) eél)(a)==a2(a) =a1ee12)(a)-a2eé3)(a)
a1[e(3)(a)+J(a)e(2)(a)]-a2eé3)(a)
= (01- a 2)e(3)(a )+a iJ(a)e(2)(a)
Now erfll) (a) for m=3,4 will require knowledge

of the second and third derivatives of e2(t) evaluated
at t=a. Differentiating (18b) with m=2, we Obtain

d[J(t)e2(t)]+

a(l) (1)

 

(2) _
(32) e2 (t) - dt (12) + b2 (t).
From (18d) we have
(33) b2(t> = éf (”(t cp(t))e1(t)e1(t)

implying that

 

(2)
bzmm =11-df $5019) e1(t)e1(t) + f(2)(t,cp(t))e1(t)ei(t).

Note that f(2)(t,m(t))

fyy(t,m(t)) for the one dimen-

sional problem. Since e1(a) = O we can conclude that

(1) _
b2 (a) — o.

40

Using (30). (31) and e2(a) = O we have

e§Z)(a) = J(a)e§1)(a) + a§1)(a)

= —alJ(a)e{2)-Q2J(a)eé3)(3)-(119;:3)(a)"a (£4) (a)

26

—a1(e{3)(a)+J(a)e{2)(a))-a2(eé4)(a)+J(a)eé3)(a))

Hence, e§Z)(a) can be expressed in terms Of derivatives

Of eo(a) by using (27b) and (28b),
e52)(a) = ai{[2J(1)(a)+(J(a))2]eéZ)(a)+J(a)eéB)(a)+eé4)(a)
+(J(a))2eéZ)(a)+J(a)eé3)(a)]-d2(eé4)(a)+J(a)eé3)(a)).

Collecting terms involving the various derivatives of

eO(a), we have
(34) e§2’(a) = (mi-a2)eé4)(a)+(2ai-a2)J(a)eé3)(a)
+2ai[J(l)(a)+(J(a))2]eéz)(a).

Differentiating (32) again, we have

d2[J(t)e2(t)]
2

 

e53)(t) = +a§2)(t)+b§2)(r).

dt
By differentiating (33) twice and using e1(a) = O, we

find that b2(2) (a) = if”) (a,cp(a))e1(_1) (a)e{1) (a). From

(2

2 )(a)=wclei4)(a)-azeés)(a) and since e2(a) = O

(30) . a

we have

41

(1) (2) (4)

(a)-a1e1 e(s)(a)

e53)(a) = 2J(1)(a)e2 (a)+J(a)e2 (a)-a

2e

(1)

+f‘2)(a,e(a))el

(1)(a).

(a)e1

Substituting for eél)(a), e52)(a), e14)(a) and eil)(a)

from equations (31), (34), (29b) and (26) respectively,

we have
e53)(a) = 2(a12 -d2)J(l)(a)eé3)(a)+2a{2)J(a)J(1)(a)eéZ)(a)
+(e12 —<12)J(a)e(4)(a)+(2c1(2)—a2 )(J(a))2e(3)(a )

+(2c3L12 [J(a)J(1)(a)+(J(a))3 ]e (2)031 )

(2)

2 [3J (a)+5J(a)J(1)(a)+(J(a))3]e(2)(a)

+a 2 [3J(1)(a)+(J(a))2]eé3)(a)+aiJ(a)eé4)(a)

l
2 eés)(a)-azeés)(a)

+0. 2 f(z) (a,cp(a))e(§2) (a)e(§2) (a).

1
Collecting on the derivatives Of e0(a), we have

(35) e§3’(a) = (ei—a2)eés’(a)+(2a§-e2)J(a)eé4’(a)
+[(5ai-2a2)J(l)(a)+(3di-a2)(J(a))2]eé3)(a)
+ai[3J(2)(a)+9J(a)J(1)(a)+3(J(a))3]eéz)(a)

+<12f(2 )(a, cp(a))e(§2) (a)ec()2) (a).

With m=3 from (18c)' we have

42

a3(t) = -a1e2(2)(t)-a2 e(3)(t)-U3 e(4)(t) so that
(36) a§p)(a) = —dle§p+2)(a)-a2e£p+3)(a)-a 03e eép+4)(a)

for any p. From (24) we have
e§1)(a) = a3(a) = -a1e e(2)(a)-OL2 e(3)(a)-OL3 e(4)(a).

Substituting (34) and (28b) for e52)(a) and e13)(a)

respectively, we have

e§1)(a) = -a1(ai-a2)eé4)(a)—d1(2di—a2)J(a)eéB)(a)

3 J(1)

1[ (a)+(J(a)) 2]e”’ (a)

-2a

J(1)

+a o (a)+(J(a)) 2]e(2)(a )

l 2[2J

e(4 )
1 2 e0

(3) (4)(a).

NJ(a)e (a)+aa (a)-OL3 e

+(1&1 2

Combining like derivates Of eO(a), this becomes

(37) e(1)(a) = (2d1a -al-d 3) e(4 )(a)+2(a1aa2-a1)J(a)eO

+[(ala -2ai)(J(a))2+2(a1a -a M)J(1)( )]e(2)(a )-

2 2

Differentiating (18b) with m=3, we have

d[J(t)e3(t)]
dt

 

e§2’(t) = +a§1)(t)+b§1)(t).

From (18d), b3(t) = é-f(3)(t,w(t))el(t)el(t)e1(t) +

f(2)(t,w(t))e1(t)e2(t) so that b§1)(a) = O as

43

a(1)

el(a) = e2(a) = 0. Substituting for (t) and using

e3(a) = O, we have

e§2)(a) = J(a)e§1)(a)-aleé3)(a)-a2e(4)(a)-a3e(5)(a).

Substituting (37). (35) and (29b) into this equation yields

8(2) (4) Ze (3)

(a)+2(al a2 -a 3)(J(a))

(a) = (Zola -ai—a3 )J(a)eO (a)

2

(l) (2)

+[(alaz-2ai)(J(a))3+2(a -a 3)J(a)J (a a)]e0 (a)

1G2

+(ala2—ai)eé5)(a)+(ala2-2a(3))J(a)eé4)(a)

+[(2a a 3 J(l)(a) (3)

1 2-50. (a )

)+(a1a2[3al)(J(a)) ]e0

3

3 (2)
1[3J

_ J2( )l(a)+9J(a)J( ’(a)+3(J(a>]e (a >

3 (2)

—a1f (a m(a))e(2’

(2) e(5 )(a )

(a)eO 2 e0

(a)+ala

OLC11C12J(a)e(4 )(a)+a1a 2[3J(1)(a)+(J(a))2 ]e(3)( a)

a2[3J(2)(a)+5J(a)J(l)(a)+(J(a))3]eéZ)(a)-a3eé5)(a).

This can be rewritten as

(2)

(38) e3 (a) = 2(ala 2 -al- (13 )e(5)

(a)+(4al a —3ai-a3 )J(a)e(4)(a )

2

+[5(a1a2-ai)J(l)(a)+(4ala -5ai)(J(a))2]eé3)(a)

2

(3)(a)+(7a a -llal)J(a)J(1)(a) +

3
+[3(alaz—al)J l 2

(2a -SGi)(J(a))3]eéZ)(a)

1&2

-Gif(2)(a.w(a))eé2)(a)e32)(a).

44

Finally for m=4 we have

eil) (a) - 9(2)

_ a4(a) = ale 3 (a)-a 2e(3)(a)— —a 3e(4)(a)- -a4 e(5)(a).

Using (38), (35) and (29b), this can be written as

ala2)e35)(a)+(3a§+a a -4a§ a2)J(a)e(4)(a >

e(1)(a.) -— 2(a% l 3

1m 1&3-
+[5(a§-aia2)J‘1’<a>+(5a§—4aia2>(J(a>>2 Je‘3’<a >

(2)

+[3(ai-aia2)J (a)+(llai-7aia2)J(a)J(l)(a)+

(Sui-2ai02)(J(a))3]eé2)(a)

4 (2)

+alf (a, m(a)) e(2)

(2)(a)

(a)eO
+(a3-aia2)eéS)(a)+(a§-2aia2)J(a)eé4)(a)

+[(2a§-5aia2)J(l)(a)+(G§-3Giaz)(3(a))2 ]e(3)(a >

-a 2a2[3J(2)(a)+9J(a)J(l )(a)+3(J(a))3 ]e(2)( a)

—aia2f(2)(a,w(a))eé2)(a)eé2)(a)

e(5)

103 60 J(a)e(4)(a)

+a (a)+aa

la 3

+a1a3[ J(l)(a)+(J(a))2 ]+(J(a))e(3)(a )

a3[3J(2)(a)+5J(a)J(1)(a)+(J(a))3]eé2)(a)

e(5)

"O‘e4o (a)

Rearranging terms, we get

 

(1) ._

(39) e4 (L

TO SU

can be used t

 

 

 

e4 Ha) reSp.

Step 1
Step 2

Step3 e(l)

45

(39) 9:1)(3) = (ZOE-Baia2+3a103+d2-a4 )eés)(a )

4 2 2 4
+(3a1-6ala +2a1d3+a2 2)e( )(a )

+[(Sa4—10aza +3d a +2a §)J(l)

12 13 (a)+

204-0.

121“

(sag—7a +a §)(J(a )) 21e ‘3’(a >

3

2 (2) 4 2
+[3(a4 l-2a1a2+ala3)J (a)+(lla1-l6a1a2 +

(1)(a)+(5a4-5a2a +01 a3)(J(a)) 3]e(2) (a)

Sala3)J(a)J l 1 3

+(d4-aid2)f(2)(a,w(a))eé2)(a)e32)(a).

1

To summarize, equations (26), (31), (37) and (39)
e‘l’ (1’(a>. e§1’(a)

can be used to determine (a), e2 and

ei1)(a)

respectively. A computationally more effective

recursive procedure for q=1 is given in Table 3 below.

TABLE 3

Step 1 Compute eél)(a), eéz)(a), eé3)(a), e34)(a). eés)(a)....

Step 2 Compute J(a). J(l)(a), J(2)(a),

e(2)

Step 3 e11’(a) = a1e O (a )
eiZ) (a) = J(a)e{1) (a)-ale (3) (a)
ef3’(a) = 2J‘1’(a)e{1)(a)+J(a)e{2’(a)-a1e(4’<a)
e{4’(a) = 3J(2)(a )e‘1’(a)+3J(1’(a)e(2)(a)+J(a)e‘3)(a )

e(S)
O“5'10 (a)

W}

 

 

Step 5 e (1

(2

e3

Step 6 (
e

4

Tab
computatlo

46

TABLE 3 continued

Step 4 e§1)(a) = —a1e{2)(a)-a2eé3)(a)

e§2)(a) = J(a)e§1)(a)—a1e{3’(a)—a2eé4)(a)

— 2J (a)e

e§3’(a) — (1) (1)(a)+J(a)e§2’(a)—a1e{4)(a)

I aze e35)(a)+f‘2)(a.m(a))e{1’<a)e{1’(a)

Step 5 e§1)(a) = ale e52)(a)—d2 e(3)(a) -a3e(4)(a)
e§2)(a) = J(a)e(l )(a)-d1e 53)(a)—a2e (4)(a)-a 3e(5)(a)
Step 6 eil)(a) = ale e§2)(a)—aze £3)(a) ~a3ei4)(a)-a4eés)(a)

Table 3 lists the equations necessary for the

(1)(a),...,eil)(

computation of e1 a) for a method which
satisfies (18) with q=1. Of course, Table 13 can be
continued. Each new error function, whose first derivative

at t=a is to be computed, requires one additional computation
at each step of the table and the addition of one more step

to the table. The only computations in Table 3 which are not
simply the result of plugging in formulae are those in steps

1 and 2. and the computation of various Frechet derivatives

as they become necessary.

 

 

Toi

I
92(1) (a), egl

in example 1

Bample 2 :

the solution

 

From equatiox

ru1e. In th

tap 1: eél

Step 2.

J(

47

To illustrate the use of Table 3 we compute e{1)(a),

(1) (1)
92 3

in example 1 using a numerical method for which q=1.

(a), e (a) and eél)(a) for the differential equation

Example 2: Suppose we are using Euler's rule to compute

the solution to

y’=y2: y(O)=.2; ogpgl.

From equation (19) we have that dk = Tii%77' for Euler's

rule. In this case Table 3 becomes

Step 1: eél)(o) = (Y(OH2 = .040000
e32)(2) = 2y(O)y(l)(O) = .016000
eé3)(0) = 2[y(O)y(2)(O)+(y(l)(0))2] = .009600
eé4)(0) = 2[y(O)y(3)(O)+3y(l)(0)y(2)(0)] = .007680

e55) (0) = 2[y(O)y(4) (O)+4y(l) (on/(3) (o)+3(y(2) (0))2]

= .007680
Step 2: J(O) = 2y(o) = .400000
J(l)(O) = 2y(1)(0) = .080000
J(2)(l) = 2y(2)(0) = .032000

 

 

Step 3: e{‘

(‘

e1

e;

e{‘

Step 4; e;

(

e2

(

‘32

Step 5. (

- e

3

(

e3

StEp 6- (

° e

4

It

180 note t
functiOn

48

Step 3: e{1)(0) = -.008000
e{2)(0) = -.008000
e{3)(0) = -.008520
e{4)(0) = -.oo9936

Step 4: e(l)(0)

2 = .002400
e52)(0) = .001810
e53)(0) = .006076 since f(2)(a,w(a)) = 2.
Ste - (l) -
p 5. e3 (0) — -.OO7950
e§2)(0) = -.oo4882
Ste 6- (l) -
p . e4 (0) — .0017780

It should be emphasized that the computations in
Steps 1 and 2 are performed by the user not by the machine.
Also note that the signs on the derivatives of the error
functions alternate in precisely the same manner as the signs

of the errors for Euler's rule.

Turning our attention to the case q=2, we assume

(18) is valid with q=2 so that (18c) and (18e) become

 

(18c) "

and

(18e) "

In
anaIOgous

abOUt the

49

m
u _ e(2k+1)
(18c) am(t) - 31:31 em“k (t t)
and
(18e)" A(z) = 2: 22k
k=oOlk

In order to construct a table for the case q=2
ana10gous to Table 3 we will need to collect more information

about the functions bm(t) for m=1,...,4. From (18d) we have

b1(t) a O
b2(t) = if f‘2’(t w<t>>e1(t>e1<t>

(40) b3(t> = if ‘3’(t w(t>)e1(t>e1(t)e1(t>+f‘2’(t.m(t)>e1(t)e2(t)
b4(t> = 24f‘4’(t.w(t)>e1(t)e1(t)e1(t>e1(t)+—f(3’<t.m(t>)

e1(t)e1(t)e2(t>+f‘2’(t.w(t>)e1<t)e3<t)

+%f(2) (t.¢P(t) )e2(t)e2 (t) .

Since our goal is to produce a table for evaluating
e(l)(a )....,eé1)(a) we examine what information will be

necessary to determine e4 a). From (18b) with m=4 we

have

(1)

e4 (a)=a4(a)
since e4(a) = b4(a) = O. From (18c)" we can conclude
e§1)(a) = -ale e§3)(a)-d2e £5)(a)-a3e{7)(a)- d4eo e(9)(a).

f

 

In genera
P
em(a)

Thus in .
to deter
(7)

91 (a)

given be

function

Zero, I;

therefOI

(3)
b2

(4)
b2

50

In general, by (18b) we see that

dp-1[J(t)em(t)]

 

82(a) = + amp-l)(t)‘*b£p-l)(t): m=1.....4.

dtp‘l

we need to know b§2)(a):

(a) we need bé4)(a); and to determine

Thus in order to determine

.15)

e§3)(a)

to determine
(7) (6) .

e1 (a) we need bl (a). The necessary computations are

given below. To simplify the notation, the arguments of the

functions are suppressed.

Since b1(t) a 0 all derivatives of bl(t) are

zero. From (40), b2(t) = %f(2)(t.¢(t))e1(t)e1(t) and

therefore:

(2)
bunt) — 14—.“ f e .42) em

2 ’ 2 dt e1 1 ee1 1
b(2)(t) _ .1. 2313.123 8 e +2(13(2))(131(1(1)+f(2)e6(2)
2 ' 2 dtz 1 1 ee1 1 ee1 1 V
e1(1)e (1)
+f( e1
3 (2) (2)
(3) _ 1. d f (2) (2)e e(1) dif 2 e(2)
b2 (t) ‘ 3' dt3 elel+3(f ee1 1 *3 dt ee1 1
+3 «113121e11)e11)+3f(2)e11)eze12)+1.})eele13)
(4) __1. d 4f(2) d 3f(2) e(1) d 2f(2) (1) (1)
b2 (t) ’ 3' dt4 e191+4 dt3 ee1 1 +6 dtz e1 e1
931121. (2) 21.11216 1613) dfg 1(2)) (1) (2)
+6 dtz elel +4 dt +12 dt e:L el
+3,1.(2)e12)e12)+11f(2)e11)e13)1:9) e14)

”1

Using the fact that e1(a) = O we can obtain

 

10(1)
(2)
b2

(41) b”)

Us
have

Q'
an
Qe

51

O‘
A
H
A
m
l
O

b(2)(a) - f(2)(a.CP(a))e(1) (a)e (1)(a )

 

 

 

2
(41) 132(3) (80 = 3[f(2) (a.CP(a))e(l) (a)e(Z) (a)
df(2) u: mm) (1) (1)
+ dtl e1 (a)el (a)]
t=a
2 (2)
(4) d f Lt.@(t)) (l) (1)
b (a) = 6 e (a)e (a)
2 dt2 t=a l l
dfu) (t 2112)) (1) (2)
+12 dtJ. e1 (a)el (a)
t=a
+3f(2)(a w(a))e(2’(a)e(2)(a )+4f(2)

(a.w(a))e{1)(a)e{a)(a).

Using the definition of b3(t) given in (40), we

 

have
(1) 1 df‘3’ 1 (3) (1) (2) (1)
b ()-a "Er-919191715 916191 ”f ) 9192
(2) (1) (2) e(1)
+ f el e2+f ele 2
2 (3) (3)
(2) l f 91.1:— em (3) e(1)9(1)
b3 (t) ‘ 6 dtz $19191+ dt e1918 1 +fee1 e1
2 (2) (2)
(3) (2) d f df (1)
“‘Ef 919191 + ‘13” “3192+2 T e1 e2
(2)
df (l)+ f(2) (2) (2) e(1) 62(1) e(2)
+2 T ele 2 e1 e2+2f +f( eel 2 .

Since e1(a) = e2(a) = O, at t=a we have

 

<
b3
(42)
t
b3

V

evaluatir

Step 1:
Step 2:

Step 3:

52

b3(1) (a) = O
(42)

b§2)(a) 22(2)(a.w(a))e{1)(a)e§1)(a).

We can now summarize the computations involved for

(l)(a).....eil)(

evaluating el a) when q=2.

TABLE 4
Step 1: Compute eél)(a), eéz)(a)....,eé8)(a), e(9)(a)....
Step 2: Compute J(a), J(1)(a),...,J(5)(a),...
Step 3: e{1)(a) = 18(3)(a)
e{2)(a) = J(a)e{1)(a)-oleé4)(a)
ef3’(a) = 2J‘1’(a)e‘“ (a)+J(a)e(2)(a)-o1e é5’(a)

e{4)(a) = 3J(2)(a)e{1)(a)+3J(1)(a)e{2)(a)+J(a)e{3)(a)

e(6)

3
efs’ (a) = k2 @110" (a)e{4‘k’ (a)—ole (7) (a)
=0

4
e{6’(a) = k2:(:)J‘k’(a)e{S‘k)(a)-a1ee(8)(a)
=0

5
ei7)(a) = k2)<:)J(k)(a)e{6-k)(a)-ole e(9)(a)
=0

 

Step 4:

Step 5;

 

AiA

A

A I

Step 4:

Step 5:

b2”) (a)

b52)(a)

bé3)(a)

13(4) (a)

eél)(a)

e(2)(a)

e53)(a)

e54)(a)

e(5)(a)

TABLE 4

O

f(2 )(a

53

continued

(1)

@(a) ))e1 (a)

)ef1)(a)

3[f(2) (a,cp(a))el(1) (a)elQ) (a)

MmNtmo1

f(2)

+ 12

dt t=a

(t ¢(t))

dt2 t=a

df‘2)(t.mlt))

dt

e{1)(a)e{1)(a)]

e{1’(a)e{1)(a)

ef1)(a)e{2)(a)

t=a

+3f‘2’(a.w(a))e{2’(a)e{2’(a)

+4f(2)(a.®(a))e{l)(a)e{3)(a)

e(3)

-Cll

J(a)e§l)

2J(l)(a)e2(l)(a)+J(a)e§2)(a)-dlel

e(7)
0‘2 e0

Iltqm

e(8)
O‘2 eo

Iltjw

e(9)
0‘92 o

(a)-a2 e

(a)_% e(4)

(a)+b§2)(a)

“@ch

(a)+b:§3) (a)

(new

(a)+b§4’(a)

(a)-o e

(5)(a)

(
2 o

(a)-o e

6)

(7)

l l

(a)

(5)(a)

(6)(a)

(a)-a1e1

(a)

54

TABLE 4 continued

Step 6: b§1)(a) = o

13:52) (a) = 2f(2) (a.cp(a))e{1) ( we?) (a)

Step 7: e§1)(a) = ole e53)(a) -o2e{5)(a)-a3 e(7)(a)
e§2) (a) = J(a)e3(l) (a)-dle 2(4) (a)- -a2 em) (a)-—a eéa) (a)

e§3)(a) = 2J(1)(a)e§1)(a)+J(a)e§2)(a)

—ale e55)(a)-a2e{7)(a)— —a3e e(9)(a)+b§2)(a)

Step 8: e§1)(a) = ale e§3’(a)- -Q2e Wé5)(a) e(7)(a)-d 4e(9)(a)
Table 4 can be continued. Each new e£1)(a) to be

computed will require two additional computations at each
existing step of the table and the addition of two new

steps to the table. One new step will contain the computations
of b£_1(a) for p=l and 2 and will preceed the step

where eéii(a) is computed. The other new step will be the

(l)(a).

computation of em

(1)

Example 3: Compute el (a)...-.eél)(

a) for the equation

y‘ = y : y(0) = .2: ogt_<_l when the trapezoid rule is

55

used as the numerical method for solving the problem.
The generating function for the coefficients in (18)

for the trapezoid rule is given by (20). Specifically,

_ _1_ -_i_ _:;l_7_._1_ -__6_2_,_1_
a1 ‘ ‘ 2' 0‘2 ‘ 120' 0‘3 ' 315 26 and 0‘4 “ 2835 28'

If the computations in Table 4 are performed for this

example the results are: ef1)(a) = .0008: eél)(a) = .000032:

6

e§1)(a) = 1.3867 x 10- : and eé1)(a) = 6.229 x 10-

As is apparent, the computation of eé1)(a) for
nu24- becomes quite involved, particularly when q=2.
However, there are applications (see Section 1 of Chapter 3)
of the pullback method with M5;4 which are more accurate

than other methods currently available.

56

Section 5. Numerical Results for Initial Value PrOblems

 

In this section the pullback method and the method
due to Lindberg [19] are compared theoretically and numerically.
Since Lindberg's approach is conceptually and notationally
quite different from ours, a discussion of his method is
included. we will confine the discussion initially to the
case M=4 and q=2 as this will suffice to point out the
differences between the methods. Unfortunately, Lindberg's
notation and that used here are in some instances in complete
Opposition to one another. The differences will be pointed

out in footnotes.

Let h:>O be the basic steplength and for k=O,l,...,4
define steplengths hk = h/2k and grids Gk % {tE=a+ihk:

A
i=0,...,2k}. In addition, let h = h/24 = h/l6*.

Assume that the numerical method being employed to

solve the problem (1) is such that the expansion

(43) Y(t§.hk) = cp(t}i()+el(t}i()}5i+e2(t]i‘)h:+e3(tis)h:+e4(t]i<)h:+o(h%0)

is valid for each i=0,l,...,2k and for each k=O,l,...,4.

Lindberg's method is as follows. Number all grid

points according to their order of occurrenceimithe finest

A
*In Lindberg's paper h is taken to be the basic steplength
and hk=2kfi‘ k=O,l,...,4. Thus the grids GR and steplengths
hk are numbered in the Opposite order.

57

grid, obtaining t =a,tl,...,t16=a+h. At each ti

0

perform as many extrapolations as possible and denote the
A
computed solution after n extrapolations as Yn+l(ti'h)'

Thus at t we can perform 4 extrapolations obtaining

16

Q a A A
Yl(t16.A). Y2(t ), Y3(t h),..., Y5(tl6.h) where

16' 16'
Y1(t16,h) is the solution computed with the numerical

. A 4
method: 1.e., Y1(ti,h) = Y(ti,h4). At t8 we can perform

3 extrapolations: at t4 and t 2 extrapolations

12'

t and t l extrapolation.

and at t , t 10 14,

2 6'

A
According to Lindberg each Yj(t,h) satisfies the

relationship

CN’C) + EEK]. Vev(t)h

v-J

A
(44) Yj(t,h)

where for each j

j-l 22p_22v
(45) Khv = H. ——3—-—-
J p=l 2 p--l

The goal here is to define Y5 at all points of the

finest grid. Initially YS is known only at the two points
A
to=a and t16=a+h=a+l6h. At each of these points form
A A A8 "10 . .
Y5(t,h)-Y4(t,h) = —Xa4e4(t)h +O(h ) by (44). Via linear

interpolation obtain an 0(hlo)* approximation to

A
- 44e4(t8)h8. Add this approximation to Y4(t8.£) to

A
obtain Y5(t8,£3+o(hlo). At the next stage form the

*This will be discussed in some detail later.

58

A A
differences Y5(t,h)—Y3(t,h) at the points t=tO,t8 and

_ , 6 , 8

t16- By (44). Y5(t,Q)—Y3(t.Q)—-A33e3(t)Q —A34e4(t)fi

+o(hlo). Using quadratic Lagrange interpolation obtain
3 . . , 6 8

an 0(9 ) approx1mation to -A33e3(t)fi -X34e4(t)fi at

t=t4 and t=t12 which is added to Y3(t,h) at these points

to obtain Y5(t,Q) for t=t and t=t However Y5(t,fi).

4 12'

t=t4,t12, is only an 0(h9) approximation to m(t). This
process is continued until Y5 has been defined at all
ti' i=O,...,16. The result, according to Lindberg is
Y5(ti,fi)=w(ti)+0(fig) for i=1,...,7 and 9,...,15 while

A A
Y5(ti,h)=m(ti)+o(hlo) for i=8 and 16.

Now the pullback interpolation method will yield
0(hlo) accuracy at all ti's and Lindberg's method yields
0499) at most ti's. Since ‘£=h/16 it looks as if Lindberg's
method is of higher accuracy. This is not so for the following
reasons: First, Lindberg presents the error analysis for
extrapolation in terms of the smallest stepsize available,
h. The '0' notation is designed to be information-
supressing and the remainder terms in extrapolation which are
said to be OKhB) are actually of the form Chag(t)
where C and B are constants and g(t) is a continuous
function of t. On a closed t interval g(t) will be
bounded and we can write Cth (t) _<_CohB = 0(h6) . We could also
write Chfig(t) = éfiag(t)5;éofi6 = 0(fi6) by defining

A A
C = 166C. The prOblem here is the fact that CO will be

59

quite large compared to CO. When the error in extrapolation
is expressed in terms of the largest or basic steplength,

h, the constant C will be smaller than 1. If this same
error is expressed in terms of the smallest stepsize, Q.

the constant 8 will be of the form (45) and is larger

than 1. In fact, when expressing the error in terms of the
smallest stepsize we are unable to prove the convergence of

the Aitken-Neville extrapolation scheme; indeed, an analysis

of (45) shows that the constants kkaw as the stepsizes
hk40.

Secondly, the magnitude of the constant 6 is
further increased in Lindberg's analysis when interpolation
is performed. Recall that Lindberg expresses the error
in interpolation in terms of the smallest stepsize also.
This means that at the first step when Lindberg is
performing linear interpolation to - 44e4(t)h8 at the

points t0 and tl6 he obtains the error

(t—to)(t-t

)
l6 (2)
2! K

e (§)h8 where t "§<(t

44 4 O\‘ Lindberg

16'
A
calls this error 0(hlo); to do this one must interpret
A
the maximum of (t-to)(t—tl6) as being ofhz). The

maximum occurs at the midpoint t8 and is easily seen to be

2 2A2
C§>2 = gr-= lézh—u To call this 0(h2) further increases

the size of the constant being suppressed. Note that in the

error analysis for the pullback method, (g)2 is inter-
preted as -1-h2 =O(h2) which causes the constant to decrease.

4

60

When Lindberg claims an error of the form 0(93)
and the pullback method claims OKhB), the inequality,
0(hB) _<_ 00/16), holds when the constants CO and 30
are taken into account. Thus OKth) in actuality is

A
smaller than 0(h9) and the pullback interpolation method

yields a more accurate solution.

This increase in accuracy is due solely to the fact
that we are using Hermite interpolation with the additional
data e£1)(a). In fact, if the pullback interpolation
method is changed so that Lagrange interpolation is used
the results are identical to Lindberg's. In this case

on?) a 0096).

Let's examine the differences between using Hermite

and Lagrange interpolation when q=2 for arbitrary M.

For j=O, we know eM(t) at two points with
accuracy oXhZ). Let PM(t) be the Hermite polynomial
constructed as outlined and let LM(t) be the Lagrange
interpolation polynomial for these two data points. We have
Pu“) = eM(t)+O(h3)+O(h2) and LM(t) = eM(t)+o(h2)+o(h2)
where the first '0' term is the error in interpolation and
the second is the error in the given data. 'Thus both
PM(t) and LM(t) are Olhz) approximations to eM(t).
When LM(t) is used numerically, as it is in Lindberg's

method, the solution at the midpoint of the interval is,

61

in most instances, less accurate than the solution which can
be obtained by just performing extrapolation (see the examples
concluding this section). This phenomenon,MflfixfliLindberg
calls attention to in his paper, does not occur when

PM(t) 15 used.

For j=1, we know e 1(t) at three points with
m-
accuracy 0(h4). In this case, PM_l(t) = eM_1(t)4-O(h4)
4 _ 4 _ 3 4
+001 ) — eM_1(t) +O(h ) and LM-l(t) - eM_1(t) +0(h ) +G(h )
= eM_l(t)+-O(h3). Here it is the error in interpolation

which determines the final accuracy so that there is a real

gain when Hermite interpolation is used.

The situation is the same when j=2, that is,

PM_2(t) = eM_t(t) +0016) +0016) while LM_2(t) = (t)

eM—Z
+ 00h5)4-o(h6). Again the error in interpolation dominates

and Hermite interpolation yields higher accuracy.

For any jyg3, it is the error to which we know

eM_j(t) that dominates, since for j;33 we have

23+l > 2j+2. The last inequality can be established by a
proof similar to that of Lemma 1. Thus the advantages of

using Hermite interpolation occur when j=O,l and 2.

We now examine some examples of results Obtained
using the pullback interpolation method when numerically
solving first order initial value problems. All numerical

computations were done on the CDC 6500 installation at

62

Michigan State University. Calculations were done in single
precision floating point arithmetic yielding l4 accurate

decimal places.

The first example in this section is a continuation

of Example 2, Section 4.

Example 4: Solve y'=y2, y(O)=.2, Qitgl using Euler's

method and pullback interpolation with M=4.

The basic stepsize is taken to be h=l and, since

M=4, the finest grid, G consists of 1? equally spaced

4!
points in [0,1]. That is,

G4 = (ti: ti=0+1 h/l6, 1=O,l,...,l6}.

The theoretical solution is y(t) = l/(S-t) and the
numerical results are reported in Table 5 below. Table 5 is
constructed to exhibit the error in the solution computed

by Euler's rule and the error in the solution computed by
pullback interpolation at each grid point of G4. In each
case the error is the computed solution minus the theoretical

solution.

63

TABLE 5

Error using Euler Error using pullback

10
ll
12
13
14
15

16

-3.15 10‘ .01 10'

—6.53 10‘ 2.97 10’
-1o.11 10‘ 8.47 10‘
—13.93 10’ 13.16 10‘
—17.99 10' 14.87 10‘
-22.31 10' 13.26 10‘
-26.92 10’ 9.17 10'
-31.82 10‘ 3.62 10’
-37.05 10’ -2.01 10'
-42.63 10' —6.58 10'
-48.58 10' -8.92 10'
-54.93 10' -8.26 10'
-71.71 10’ —4.97 10'
—68.95 10‘ -l.26 10‘
-76.68 10' - .54 10‘
—84.96 10’ -3.30 10'

 

64

The results given in Table 5 can be quickly summarized
by noting that Euler's method can guarantee only two accurate

digits at all grid points of G while the pullback method

4

can guarantee 5 accurate digits.

Pullback interpolation based on Euler's rule with
M=4 was also used to compute the solution to yH=y2,
y(0) = .2 with basic stepsizes of h=l/2 and h=l/4.
In the first case the solution is computed on [0,1/2]
and in the second case the solution is computed on [0,1/4].
The results are summarized in Table 6 below. The first two
columns are the results for h=1/2 and the last two are
the results for h=l/4. Table 6 is arranged to exhibit the
correspondence between grid points for the two different

mesh sizes.

 

65

TABLE 6
h=bQ h=y4
i Error using pullback i. Error using pullback
0 0 0 0
1 .15 x 10'10
1 1.96 x 10"10 2 11.93 x 10’10
3 33.93 x 10‘10
2 189.25 x 10"10 4 50.92 x 10'10
5 55.90 x 10‘10
3 532.79 x 10"10 6 147.10 x 10'10
7 27.94 x 10'10
4 817.32 x 10"10 8 3.45 x 10"10
9 21.15 x 10"10
5 905.78 x 10‘10 10 40.66 x 10"10
11 50.02 x 10"10
6 777.95 x 10"10 12 45.84 x 10'10
13 29.33 x 10‘10
7 486.78 x 10’10 14 9.12 x 10’10
15 1.17 x 10‘10
8 110.79 x 10"10 16 .40 x 10'10
9 269.05 x 10"10
10 571.77 x 10"10
11 718.91 x 10"10
12 657.64 x 10'10
13 408.36 x 10’10
14 108.33 x 10"10
15 26.89 x 10'10
16 32.13 x 10-10

66

As can be seen by examing Table 5 and Table 6 the
smaller the basic steplength is the more accurate are the
computed solutions. This is in agreement with results for

exprapolation (see Lambert [18] and Gragg [11,12]).

Example 5: Solve y'=y, y(O)=l, ogtgl using the trapezoid
rule, and pullback interpolation with M:4. Using a basic
steplength of h=l we again obtain the solution at 17
equally spaced points in [0,1]. The theoretical solution

of this equation is y(t)=et and the error reported is the
computed solution minus the theoretical solution. The results
of Example 5 and the next example will be presented simul-

taneously in Table 7.

Example 6: Solve y'=-sin t, y(O)=1, 03631 using the
trapezoid rule and pullback interpolation with M=4. The
basic steplength is taken to be h=l and the grids are the
same as in Example 5. The theoretical solution is given by

y(t)=cos(t).

67

TABLE 7

Error in pullback Error in pullback

N F4 O +4

W

10
11
12
13
14
15
16

for y'=y with h=l for y'=-sin t with h=l
0 0

.12 x 10-10 - x 10-13

.48 x 10"10 — x 10-13
1.05 x 10"10 — x 10'“13
1.70 x 10"10 - x 10'13
2.30 x 10’10 - x 10"13
2.74 x 10‘10 - x 10’13
2.94 x 10.10 O. x 10'-13
2.90 x 10"10 - x 10"13
2.65 x 10'10 - x 10"13
2.34 x 10.10 - x 10-13
2.11 x 10’10 x 10"13
2.16 x 10'10 1. x 10'13
2.63 x 10"10 2. x 10"13
3.57 x 10"10 3. x 10'13
4.86 x 10-10 2. x 10”13
6.13 x 10'10 -2. x 10'13

The equations in examples 5 and 6 were solved in the

68

same manner with a basic stepsize of h=l/2 to obtain

solutions at 1? equally spaced points in

[0.1/2].

The

fireesstnlts of these computations are presented in Table 8.

H0?”

10
ll
12
13
14
15
16

TABLE 8

Error in pullback
for y'=y with h=l/2

N

@014:-

O

X

10-13

10'13

10-13

10-13

10-13

10-13

10—13

10—13

10-13

10-13

10-13

10-13

10-13

10"13

10—13

10—13

Error in pullback
=—sin t with h=l/2

for y'

0

X

10-14

10-14

10-14

10-14

10'14

10-14

10"14

10—14

10-14

10-14

10-14

10-14

10-14

10-14

10—14

10-14

 

69

As the results presented in Table 7 and Table 8
so vividly indicate, pullback interpolation will yield

uni form accuaracy at all grid points of the finest grid.

Lindberg in [19] presents numerical results

obtained when solving y'=y, y(O)=l, on [0,1] using the

trapezoid rule and his interpolation method. The largest

error Lindberg obtains is at t5 and has magnitude l6xlO-lo.

Examining Table 7 we see that the largest error produced for

this equation when using the pullback method is the error
in extrapolation at the endpoint, tl6=1. The magnitude of

this error is 6.13xlO-lo.

In the next series of examples extrapolation,

Lindberg's method and pullback interpolation are compared

numerically.

Example 7: Solve y'=y2, y(O)=l, ogtgl using the trapezoid

rule . Choosing a basic steplength of h=l and M=3, the

solution is computed first by using extrapolation as often
as is possible at each grid point of the finest grid,

G3 = {ti: ti=0+i h/8, i=O,l,...,8}. Secondly, the solution

is Computed on G3 by using Lindberg's method as described

earlier in this section. Lastly, the solution is computed

on G3 using pullback interpolation. All errors are
given as the computed solution minus the theoretical solution.

The results are presented in Table 9.

70

TABLE 9

Error using
pullback

Error using
Lindberg's method

Error in the best

i. extrapolated value

(1 -0 0 0

1 1.7 x 10'.6 .8 x 10-10 - .4 x 10-10
2 -1.0 x 10’8 -26.7 x 10'10 —3.0 x 10’10
3 5 9 x 10"6 -30.3 x 10‘10* —4.7 x 10’10
4 3.7 x 10'10 - 5.9 x 10'10 -1.9 x 10'10
5 1.2 x 10'5 22.0 x 10-10 3.7 x 10'10
6 -5.3 x 10‘8 23.9 x 10'10 6.2 x 10'10*
7 1.9 x 10-5* - 4.0 x 10—10 .8 x 10-10
8 -1.5 x 10‘10 - 1.5 x 10'10 —1.5 x 10"10

*greatest absolute error for the method

Example 8: Solve y'=-sin t. Y(0)=1. 053131 using the

trapezoid rule. We again take h=l and M=3 and compute

Solutions using extrapolation, Lindberg's method, and

IplllJLIxack interpolation. The results are given in Table 10-

71

TABLE 10

Error in the best Error using Error using

i extrapolated value Lindberg '3 method pullback

0 0 0 0

1 1.02 x 10"5 178.36 x 10’10 9.46 x 10’10
2 - 4.22 x 10"8 37.54 x 10'10 6.33 x10-10
3 9.05 x 10"5 - 47.24 x 10"10 — 2.47 x 10'10
4 9.96 x 10'10 - 9.23 x 10’10 .37 x 10'10
5 24.62 x 10‘5 31.70 x 10-10 4.92 x 10‘10
6 —36.48 x 10‘8 - 63.15 x 10'10 —17.54 x 10"10
7 46.76 x 10‘5* -239.94 x 10‘10* -62.65 x 10‘10*
8 .96 x 10‘10 - .96 x 10'10 - .96 x 10'10
*greatest absolute error for the method

Comparing the results in Example 7 and Example 8,

we see that the solutions obtained with both Lindberg's

method and the pullback interpolation method,

mediate grid points

the best extrapolated values.

t.,
1

i=1,

...,7,

at the inter-

are more accurate than

The pullback interpolation

method is the most accurate of all, being almost a full

decimal place better than Lindberg's method at all inter-

mediate grid points.

As mentioned before,

one drawback to Lindberg's

method is that his results at the midpoint of the interval

are no better,

and in fact often worse,

than the results

72

obtained by just performing extrapolation. This phenomenon

does not occur when the pullback interpolation method is

used, as a comparison of the errors at t4 in Table 9 and

Table 10 clearly demonstrates.

Examining Table 10 we note that the pullback method

is dramatically less accurate at t7 than it is at all

other grid points. In an attempt to explain this behavior we

computed the following example.

Example 9: Solve y'=-sin(t), y(0)=l on Qgtgl/Z. The

only difference between this example and Example 8 is that we

now take h=l/2 as the basic stepsize. The results are

presented as

TABLE 1 1

Error in the best Error in Error in
i extrapolated value Lindberg's method pullback
O O O O
1 6.36 x 10‘7 71.71 x 10'12 3.89 x 10‘12
2 —6 62 x 10"10 15.02 10'12 2.55 x 10'12
3 5.71 x 10"6 —19.21 10'12 —1.15 x 10"12
4 -3.92 x 10"12 — 3.87 10"12 .03 x 10'12
5 1.58 x 10‘5 12.74 10"12 1.99 x 10'12
6 -5 89 x 10‘9 —25.85 10‘12 —7.54 x 10"12
7 3.07 x 10'5* —98.l4 10'12* -26.90 x 10’12*
8 —1.1 x 10‘13 .11 10"12 - 11 x 10"12

*greatest absolute error for the method

73

Once again t7 is the least accurate solution for

all three methods. The only reasonable explanation which
suggests itself is that the solution (xx; t is changing
very rapidly (relative to its behavior on the rest of the
interval) between t7 and t8. This behavior is compen-

sated for by extrapolation at t8, but since no extrapolation

is done at t no correction is possible there.

7

From these examples it appears that pullback

 

-_-7._s
I
II

interpolation is quite sensitive to the accuracy to which
we know the solution at the intermediate grid points and
the type of equation we are solving. Thus it appears that
it is the error in the original data and not the error in
interpolation that is determining the overall accuracy of

the method.

We next examine what happens when we solve the same
equation on intervals [a,b], [b,c] and [c,d] using the

last computed solution as initial data for the next interval.

Example 10: Solve y'=y2, y(0)=l, 03tg3/2 ‘by computing
the solution on intervals of length L/2 and using the com-
puted solution at the endpoint of the previous interval as
the initial data for the equation on the next interval.

The method of solution on each interval will be the trapezoid
rule and pullback interpolation with M=3. The results are

presented as Table 12.

74

TABLE 12
[0.1/2] [1/2.1] [1.3/2]

i Error i Error i Error

0 0 0 0 0 0

1 - 01 x 10’12 1 56 x 10‘12 1 1.13 x 10‘12
2 - 55 x 10‘12 2 -1 20 x 10'12 2 -4 39 x 10'12
3 - 74 x 10"12 3 -1 95 x 10"12 3 —6 93 x 10"12
4 .13 x 10"12 4 .27 x 10'12 4 - 64 x 10’12
5 1.59 x 10’12 5 4.19 x 10'12 5 11.24 x 10’12
6 2.21 x 10'12* 6 5.85 x 10‘12* 6 16.07 x 10'12*
7 94 x 10’12 7 2 39 x 10"12 7 5.04 x 10"12
8 56 x 10'12 8 1 13 x 10‘12 8 .84 x 10"12

*greatest absolute error

Some deterioration in the accuracy of the computed
solution for larger t can be observed in Table 12. This
can be partly attributed to an accumulation of roundoff
errors and partly to the sensitivity of the pullback method

to the accuracy of the data it receives.

If we compare the first two columns of errors in
Table 12 to the results given in Table 9, we see that in
terms of greatest absolute error it is better to take a
smaller basic steplength and solve the prOblem several times
in succession. However, we should point out that this can

be quite expensive computationally.

 

75

Also, note that the largest absolute error always
occurs at the same relative position of the three grids,
namely at t6. This is further support for our contention
that the nature of the equation and the accuracy of the data
points used when interpolating are the factors which control

the overall accuracy of the pullback method.

In each example computed using the trapezoid rule the
largest error in the pullback method has occurred to the right
of the midpoint. This is no coincidence. The trapezoid rule
yields more accurate solutions at grid points nearer the init-
.ial point and consequently extrapolation will be more accurate
in the first half of the interval. Also, when we perform
Hermite interpolation we have an extra data point at the
initial grid point. All of these factors combined make it
reasonable to expect that the largest errors will be produced
to the right of the midpoint. Thus, a solution computed
using pullback interpolation should be most reliable in the

first half of the interval on which the solution is computed.

Note that the largest error in Lindberg's method
can occur in the first half of the interval, as an examination
of Table 9 reveals. This is due to the fact that it is the
error in interpolation which determines the overall accuracy

of Lindberg's method.

In summary, pullback interpolation coupled with the

trapezoid rule will yield highly accurate solutions at all

 

76

grid points. In fact, if the solutions obtained with the
trapezoid rule and extrapolation are sufficiently accurate,
pullback interpolation will yield uniform accuracy at all

grid points.

Pullback interpolation coupled with Euler's rule
is not recommended as a viable solution technique. Euler's
rule is simply not accurate enough to enable pullback

interpolation to operate effectively.

CHAPTER II

TWO POINT BOUNDARY VALUE PROBLEMS

Section 1. The Problem and Its Discretization

 

In this chapter we will consider two point boundary

value problems of the form

x"(t)-f(t,x(t) ,x' (t)) = 0

(l)
x(a) = A, x(b) = B.

1
We assume that f(t,x,x') 6C [[a,b] x(-”.”) X(-”.”)].

f(t,y,z) is uniformly Lipschitz continuous in y and z,

0 < e < g;- and [€515;K where K is a constant. Under

these assumptions (see Keller [17]) problem (1) has a

unique solution which we will denote by 0(t).
The continuous problem (1) will be denoted by
(1') F(x) = 0.

The operator F(x) 2 x"(t) — f(t,x(t),x'(t)) maps the
Banach space of twice continuously differentiable functions

defined on [a,b] into C[a,b].

To obtain a discrete version of (1), let n;12 be

any natural number, define h=(b-a)/n and form the uniform

mesh [tk=a+ih}ril=O in [a,b]. The discrete problem is then

given by

77

 

78

 

 

X. -2X.+X. X. -X.
144. 5. l l _ f(t.,X., 1+1 1—1) _ 0, 1:1, ,n-l,
h 1 1 2h
(2)
X0 = A, Xn = B,
where we have introduced the notation Xi==X(ti). Problem

(2) may be thought of as a nonlinear system of equations in

En”l with the unknown being the vector (Xl'°"'Xn-1)T'

The solution to (2) will be denoted by X(h) and we

introduce the Operator notation
(2') Fh(X) = 0
for problem (2).

In order to set up the correspondence between the
continuous and discrete prOblems we will need to define a
space discretization Operator wh' Thus, let z(t) be an

arbitrary function defined on [a,b] and define ”h acting

on 2 by

mhz mh(2(t)) = (2(t1).....2(tn_1

E (21,...,zn_l

Note that mnz is the vector in En.l whose components

are obtained by evaluating 2 at the grid points.

Problem (1) and the discretization (2) have been

studied by Pereyra [24,27]. Stetter [31] and Pereyra [26,28]

ihave also studied the special case of (1) when x' is

not present in the equation.

 

79

To formalize what is meant by convergence of the
discrete solution X(h) to the continuous solution 0(t)
we have the following definition (see Lambert [18] and

Pereyra [28]).

Definition 1: We shall say that the discrete solution X(h)

 

converges discretely to the theoretical solution 0(t) if

and only if

(3) iig\\x(h)—%m][(h) = o,
where “'“(h) is the maximum norm on En-l.

The subscript (h) will be omitted from the norms
throughout the remainder of this chapter. We should point
out that Pereyra in [27] utilizes a different norm than
the one we have used here. while in his earlier work [24] he

uses the maximum norm.

Typically, discrete convergence depends on the two
properties of consistency and stability of the discrete Operator
Eh in (2'). In formulating the definitions of these concepts
xMe have followed the approach of Pereyra [25.28].

Definition 2: The Operator Fh is said to be consistent of

<arder p:>0 if, for each solution 0(t) of (l) and for all

h_<_ ho we have

(4) [[13th n = 0011”) .

 

80

Definition 3: The operator Fh is stable if for any pair

 

of discrete functions U and W and for all hfiho there

exists a constant CI>0, independent of h, such that

(5) \[U-WH _<_ 0]]Fh (U)-Fh (w) \\.

The standard result for discrete Operators, that
consistency and stability imply discrete convergence, is
valid here. This theorem, in the context of linear multistep
methods for ordinary differential equations, is originally
due to Dahlquist [8,9]. For completeness we present Pereyra's

statement of this result and briefly summarize his proof.

Lemma 1 (Pereyra [ ]): If Fh is stable then it is locally
invertible around mhm and the inverse mapping is locally

Lipschitz continuous for all hiiho:

Proof: Let Bh = B(mhm,p), be the Open ball of radius p
centered at whm, where p::0 is independent of h. If

U,W¢=B then stability implies that Fh is one—to-one

h
on Bh’ for otherwise we can violate (5). Therefore
thBh4Fh(Bh) is a bijection, implying that Fgl exists on
Fh(Bh). If X,Y€EFh(Bh) then we can write (5) as

\\F;11(X) 41:1 (Y) ]]_<_C]\X-Y]]. C]

Theorentld (Pereyra.[25 28]): Assume (1) has a unique solution

 

w and Ph is stable on Bh E 8(whm,p) for all hi3h0°

 

81

Let Fh be consistent of order p with P. Then there

exists an h02>0 such that:
(a) thgfio there exists a unique solution X(h)
for the discrete problem Fh(X) = 0, and
(b) the solution satisfies
(6) “X(h)-whch = 0(hp).

Equation (6) can be summarized by saying the solutions are

(discretely) convergent of order p.

Proof: By Lemma 1, Fh is a homeomorphism between its domain
Eh and its range Rh 5 Fh(Bh) thgho. Brouwer s Invariance
of Domain Theorem implies Fh

onto the interior of- Rh and the boundary onto the boundary.

maps the interior of Bh

If v {—88 stability implies that g: []Fh(V)-Fh(whtp) H

h,

and therefore by letting V vary over 68 we see that

h
B(Fh(mhcp),%) ash.

Consistency implies “Fh(whw)“ = 0(hp) and therefore

“Fh(mh®)u * O as h"0- Thus, there exists h such that

oého
R ./ fl .

\\Fh(whcp)\\<c VhJSO and O€B(Fh(whcp). C). But Fh ls

one-to-one and onto and therefore there exists a unique

X(h) CB such that Fh(X(h)) = 0.

h
Stability implies that
“X(h)—when“ g C\\Fh(X(h))-Fh(uhcp) \\ = C]]Fh(whco) \] = 0(hp). [:1

Thus, in order to establish the discrete convergence

of the numerical solution to the theoretical solution it suffices

 

82

to show the numerical method is both consistent and

stable.

The usual technique for proving that Eh is
consistent is to examine the local truncation error at each
grid point ti' The local truncation error Th(ti) is

obtained by applying F to the discretization of the theor-

 

 

h
etical solution and evaluating the result at ti' Formally
(7) _ _
_ m<ti+1) 24(ti)+w(ti_l) - f(t (t ) w(ti+1) w(ti_
‘ 2 i'cp i ' 2h
h
for i=1,...,n-l.

If we assume that f(t,y,z) has M total derivatives
with respect to t on [a,b], then by expanding Fh(wh®)(ti)
in a Taylor series about 0(ti) it can be shown that the

local truncation error ¢h(ti) satisfies

M

 

_ 2k 2 (2k+2) _ ,
(8) ¢h(ti) — F(c;>)]t + 23h [(2k+2): 4 (ti) 92k( ))
1 k—l
+ O(h2M+2),
where ti is any grid point in [a,b]. The functions
92k are obtained from
M 2k _ M. :1 3kf ,
(9) 23h 92k(°) ' 'ET' k }
k=1 k=l ' az t.,cp(t.).cp'(t.))
1 1 1
M 2j .
h (23+l) k
{jEi(2j+l): 9 (ti)] '

by rearranging the right hand side in powers of h.

)

 

83

For instance ,

92“) = 3‘L fig] (ti,cp(ti).cp' (tin'pm‘ti’
and
94") = EL ‘65:] (ti,cp(ti).cp' (tin'pm (ti)
+ ~51- i-f- (343- cp”(ti))2.
az (ti,cp(ti).cp'(ti))

 

Because equations (8) and (9) are well-known we have

not given the details for constructing them. The reader

interested in more detail is referred to [27] or the work

in Section 2.1 of Chapter 3.

By observing that f(cp) :- 0, equaticn (8) implies
that fh(x) is consistent of order 2.

In the next section we establish the stability of the

finite difference operator Fh'

 

84

Section 2. Stability

In this section we will investigate the stability
of the discrete scheme (2) . The fact that (2) is stable in
the uniform norm has been proven by Pereyra [24]. The proof

given by Pereyra uses the theory of monotone Operators and is

 

quite technical in nature. We present here a simple and I
direct proof of the fact that Fh is stable.
F
In what follows we will be considering vectors V a
which prOperly belong to En+l. However we will restrict

these vectors to the (n—l) - dimensional subspace consisting
of all vectors in Eml whose first and last components are
given by the boundary conditions in (2) . That is, all vec-

tors will have the form V = (A,Vl.....Vn_1.B)T. and we

will regard Fh as an Operator on V = (V1.....V _1) 6E

Let U and W be any two vector in E . For

2 S... i Sn-Z the ith component of Fh(U)-Fh(W) is given by

(10) [Fh(U)-Fh(W)]i

 

 

 

 

= h'2(Ui+l-2Ui+ui_l) — f(ti,Ui.Ui+£;Yi'l)
- h'2(wi+l-2wi+wi_l) + f(ti,wi,wi+£;Yi'l)
= h‘2[ (Ui+1'wi+1) ’2 (Ui'wi)+(Ui-1‘Wi—1)

_ [f(ti'Ui'Ui+12:Ji-l)_f(ti’wi'wi+12-:i-l)J

85

Using the Mean Value Theorem for continuous functions

of? two variables (see Widder [33]) and regarding

 

 

 

 

V‘ -V° V. -V-
1+1 1—1 . 1+1 1-1
f(tii'vi' 2h ) as a function of Vi and 2h ,
we can obtain
U0 -U0 w. -Wa
1+1 1—1 1+1 1-1 I

 

= fy(ti,ai.fii)(Ui-wi)+fz(ti,ai,Bi)-

 

 

 

 

 

 

i:
(”i+1'Ui-1 _ Wi+1’Wi—1) I)
2h 2h
whe1:e3
1 = -
( 2a) ai wi + 8i(Ui wi)
and
(121)) B _ wi+1'wi—1 + e (”i+1'Ui-1 _ Wi+1'Wi-1)
i ‘ 2h i 2h 2h

With 0/ 9 <1.
Substituting (11) into (10) we have

(13) [Fhun-Fhfln]i

= h"2 (U. -w 2

1-1 i-l )

)-2h'2(Ui—wi)+h' (U

i+1-Wi+l
1
' fy(ti'o'i'BiHUi‘WiH2h fz(ti’ai'Bi)(Ui-l-Wi-l)
—-:-L-f (t a BHU -W)
2h 2 i' i' i i+1 i+1

h

h‘2[(1+-2—

fz(ti'ai'6i))(Ui-l-wi—l)

86

2

h
+(1"2'fz(ti'0‘i'f5i”(”i+1’wi+1)]°

For i=1 and i=n-l we get similar expressions except the

terms multiplying (UO-WO) and (Un-Wh) are identically

zero in the respective equations because UO 5 W0 and
Un E Wn 'by our earlier convention.
If we write (13) for i=1,...,n-l as a matrix system

we obta in
(14) Fh(U)-Fh(W) = Mh(U.W)(U-W).

VVhere ”H1ULVW is an (n-l) x(n-l) matrix whose rows are

given by
[ (U W)] =h-2(-2-h2f (t O B ) 1----1| f (t on B ) 0 0)
M11' 1 yl'l'l' 221'1'1""°"
[Mh(U,W)]i-—h (0,...,0,l+2 fz(ti,ai,Bi), 2 h fy(ti,ai,Bi),
l-h-f (t a B ) O O)
2 z i' i' i ' ”"'
for i=2,...,n-2, and

'2(O,...,O,1+l‘- f

[Mh(U"m]n-l=h 2 z(tn-l'an-l'Bn-l)’

2
’Z-h fy(tn-l'an-1'Bn-l) ) I

With the nonzero entries for the ith row, 2_;i;;n-2, occuring
in the (i-1)st, ith and (i+l)st positions. The arguments

Qi and Bi are given by (12).

 

87

Taking norms in (14) we have

\\Fh<u>-Fh<w> u = Wm (M) u .

Multiplying both sides of this equation by W, using

the linearity of Mh(U,W), and noting that

 

 

 

fl“ has norm one, we have for all U7=’W
[\Fh (u) -Fh (w) H \[Mh(U,W) (U—W) \\ !.
])U-WH = “II-W)] _
= “Mh(U,W)“%E¥““ El
2 inf (U,W) Z .
\\Z\\=1 “Mb \\

Thus, if “Mh (U,W) Z \] is bounded below by a constant C,

independent of h, for all vectors Z belonging to the

unit ball in En—l we can write

[\Fh (U)—Fh (w) ]] 2 c\\U-w\\.
which proves that Fh is stable.

To establish that ][Mh(U,W)Z[] is bounded below we

PrOceed as follows:

let h0 = 2/K and define

a . =
l

ro|:3‘

fz(ti.ai.Bi)

arua
b. =
1 fy(ti,ai.ffii)

for each i=1,...,n-1.

88

For any h (ho, because of our hypotheses on fy and fz,

and Bi>e>0 for all i, l_<_ign-1.

we see that [a1] <1

In this notation the matrix Mh(U,W) is given by

 

_ '2 2 .
[NLh(UIW)]l, - h (-2-h blll-alloooo-oo)o P
— -2 2 P. -
[M-h(UIw)]i. — h (OI-o-oool+aio’2‘h bi,1-ai,O,...,O) '
r7
for 2;i_<_n-2; and 1"
[ (U M] = h‘2(0 0 l+a -2-h2b ) J
Mh ' n-1- '°°°' ' n-1' n-1 °

Now let Z = (Zl"°"zn-1) be any unit vector

in En"1 with respect to the maximum norm on En-l. This
means that at least one component of Z has absolute value

one and no component of Z has absolute value larger than
one.
Suppose Z1 = 1, then the first component of

Mn (U,W)Z has absolute value
-2 2 -2 2
h ]-2-h b14(1-a1)zz];3h (\2+h bl\-\l-a1] ]z21)

:>h'2(2+h2b "2)Zz))

l

211-2(24-h2b1 -2)

89

We have used the facts that b11>0. [a1\<:1 and \zz]:gl

in making the above estimates.

Similarly.1Lf Zn-l = 1 then the last component of

Mh(U,W)Z has absolute value larger than 6.

Suppose Zi = l where 2_<_i_<_n-2, then the ith

component of Mh(U,W)Z has absolute value

-2-h2b.+(l-a.)Z.
1 1

h—2](l+ai)zi_ 1+1)

1

-2 2
4,h (2+h bi-](1+ai)Zi_l+(l-ai)zi+l[)

\ -2 2
4-h (2+h bi‘()Zi—1+Zi+1)+)ai) )Zi—l‘zi+1‘))

h‘2(2+h2

\ ,
IV

bi’()Zi-1+Zi+1)+)Zi-1’Zi+1)))

h'2(2+h2

IV

bi-Z)

The next to last inequality is valid for the following reason.

Consider the expression

[0+B]-+]0-8] with \0]5;1, \B]g;l.
Squaring we obtain
(15) (0+8)2+2]02—82]+(0-8)2 = 202+282+2\02-821.

If a2 = 82 then this is trivially 3.4 so without loss of

generality we can assume 02 > 82. In this case equation (15)

becomes

 

9O

(]a+B]+]a-B])2 = 2a2+282+202-282 = 402 g_4.

Thus for any two real numbers a and B .with \a]5;l

and [B] :1 we have
[9831+ ]a—B] g 2.

Therefore, for any unit vector Z, some component
of DQIULVHZ has absolute value larger than c. This implies
that UMh(U,W)Z“Z>e and consequently
inf \‘1 (u w)z\\ \ e
uz\\=1‘”h ' 1- '
With this our proof of stability is complete. In

order to obtain the desired result we quite explicitly

assumed that 0<<e<iggu The formal hypothesis given by
Pereyra in [24] is that 0 g_%§u However, he implicitly
assumes that 35- is bounded away from zero in his proof

of stability. Thus our proof is as general as Pereyra's,

and considerably more elementary.

Since the Operator F is both consistent and stable

h
we can apply Theorem 1 to conclude that the discrete scheme
(2) has a unique solution, X(h), which converges discretely

to the theoretical solution, w(t), of problem (1).

 

91

Section 3. The Numerical Method

 

In this section we will summarize the results of
Pereyra [24.27] and Stetter [31] concerning the global
discretization error for the method (2). Using the results
obtained by these authors we will develOp a pullback inter-
polation scheme for solving boundary value problems of the

form (1).

The global discretization error for the numerical

method (2) is defined to be

X (h) - mhcp.

Both Pereyra [24.27] and Stetter [31] have presented a general

theorem which, when applied to our problem, becomes

Theorem 2. (Pereyra, Stetter):

 

Let F and Fh have M+l continuous Frechet deri-
vatives. Then for sufficiently small h the glObal dis-

cretization error satisfies

M
(16) X(h)-whet) = ZIhZRmhekun +0(h2M+2).
k:

The functions ek(t) are independent of h and satisfy the

linear two point boundary value problem

e];(t)-fz(t.cp(t) .cp' (t) )e];(t)—fy(t,cp(t) ,cp' (t) )ek(t) =bk(° ) ,
(17)

ek(a) = 0, ek(b) = O.

The functions bk depend on previous error functions

 

92

. . . . /
e1,...,ek 1’ various derivatives of m, and various Frechet

derivatives of F.

We will prove an analogous theorem in Chapter 3 and
the reader interested in details can refer to either Section 2.3

of the next chapter or to the aforementioned papers.

To set up the numerical method, define stepsizes

b—a

 

hk = —k:1' for k=O,l,...,M and construct the uniform grids
2
k , ._ k+l -
Gk — [ti—a+1hk.1—0,l,...,2 ] c [a,b].
Using the numerical method (2), compute a solution vector

X(hk) on each grid Gk'

The computation of X(hk) for each k requires us
to solve a system of 2k+1—l equations. These equations
will be nonlinear whenever f is nonlinear in either x or
x'. In the nonlinear case we must use a root finding procedure
to solve the system (e.g. Newton's method). In the linear
case, the matrix is tridiagonal and can be easily solved using

an LU decomposition and back substitution (see Isaacson

and Keller [15]).

Once the solution vectors X(hk) for k=O,l,...,M
have been obtained, extrapolation can be performed to Obtain
gM+2) approximation to 0(t) at t==égéa This has

been studied by Pereyra. As was the case for initial value

an 0(h

problems, extrapolation at other grid points does not yield

comparable accuracy.

93

The task of implementing a pullback interpolation

2M+2)
0

GM is simplified in this case. This is due to the fact that

from (17) we now have two pieces of information available

scheme for obtaining 0(h accuracy at all grid points of

to us concerning each of the error functions ek, namely,

ek(a) = ek(b) = 0. Recall that for initial value problems I.—
only one piece of information was readily available concerning
ek. Most of the effort involved in implementing pullback t

interpolation for initial value problems was expended in ;[

 

obtaining an approximation to efi at the initial point.
However, by using the boundary conditions ek(a) = 0 and
ek(b) = 0, we will have enough data points to enable Lagrange
interpolation to yield accuracy comparable to that of our data

(see the discussion in Section 5 of Chapter 1).

Thus a pullback interpolation scheme based on
Lagrange interpolation can be devised for solving boundary

value problems of the form (2).

b-a

At tw=—§—-Ea+ho we solve an (M+l) X(M+l) system
to obtain the solution vector (m(t),el(t),...,eM(t)) ‘with
2M+2

accuracy (0(h ), 0(th),...,0(hg)). Thus we know

0
eM(a) = eM(b) = 0 exactly and eM(a+hO) with accuracy
0(hg). Construct LM(s), the second degree Lagrange inter-
polating polynomial to the above data. Then

_ 2 3 . . . .
LM(s) — eM(s)4-o(ho)4-o(ho). The first 0. term is the

error in the initial data and the second is the error in

interpolation.

94

Notice that we are using a sharper order estimate
for the error of interpolation than was utilized in Chapter 1.
It is well known that when performing interpolation on equally
spaced pointstflmzerror is of order a power of the stepsize
(see Isaacson and Keller [15]). However, in Chapter 1, when
interpolating on the three points a, a-tgy and a+h we
referred to the error as being 0(h3). Actually. this error

is 0((%)3)<;0(h3), so we were quite generous in our error

estimates.

In order to show that pullback interpolation for solving

2M+2
0

grid points of the finest grid, we must use the sharper

boundary value problems yields 0(h ) accuracy at all
estimate for the error in interpolation when interpolating

to e When interpolating to the other error functions

M.
we can again be 'generous' and interpret the error as being
in terms of hO which will be larger than all the other step-

sizes.

The details for constructing the pullback inter-
polation scheme for boundary value prOblems are nearly iden-
tical to those given in Chapter 1. The exceptions are
knowledge about e}; is replaced by ek(b) = O and we have
one more data point each time we interpolate. Therefore.

2M+2) accuracy at

pullback interpolation will yield 0(h
all grid points of the finest grid when used in conjunction
with the numerical method (2) for solving two point boundary

value problems of the form (1).

 

CHAPTER III

THE NUMERICAL SOLUTION OF DIFFERENCE DIFFERENTIAL
EQUATIONS WITH CONSTANT RETARDATION

Section 1. First Order Equations
In this section we will be concerned with difference

differential equations of the form
(1) x(t) = f(t,x(t),x(t-r)), r>0, t>0.

Since r >0, (1) is an equation of retarded type. For
brevity, we will sometimes refer to (l) as a delay differ-

ential equation.

This equation and variants of it occur frequently in
mathematics. For instance, Lord Cherwell is credited by

Wright [34] for using the equation
x(t) = -ox(t-l)(l+x(t))

in the study of the application of probability methods to
the theory of asymptotic prime number densities. Cunningham
[ 7] uses a variant of this equation as a growth model to
describe a fluctuating pOpulation of organisms under certain
conditions. Cooke and Yorke [44] use a variant of (l) as

an economic model for the growth of capital stock.

Because of the presence of the delay term x(t—r)
in (1), it is necessary to know the solution at time t-r
to Obtain a solution at time t. Therefore to solve (l) on

the interval (0,r], it is necessary to specify the solution

95

 

96

on [-r,0] by means of an initial function, 0(t). A
solution of (l) with initial value 0 is a continuous function
which agrees with w(t) for t E[—r,0] and satisfies (1)

for t;>0. The solution for tj>0 will also be denoted by

w(t). I

The usual method of obtaining the theoretical sol-

ution to (l) is called "the method of steps". A continuous

 

initial function m(t) is specified on [-r,0] and the i]

initial value problem

int) f(t.X(t).CO(t-r))

(2)
x(O)

49(0). ogtgr

is solved. Denoting the solution to (2) by ¢(t), the process
is repeated on intervals of length r to obtain a solution,

m(t), defined on [-r,.).

In general, the solution 0(t) will have jump
discontinuities in various derivatives at the multiples of r,
zero included. In terms of smoothness of the solution,
the worst possible situation is that a jump discontinuity
will occur in the (k+l)st derivative of 0(t) at the pohit
kr for k=O,l,... . By adding conditions on 0(t) for
t.e[-r,0] these jump discontinuities can be avoided to
some extent. For instance, if we require that the left hand

derivative of w(t) at t=0, w'(0-), exists and satisfies

97

$‘(0-) = f(0.X(0).w(-r)):

then the jump discontinuity at kr is in the (k+2)nd
derivative for k=O,l,... . By imposing similar conditions
on the initial function we can force the jump discontinuities

to occur at whatever derivative of the solution we desire.

For f to have n continuous total derivatives with
respect to t on (O,r] it is necessary that m(t) be n
times continuously differentiable on [—r,0]. In this case
the solution will be n+1 times continuously differentiable
on (O,r]. However» without the added conditions discussed
above on the left hand derivatives of 0(t) at t=0 there
will still be a discontinuity in the first and all higher

order derivatives of the solution at t=0.

Note that as we proceed to the right using the method
of steps to solve (l) we gain differentiability of the
solution provided we do not lose any differentiability of f.
That is, if f is n times continuously differentiable with
respect to each of its arguments on L—a,m) and the initial
function is n times continuously differentiable on [-r,0],
then the solution to (1) will be n+k times continuously

differentiable for t €((k-l)r,kr], k=l,2,...

In order to insure the existence of a unique solution,
0(t), to (1) which is continuous on [-r,W) we will make the

following assumptions (see Hale [13]):

 

98

(i) the initial function m(t) is continuous on [-r ,0];

(ii) f(t,u,v) is a continuous function of each of its

arguments:

(iii) f(t,u,v) satisfies uniform Lipschitz conditions

with respect to u and v.

Under these hypotheses we can also show that the sol-

ution to (1) depends continuously on the initial function

0(t) (see Hale [13]).

In order to obtain a numerical solution of (l) we
shall require more stringent conditions on f and the initial
function. These conditions will be stated as they become

necessary.

To solve (l) numerically, given an initial function
m(t) on [-r,0], define F(t,x(t)) s f(t,x(t),m(t—r)) and

numerically solve the initial value prOblem

x(t)

F(t.X(t))
(3)

MO) 42(0). ogtgr.

To solve (3) numerically we select a basic stepsize
h in such a manner that r is an integer multiple of h:
i.e., r = Nh for some integer N:>O. The reason for this

restriction on h will become apparent later.

 

99

As in Chapter 1 we define stepsizes hk = h/2k
. _ k_. .._ k
and grids Gk — [ti—lhk. i—O,l,...,2 ] for k=O,l,...,M.
We assume that the numerical method employed to solve (3)

is such that an asymptotic expansion of the form

k _ k M k jq (M+1)q
(4) X(ti'hk) — 428:1) + 3:31 ej(ti)hk +001k )

. . k k .
is valid at each ti for every hk. Here, X(ti'hk) is
the numerical solution to (3) at the grid point t: obtained

with stepsize hk and 0(tt) is the solution to (l) at t:.

The existence of an expansion (4) will, in general,
require the existence of (M+l)q continuous derivatives of
the solution w(t). This wil be assured if we assume that
f(t,u,v) has (M+1)q-l continuous derivatives with respect
to each of its arguments and that the initial function ¢(t)

is (M+1)q-l times continuously differentiable on [-r,0].

If the numerical method used to solve (3) is such that

q=1 or 2, then we can proceed as in Chapter 1 to obtain a

numerical solution X(t) which satisfies
X(t) = cp(t) + 0(‘h(M+l)q).

This relationship will be valid for all t belonging to the

finest grid, GM'

The presence of jump discontinuities in derivatives

of the solution at multiples of the delay r does not affect

 

100

the applicability of the pullback interpolation method.

The reason for this is that the pullback interpolation method
employs derivatives that are computed using the differential
equation at the initial point. Since the differential
equation is valid only to the right of the initial point all .
derivatives considered in Chapter 1 are right hand derivatives l—J
at the initial point. If f(t,x(t),x(t-r)) has M total

derivatives with respect to t and ®(t) has M right —_45

 

hand derivatives at t==—r, the solution of our delay differ— Li
ential equation will have M right hand derivatives at zero

and therefore at all other multiples of r also. These

derivatives can be computed by differentiating the differen-

tial equation in the same manner as was done in Chapter 1.

Repeating the solution procedure for the problem,

x(t) F(t,x(t))

X(h)

XUU. hgtgzm

(M+1)q)

we can obtain an 0(h solution to (l) at the points

{h+ihM: i=0, 1, . . . , 2M] .

If we repeat the above procedure N times we will

have a computed solution, X(t), satisfying
(5) X(t) = wt) +O(h(M+l)q)

for all t such that

101

t.€[jh+ihM:i=O,l,...,2M, j=0,1,...,N-l]=[ihM:i=O, ,...,N2M}

C [O,r].

Let's examine what happens when we now try to obtain
a solution to (l) on [r.2r]. Suppose t is any grid
point in [r,2r]; because GkCGM Vk, t can be represented
as t=r+ihM where 05;13;N2M. Since t.€[r,2r], t-r €[O,r]
and in order to solve (l) we must be able to evaluate
0(t-r). Using the above representation of t, we can write
M

t-r = r+ihM-r = ihM with O_<_igN2 . This is a grid point

in [O,r] and by (5) we have

(6) w(t-r) = x(ihM) + 0(h(M+l)q).

It is easy to see that there is a 1—1 correspon-
dence between grid points in [r,2r] and grid points in
[O,r]. This fact and equation (6) emphasize the importance
of the pullback method for obtaining numerical solutions to
difference differential equations. No numerical method is
any more accurate than the data it receives and the solution
of (l) on [r,2r] requires the solution on [O,r] as data.
The uniform accuracy guaranteed by (6) means that the entire
process used to obtain the solution on [O,r] can be repeated

to obtain a solution on [r,2r].

Because normal extrapolation without pullback
interpolation produces a solution whose accuracy varies

from grid point to grid point, it is not an effective pro-

cedure for this problem.

'm —-‘
.

 

 

102

Extrapolation has been used to obtain solutions to
problems of the form (1). Feldstein [10] has developed
two algorithms based on Euler's rule that have (partial)

asymptotic expansions of the form
(7) X(t) = co(t) + e(t)h + 0(h2).

An expansion of the form (7) justifies the use of one

extrapolation.

Cooke and List [3 ] have develOped an algorithm
which they call the "modified Euler-Feldstein" algorithm
which can be used to solve delay equations with non-constant
retardation. They indicate that a sketch of a proof that
their method has an eXpansion of the form (7) has been
supplied by John Hutchison. However, neither the proof

nor its sketch is contained in [3 ].

When using the "Euler—Feldstein" or the "modified
Euler—Feldstein" algorithm the basic stepsize h is taken
to be quite small (usually r/l6 or smaller) and the
prOblem is solved on several grids with stepsizes h/2k.
Extrapolation is performed at all multiples of the basic
stepsize. The retarded term x(t-r) is obtained by using
the solution in the previous interval that was computed

using the same stepsize with which one is now working.

This procedure appears to work well in practice

but is computationally quite expensive since the basic

 

103

stepsize is so small. Also, the existence of an eXpansion
of the form (7) is not sufficient justification for using

more than one extrapolation.

Of course, systems of delay equations can be solved
numerically by generalizing the above procedure in the I

obvious manner. A more interesting prOblem is the equation

x(t) = f(t,x(t),x(t-rl),x(t-r2), ...,x(t-rn)) .

 

where riZ>O is a rational number for i=1,...,n. The

same numerical procedure will work for this problem provided
the basic stepsize is chosen in a manner which insures that
the information required by the delayed terms is available.

Without loss of generality, assume the delays r. are

i
ordered : ogrlgr 2 g . . . . g rn. Let d be the least common
multiple of the denominators of the ri for i=1,...,n,

then with h=l/d this equation may be solved in the same
manner as (1). This choice of the basic stepsize insures
that each ri will be an integer multiple of the smallest
stepsize employed and therefore t-ri will be a grid point
at which the solution is known Vi whenever t itself is

a grid point.

We close this section with two numerical examples.

Example 1: Solve x(t) = —x(t-l), 05;t5;3, numer-

ically for the initial function m(t) = t2 on [-1,0].

104

Here, the delay is r=l and the theoretical solution is

 

 

 

 

given by
r 3
-Lt-l3) -1 ogtgl
(t-2)4+4t-9
mu)=< 12 igtgz
-(t-3)5-10t2+65t-96 , /
k 60 Zet\3

Note that 0(t) is given by a polynomial on each subinterval.

To compute a numerical solution on [0,3], we first
solve the equation numerically on [0.1] using the trape-
zoid rule, extrapolation and pullback interpolation with

Ms3. The finest grid G contains nine equally spaced

3
points in [0,1]. Denote the computed solution at these
nine points by Xl(t). Using X1(t) as the (discrete)
initial function we next solve the same equation on [1,2]

in the same manner to Obtain X2(t). This procedure is

repeated a third time to obtain X3(t) on [2,3].

The computations of ei(t), eé(t) and eé(t) at
the initial points 0,1 and 2 are greatly simplified
for this problem. One factor contributing to this simpli-
fication is that the Jacobian J(t). is identically zero.
The second factor is that all derivatives of the solution
m(t) at the points 0,1 and 2 can be eXpressed in terms

of the derivatives of the initial function cp(t)=t2 at

 

105

=-1 and the computed solution at smaller multiples of

r=l. In fact, by noting that for any k=O,l,2,... and for

any j.>k we have x(j)(kr) = x(j-1)((k-l)r) it is easily

seen that x(3)(kr) = x(J-k-1)(-r). While, if jjgk, then

x(j)(kr) = x((k-j)r).

The numerical results are presented as Table L3. The
error given is the computed solution minus the theoretical

solution. Each column of Table L3presents the errors at

 

t t . . t contained

0' 1’ ' ’ 8
in the interval indicated by the heading of the column.

the nine equally spaced grid points

TABLE 13
i Error on [0,1] Error on [1.2] Error on [2.3]
0 0.0 0.0 x 10"13 -1.5 x 10'13
1 0.0 x 10'13 - .3 X 10"13 0.0 x 10"13
2 0.0 x 10'13 - .4 x 10'13 .2 x 10"13
3 0.0 x 10"13 - .5 x 10"13 .3 x 10"13
4 0.0 x 10’13 - .6 x 10'13 .4 x 10"13
5 0.0 x 10"13 - .8 x 10"13 .5 x 10'13
6 ' 0.0 x 10‘13 -1.0 x 10"13 .6 x 10"13
7 0.0 x 10'13 —1.3 x 10'13 .7 x 10"13
8 0.0 x 10'13 -1.5 x 10"13 .8 x 10'“13

due to the limitations of machine accuracy and the particular

The errors reported

in Table L3are almost entirely

106

routine used to compute the Hermite interpolating polynomials.
Only the first twelve digits can be absolutely guaranteed

to be accurate and to twelve decimal places all errors are
zero: This is not particularly surprising, since the trape-
zoid rule and extrapolation yield extremely good results I
when the theoretical solution is a polynomial, as it is in —1

this case.

 

Example 2: Solve x(t) = -x(t-l), 05;t§;3 numerically for ;i
the initial function 0(t) = et on [—l,O]. The differen-
tial equation is the same as that in Example 1. The theor-

etical solution for this initial function is given by

 

{-e(t"1)+1+e":L 0_<_t-_<_1

wt) = ( e‘t'2)-(l+e‘l) (t-l) 1_<_tg_2
(t—3) -1 t2 -1

k—e +(l+e )T-2(1+e )(t-l) 2_<_t_<_3

In this case the solutions are not polynomials.

The same numerical method was used to compute the
solution as was used in Example 1. The results are presented

as Table 14.

107

TABLE 14
1 Error on [0,1] Error on [1,2] Error on [2,3]
0 0.0 -1.2 x 10‘10 - 7 2 x 10"10
1 -5.3 x 10"10 4.5 x 10’10 -12 8 x 10'10
2 -3.2 x 10-10 3.0 x 10’10 -11 2 x 10'10
3 2.3 x 10’10 -2.6 x 10"10 - 5 6 x 10"10
4 .7 x 10"10 -1.2 x 10 1° - 7 0 x 10 10
5 -2.4 x 10"10 2.1 x 10"10 -10.2 x 10"10
6 12.0 x 10"10 —12.7 x 10'10 4.7 x 10"10
7 42.4 x 10—10 -46.7 x 10 10 38 9 x 10"10
8 —1.2 x 10‘10 - 7.2 x 10'10 .3 x 10"10

Since the theoretical solution is not a polynomial,
the computed solution is not as accurate as the computed
solution in Example 1. However, we do Obtain eight reliable
decimal places at all grid points, which is comparable to
the accuracy that extrapolation yields at the endpoints

1, 2 and 3.

Also, it should be noted how stable the numerical
results are. There is very little deterioration in the
accuracy of the computed solution as we move to the right.
Thus, extrapolation and pullback interpolation appears to

be a viable solution technique for first order delay equations.

 

108

Section 2. Second Order Equations
For the remainder of this chapter we will be working
with the second order difference differential equation of

retarded type
(8) x(t) + f(t,x(t),x(t-r), x(t),x(t-r)) = 0. r50. t>0.

By writing (8) as a system of two first order equations
and imposing conditions (i), (ii) and (iii) of Section 1
or the more general conditions of Sansone [30] on this system
we can insure the existence of a unique solution to (8) which
depends continuously on the initial data (see Hale [13] and

Norkin [23]).

In either case, f is assumed to be a continuous
function of its arguments. For our purposes it is enough to
assume the existence of a unique solution to (8) for tI>O
whenever we are given a continuously differentiable function

m(t) on [-r,O].

To solve (8) theoretically, the method of steps can
again be used. The analysis of the behavior of the solution
to (8) is completely analogous to that given for first order
equations in Section 1 of this chapter. Accordingly, we will
not go into this in detail. When specific aspects of the

theoretical solution become important they will be mentioned.

If (8) is written as a system we can of course use

the method discussed in Section 1 to solve this system.

 

“use

109

The rest of this chapter will be devoted to an investigation
of a direct method for solving (8) which does not involve

reducing it to a system of first order equations.

 

We introduce the notation I
‘1‘?!
(8') F(x) = O
for the continuous Operator defined by (8) and we will refer .5
to (8) as the continuous problem. ;[

For ease in referring to partial derivatives we will

refer to equation (8) as
x(t) + f(t,u,v,y,z) = O.

In setting up a discrete analog of (8) we actually
get two slightly different problems; one when we are trying
to obtain the solution on [O,r] and the second when we seek
a solution on an interval of the form [(L—l)r,Lr] with
L:>l. The version of the problem for the interval [O,r]
is a special case of the version for the general interval
[(L-1)r,Lr] and will be explained in the subsection on
implementing the algorithm (see section 2.4). Thus until
further notice we will concern ourselves with the discrete

analog of (8) on the interval [(L-l)r,Lr] with L:>l.

Let R2>l be any given natural number and define

. _ . R
h — r/R. Form the uniform mesh [ti-(L-l)r+ih]i=_(R+1).

110

We assume the solution has already been obtained at the
grid points ti for '=-(R+l),-R,...,O and we denote this

solution by w(ti).

Setting Xi = X(ti). the discrete problem corresponding

to (8) is given by

 

 

Xi+1“2Xi+Xi-1 +f(t x x Xi+1’Xi-1 Xi+1-R‘Xi—1-R) _ 0
h2 i' i' i-R' 2h ' 2h ‘ '
(9) i=0,1,...,R-17
X . = t .), '=O,l,...,R+1.
_3 w( _3 J

The discrete problem (9) will be denoted by
(9 ) Fh(X) = 0.

Since )g_. is known for j=O,...,R+l, this may be thought
of as a (in general) nonlinear system of equations in ER

with the unknown being the vector (X1,...,XR)T. The solution

to (9) will be denoted by X(h).

Since each equation in (9) taken in the order
i=O,l,...,R-l introduces only one unknown, we can solve the
system uniquely by forward substitution. In general, the non-
linear nature of f will necessitate the use of a root finding

procedure (e.g. Newton's method) to solve each equation.

We will adopt the same space discretization operator
wh that was utilized in Chapter 2. That is, for any

continuous function z(t) on [—r.”)

 

111

T
(2(t1).....2(tR))

:5
N
A
I:
II

)T.

(21,...,zR

We will continue to use the definitions of discrete

convergence, consistency and stability (definitions 1, 2

 

and 3 respectively) that were adOpted in Section 1 of fl_
Chapter 2.
The proofs of Lemma 1 and Theorem 1 in Section 1 :3”!

 

of Chapter 2 are valid in this setting. In fact, Pereyra
mentions in [28] that these results are not truly tied up

to the two—point boundary value problem.

In order to establish the existence and discrete
convergence of a unique solution X(h) to (9) we will

again show that (9) is consistent and stable.

In order to prove that Fh is stable we have to make

some restrictive assumptions about the equation (8). These
assumptions are necessary for our proof of stability to be

valid. However, we believe that F is stable for a much

h

larger class of functions f(t,u,v,y,z) than what we are

able to establish.

 

112

Section 2.1. Consistengy,

 

Let 0(t) denote the theoretical solution to the
continuous prOblem, (8). If we apply the operator Fh
defined by (9) to the discretization of cp(t) we obtain

what is known as the local truncation error, T

 

 

 

h u.
Formally q._J
Th(ti) E Fhmhw) (ti)
9(t. )-zw(t.)+w(t. ) . ’1
= 1+1 1 1-1 ‘3
cp(ti+1)""’(ti--1) m(ti+l-R)-m(ti—l-R))
2h ’ 2h
where ti is any grid point i=1,...,R.

If Th(ti) = 0(hp) Vi then we obviously have
that our method is consistent of order p. Instead of
investigating (10) directly we let x(t) be any sufficiently
differentiable function and investigate Fh(“h¥)(ti)' If
we make the "localizing" assumption that no previous
truncation errors have been made, we may use Taylor's Theorem

repeatedly to obtain an asymptotic expansion for

Fh(mhx) (ti) .

Let ti be any grid point, i=1,...,R. If

x(t) has Ml+2 derivatives at ti we can write

2
(11) x(ti+1) = x(ti)+hx'(ti)-+%Tx"(ti)
M +1
1 j . M+2
+2 l'1.---.-x(3)(t.)+c}(h l )
j=3 j. i

113

 

 

and
1 hz
( 2) X(ti-l) = x(ti)—hx (ti) +-2—:-x"(ti)
M +1 . .
l _ j j . M +2
+ Z) _L_1.l'__11_ x(j) (t.)+0(h 1 ).
._ j. 1
3—3
Using (11) and (12) we have
(13) x(ti+l)-2x(ti)+x(ti_l)
2
11 M
[71’] 2k 2[¥l]+2
=h ”[h x"(t. HZ gT<-,-x(2k)(t.)+o(h 2 ))
. i
k= 2
K 2k-2
= x"(t. )+2 2255-" (2") (t. )+G(h2K
k= 2 '
M1
where [-] is the greatest integer function and K = [3?].
Expanding f(°,',°,y,°) in a Taylor series about
x(t. )-x(t. )
X'(ti) where y = 1+12h 1-1 we obtain

(14) f(‘p'o'pYo°) = f(‘l'1'1X'(ti)l°)+fy(°I°I'oX'(ti)l')(Y-x.(ti))

 

(y—x'(ti))2
+fyy(°’°'°’x'(ti)'°) 2:
M2 j .
+Z___a.—f (y-x'(ti))]

J- J u
j=3 BY [(‘I'I'IX (ti-)0.)

M2+l

+O[(Y-X'(ti)) ]o

'where M2+l is the number of partial derivatives of

f(t,u,v,y,z) with respect to y which exist.

[_flfll'

 

W;__‘“““‘

114

Using (11), (12) and the definition of y, we obtain

K h2j+l

1
= 2h[2hX' (t. 1)+2jE£ 723:1TT (ti )+O(h

(15) y-x'(ti) — 2h {x(tfl

2K+1 .
M} (ti)

K 2j

2 h
j=1 (2j+l):

x(2j+l)(ti)+o(h2K).

Substituting (15) in (14), we have

f(".'.'Y'.) = f(‘o'0'0X.(ti)o')+fy(°r°o'IX.(ti)I').

K 21 -
h (23+1) 2K
[1:31 (2‘j+1")': X (tiHMh H
K 23'
+%_f (°,°,°,X'(ti)p°)[ Z (2?+1):X

(21+1)
(t. )+0(h
yy j=1 1 '

2
+21—1—5f 1.
. .l.IX'(ti)l.)

(2j+1) 2k)}k

(ti)+0(h

2j . M +1

K 2j . M +1
h 2K 2
BUt ”({jE] (2j+1): X )1 )

M +1
0( 0(h 2)+0(h2K) } 2 } = 04h ), and for any i=1,...,M2.

we have, from the binomial theorem,

 

2K)]2

115

K 2j 0 K 2j o
h (23+1) 2K 2 _ h (23+1) z
{jEi (2jt1): x (tiHMh )} ‘ (jEQ (2j+1): X (ti))
+0(h2K) .

Consequently we can write

(16) f(.’.'.,y'.) = f(.,.,.,x'(ti),.)+fy(.,.,.,x'(ti),.) -
K 2j .

h (23+1)
jg: (2j+1): x (ti)

M k
k:
k=2 BY (°'°'.'X'(ti)l.)

N

K 2j
Z: h

(2j+1) k
[j=1(2j+1): x (ti)}

2M2+2

+0(h ) + 0(h2K) .

Next, expand f(.,.,.,x'(ti).z) and its partial derivatives
with respect to y as functions of z in Taylor series about

x' (ti-r) and then set 2 = 1+1-R 2h l-l—R to obtain

(17) f(‘o'I'IX'(ti)Iz) = f(‘l’r .IX'(ti)IX'(ti-r))

+ fz('a '1 .IX' (ti)oX' (ti-1.)) (z-X. (ti-r))

M

 

3 L
+ 21.17: 5
i=2 ' 02 (-,-,-,x'(ti),x'(ti-r))

(z-x'(ti-r))‘

M3+1
+O([z-x'(ti-r)] ).

 

 

116

In the same manner, for m=l,...,M2, we can obtain

m emf
(18) = m
Oy ('o’p'ax.(ti)oz) BY ('I'o°oX'(ti)IX.(ti-r))

 

am+lf
+ m (z-X'(ti-r))

BzBY [(4 w n XTti).X'(ti-r))

M

3 1 am+£f ' 1
+52 7::- azlaym‘(.,.'.,xu(ti)'xn(ti_r))(z-X (ti-r))
M3+l

+G([z-x'(ti-r)] );

where M3+1 is the number of partial derivatives of f(t,u,v,y,z)

with respect to 2 which exist.

Since ti-r E ti-R' an expansion similar to (15)

 

O
is valid for z-x'(ti-r). Introducing the notation gyf
for the function f we can, in direct analog to (16), rewrite

(l7) and (18) as

m

a f amt
(19) — I = —__ u I
aym ('I'I°IX (ti)lz) Oym ('o'o°ox (ti)'x (ti-r))
M
3 m+£
+ Z){i%.ilisl%] 1°
b=l ‘ 62 By (°,-,-,x'(ti),x'(ti-r))
K 2j .
h (23+l) z
{3.31 72—3757 X (ti-r”
2M +2
+001 3 ) + 0(h2K).

for m=O,1,...,M2.

 

 

(20)

117

Combining equations (16) and (19), we have

f(',',°,Y,Z) = f(-,-,-,x'(ti),x'(ti-r))

M l. K 2

 

 

 

 

 

3 j .
1 a f - h (21+l) z
23 --.- { —---—rx (t.—r))

L11 1” 02‘ ('I'l'lx.(ti)lx.(ti-r))j=l(23+l). 1 1
2M +2
0(h 3 > + 0(h2K) 1

K 23 . I :
5f h (23+1)
fi('l'o'lx'(ti)lx'(ti-r))j§1(2j+l): X (ti) ‘3!

M3 1 1+1.
( Z:[7F'0 f 1'
”=1 az‘ay ( .-.-.x-(ti>.x'(ti-r>)
K 2j . K 2j -
h (23+l) _ L . h (23+l)
[E3 (2j+1): X (ti r” 7 Z: (2j+1): X (ti)
3-1 j—l
2M +4
0(h 3 ) + 001””)
M M
2: lug i akuf -
. [ . ]
k=2 1" L20 ‘° 1 X]
Oz BY (.I.I.IX'(ti)IX'(ti-r))
K 2j . K 2j .
h (23+l) L h (23+l) k
[‘13 - .X (t.-r)])'(Z————-X (12.)) }
j=1z23+1)' 1 j=l(2k+l): 1
2M +2 2M +2
0(h 2 ) + 0(h 2 ) + 0(h2K).

Combining equations (13) and (20), and evaluating

all derivatives at (ti,x(ti). x(ti-r),x'(ti),x'(ti—r)) we

can write Fh(uhx)(ti) as

118

X(t )- 2 (ti )+x(t.

 

 

 

 

Fh(mhx)(ti) a “1 hz “’1 +f(t. xi(t ) x(ti_ R)
X(ti+1)‘X(ti—i) X(ti+1- R) X(ti- 1-R ))
2h ' 2h
2k- 2
= x"(ti)+21:82};---2--———k).x(‘7‘}‘)(ti)+c3(112K
+ f(tilx(ti)ox(ti‘r)IX'(ti)oX'(ti-r))

M

3 L 2] .

;L_§_£ h (23+1) _ L
+ 1:31 L: 321 (321 (2j+1)’ (ti X”
M M3 .
2 k+ L 23 .
l 0 f h (23+l)
+ Z —.-{[ Z -——-(Z . .x (t.-r))]
k=1 k' [=0 azi'ayk j=1 (23+1)
K 2j
h (23+1)
[:13 (2j+1)’ (t H )
j—l
2M +2 2M +2
+ 0(h 2 ) +0(h 3 ) +0<h2K
Define M = min(K-1,M2,M3); then, reindexing, we have
2k- 2
(21) Fh(th) (ti) = 19(th )+2 3212—157" x‘ZX’ (ti)
M L M 2j

J. a f h (23+1) L
+ Z —T— (Z ' I X (t 17))

L=1 L 62L :1 (23+1) i
M M k+L M 21'

l 0 f h (23+1) L
+z—.—{[2 (23 . . (t rHJ
k=1 k [:0 azLayk j 1 (234-1

M 2j .

h (23+l) k

[151W X (tl)] }

 

119

If we regroup (21) in powers of h we obtain

M O
I _ 2] (’) 2M+2
(21 ) Fh(whx) (ti) — F(x) (ti)+j_2;3l gzj x,ti+0(h )

I

where the arguments of the functions gzj(-) involve various

derivatives of x evaluated at ti and ti-r, and various

partial derivatives of f(t,u,v,y,z) with respect to y

and/or 2 evaluated at (ti,x(ti-r),x'(ti),x'(ti—r)). For

instance

 

92") = 113x”) (ti) +33:— X(B) (ti-r) +3.5. X(3) (ti))
and
2
g4(°) = 'é'ZTX(6) (ti) +31..— 1 2 5; (X(3)(ti-r))2
° ° (3:) 82

To obtain the local truncation error for the discrete

prOblem (9), by (10) we need only take x(t) to be the

theoretical solution m(t) of problem (8) and use either

equation (21) or (21'). Thus,

M .
_ 2] . 2M+2
(22) Th(ti) — F(Cp) (ti) +j§lh 92j( ) cp,t.+0(h )

from equation (21').

To establish the consistency of our discrete Operator

we need only note that F(w) E 0. Consequently

 

 

120

th(ti) = 0(h2) for every grid point ti, i=1,...,R.

our discrete Operator is consistent of order 2.

Thus

 

 

 

121

Section 2.2. Stability

 

 

In this section we investigate the stability of the

discrete scheme (9). The first part of the analysis is done

 

for the general method and is analogous to the work done in

Section 2 of Chapter 2.

3‘2

In investigating the stability of the discrete

problem (9), it is convenient to view Fh as a function .

from ER to ER. The domain of Fh will consist of vectors 5] “a

_ T R
V — (V1,...,VR) (SE . Of course, Fh

the initial conditions

 

continues to have

v_j = o(t_j), j=O,l,...,R+1.

When considering the jth component of Fh(V) we will
suppress the (constant) first, third, and fifth arguments of

f(t,u,v,y,z) and write

 

f(t V V 31:51.2. Vj-R-Vj-z-R)
j-l' j-l' j-l-R' 2h ' 2h
v.-v._2
= f('rv' I'I _J___J____’ ')
j-l 2h

Note that the jth component of Fh(V) is obtained from the

(j-l)st equation in (9) for j=l,...,R.

Let U,W ‘be any two vectors in ER with Uaiw.

l”

The cases j=1 and j=2 are considered individually

h

later. For 3gj=_<_R the jt component of Fh(U)-Fh(W) is

122

given by

(23) [Fh(Y)-Fh(W)]j

= h-2(U -2U +U )+f(- U - 21:31:; -)
j_2 j_l j I j_lI I 2h I
_h-2 (w -2W +w )+f ( , w , Yfli .)
j_2 j-1 j I j_lI I 2h I

_ -2
— h [(U. -W._2)-2(Uj_1-Wj_l)+(Uj-Wj)]

 

3-2 J
U.—U._2 w.—w._2
+f(°IUj_lI°I 2h I.)-f(.le_ll.l 2h I.)°

Proceeding as in Chapter 2, we use the Mean Value
Theorem for continuous functions of two variables (Widder [33])

to obtain

(24) f(.yu_ ";L_j:g..)_f(.fiq. . _j__i:Z .)

J-l’ 2h I 3.1! I 2h I
= fu( Iajl IBjI ')(Uj_1-WJ_1)
+f ( a , B ,)(Eifli_Zj_—XJ:£)
y ' j' ’ 3’ 2h 2h '

where
25 . = w. + . U. —W
( a) a] 3-1 83( 3_1 _1)
and

W.-W._-2 U.-U._2 W.-W._2
(25b) 6 = _J__J_+ 9.(_J___L.___.J_.J—)'
3 2h. 3 2h 2h

Substituting (24) into (23) we have for

nggn

123

(26) [15h (U)—Fh (w) ]j

= h—2

1 1
+ f ( Ia'l 06°: )(Uj-Wj) _2th(.'O’j'.'Bj'°)'

(Uj_2-Wj_2)
2h

h

a h“2[ (1--2-fy(-.a '°'Bj’ '” (”j-2‘Wj-2)

J'

)

+ (‘2+h2fu(' Iajl 'fijl ' )) (Uj_1-w

j-l

h

+ (1 +§fy(-,a ..8j. .)) (Uj-Wj)].

1'
For j=1 and 2 we can reason in the same manner

to obtain

(27) [icm—Fhm 11 h’2 (1 +§fy(o.a1.-.81..)) (Ul-Wl)

and

-2 2
(28) [Fh(U)-Fh(W)]2 h [(-2+h fu(-,02,-,f32,-)(U1-W

l)

h
+(1 +§fy(.0a2I 'IBZI ’) (UZ‘W2)].

The functions 01,81,02, and 62 are given by

(29a) 91 = ”(to) 5 WC.
W -W U -W
_ 1 -1 1 1
(30a) (12 = Wl + 92 (Ul—Wl).

 

 

 

124

 

and
W -W U -W
_ 2 0 2 2

with O< 81<l and O<Bz<l.

Writing (26), (27) and (28) as a matrix system we

Obtain

(31) Fh(U) - Fh(W) = Mh(U.W) (U-W)

 

where Mh(U,W) is an R.xR matrix whose rows are given by:

[Mh(U.W) ]1. =h'2(1 +3§fy(-.01.-Bl.-),0,0,...,0);

(32) [M}1(UIW) ]2. =h-2 (-2+h2fu(°Ia2I 'BZI ') I1 +-}zlfy(.la20 .le ’)I

_ -2 .9 , , ,
[%(U'W) ]j. h (0' O 0 0,0,1 2fY( 'aj’ Bj' )'
2 h
"2+h fu( Iajl IBjI')I1+2fy('IajI°IBjI‘)I

O,ooo,0)

for 3gj_<_R:

with the non-zero entries for j_23 occurring in the (j-2)nd,

(j-l)st and jth positions. The arguments aj and Bj for

j=l,...,R are given by the apprOpriate equation (25). (29)

or (30). Note that since R = r/h the dimensionality of

Mh(U,W) depends on h.

 

 

125

Since Mh(U,W) is not diagonally dominant, as it
was forboundary value problems, we are not able to prove
directly that Mh(U,W) is bounded below. However, we do

have the following equivalent formulation of the concept of

 

stability. "I

Suppose Mh(U,W) has an inverse Mil(U,W) that is

bounded above in norm, independent of h. Then we can write y

(3]J as 9] ‘fi

U-w = Mgl (U,W) (Fhun -Fh(W) >.

 

Taking norms and using the fact that the norm of a product

is less than or equal to the product of the norms, we have

“U-W“ S \WIWM) \) \\Fh(U)-Fh(W) \\ S C\\Fh(U)-Fh(W) \\

 

where C is the bound on “MQI(U,W)“. This inequality is
precisely the definition of stability given in Section 1 of

Chapter 2.

The norm inequality we have used is valid as long as
we use the matrix norm induced by our vector norm. Since we
are working with the maximum vector norm on ER. the induced
matrix norm is the maximum absolute row sum (see Isaacson and

Keller [15]).

We are also unable to prove directly that
“M;1(U,W)“ is uniformly bounded above. However, we are able

to establish this fact with some additional hypotheses on f.

126

For notational convenience we will assume that r=l.
This is a rescaling of the delay interval and has no effect
on any of our results. We will also write the dimensionality

of the matrix system (34) as n ER=l/h.

The special case for which we will prove stability F1

is given by the following assumptions:

 

ii = ...;

By ‘ 0' +.
(33) .

—2 < -K<-§-—E<—O<O,

where 0 (EC are positive constants.

Write bi = fu(°”1i'°'Bi'°)' Then the assumptions

(33) imply that bi is negative for all i: in fact,
(34) —K < bi < -0.

Rewriting the matrix Mh a Mh(U,W) given in (32)
in terms of bi and n and imposing the assumptions (33),

we obtain

 

127

 

 

 

 

”- 1 . o . o , o . 07
b2
-2 +7, 1 I O I O p o
n
b3
]- I "’2 +-—_l ]- I O I
2
n
_ 2 b
M'h-n O I l I ‘2 +—'4-I l I
2
n
l, ’2 + 2‘1, 1! O
n
bn
b O O, 1' -2 +7; 1 o
n .J
Define
. -_1_ .
O I 4! O I O o O
-l
O , 0 , Z, O
1 -1 °
Sh- 2 O I O I O I Z
n
° 0
-1
O I O I z
n n-1 n-2 n-3 2 l
L_4n‘2 4n2 4n 4:12 4n2 4n2j

 

 

 

 

then
r -1---—b‘°‘
2 4n2
-_1.
4
O
Sth =
O
‘ 4n4
write
(35) Sth =
with
F1 b2
3*“?
4n
1_ .
4
Th = O I
—(n--l)b2
. 4n4
where I is the

h

128

 

z I O ‘ ‘ O .
L331. -l
2 4n2 4 .
-l ° _1_-_.b_n_
4 2 4n2
n- b
( 2) 3 . 2bn-l bn
4n4 4n4 4n4 '
.. Th
1
Z I O o o o o o o
1.33.. l
2 4n2 4
. l l+iri .
( 2 b 4 2 4n2
n- ) 3 -2b -1 -bn 1
4 p o I _—_' -
4n 4n4 4114
n xn identity matrix.

I
bud O

" l
5 0‘
A)

 

 

 

 

129

Because of (34), for n sufficiently large (h

sufficiently small) we have

 

  

0 i’“ 1 6'33
4n 4n
for i=2,...,n. Therefore the first through (n-l)st rows
of Th all have absolute row sum less than l--2:-.
4n

All entries of the last row of T are positive,

h
so the absolute row sum of the nth row is
-l
1 1 n
1--—-———-—— Zib .
4n2 4n i=1 “+1'1
n-l
l l
=l———+-—-—Zi\b .\
4n2 4n4 i=1 n+1-1
n-l
_<_l--l§+'—lZ'K Zi
4n 4n i=1
=1____]_.§_+Knn:ll).
4n 8n
Let 6 = 2-KZ>O, then
4n 8n 8n
>—l-2-(2-K)
8n
=_§._
8n2
Therefore,
n-l
1 1 - 1 Kala-ii
1--2--—-z>?lbn+1-i.<.1- --2-- 2
4n 4n i=1 4n 8n
6

 

130

Thus,

5
\T H _<_ max (l--g-—, 1—-—-)
\ b 4n2 8n2

= 1..

 

- 6
min (O,-0.
4n2 2

Consequently, for all n “Th“ < 1. This implies
by a well known theorem in linear algebra (see Varga [32]

or Isaacson and Keller [15]) that (Ih"Th)-1 exists and

moreover

(36) (I -T)-1 < 1 .
)) h h “_m

Sh is easily seen to be invertible by observing that
it is row equivalent to Ih' Therefore Mh must be nonsingular

and from (35) we have

-1 _ -l
(37) “11 — (Ih-Th) Sh'
Taking norms in (37) and using the norm inequality

(36), we can write

(38) \uw;1\\_<_ “uh-Thrlu - ushu
ushu

1- Th

 

. 1 . 5
Since “Th“-S 1--4n2 min (0,3), we have

1 .. “Th“ 2 1 _ (1 - if, min (0.31))

 

min (0,59) .
4n

.1.

131

 

 

 

But
n
1 l .
[\S H = max —. --- 21)
b 4n2 4n4 i=1
= max ( l 1 njn+l))
4n2 4n4 2
= -l’- max (1, 95%)
4n
= __1_
4n2
Finally, using the above estimates in (38), we have
1
"‘2'
-1 4n
HM}. H _<_ 1 , 5 = —L3-.
7111111 (0:3) min(0,2)

4n
Since this is independent of n, and therefore of h, we can
conclude that “MRI“ g_C: where C is a finite constant

independent of h.

Thus F is stable and Theorem 1 of Chapter 2 allows

h
us to conclude that (9) has a unique solution X(h) which

converges discretely to w(t).

Clearly, the assumptions (33) are required by our
method of proof. We believe that there is nothing inherent
to the Operator Fh
but at this time we are unable to construct one.

In the next section we will return to the study of
the general operator Fh' without the assumptions (33), and

examine its global discretization error.

- 4.32:“.

 

which would prevent a more general proof,

 

132

Section 2.3. The Global Discretization Error

 

In this section we investigate E(h)!EX(h)—mhw,
the global discretization error for our numerical method.

Assuming the discrete convergence of our method, we know

 

that E(h)-90 as h-+O. We seek more information as to the
nature of E(h) as a function of h. Ideally we would like
E(h) to have an asymptotic expansion in even powers of h.

That is, we would like to show that Vt E [ (L-l)r,Lr] .

 

M
(39) E(h) = %(k21ek(t)h2k) +0<h2M+2)

with the functions ek(t) independent of h. Equivalently,
‘we want to show how to determine ek(t) for k=l,2,...,M and

t.€[(Lrl)r,Lr] so that

M
(40) S(h) a E(h)-wh( Zlek(t)h2k) = 0012””).
k:

If we write

M
u(h) = Zek(t)h2k
k=1

then we can rewrite (40) as
S(h) = E(h) - mh(u(h)).
Let I = Fh(wh(w+u(h))), then since Fh(X(h)) = O we have
\\I\\ = \\Fh(wh(q>+u(h)>) - Fh(x(h))\\.

Using the linearity of wh and the stability of Fh' we

can write

133

“I“ 2215 “X00 - whcp- whum \\
=-}; “E(h) - mhuon \\
= $— nsm “-

Consequently, if we can show how to choose ek(t)

2K+2)

for k=l,...M so that “I“ = 01h then (40) will be

valid and we will have the desired asymptotic expansion.

Theorem 1. Let Fh(V) have Mi-l continuous FrEChet deri-
vatives with respect to V. Then,

1% thwhek(t) +O(h2M+2)

X(h) - cp =
m“ k=l

where the vectors whek(t) are solutions of

‘41) WW) whek = My

The vectors uth will be defined in the proof.

Proof. Expand I = Fh(wh(¢+u(h))) in a Taylor series

about the vector ufim to get

M
(42) I = FhwhmHFfiNth) whu(h)+k§2 El;- Erik) (whee) (whufll) ]k+RM+1°

Here Fék)(uhw) is the kth Fréchet derivative of Rh
evaluated at whm. The remainder term RM+1 involves
[[J,(h)]M+1 which is o(h2M+2) and therefore R is

M+1

2M+2

04h ).

 

 

 

134

Fh(whm) is the vector of local truncation errors

for the discrete Operator Fh‘ By equation (22) each component

of Fh(whm) is given by

M .
_ 2] . 2M+2
Th(ti) — F(o)(ti) + jEfih g2j( )‘¢:t1+0(h

).

Using the space discretization wh' we can write

(43) Fh(thP) = wh'rh
M .
23 2M+2
Map)??? wh(gzj(o)1.p>+ om >.

2M+2

Here we have used the linearity of wh' Also, the 0(h )

term in (43) is understood to be a vector each of whose
2M+2)

components is 0(h

O and substituting (43) into

Noting that F(o)

(42), we obtain

M .
(44) I = j§1h23%92j("‘co + Filmhcpmhum

(R)

5‘15- Fh (mhco)[whu(h) ]k + 003””).

P43

+
k

2

Using the definition of g(h), it is immediate that

2

the coefficient of h in (44) is

whgz + Fi‘whw’ whey
In order to eliminate this term from the expansion (44),

define thl = —whgz and take “hel to be the solution of

Ffimhq’) whel = th1°

 

135

With whel so determined we can write (44) as
M .
I: 73112] wh92j+F' (whtp)wh(j:322jhe
j=2

1.1“); (whcp)[whu(h) ]k + 0mm”).

3...»
ll 3
;:

Collecting the terms involving h4 we have

F(2)
1'14 ”01194 + F' (mhfo)h4whezi- (mhcp)h2 L“r1e11"2""t~1€31'
Define B = g F(2 2)( o) e e and obtain e
“’h 2 “*n 4" mh “’h 1%1wh 2
by solving
F111(0‘)hq))u’)heZ = th2
Substituting whez in the expression for I we can eliminate

the terms involving h4.

Suppose whe1,whe2,...,wheJ have all been determined
as solutions to equations of the form (41), so that the

expression for I is

2j

I = 23 h g + F'( m) (h)
j- —J+l wh 2j( .))m wh ”h”
(k) k
+ Z) . F ( w)[ (h)
k=—J+1k1"h wh whu ]
JM
+ 2 62k( )h2k + 0012“”) .
k=J+l

The arguments of the functions sz(-) involve various error

. . I . .
vectors whej for J=l..-..J and various Frechet derivatives

 

 

 

136

Féj)(wh¢) for j=2,...,J. If

tfiQ

%?Féj) (thP) [ mth) ]j

j 2
is expressed as a power series in h2 then for k=J+l,...,JM
the function G2k(') can be obtained as the coefficient of

th in this expansion.

Collecting terms that include h2J+2 we have

2J+2

2J+2
"1192323r1 +G2J+2 ( )h

+F'(““hcp)""hea'H'

Because the function G2J+2 includes only the error vectors
whel""’wheJ as arguments of the various Frechet derivatives

it contains, we are able to define thJ+l as

thJ+1 = ’“h92J+2"G2J+2'

Therefore, we can determine wheJ+l as the solution to

I (1) =
Fh( hm)wheJ+l thJ+l‘
Hence by induction Theorem 1 is valid. D

Thus the global discretization error will have an
expansion of the form (39) provided Fh(V) has M+l continuous

Fréchet derivatives and we can solve equations of the form

(41) .

In the finite dimensional case, the Fréchet derivative
of a discrete differential Operator Fh(V) is just the

JacObian matrix of the operator considered as a vector-valued

 

137

function of the R—dimensional vector V = (V1,...,VR)T. For
our particular difference operator defined by (9) the jth
row of the JacObian matrix is obtained by computing the

partials with respect to the components of V of equation

(9) with i=j—l.

We will denote the jth row of Ffi(V) by [Ffi(V)]j .

Introducing the notation

fk(v) _ f( v v Vk+l-vk-1 Vk+1—R’Vk-1—R)
‘ tk' k' k-R’ 2h ' 2h

 

 

and using it for the partial derivatives of f(t,u,v,y,z),

we can write Fh(whw) as

[Ffi(whw)]1, = h'2(1-+§f;(whw).o,...,0);

[Ffi(wh@)]2. h-2(-2+h2f:(u%¢0.l‘tgf;(wh¢).0....,O);

and for 3_<_ng.

I _ -2 b. j-l _ 2 j-l
.11 j—l
1+2fy (whcp).o,...,o)
with the non zero entries in the jth row for 3_<_ng

occurring in the (j—2)nd, (j-1)st and jth positions.

Since Fh(whw) is lower triangular, it will be
nonsingular provided none of the diagonal entries vanish.

For any j the (j,j)th entry of Ffi(wh¢) is

138

This will be non zero provided

:w'lkJ

(45> fym ,1

If we assume that fy(t,u,v,y,z) is bounded, then as h-+O

the left hand side of (45) is bounded while the right hand side
diverges to -w. Thus for sufficiently small h, (45) will

be valid and Ffi(whw) will be nonsingular. Consequently,

the equation (41) will be uniquely solvable.
Now let us consider the system of equations
(46) Fh(u%¢)ufiek = O

in more detail. The jth equation in (46) is given by

2

h-2 (1 - 121—£34 (thP) )ek (tj_2) +h‘2 (-2+h fg‘l (when) )ek (tj_1)

-2 :1. j—1 _

+h (1+2fY (thpHekuj) — o.

Equivalently, we can write this as
ek(tjf2)-2ek(tj_l)+ek(tj)

j-l
(47) h2 + fy (wh

 

cp) ek (13)—e? (tj-z)
2h

j-l _
+ fu (mhcp) ek (tj-l) — O.
This equation is the discrete analog of
(48) e]g(t) + fy(cp)e};(t) + fu(CP)ek(t) = 0
at t=tj_1. where fy(w) = fy(t,m(t).¢(t—r),¢'(t),w (t—r)).

Equation (48) is the linear variational equation associated

 

 

 

139

with the continuous Operator F(x) = O and is often denoted

by F'(Cp)ek = 0.

By examining the first and second equations in (46)
we see that ek(t_l) = ek(to) = O. In fact, since
E(h) E X(h)-whm and we have taken X_j = w(t_j) for all '1
1

j=0,1,...,R+l we must have ek(t_j)==o Vk and Vj==O,...,R+l. “...

One method for obtaining (48) is to let h-OO in V

 

(47). If we do this, the points t_j become dense in the b1
interval [(L-2)r,(L—l)r] and we therefore must have

0

ek(t) 50 on [(L-2)r,(L-l)r]. This implies that ei(t)
on [(L—2)r,(L-l)r]. Since the solution to (48) will have a
continuous first derivative, we see that the discrete problem
(47) with the initial conditions ek(t_1) = ek(to) = O corres-

ponds to the continuous prOblem (48) with the initial con-

ditions ek(t0) = ei(to) = 0.

Consequently, the vectors “hek in (41) are actually
discretizations of differentiable functions ek(t) which

satisfy

(48') e§(t)+fy(m)ei(t)+fu(w)ek(t) = Bk(')

ek(to) = e];(t0) = 0.

To summarize, we have established that the global

discretization error has an asymptotic expansion in even powers

of h:

140

M
X(h) (ti) = cp<ti> +k§16k<ti)h2k + 003””).

Moreover, the functions ek(t) are independent of h and

satisfy ek(to) = ek(to) = O.

In the next section we will discuss the implementation

of this method and the modifications necessary to obtain a .1

solution on [O,r].

 

141

Section 2.4. Implementing the Second Order Method

In order to obtain a numerical solution to

(8) £(t)+f(t.x(t).x(t-r).x(t).k(t-r)) = o, r:>o

on [O,r], some modifications of the discrete method are

 

needed. The first modification is based on the fact that E]__F
on [O,r] we have more information available to us concerning 1

the behavior of the theoretical solution than we do later.

 

To obtain the theoretical solution to (8) on [O,r] {3
we start with a given initial function m(t) whose derivative
m'(t) is known on [—r,O]. The extra information we have
in this case is knowledge of ¢'(t). This is incorporated

into the discrete version of (8) as added initial conditions.

Let R:>O be any given natural number and define
h==r/R. Construct the uniform grid {ti=0+ih)§=—R' Since
we are given ¢(t) and ¢'(t) for t.€[-r,0] 'we may
define a discrete version of (8) for t.€[O,r] by

‘ZX'+Xi—1 Xi+i'Xi-1

X.
1+ 1 ' _ ._ - .
2 + f(ti'xi'wi-‘R' 2h 'mi-R) _ 00 1-0. 1' o o O'R 1,

1

 

 

GPi-R = cp(ti-R) E Cp(ti'r)'

' ' t - -= - .
w (ti-R) a o ( i r), 1 O,l,...,R 1.

6
I
71
ll

X
l

_1 - w(t_1) a w(4h).

X
I

O — Cp(to) E cp(O).

142

The other modification that we must make is based on
the fact that the theoretical solution to (8) will, in general,
not have a continuous second derivative at the origin.

X'+1‘2X1+Xi-1

Since the expression 2 for i=0 in
h

(49) is an approximation to the second derivative of the sol-
ution at t=0 ‘we must take steps to insure that the second
derivative exists and is continuous at t=0. However, we
want to accomplish this in such a manner that the initial

information X X_1, m._R and w i-R for i=O,l,...,R—l

O' i

is not altered.

One way of accomplishing this is to modify the given
initial function m(t). Select. e<:O such that -h<:e<(0
and define a new initial function $(t) by using an inter-
polating polynomial in place of m(t) on [c,O]. More
precisely, let H(t) be the fourth degree Hermite interpolating
polynomial which satisfies P(e) = m(e). P'(e) = m'(e).

P(O) = cp(0). P'(O) = CP'(O) and P"(O) = f(O.CP(O).cp(-r).

¢'(O).w'(-r)). Define a new initial function E(t) on
[-r.0] by

CP(t) —rg_tge
(50) E(t) =

HR) 6 gtgo.

By construction E(t) will be continuously differentiable
on [—r,O] and since $"(o) = f(o,o(o),o(—r),w'(o),o'(-r)). the

second derivative of the solution obtained with the initial

 

 

il

143

function Q will be continuous at t=0. Now $(t) agrees
with w(t) for -r5;t5;e and also for t=O. This includes
all points in [-r,O] at which initial information is needed
in (49). Since, as mentioned earlier, the solution to (8)
depends continuously on the initial function, the use of

E(t) in place of m(t) will not radically alter the character
of the solution to (8). It will however make (49) a reason-

able discretization for (8).

Note that, as we move to the right solving (8) with
the method of steps, we gain differentiability of the solution
(see the discussion for first order equations in Section 1
of this chapter). Thus the solution to (8) will have a
continuous second derivative for all t:>O and the approx-

. . Xi+1‘2Xi+Xi—1 . . . .
imation 2 w111 not cause any difficulties at

h

other multiples of r.

The proof that the discrete prOblem (49) is consistent
is a special case of Section 2.1. The only difference is that
now we are using the actual value of ¢'(t-r) instead of an
approximation. This will simplify the expression for the
local truncation error but does not alter the fact that the
truncation error has an asymptotic expansion in even powers
of h. Also, note that the proof of stability with the added
hypotheses (33), and the analysis of the gldbal discretization

error given for (9) carry over directly to (49).

 

 

:2:
.

144

Using the appropriate discrete scheme, (9) or (49), we
can compute a solution to (8) whose global discretization
error has an asymptotic expansion in even powers of h.

Since the coefficient error functions in this expansion are

independent of h. extrapolation can be performed to obtain

 

a more accurate solution.

 

Let R be any natural number and let h0 = r/R be
the basic steplength. For each k=l,...,M define
= h /2k = -£—- Define grids b
bk 0 k ' Gk Y
2 R
Gk -—- {.t§:t]i‘=o+ihk,i=-2kR,-2kR+1, . . .,o, 1, . . .,2kR] for k=0, 1, . . .M.

As before, the grids Gk are nested; Gk<:Gk+l'

Given an initial function m(t) such that w(t)
and m'(t) are continuous on [—r,O], select 3:10 such
that -hM<’e<ZO. With this 6 define a new initial function,
again denoted by w, as in (50).

On each grid Gk compute the solution X(hk) to
(49). Denoting the solution to (8) on [O,r] also by w(t)

we know from our previous work that

2M+2 k

“23' k
= h ej(ti)+0(hk ), for i=1,...,2 R.

_ k
(51) [X(hk) ]i — wti) + '23
j 1
If R=l, then we can use pullback interpolation as

described in Chapter 1 to obtain a solution, X(h), which

 

satisfies

145

2M+2

(52) [X(h)]i = cp(ti) + who

)

for each. t.i EGM. Actually the pullback interpolation method
does not involve as much work in this case as it does for first
order initial value problems. This is so because from (48')

we have ej(O) = O for j=l,...,M. Therefore we do not have

to go through all the work of constructing a table analogous

to Table 4 of Chapter I to obtain e3(O).

If R;;2, then we can obtain uniform accuracy at

all grid points of G without using any information about

M
e5(O). The reason for this is that with R::2 we have

enough data points at each stage to do Lagrange interpolation
with accuracy comparable to that of the data points. That is,
when performing Lagrange interpolation to accomplish the
pullback with R¢;2, it is the accuracy to which we know the
error functions that dominates (see the discussion in Section 5
of Chapter 1). Thus we do not need to do Hermite interpolation

in this case and therefore the information about ej(O) is

superfluous.

Once the solution X(h), satisfying (52),is obtained
we can use the discrete method (9) to Obtain a numerical
solution on [r,2r]. This process, theoretically at least,
can be repeated indefinitely on intervals of length r, to

obtain a solution which is 0(h2M+2)

O

 

 

146

It should be noted however, that when using (9) to
solve (8) on intervals of the form [(L-l)r,Lr], L;gl, the
initial data is taken from the computed solution on the
previous interval. Since the initial data will not be precise

2M+2)) an accumulation

(the error in the initial data is 0(h
or errors is unavoidable. Of course, this also occurs when

any numerical method is used on successive intervals.

A variety of other schemes, based on (9), (48) and
pullback interpolation, to solve (8) are possible. We could,
for instance, use (49) and pullback interpolation on [O,ho]
to obtain a solution. Then by modifying the initial con—
citions in (49) we could repeat this procedure on [ho,2ho].
Repeating this R times we would have the solution on all

of [O,r]. We could then follow the same format, using (9)

in place of (49), to obtain a solution on [r,2r], etc.

We do not make a comparison of these schemes here.
In the next section we will merely exhibit some numerical

examples computed using the first described solution tech-

nique.

 

 

147

 

 

Section 2.5. Numerical Results for Second Order Equations
This section consists of two examples. The first of

these is

Example 3: Solve x(t) = X(t-l), Oj;t3Ll numerically for

the initial function m(t) = 12(t+l)2-+24 on [-l,O]. The
method of solution is the discretization (49), extrapolation
and pullback interpolation with four grids. The steplengths

employed are hk = 1/2k for k=O,...,3 and the grids are

_ k _ . . ._ k
Gk _ {ti _ o-+ihk. l—O,l,...,2 }

for k=O,...,3.

Note that the initial function ¢(t) satisfies
m(-l) = m"(O-), so the second derivative of the solution is
continuous at the origin. Thus, the initial function does

not need to be modified as in (50) for this example.

The theoretical solution is given by

m(t) = t4.+12t2-+24t-+36

for t €[O,l]. The numerical results presented in Table 15
are for the nine equally spaced points in [0,1] which
belong to the finest grid G3. The error reported is once

again the computed solution minus the theoretical solution.

 

 

 

148

TABLE 15
i Error
o 0.0
1 0.0 x 10‘12
2 -.23 x 10"12
3 -.23 x 10'12
4 -.23 x 10"12
5 - 23 x 10'"12
6 -.23 x 10'12
7 0.0 x 10'12
8 0.0 x 10"12

The high accuracy of the numerical solution in this
example is due to the fact that the theoretical solution is
a low degree polynomial. Thus, extrapolation and consequently
pullback interpolation are quite accurate in this case.
This accuracy is not to be expected for more general equations

as the following example illustrates.

Example 4: Solve x(t) + %x(t) - %x(t—w) = O on [O,w]

for the initial function cp(t) = l-sin t, —7Tgt_<_O. This
example is from Norkin [23] and the theoretical solution is
also given by cp(t) = l-sin t, ogtgw. Thus, the second
derivative of the solution is continuous at the origin and

no modification of the initial function is required.

 

 

149

The method of solution is to use (49), extrapolation
and pullback interpolation with four grids to Obtain a num-
erical solution on [O,W/4]. Using this numerical solution
as initial data, the discretization (9), extrapolation and

pullback interpolation with four grids is employed to obtain

a solution on [%3%J. This process is repeated twice more
to Obtain numerical solutions on [%u%gj and [%¥,w]. On

each interval of length w/4 the stepsizes used are

hk = Trk for k=O,...,3. Thus, on each Of the four sub—
4-2
intervals the finest grid contains nine equally spaced points

 

which we number as t .,t

O"° 8'
The accuracies Of the numerical results are reported

as Tablele‘with the error being computed in the usual manner.

 

 

 

150

 

 

 

TABLE 16
Error on [O,TT/4] Error on [IT/4,7T/2]
0.0 .5 x 10'8
15.9 x 10‘8 - 1.0 x 10"6
37.6 x 10‘8 — 4 5 x 10'6 ‘a
30.2 X lo-8 - 6.4 X 10"6 l~—E"
— 1.6 x 10-8 - 6.6 x 10"6 ,
-33.2 x 10"8 — 6.7 x 10"6 ~—%
~40 O x 10"8 - 8.4 x 10'6 uj
-l7.l x 10'8 -ll.5 x 10‘6
.5 x 10"8 -11.9 x 10—6
Error on [Tr/2,3Tr/4] Error on [3TT/4,TT]
-11.9 x 10'6 -27.2 x 10‘6
-l3.7 x 10'6 -28.6 x 10'6
-18.7 x 10‘6 —31.9 x 10'6
-21.2 x 10"6 -33.2 x 10‘6
-21.0 x 10'6 -32.4 x 10‘6
-2o.6 x 10"6 -31.4 x 10‘6
—22.6 x 10"6 -32.0 x 10'6
-26.6 x 10‘6 -34.4 x 10'6
-27.2 x 10"6 -34.1 x 10’6

151

The accuracy Of the computed solution is considerably
higher on the first subinterval than on the last three sub-
intervals. Also the errors are steadily increasing at a
very low rate over the last three subintervals. This suggests

an accumulation of truncation errors.

This accumulation is to be expected for the following
reason. For a linear equation, the discretizations (9) or
(49) are similar to a well known linear multistep method
for solving second order initial value problems (see Lambert
[18]). The only difference is the manner in which the
discretizations are initialized. The stability theory for
linear multistep methods for solving second order equations
predicts an accumulation Of truncation errors. Thus, based
on comparisons with similar methods, our numerical results

are consistent with how our method should behave.

Note that the equation x(t) + %x(t) - %X(t-W) = O
is not covered by our proof Of stability. Yet the numerical
results seem to indicate that the method is stable when
applied to this equation. Consequently, we again reiterate
our belief that stability is valid for a much larger class

Of functions than our proof indicates.

 

 

 

BIBLIOGRAPHY

 

10.

ll.

BIBLIOGRAPHY

Aitken, A.C., On interpolation by iteration Of
prOportional parts, Proc. Edinburgh Math. Soc.,
2, 56-76 (1932).

Bj5rck, A., and V. Pereyra, Solutions of Vandermonde
systems of equations, Math. Comp., 24, 893-903
(1970).

Cooke, K.L., and S.E. List, The numerical solution of
integrO-differential equations with retardation,
Department Of Electrical Engineering, University Of
Southern California, LOs Angeles, Technical Report
NO. 72-4 (1972).

Cooke, K.L. and J.A. Yorke, Equations modelling pOpulation
growth, economic growth and gonorrhea epidemiology,
preprint.

Corliss, G.F., Parallel Rootfinding Algorithms, Ph.D.
thesis, Michigan State University, East Lansing,
Michigan (1974).

Cullen, C.G., Matrices and Linear Transformations,
Addison-Wesley Publishing Company, Reading, Mass.
(1967).

Cunningham, W.J., A non-linear differential-difference
equation Of growth, Proc. Natl. Acad. Sci., U.S.,
40, 709-713 (1954).

Dahlquist, G., Convergence and stability in the numerical
integration of ordinary differential equations,
Kungl. Tekniska Hogskolans Handlingar NO. 130 (1959).

Dahlquist, G., A special stability prOblem for linear
multistep methods, BIT, 3, 27-43 (1963).

Feldstein, A., Discretization Methods for Retarded Ordinary

Differential Equations, Ph.D. thesis, University Of
California, Los Angeles (1964).

 

Gragg. W.B., Repeated Extrapolation to the Limit in the
Numerical Solution of Ordinary Differential Equations,
Ph.D. thesis, University Of California, Los Angeles
(1964).

152

 

 

 

12.

13.

14.

15.

l6.

17.

18.

19.

20.

153

Gragg, W.B., On extrapolation algorithms for ordinary
initial value problems, SIAM J. Numer. Anal., 2,
384-403 (1965).

Hale, J., Functional Differential Equations, Springer-
Verlag, New York, N.Y. (1971).

Huygens, C., De Circuli Magnitudine Inventa, Leiden
(1654).

 

Isaacson, E., and H.B. Keller, Analysis of Numerical
Methods, John Wiley and Sons, Inc., New York, N.Y
(1966).

 

Joyce, D.C., Survey Of extrapolation processes in
numerical analysis, SIAM Review, 13, 435-490
(1971).

Keller, H.B., Numerical Methods for Two-Point Boundarye

 

 

 

Value Problems, Blaisdell Publishing Company,
waltham, Massachusetts (1968).

Lambert, J.D., Computational Methods in Ordinary Differ-

 

ential Equations, John Wiley and Sons, London (1973).

Lindberg, G., A simple interpolation algorithm for
improvement of the numerical solution of a differ-
ential equation, SIAM J. Numer. Anal., 9, 662-668
(1972).

Milne, W.B., and R.R. Reynolds, Stability of a numerical
solution Of differential equations, J. Assoc. Comput.
Mach., 6, 193-203 (1959); Part II, J. Assoc. Comput.
Mach., 7, 46-56 (1960).

Milne, W.B., and R.R. Reynolds, Fifth-order methods for
the numerical solution Of ordinary differential
equations, J.Assoc. Comput. Mach., 9, 64-70 (1962).

Neville, E.H., Iterative interpolation, J. Ind. Math.
Soc., 20, 87—120 (1934).

Norkin, S.B., Differential Equations of the Second Order
with Retarded Argument, American Mathematical Society.
Providence, Rhode Island (1972).

Pareyra, V., On improving an approximate solution Of a
functional equation by deferred corrections, Numer.
Math.. 8. 376—391 (1966).

 

25.
26.

27.

28.

29.

30.

31.

32.

33.

34.

154

Pereyra, V., Iterated deferred corrections for nonlinear
Operator equations, Numer. Math., 10, 316-323 (1967).

Pereyra, V., Accelerating the convergence Of discretization
algorithms, SIAM J. Numer. Anal., 4, 508-533 (1967).

Pereyra, V., Iterated deferred corrections for nonlinear
boundary value problems, Numer. Math., 11, 111-125
(1968).

Pereyra, V., High order finite difference solution Of
differential equations, Seminar Notes, Stanford
University, Palo Alto, California (1973).

Richardson, L.F., The deferred approach to the limit,
I-single lattice, Trans. Roy. Soc. London, 226,
299-349 (1927).

Sansone, G., Teorema di esistenza di soluzioni per un
sistema di equazioni funzionali differenziali, Ann.
Mat. Pura Appl., 39, 65-67 (1955).

Stetter, H.J., Asymptotic expansions for the error Of
discretization algorithms for non-linear functional
equations, Numer. Math., 7, 18-31 (1965).

Varga, R.S., Matrix Iterative Analysis, Prentice-Hall,
Inc., Englewood Cliffs, New Jersey (1962).

Widder, D.V., Advanced Calculus, Prentice—Hall, Inc.,
Englewood Cliffs, New Jersey (1961).

Wright, E.M., A non-linear difference—differential equation,
J. reine angew. Math., 194, 66-87 (1955).

 

 

SECURITY CLASSIFICATION OF THIS PAGE (When Dete Entered)

 

REPORT DOCUMENTATION PAGE BEFTDEEDC‘SSEEEEESNSORM

 

‘0

REPORT NUMBER 2. GOVT ACCESSION NO. 3. RECIPIENT'S CATALOG NUMBER

 

 

C.

SOLUTIONS TO DIFFERENTIAL EQUATIONS USING
EXTRAPOLATION AND INTERPOLATION.

TITLE (Md Subtltle) [INI FORMLY ACCURATE NUDERICAL 5. TYPE OF REPORT 6 PERIOD COVERED

 

6. PERFORMING ORG. REPORT NUMBER

 

 

7. AUTHOR(3)R]° Chard Allan Rogers 8. CONTRACT OR GRANT NUMBER“)
AFOSR-72-227l
9. PERFORMING ORGANIZATION NAME AND ADDRESS IO. PROGRAM ELEMENT. PROJECT, TASK

, , . . AREA 6 WORK UNIT NUMBERS
Michigan State UniverSity

Department Of Mathematics
East Lansing, Michigan 48824

 

 

 

 

 

 

 

 

 

 

H. CONTROLLING OFFICE NAME AND ADDRESS I2. REPORT DATE
Air Force Office of Scientific Research August 1974
1400 Wilson Blvd. I3. NUMBER or PAGES
Arlington, Virginia 22209 154
A. MONITORING AGENCY NAME 6 ADDRESSU! dIllerent from Controlling OHI'Ce) 15. SECURITY CLASS. (0! this report)
UNCLASSIFIED
ISe. DECLASSIPICATION/DOWNORADING
SCHEDULE
‘6. DISTRIBUTION STATEMENT (0! this Report)
Approved for public release: distribution unlimited.
I7. DISTRIBUTION STATEMENT (o! the abstract entered In Block 20, H dlflerent Irom Report‘-
13. SUPPLEMENTARY NOTES

 

KEY WORDS (Continue on reverse side i! necessary and Identity by block number)

 

 

20.

ABSTRACT (Continue on reverse side I! necessary and Identtly by block number)

In this work we are concerned with numerical methods for
solving ordinary differential equations. We consider those
methods that have asymptotic error expansions involving all
powers Of hq, where h is the steplength and q is a fixed

integer. The process Of extrapolation can be employed with such

methods to Obtain highly accurate solutions at grid points
belonging to the coarsest mesh. In Chapter I we develOp the

 

 

DD ,FORM 1473 EDITION OF 1 NOV 65 Is OBSOLETE

JAN 73

 

SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered)

1..

 

SECURITY CLASSIFICATION OF THIS PAGE(When Dare Entered)

 

 

"pullback interpolation method". This method combines extra-
polation with Hermite interpolation Of the coefficient functions
for the asymptotic error expansion to produce a highly accurate
solution at all grid points Of the finest mesh. When q is

l or 2 pullback interpolation yields uniform accuracy at all
grid points Of the finest mesh.

In Chapter II the pullback interpolation method is mod-
ified so as to be applicable to boundary value problems. In
addition, an elementary proof Of the stability Of V. Pereyra's
finite difference scheme for solving two point boundary value
problems is given.

In Chapter III we consider difference differential equations
with constant retardation. The methods Of Chapter I are shown
tO be applicable to first order delay equations. Because Of
the presence Of the delay term, the uniform accuracy Obtained
through pullback interpolation is indispensible for these
problems.

A finite difference scheme for directly solving second
order delay equations is constructed and analyzed in Chapter III.
The global discretization error is shown to have an asymptotic
error expansion in even powers Of the steplength.

 

 

 

 

SECURITY CLASSIFICATION OF THIS PAGE(When Data Entered)

 

 

”IIIIIIIIIIIIIIIITE IIIIIIIIIIIIIUIBIIAIIR'ES
31293 03168 8371