Boundary value problems and dichotomic stability

Boundary value problems and dichotomic stability

Citation for published version (APA):

England, R., & Mattheij, R. M. M. (1988). Boundary value problems and dichotomic stability. SIAM Journal on Numerical Analysis, 25(5), 1037-1054. https://doi.org/10.1137/0725060

DOI:

10.1137/0725060

Document status and date: Published: 01/01/1988 Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

BOUNDARY VALUE

PROBLEMS

AND DICHOTOMIC

STABILITY*

R. ENGLANDt AND R. M. M. MATTHEIJ$

Abstract. Sincethe conditioning ofaboundaryvalue problem(BVP)isclosely relatedtothe existence ofadichotomic fundamental solution (i.e.,whereonesetofmodes is increasing andacomplementaryset

is decreasing),it isimportant tohave discretization methods that conserve this dichotomyproperty. The conditions thisimposesonsuchamethodareinvestigated in this paper.

Theyareworkedoutinmoredetail for scalar second-orderequations(thecentraldifferencescheme), andfor linear first-order systemsaswell;for the latter type both one-step methods(including collocation) and multistep methods (thosethat may be used in multipleshooting)areexamined.

Keywords, boundaryvalueproblems,dichotomy, stability AMS (MOS)subject classifications. 65L05, 65L10, 65L15

1. Introduction. Inthe studyofboundaryvalueproblems

(BVP)

stability notions, describingandinterpreting theeffectsof

(small)

localperturbations, playanimportant role.Forinitialvalueproblems

(IVP)

thestability theoryof numerical methods is very

welldeveloped

(cf. [3], [5], [14],

[24]);

for BVPthisstability question is much more

complicatedand less developed, althoughthere hasbeenrapid progress,inparticular, forsingularlyperturbed problems

(cf. [1], [13], [21], [25],

[28]).

The basic difficulty for

BVPs

is that the global errors depend onthe data of the entireinterval on which theproblemis defined,whereas they only dependon data of the past interval in IVPs. Nevertheless thereis much similarity.

Indeed,

thinking ofa

linearproblem, where the solution space of the homogeneous equations can besplit

into subspaces of decaying modes onthe one hand andgrowing modes on the

other,

it is knownthataconditionnumbermainly dependsontheboundaryconditions

(BCs)

imposed;these should be such that the "initial conditions" controlthe decaying

modes,

andthe "terminal conditions" thegrowing modes

18].

Thisquestioniscloselyrelated to the notions of conditional stability

(cf.

[23])

and dichotomy

(cf.

[4]) (see

also

[2],

[15], [17],

[20]).

The latter concept will also play animportantrole here. In principle

dichotomy denotes a splitting of solution spaces into subspaces of solutions with a

markedly differentgrowth behaviour, likeincreasing o decreasing, increasing faster thanacertainexponentialrateoincreasingslowerascomparedwith this rate.

Recent

results show that inawell-conditioned

BVP,

theordinarydifferentialequations

(ODEs)

should be dichotomic in the sense that there is a splitting into solutions that do not increase significantly on the one hand and do not decrease significantlyonthe other

(see 11]).

Sincethe BCs also control the modes ofthediscretized problem,it isclear that it makes sense to aim at discretizations that produce a decaying (growing) approximate mode if the corresponding continuous mode is decaying (growing), in particular, for singularperturbation problems.

Among

existing BVP algorithms, it seems that multiple shooting

("stabilized

march")

types of methods have hardlybeen investigated from this point ofview, in contrastto global approaches such as collocation

(cf. [1],

[25]).

In [1],

for example,

Receivedbythe editors March 1, 1984;acceptedfor publication (inrevisedform) August11, 1987.

This workwas supportedin part byagrant from the Netherlands Organization for the Advancement of PureResearch(ZWO).

Instituto deInvestigacionesenMatemfiticasAplicadasyenSistemas,Universidad National Aut6noma deM6xico, ApdoPostal 20-726, 01000 Mexico D. F. Present address, Facultyof Mathematics, TheOpen University, WaltonHall,MiltonKeynesMK76AA,United Kingdom.

$Faculteit Wiskunde en Informatica, Technische Universiteit Eindhoven, Postbus 513, 5600 MB Eindhoven, the Netherlands.

1037

afairly detailedinvestigation ismade of thedampingof fastcomponentsin theproper

direction.

As

we shall show in this paper a similarstudy can also bemade for other

methods,

including multiple shooting.

A

distinction must be made between "global" and "local" discretizations.

A

recurrence relation for approximateoutput

values,

suchasisfound inmultiple shooting orcondensedcollocation atthematchingpoints, is called theglobal

_difference

equation. The actual discretizationmethod that is definedonthe finergrid, e.g., bya

Runge-Kutta

method in multiple shooting, is called the discretization method. Unless the finer grid coincides with the coarser, the behaviour of particular interest is the growth of the modes of the global difference equation;

however,

this partly reflects that for the discretization scheme from which it results by internal condensation. The approach

will be based upon investigating increments of the global difference equation for suitable modelproblems

(as

is donein

IVP

stability

theory).

2. Dichotomic stability.

As

remarked in 1 a satisfactory numerical method for

solving

BVPs

should approximate bothdecreasingandincreasingmodesproperly.The word "accurately" is deliberately avoided, because we are rarely interested in fast modes outside narrow

(boundary)

layers. In such layers accuracy may be desirable,

but outsidethem therequiredsolution willoftenbequitesmoothsothat we may wish to uselarger stepsizes. This ispreciselythe

BVP

analogue of what is called "stiffness" in

IVPs (cf. [5],

[24]). Note,

however,

that the related stability question for BVPs is different. Indeed for stiff

IVPs

it is desirable that subdominant modes and spurious

modes damp out as increases, something that does not make sense for BVPs. It is then of interest that these modes which make little or no contribution to the true solution damp out inthe appropriate direction.

As

follows from what has been said in 1, the nature of the BCs appears to be most important

here,

as they somehow have to control dominant, subdominant, and any possible spurious solutions. Before

amoreprecise definition is given of what should be called

stable,

an example is first described.

Example

2.1.

Let

y"=

A2y+f(t)

(A

a

constant)

and

y(0)=y(1)=

1. Divide

[0,

1]

into N subintervals of length h.

Denote

_f

=f(ih),

Yi the approximant to y(ih) and

a

h2A 2.

Consider the following discretization:

(2.2)

-Yi+3+4yi+z-5yi+l+(2-c)yi=hZf,

O<=i<-N-3.

Thisclearly giveslocal discretization errorsof

O(ha).

Since it isathird-order difference equation, but has onlytwo

BCs,

it is necessary to provide an extra BC in order to solve fory uniquely.

Suppose

the relation

(2.3)

Yi+-(2

+

a)y

+

Y-I

h2f

isused to express either y interms of_Yoand _{Y2 or YN-1} in terms ofYN-2and YN.

(It

also has local discretization errors of

O(h4).)

This "trick" virtually creates an extra

BC at 0 or 1, respectively. The first choice leads to the system

(2.4)

4-5/(2+a)

-1

-5+(2- a)/(2+ a)

4 2-c -5 -1 4 -1 4 -1 -5 4

5fl

_)-(2-a)+

5

h2

_f-2+

a/ 2+a

h2fl

₁₊₂₊

hA

hf_

h2fN_3

+

1

By

choosing -A2as the forcing term

f(t)

the exact solution to the original BVP

is

y(t)=

1 for

t[0, 1];

of course both

(2.2)

and

(2.3)should

integrate this solution

exactlysothat the onlyerrors tobeexpected arethose due torounding. In Table 2.1 the errors

lYi-y(ih)l

are given some values of where A2

and h were chosen as 800 and

1/80,

respectively. The computationswereperformedon an

IBM/MVS

4331 with double precision (relative machine error

10-6).

The second choice, that is, eliminating yu_,leads to the system

(2.5)

-5 4 2-a -5 -1 4 -1 2--a -5 4 -1 -5

4-1/(2+a)

2-a

-5+4/(2+

a).

A

table corresponding to Table 2.1 only shows error <2x 10

-15.

Therefore it would appear that

(2.5)

is stable and

(2.4)

is unstable

A

simple explanation can be given as

follows the characteristicpolynomial X ofthe homogeneous part of

(2.2)

is

(26)

x(r)

r --4r2

+

5r- 2

+

a. TABLE 2.1 Yi-y(ti) 0 0 0 10 0.125 2.00 10-14 20 0.25 1.11 10-12 30 0.375 6.41 10-ll 40 0.5 3.30 10.9 50 0.625 1.88 10-7 60 0.75 1.10 x 10--5 70 0.875 6.40 10-1 80 0

For h

o

and A2= 800

(2-a

),

X has the zeros

__1 ₅

(2.7)

that is, two "unstable" roots

_{(lr l>}

a)

and one "stable" root

(Irjl<

1).

Hence the

homogeneouspart of

(2.2)

has twoincreasingbasis solutions andone decreasingbasis solution. From an analysis such as that givenin

[18]

it follows that it is necessary to have two terminal conditions and one initial condition. Therefore

(2.4)

cannot be

stable, whereas

(2.5)

is likely to be stable.

By

comparing the errors in Table 2.1 at

points with distance

10h,

it becomes clear that the solution

{rl}

which is second in

growth,is responsible for the errorgrowth

(a

factor 58 after each 10

steps). It

turns

out that such a dichotomy

(cf.

[4])

also holds for c

(0, 12];

hence the stability of

(2.5)

is uniform in h ifh

=<

,

althoughfor values ofc

->

7

the situationbecomes more

complicated.

To

examine it, Fig.2.1 shows agraphof

x(r)

forc--0.

As

canbe seen there is adouble root and athird root 2 and there isalocal minimum

_of-7

at For any nonzero value of ce, it is easy to check the roots by shifting this graph up a

distancece.

Hence

aslongasc

<

₂

there arethree real

(and

positive) roots,onesmaller

than, and the other two larger than one. For

4/27

<

ce

<

12 there is a real root of modulus less than and two complex conjugateroots of modulusgreaterthan 1.

For

c

>

12,

however,

the smallest root becomes smaller than-1.

Therefore,

ifh andA are

suchthat

h2A2>

12, both methods are expectedto beunstable.

,/

FIG.2.1

For

thisproblem, thehomogeneousdifferential equation has two

modes,

agrowing

mode e

’

controlled by a terminal condition

y(1)=

1, and a decaying mode e-At

controlled by an initial condition

y(0)=

1. Thediscretization

(2.2)

has three

modes,

one spurious modebeing added by the use ofa difference equation of higher order than the differential equation. While

heA2<

12 it has two growing modes and one

decaying mode. The initial condition controls the decaying

mode;

the final condition controls the dominant (spurious) growing mode; and to obtain a stable numerical

solution, it is necessary to impose the extraboundary condition at the end so as to

control the secondgrowingmode. When

hA

>

12 theoriginal decayingmode becomes

an oscillatory growing mode and the given initial condition canno longercontrol it.

It

is now appropriateto give the following definition.

DEFINITION2.8. Let@ denoteaclass of linearBVPsonaninterval

[a,/3 ],

where:

(a)

the nth order

(n->2)

ODE, Dy=f(t),

is such that the homogeneous equation,

Dy

=0, has k independent solutions

gj(t)

(j 1,...,

k)

which grow in magnitude,

and

(n-k)

independent solutions

gj(t)

(j= k+ 1,...,

n)

with decaying magnitude;

and

(b)

the BC can be written as

Py(a)=p

((n-k)

equations) and

Qy()=q (k

equations).

Let Mh

be a consistent discretization method defined on anequally spaced grid

withmesh size

h,

givinganmthorder

(rn

_->

_n)

lineardifferenceequation

Dh(y

h)

=fh(t)

and BC

Phyh(a)=ph

((m-l)

equations

say),

Qhyh()=qh

(1

equations), where

yh

is defined onthe grid. Let the basis solutions of the discrete problem

Dh(y)

=0, be denoted by

r(ih),

j 1,

,

m,and ordered in such away that

r(ih)

gj(

ih

+

O(

h

),

j 1,...,n. Then Mh

is said to be dichotomieally stable for @ iffor each BVP

of this subset and its discretization and each mesh size h there holds"

(i)

Ir(ih)l

is an increasing function of (j 1,...,

k);

(ii)

]r(ih)]

is a decreasingfunction of (j= k+ 1,...,

n);

(iii)

Of the remaining

(m-n)

solutions of the discrete problem,

)(ih),

j=

n+

1,..., m,

(l-k)

are such that

]r(ih)

increases and

(m-l-n+k)

are such that

]r(ih)l

decreases with i. [3

Remark 2.9. Definition 2.8 requires that to the spurious modes

(only

if m

> n)

there correspondthepropertype of

BCs (e.g.,

initial conditions fordecaying

ones).

It

does not require,however,that these BCs actuallycontrolthe discrete solutionsatall.

Includingthisrequirementwould make it necessaryto introduceacondition number orsomethresholdforit,notonlyfor the discretization but alsofor theoriginal problem.

Since the actual determination of the discrete

BCs

might be uncoupled from the determinationof the difference equation,ourdefinition is still ameaningful one, asit

gives a necessary condition for a stable or well-conditioned discrete problem.

Thus,

dichotomic stability is by no means a sufficient condition fora proper discretization.

Quite apart from accuracy criteria, the conditioning of a problem

(continuous

or

discrete)

dependsonthe boundary conditions,which are not considered in Definition 2.8.

Inwhat

follows,

attention willbe mainly restricted to difference equations which are of the same order as the differential equations. The following simpler notion,

appliedto an equation with constant coefficients, maythenbe used.

DEFINITION2.10. Letthe basis solutions of thehomogeneouspart ofannthorder linear differential equation have components proportional to

e.;

(j 1,...,

n).

Let Mh

be a discretization method giving an nth order linear difference equation with

correspondingbasissolutions of the discreteproblem having componentswithgrowth

{(r)i},

wherer-h

+

_hA

+ O(h2).

ThenMh

is

di(chotomically)

stableon aregionRcC

if and only if, for

hAj.

R,

(i)

Re

(hAj)_-<0:=Jrl=

<_1,

(ii)

Re(hA)_->0=Jr[_

>-1. [3

Obviously, if m n in Definition 2.8 condition

(iii)

is

redundant,

while, for the

constantcoefficient

_ase,

Definition

2.8(i)

and (ii) are equivalentto Definition2.10(i)

and

(ii).

The class @ then consists of BVPs for which the homogeneous part of the

ODEs has basic solutions that foragiven h satisfy Definition

2.10(i)

and (ii). Inthis sense Definition 2.10 is consistent with Definition 2.8. Of course, Definition 2.10 resembles the more familiar IVP notions of stability such as absolute stability and

A(a)

stability

[5],

[29].

Example 2.11. If M is the forward Euler method, then R

c

_IRe

(z)>=O}U{z Cllz+ 11<_-

1}.

3. Aninvestigation of the di-stability ofsomemethods. Inthis sectionthe di-stability ofa numberof methods is investigated to demonstrate the use of the concept, rather

than to give a detailed discussion of all its aspects. In 3.1 consideration is givento the classical central difference scheme for scalar second-order ODEs. Then in 3.2 one-step difference equations for systems of first-order differential equations are examined.Finally,in 3.3amorespecialone-stepequationisconsidered,viz. multiple

shooting relations, and attention is paid to questions to selection of appropriate

discretization schemes forintegration overa shooting interval.

3.1. The central difference scheme. Consider the scalar second-orderODE

(3.1)

y"

+

py’

+

qy O.

Perhaps one of the oldest numerical discretizationmethods uses central differences

(3o2)

1

+

h

y/h+l

+

(--2

+

qh2)yh

+

1 p 12i_

=o.

As

may be seen the basis solutions of

(3.1)

are e

’’,

e2t with

(3.3)

A1, /2 p-1- _{p -4q,}

2

and those of

(3.2)

are

{(rl)h

i}, {(r2)i}

with

-q

_h2

+/--

h,/p2-4q

+ q2h2

1

+-

h

(3.4)

r

h,r2

h=

2 2

Note that Definition 2.10 may be used to investigate the di-stability region R. The

expressions in

(3.3)

and

(3.4)

look ratheruninviting to use for exploring this R. It is therefore possible toproceed as follows"

Write forshort

-2+qh2

1-(p/2)h

(3.5)

Kh

L

l+(p/2)h’

l+(p/Z)h

Itis easyto seethat

(3.6a)

(3.6b)

(3.6c)

h

Ir,

I,

Irhl

areboth <1 iff

Igl

1<Lh_-<: 1, h

[rl[, [rh[

are both >-1 iff

[Kh[+Lh<=-I

or

[Kh[-Lh<=I<=L

,

h

oneof

[rl]

Irl

is ->-1 and the other_-<1 iff

]Kh[>=

]Lh+ 1[.

h h

In order to find outthe values of p, q for which r and r2 are smaller or larger

than 1 in modulus, it is necessaryto investigate these inequalities separately for the case

+ph/2>O

and l+ph/2<O

(cf.

(3.5)).

In Fig. 3.1 the various regions in which

(3.2)

has similar dichotomy properties are

sketched,

using as abscissa qh2

and as ordinateph.

As

may beseen,the lines qh2=

_O,

_{qh2= 4 andph 0form the boundaries of these}

regions:

(3.7a)

_r

l, [rhl

are both 1 iff 0<

_qh2<4,

_ph>

_O,

(3.7b)

areboth >1 iff 0<_qh

_2<4,

_{ph <0,}

(3.7c)

h

one of

Irl

l,

[r[

is >-1 and the other is -<1 iff

]qh

-

₂₁

>- 2.

ph:2

ph=O

b

ph -2

ph

IX

qh

qh=O

qh:Z,

FIG.3.1

In greater detail the regions may be broken down as shown into

Ia"

_(qh2

_< 0,= ph>

_2)"

_r

_h>l=

_{1, 0 >--}

_rz

>_-l, h>₁

_0<rzh<l

Ib:

(qhZ=<0,]ph]<2)" r=

Ic"

_(qh20,

_ph

_<-2)"

rlh_-<-1, 0=<r2h<_-1, h<₁

_0>rzh>--I

II"

([qh

_2-21

<=2, ph>-_

_2)"

0_-<r, h

r2h<l

III"

(0

=<

qh2

<-_{2, O<-ph<=2,}

(ph)2+(qh2-2)2->4)

0<-rl, h

rzh>_l

IV:

(2=<qh2_-<4,

O<=ph<=2,

(ph)Z+(qh2-2)2->4)

O>=r,

h

V:

(ph->O,

(ph)2+(qh2-2)2<4)

]rh]=]r21<=l,

complex conjugates,

h h

VI"

(ph<=O,

(ph)2+(qh2-2)2<4)

]r[=[r2]_->l,

complex conjugates, VII"

(0=<

qh2=<

2, 0->ph>-2,

(ph)Z+(qh2-2)2->4):

rh,

r2= h

r2h<_l

VIII"

(2=<qh2=<4,

0=>.ph>-2,

(ph)Z+(qh2-2)2>=4)"

rl, h>₁ IX:

([qh

2-

₂₁

_<=

_{2, ph}

_<

_-2):

_r

h _-<-1, /’2

h<l,

_rzh<--I

Xa"

_(qh2->4,

ph->2)"

0=<

_r

h> _1,

_h<--I

Xb"

_{(qh2>-4,]ph]<2)"}

0->r=-r2= Xc"

_(qh2>4,

ph<-2)"

0>rh> 1, /’2---

h>l

It is interestingto compare Fig. 3.1 with Fig. 3.2, where regions are shown in which

p

FIG. 3.2

the original differential equation

(3.1)

has similar dichotomy properties.

As

before,

the lines q 0 and p 0form the boundaries of the regions:

(3.8a)

Re

(A1),

Re

(A2)

areboth =<0 iff q_->_0, p=>0

(B, C),

(3.8b)

Re (1),

Re

(2)

are both _>-0 _{iff q>=0, p-<0}

_{(D, E),}

(3.8c)

one of Re

(1), Re

(2)

is ->0 and the other is <-0 iff q<-0

(A).

Intheregions

C,

D

shown,

where4q->

p2,

A1,/2arecomplex conjugates,butelsewhere

they are real. The scales in Fig. 3.2 are somewhat arbitrary and the position of the

separatrix

p2=

4q does notchangeif

qh

2

andph areusedinsteadof q and p as abscissa and ordinate. The two figures arethus directlycomparable.

Ifq

<

0(region

A),

the differentialequation

(3.1)

hastwomodes: agrowingmode with

>

0andadecayingmode with

<

0.Forany positive value of

h,

the discretization method

(3.2)

also has two modes: agrowing mode with

Irl>

1 and a decaying mode with

]r

<

1 (regions I

a,

I

b,

IC).

Thus for all values of p and

h,

the method is di-stable for q

<

0. Results will however deteriorate for

Iphl

much greater than 2 as in regions

Ia, I where one of the two modes becomes oscillatory

(r

< 0),

althoughit still decays

inthe appropriate direction.

If q>0, p>0 (regions

B,

C)

the differential equation

(3.1)

has two decaying modes with Re

()<

0. If qh

<-4,

the discretization method also has two decaying

modes with

_Irl

<

(regions

II,

III, IV, V).

Similarly, if q

>

0, p

<

0 (regions

D, E)

the differential equation

(3.1)

has two growing modes with

Re

(A)>

0. If

qh<-4,

the discretizationmethod also has twogrowing modes with

_Irl

>

(regions

VI, VII, VIII,

IX).

Thus for all values ofp, the method is di-stable for

qh<=

4. Once again, results will deteriorate for

Iphl

greater than about 2, as in regions

II,

IX where one ofthe

twomodesbecomesoscillatory,whilethoseof the differentialequationdo not oscillate in regions

B,

E.

Clearly the method is not di-stable for

qh2>

4, as the differential equation still hastwo decaying modes if p

>

0, ortwogrowingmodes if p

<

0, but the discretization method has one decayingmode and one growingmode inregions X

a,

X

,

X

.

3.2.

One-step

difference equations for systems.

Most

BVP algorithms for systems of first-order differential equations are based on finding (directly or indirectly) a

one-step recurrence relation for approximate solutions onacertaingrid. Togetherwith

the BC this then leads to a linearsystem from which the (approximate) solution can be found. Examples are the Box scheme

[12],

higher-order difference schemes

[16],

collocation

[1], [22],

and multiple shooting

[6],

[19].

Inthis section consideration is

givento one-step difference equations that arise from one-step discretizationschemes,

such as

Runge-Kutta

or one-step Obrechkoff formulae. It is well known that a large

class of collocationmethods can be interpretedasimplicit

Runge-Kutta

methods

[27].

These methods have beenextensively investigatedfor theirIVP stability properties

(cf.

[1], [3],

[26]).

A

general

Runge-Kutta

method for the system of differentialequations

y’--f(t,

y)

may be written

(3.9a)

q

y)

+

h

_E

atf(

t,

+ cth

qt),

k 1,2,.

.,

m,

/=1

(3.9b)

y/h+l

yh h

₂

blf(

ti

+

clh, q,)

!=1

where the coefficients akl

(k,l=

1,2,...,

m),

bt

(1

1,2,...,

m)

define the

method,

Ck

Y=I

akt, k 1,2,’’’, m, and forconsistency

YI

bl

1. Clearly anypermutation of the rows ofthe matrix

A

{akl}

applied alsotoits columns and to theelements of the vectors b=

{b},

e=

{Ok},

will not change the method. When written in the form

(3.9),

suchamethod expresses

y+l

in terms of

y,

but for

BVPs

it isofequalinterest

to express

y

in terms of

y+,

as thoughthe integration were to beperformed inthe

opposite direction. Thus

(3.10a)

qk=Y)+,

-h

₂

(bl--akt)f(ti+l--(1--Cl)h, ql),

k= 1,2,"" ",m,

/=1

(3.10b)

y=y+-h

btf(ti+-(1-c)h,

/=1

and ifthe coefficients ak,

bt,

ck representamethod applied in the forward direction, the samemethod applied inthe backward direction is represented bythe coefficients

b, agl

),

hi,

(1

Ok).

DEFINITION 3.11.

A

Runge-Kutta

method will be called symmetric if it is the

samemethod whether applied in a forward or backward direction.

Thus,

letting e be the vector of rn unit

elements,

the method

(3.9)

is symmetric

ifthere is apermutation matrix P such that

(i)

ebr

a

PAP

,

(ii) b

Pb,

(iii)

e- e Pc.

As

defined in Definition 2.10, the concept of di-stability is given in terms of constantcoefficienthomogeneouslinear differential equations. The mostgeneral first-order system of this typemaybe written

y’=

Jy

whereJisaconstantmatrix,but since the method

(3.9)

onlyusesthe function

f(t,

y)

inalinearmanner,it is notaffectedby a linear change of variables which causes J to undergo a similarity transformation.

For practical purposes, it is sufficient to consider

J

a complex diagonal matrix, and

thus to examine the effect of the method on one single differential equation ofthe

form

y’=

Ay, as is normally done for linear stability analysis of methods for initial value problems.

When method

(3.9)

is appliedto thisstandard test equation, it is found that

(3.12)

yh+l

yh

r(hl

where

r(z)

is a rational function ofdegree notexceeding m in the numerator andthe

denominator. Explicitly

(cf. [3]),

(3.13)

r(z)

+

b

z(I-Az)-le

the denominatorbeingthe determinant ofI-

Az.

For

a symmetric

Runge-Kutta

method,

and indeed fora largerclass of

Runge-Kutta

methods which could be called linearly symmetric, the rational growth factor

r(z),

which characterizesthe fundamentalsolutions ofthedifference equation,is such

that

(3.14)

r(-z)=l/r(z).

For

one-stepObrechkoff formulae

(cf. [14])

tfa

(3.15)

y?+

yh +

₂

h \k.O/fi+l"1

]l/yl

1))

/=1

where

yt

is the lth derivative of

y(t)

obtained by repeated differentiation of the differential equation

y’=f(t, y)

at (t,y)=(ti,

y)).

Applying such a method to the standard test equation

y’=

&y,

it is again foundthat

(3.12)

is satisfied, where now

(3.16)

r(z)=

1+

_{Y ]3,,z’}

_{1-2 ]3o,Z’}

/=1 /=1

Thus such a method may be defined to be symmetric when

/31=(-1)-1/3o

(/=

1,2,...,

m)

which is equivalentto

(3.14).

There is clearlyaone-to-one relationshipbetween one-step Obrechkoff formulae and theirgrowth factors

r(z).

However,

ingeneral, many different one-step methods

may giverisetothe same rationalgrowthfactor

r(z).

Nevertheless,

their linear stability

properties suchasabsolutestability,A()-stability,and also di-stability will beentirely

determinedbythe function

r(z).

Thus,

by using Definition 2.10,aone-step method is dichotomically stable on a region Rc_{C if and onlyiffor all} z

R,

its growth factor

r(z)

satisfies

(3.17a)

Re

(z)O]r(z)

1,

(3.17b)

Re

(z)

Olr(z)l

1.

THZORZM 3.18.

A

one-step method that is symmetric and

A()-stable

is di-stable

on the region

R

:=

{z

C

_lRe

(z)l

sin

lira

(z)]

cos

} (0<

N

/2).

Proof

Sincethe method is

A()-stable,

the growth factor

r(z)

satisfies

Ir(z)N

1

on

R-:={zeCl-Re(z)

sinllm(z)cos}.

Since it is symmetric

r(z)l=

1/lr(-z)]

on

R+:

{zClRe

(z)sin

_lIm

(z)[

cos

}.

Thus

(3.17)

is satisfied on

R=R-UR

+.

Sowe may conclude thatastabilitypropertyinthe left halfplane plussymmetry

as investigated in

[1], [25]

gives indeed asufficient criterion for di-stability.

Example 3.19. The implicit midpoint rule

(which

is also the Box

scheme,

or a

collocationmethodat asingle Gauss quadrature point), orthe trapezoidalrule

(which

isacollocationmethodat twoLobattopoints,orasimpleone-step Obrechkoff

formula),

areboth implicit

Runge-Kutta

methods withgrowth factor 1

+5z

(3.20)

r(z)

These methods are A-stable

(A()-stable

with

/2)

and symmetric. Thereforethey

are di-stable on the whole of the complex plane C.

Example

3.21.

Any

one-step method whosegrowth factor

r(z)

is adiagonal Pad6

approximant to the exponential e is A-stable

[3].

These Pad6 approximants also satisfy

(3.14),

andsosuch methodsaredi-stable onthe whole of thecomplex plane C.

Since asymmetricmethod hasagrowthfactor which satisfies

(3.14)

it also satisfies

lr(z)l

2=

_{r(z)?(z)= r(z)r()= r(z)r(-z)=}

_{1 whenever}

_{Re (z)=}

_{0. Itis thustemptingto}

suppose that any symmetricmethod mightbe di-stable on the whole ofthe complex plane C. The following counterexample shows that this is not the case.

Example3.22.

A

one-stepmethod

(Runge-Kutta

orone-step Obrechkoff

formula)

whose rational growth factor is

4 1

+5z

+z

(3.23)

r(z)

4

Z

+

5Z

is symmetric and

A()-stable

with

/6,

butis not

A-stable,

andhence not di-stable

on the whole of the complex plane C. Itis in fact di-stable onthe region R :=

{z=

x

+

iy3y2-xZ+(4x4+

9}

which containsthe smallerregion

/

:=

{z

C

IRe

(z)l

=>

_IIm

(z)l

}.

General one-stepmethodswillhave their mesh size restricted forstabilityreasons. This is well known forexplicit formulae applied to initial value problems, in which case it is necessarythat, for each eigenvalue ) ofthe Jacobian matrix

J,

the product hA must lie in the absolute stability region

R:-{zC[]r(2)l<l}.

When applied to boundaryvalueproblems,implicit formulae also suffer from suchastability restriction

upontheir stepsize, asthe growing modesmust beproperly represented.

However,

a symmetricmethod which is alsoA-stablewill not suffer from any mesh size restriction

owing tostability, as it is di-stable onthe whole of the complex plane C.

A

symmetric

A(c)-stable

method will beequallyefficientprovided theeigenvalues of the Jacobian matrixJlie withinthe appropriate sectors of di-stability of thecomplex plane,although

these statements areto some extent dependentuponthe assumption that the matrix

J

is constant.Theboundaryconditions could also disturb the overall stabilityof such a discretization, as is well known in the case of the implicit midpoint rule and the trapezoidalrule,whichsuffer fromthesamekindof oscillationasthe central difference scheme for second-order equations.

3.3.

_Mtfltille

shooting.

As

illustrated in Example 2.1, the use of a difference equation of order higher than that of the differential system imposes the need for

additionalboundary conditions. Such a situation most naturally arises whenusing a

central difference scheme ofhigherorder than the differential equation, e.g., fourth-order(five-point) central differences forasecond-order differentialequation.

In

those

cases, an analysis similar to that performed in Example 2.1 should show where additionalboundary conditions are needed.

However,

when treating systems of first-order differential equations, shootingor multiple shooting areverynatural approaches. Onewayof viewing multiple shooting

is to think ofa basic discretization method

(one-step

or multistep) but to eliminate

from the algebraic equations the solution values internal to each shooting interval

(internal

condensation)

orall the internal variables in the case ofsimple shooting.

In

practice the difference equations are set up sequentially, and internal variables elimi-natedas soon astheyarenotneeded. The step size may also be determined and varied

dynamically.

If the basic discretization method is a one-step

scheme,

it is important that it should be dichotomically stable for the problem to which it is applied. Then any decaying (growing) mode of the differential system generates a decaying (growing)

mode of the basic difference equation, which is controlled by an initial

(terminal)

condition. Theonlydifference from direct solutionof theone-step differenceequation

(apart

from areduction in storage requirements) is thatsequential block elimination of internal variables may be an unstable process, and lead to swamping ofdecaying

modesbyrounding errors inthe growing modes.This is oneoftheprincipalmotivations behind multiple shooting, in which, unlike simple shooting, the sequential block elimination is not carried too far before a new uncontaminated set of fundamental solutions isagainrestored. The globaldifference equation, which relates values atthe ends of the shooting intervals, mustthenbe

solved,

by somestable recursion process, involving decoupling ofdecaying and growing modes andnot by direct block elimi-nation.

Since, in eachshooting interval,we aresolvinganumber of initial valueproblems,

it isnatural to considermultistep

schemes,

instead of one-step schemes,for use in the basic discretization.

As

always, the use of a multistep discretization method for a

system of first-order differential equations introduces spurious

modes,

and imposes

theneed for additionalboundaryconditions.

In

general,sincethediscretization should bedichotomically stableforthe problemto be

solved,

Definition 2.8providesacriterion as to how many additional conditions should be imposed at the beginning or at the

end ofthe interval--an initial conditionfor each spurious decaying modeand a terminal conditionforeach growing mode.

However,

inamultiple shooting context,amultistepdiscretizationschemewould be usedto solve initial value problems, and sequentially eliminate internalvariables.

Thus,

for practical reasons the additional boundary conditions should all be extra startingvalues, which are always required in conjunctionwithmultistep methods for initial value problems. These extra initial conditions could be generated by some one-step method of high

order,

or by low-order methods ofthe multistep family at

small step size as is done in automatic, variable-order, initial value integrators. The important point is thatthe multistep discretizationmethod togetherwiththe starting procedure

(which

provides theadditional initial

conditions)

should be dichotomically stablefortheproblemtobesolved. This means thatanyspurious mode of themultistep method should be ofthe decaying type.

These considerations leadustoconsider new stabilitypropertiesof discretization

formulae for solvinginitialvalueproblems.

In

the context of stiff initial valueproblems, the need to represent decaying modes by decaying approximations has lead to the definition of

A(a)-stability

(cf.

[5],

[29]).

Clearly,when considering multiple shootingforboundaryvalue problems, there

will normally be growing modes present,and it wouldbemeaninglessto require that a discretization represent them by decaying approximations, or be absolutely stable

forsome hA with positive real part. The first idea mightbe to consider some form of relativestabilityto ensure that numerical approximations do notgrowfaster than the true growingmodes. This leads tothe following definition ofR(fl)-stability.

DEFINITION 3.24. If a discretization

M

h

applied to the test equation

y’=

Ay generates a difference equation with basis solutions

_{(r)i},

j 1,2,...,

k,

it will be

called

R()-stable

if

Ir l<le l

(j-1,2,,,,

,k)

for al step sizes

h,

and A C such that hA

R+:={z6C[Re(z)

sin>JIn(z)Jcos}.

[3

This may be interpreted by saying that R/

lies within the white fingers ofthe order star ofthe method

(cf.

[26]).

Definition3.24does not contributetowards the dichotomic stability of the discretiz-ation.One important requirement fordichotomicstabilityisthat growingmodes should be represented by growing approximations, or that the discretization be absolutely

unstable forhA withpositiverealpart.This idealeadstothedefinitionofthefollowing

desirable property.

DEFINITION 3.25. If a discretization

M

h

applied to the test equation

y’=

Ay generates a difference equation with fundamental solutions

{(rh)i},

j 1,2,’" ",

k,

it

h

will be called

A(3,)-unstable

if oneofthe

r;,

say rlh satisfies

Irl

[>

1 for all stepsizes

h,

and A C such that hA

g+:={zC[Re(z)

sin

y>[lm

(z)[

cos

y}.

[3

Most

methods for stiff initial value problems are not A(y)-unstable, as they

concentrate on decaying

modes,

and often require that

_liml

_{I_}

_[r.[<l.

(j= 1,2,

,

k).

However,

at leastfor boundary value problems, it does seem that

A(y)-instabilityis a desirableproperty.

Nevertheless,

it is still not sufficient to determine a regionof dichotomic stabilityofthe global multiple shooting equations.

For

this it is

necessary

that,

while the dominant approximate mode grows, the spurious modes of the formula decay. Thus there should be a region R+ in which

[r[=>l

and

[r[_-

<_1,

j 2, 3,

,

k. Ifthe discretization satisfiestheseinequalitiesin

R

/

of Definition3.25,

and is also

A(

c

)-stable,

then the global equation is dichotomically stable on R

+,

togetherwiththe relevant region inthe left-hand side of the complex plane.

Example3.26. The methodsgiveninExample3.19 are

A(

c

_)-stable

withc

r/2,

and also

A(

y)

-unstable with 3,

7r/2,

since there is only one fundamental solution

{(r)i}={r(z)i},

and

_]r(z)]<l

when

Re(z)<0

while

]r(z)]>l

when

Re(z)>0.

However,

these methods are not R

(/3 )-stable

since

]r(z)]-->

oe when z-->2.

Example 3.27. The following one-stepmethods have growth factor

1

+

1/2z

+

z

(3.28)

r(z)

-1/2z+z

which is the Pad6

(2.2)

approximantto the exponential ez"

(i) A

one-step Obrechkoff formula:

1

,).

y/h+l

yh

_{+_}

h(y’i+

+

Y’i)

--

h2(y+l

_Y

(ii) The two-stage implicit

Runge-Kutta

method using Gaussian quadrature

points:

q Y

+h[f(ti+(-?)h,

ql)+(-7)f(ti+(l+"/--2

-7-

q2)]

zo]

h,

q:z=Yhi+h[(+-36

)f(ti+(--36 )h,

ql)

-4

y)+l

y)+h

[f(ti+

(-)h,

q,)

+

f(ti+

(+)h,

q)

].

(iii)

The three-stage implicit

Runge-Kutta

method using Lobatto quadrature points

[10]:

q=y)+h

5

_y:+f

_t+h,

qe

2 1 1

y)+=y)+h

_{y:+f t+h,}

qe

_+Y:+l

These methods are

A()-stable

with

=/2,

and

A(7)-unstable

with T=

/2.

In

addition theyare

R()-stable

with 16.9

.

In contrasttothe previous example,the

positive real axis is in awhite finger ofthe order star.

These exampleswere allof one-step methods withonlyonefundamentalsolution. The next examples are of genuine multistep methods.

Example3.29. Thefollowingtwo-step

methods,

when appliedtothe testequation

r,),(r)"

y’

Ay,

generateadifferenceequationwiththesamepairof basis solutions

(i) A

two-step Obrechkoff method

(Enright’s

two-stepsecond derivative

method)

[9]:

(294b

_{+5 -1}

)_1,2

y)+l

y)

+

_h=

y:+l

y

y:-

,

y+,.

1050 R. ENGLAND AND R. M. M. MATTHEIJ

(ii) A

hybrid implicit stiffly stable method

[7]"

21 3 1 3 q2

=

Yhi+l

+-

Yhi

--

Yhi-,---

hy’i+, 33

yh

7 7 1 5

--32

y+h

)1]

y)+=y)+h

_y:+f

_t+h,

q2

+gy:+

Forboth methods

(3.30a)

[

5

(3

43

A

1

)’/2]/[

(

29 1

)]

r=

l+hA+ l+-hA+h

----h3A

2

1-hA+-hZA

12 4 192 96 48 8

h 5 3 43 29

r2 1

+

hA 1

+-

hA

+ hA

---

h3A

2 1-hA

+- h2A

2

4 192 96 48 8

These methods aretherefore

A()-stable

with

/2,

and

R()-stable

with 17.4

as may be seen by determining the order star.

However,

they are not A(y)-unstable

since rh

,

r2h 0when hA

.

Iftheyareused with too largea step size h then allthe

basis solutions of the difference equation will decay, while some of those for the

original differential equation will grow.

Thus,

for the global difference equation, the region ofdi-stabilityisbounded in the positive half of thecomplex plane,andthe step size will be limited for stability reasons, as is that for the central difference method

(3.2)

when the coefficient q

>

0.

Example 3.31. Consider the three-step Obrechkoff method

(

2471 1

)19

y

y)+=y)+h

79

_{Y:+I+Y:-y:-I+Y:-2}

_-h2(y+,-

).

When applied to the test equation

y’=

Ay,

this method generates three fundamental solutions

{(

rhl)},

{(r)},

{(r)};

r

approximates eh for

_Ih]

small.

In

fact

r

eh

h h

O(h6A 6)

as hA 0. Thespurious roots

r),

r3 are both zero for hA 0.Also

r

1 and

h

r

0 as hA

.

It

may be checked by determining the absolute stability region

F2

and orderstar, that themethodis

A()-stable

with 88.3

,

and A(y)-unstablewith

=

/2.

Itis also

R()-stable

with

/12 (15 )

the largestpossible valueof for

a discretization method with fifth-order accuracy. Fuhermore the spurious roots

h h

satisfy

_]r]<l,

_]r2<l

for Re

(hA)0 (and

indeed for all values of

hA),

and so the

globalmultipleshooting equationshaveanunboundedregionofdi-stability containing

R-U

R+ where

R-:=

{z

c

_l-Re

(z)

sin88.3>

_Jim

(z)

cos

88.3},

e

+_:=

{z

C

]Re

(z)

>

0}.

Withthis

method,

the step size will not be limited forstability

reasons,

unless the Jacobian matrix J has eigenvalues very close

(within

2

)

to the imaginary axis, or unless stability is disturbed by either the variation of

J,

or by boundary conditions which do not actually control the modes of the discrete solution.

Very

few

(if any)

multistep methods of this typeareto be found inthe literature.

A

paper hasbeen published

[8],

which specifies precise familiesof such methods; we hope to usethese forsolving boundary valueproblems in the future.

4. Examples. This initial paper introduces the concept ofdichotomic stability, and its importance for numerical methods for the solution ofboundaryvalueproblems.

It

is notpossible,inalimitedspace,togive anexhaustiveanalysis oftheimplications

of the concept, oracomplete description ofa multiple shootingmethod which takes this analysis into account.

For

otherpapersthat extend the analysis andimplementation

aspects, see, e.g.,

[8].

Nevertheless,

this paper would not becomplete withoutan exampleto illustrate

(i) the necessity, for certain "stiff"problems, ofalargeregion of di-stability, and

(ii)

the failureof other

methods,

e.g.,A-stable

methods,

for suchproblems.Implicitly, the importance of dichotomic stability forglobal methods

(in

particular

collocation)

was demonstrated in

[1]. A

shootingapproach is considered here.

In

subsequent papers,

details are given as to how an integrator should select appropriate step sizes, which are commensuratewith theactivityofadesiredparticularsolution.

Here,

it isassumed thatasuitable constantstep sizecanbe

found,

and bounds upon its valueareconsidered.

Example 4.1. Considerthe following pairof ordinary differential equations:

(4.2a)

y,=[0

A]

[(l-A)

e

]

A 0

y+

(l-A)

e

(where

A is areal positive

constant),

withthe boundary conditions

(4.2b)

[10

00]Y(0)+[: 0]

y(1)=I

1]

1 e

The solution of

(4.2)

is

y(t)=(et,

e’)

,

while the basis solutions of the homogeneous partform a fundamental set

(t)

[

e-t

e ’l

_e-At

earl"

Problem

(4.2)

is well conditioned, with a condition number independent of A

(cf.

[18]),

and itmightbeexpectedthat for anygiventolerance

(TOL),

there should exist

amaximum step sizeh (independentofA such that theglobalerrorisboundedby TOL.

Considerthe use of the Backward Euler

method,

and a large value of A. Since

the step size h should notdepend upon

A,

the product Ah mayalsobeverylarge. The discretization of

(4.2a)

is

(4.3)

Y/h+l

Ah 0

[Y)

+

hfih+l]

where

(1

A)

eik

]

f)=

(l-ae

and ifAh is largethe discrete fundamental solution set is approximately

(4.4)

_(Ah)_

(_Ah)_

One ofthe difficulties ofshooting is clearly

absent,

as there isno growth of either of the discrete fundamental

modes,

which might swamp other modes or the particular

solution.

However,

at t--ih 1, and as h

,

[[cbhll-*0,

and it can be seenthatthe secondof the boundaryconditions

(4.2b)

fails to control either ofthe discrete funda-mental

modes,

unlessthestep size his alsoreduced sothat hhremains oforderunity.

For

largevalues of

hh,

the discrete problem is ill conditioned, and so the step size h

is highlydependent upon h.

In

thisexample,theeigenvalues ofthe Jacobian matrixJwere+/-h. _{Methods with}

A-stability, such as the Backward Euler

method,

are veryuseful forstiff initial value

problems, but their treatment ofeigenvalues with positive real parts is verybad. This is also true of any one-stepmethod whose growth factor

r(z)

is a Pad6 approximant

to the exponential e with the degree of the numerator smaller than that of the denominator.

For

large positive eigenvalues h, lackof di-stability forces the step size down in order to recover a well-conditioned discreteproblem, inmuch the same way as lack ofstiff stability does forproblems withlarge negative eigenvalues.

Example4.5. Consider the same continuousproblem

(4.2);

but now it is discretized withthe trapezoidal rule

(3.19)

which is di-stable on the whole of the complexplane C. The discretization is

]

-

[[’]-

1 hh 1

o

1

01] [0

A0h]

1

0]

[0

Ah

]]yi+hgi

_}

Ah

where

g)=1/2(f)+fL1),

and if lh is large the discrete fundamental solution set is

approximately

(4.7)

,h[

(--1)

(-1)

]

-(-1)

(-1)’

For

such largevalues of

Ah,

these discrete fundamentalmodes do not provide agood approximationto the continuous fundamental

modes,

but they are good enough for

theboundaryconditions to control them.

It

is straightforwardto seethatthe resulting discreteproblemis well conditioned. Givenany tolerance

(TOL),

the maximum step

size

h,

suchthattheglobalerror isbounded by

TOL,

depends only upontheparticular

solution of

(4.2),

and not upon A.

Furthermore,

there is no difficulty associated with the solution of the discrete problem bysimple shooting.

It

maybenotedthat dichotomicstabilityis not a sufficientcondition toguarantee

well conditioningof the discreteproblemalthough it is anecessary condition, atleast when the eigenvalues A of the Jacobian matrix become large, whether with positive

ornegative realparts. Ifboundaryconditions

(4.2b)

werereplaced with

[1

00]

[0 00]

[1]

(4.8)

0

y(0)+

1

y(1)

e

both the solution, and the conditioning of the continuous problem would remain

unchanged.

However,

the boundary conditions

(4.8)

would not actually control the discrete fundamental modes of

(4.7),

andthe resulting discreteproblem would be ill

conditioned, inspite ofthe di-stabilityof thetrapezoidalrule. The difference between

(4.2b)

and

(4.8)

is similar to that occurring with the central difference scheme when the number N ofsubintervals is changed from an oddto an even value.

5. Conclusion.

A

problem with fast growing and decaying fundamental modes

maybe,ifappropriate boundaryconditions aregiven, averywell-posed problem,and this would be reflected in a small condition number

[18].

If a discretization is used which does not approximate the fundamental modes of the continuous problem by

similar

(growing

or

decaying)

modesof the discreteproblem,thenboundaryconditions that control the continuous modes cannot control the discrete modes. The resulting

discrete problemwould then be ill conditioned, with a large condition number

[11],

and small discretization errors would give rise to large global errors. Even if care is taken in the boundary layers, with accurate approximations obtained by using

sufficiently smallsteps, theuseoflargersteps outside theboundary layerscouldchange

the nature of some fundamentalmodes, fromgrowing to decaying, orvice versa.The

resulting discrete problem would not have the correct dichotomy, and large errors could result.

Dichotomic stability, as defined in Definition 2.8 or 2.10, is a property of the global discretization method for a boundary value problem.

It

guarantees that the fundamental modes ofthe continuous problem are approximated bythepropertype

(growing or

decaying)

ofmode forthe discreteproblem. It appears to be anecessary

condition to ensure that the conditioning of the discrete problem is not worse than that of the continuous problem. Itis not, in itself, a sufficientcondition, withoutany

consideration ofthe boundaryconditions. Itis notimpossible foradi-stable discretiz-ation to produce an ill-conditioned discrete problem from a well-conditioned

con-tinuousproblem.Thiscanhappenifthe fundamental

modes,

while having the correct type

(growth

or

decay)

ofbehaviour, areneverthelessdistorted

(in

the n-space ofthe dependentvariables

y)

insuchaway thatthe

BCs,

whilebeingcorrectin

number,

do not actually control the discrete modes atthe correct end of the interval.

Inthe context ofmarching, or multiple shooting, type

methods,

the requirement

for dichotomic stability of the global discretization impliesthe sameproperty for the basic discretization

scheme,

which is used in a marching mode to solve

IVPs.

A(c)-stabilityandA(y)-instability

(3.25)

thenjointlycontributetowardanunboundedregion

ofdistability,but arenotofthemselves sufficient,since A(y)-instabilitydoes notplace

any condition upon the spurious modes

(if any)

ofthe basic, discretization.

R(/3)-stability is also defined

(3.24),

but is of interest principally for

IVPs,

and does not contribute towards the di-stability ofa discretization for a BVP.

Aeknowlelgments. The authors thank the Instituto de Investigaciones en

Matemticas Aplicadas y en Sistemas, Universidad Nacional Aut6noma de M6xico,

and Mathematisch Instituut Katholieke Universiteit Nijmegen, whose funding made possible the contacts which led to this paper. Theyare also grateful to

ZWO,

which funded alongervisitby the firstauthorto the

Netherlands,

thus enabling completion

of the work.

REFERENCES

[1] U. ASCHERAND R. WEISS,Collocation_forsingular perturbationproblemsI,SIAMJ. Numer. Anal.,

20(1983),pp. 537-557.

[2] C.DEBOOR,F.DEHOOG,ANDH. KELLER, The stability_ofonestep schemes_{for first-order}two-point boundaryvalueproblems, SIAM J. Numer. Anal.,20(1983),pp. 1139-1146.

[3] J. C. BUTCHER, Runge-Kutta methods, in Modern Numerical Methods for Ordinary Differential Equations,G.Hall andJ. M. Watt eds.,ClarendonPress,Oxford, 1976, Chaps. 5, 10.

[4] W. A. COPPEL, Dichotomiesin stability theory, Lecture Notes in Mathematics 629, Springer-Verlag, Berlin, 1978.

[5] G. G. DAHLQUIST, Aspecial stabilityproblem_forlinear multistepmethods, BIT,3 (1963),pp. 17-25. [6] R. ENGLAND, Aprogramforthe solution_ofboundaryvalueproblemsforsystems_ofordinarydifferential

equations,Culham Laboratory Report PDN3/73, Culham, U.K., 1976.

[7] Some hybrid implicit stiffly stable methodsforordinarydifferentialequations,Numerical Analysis

Proc.,Cocoyoc, Morelos, J. P. Hennart,ed.; Lecture NotesinMathematics909, Springer-Verlag, Berlin, 1982, pp. 147-158.

[8] R.ENGLANDANDR. M. M.MATTHEIJ,Discretizationswithdichotomic stabilityfortwo-pointboundary valueproblems,inNumericalBoundaryValueODE’s,U.M.Ascher andR. D. Russell, eds.,Series

ProgressinScientificComputing, Vol.5, Birkh/iuser, Boston(1985),pp. 91-106.

[9] W. H. ENRIGHT, Second derivative multistep methodsfor stiffordinarydifferentialequations,SIAMJ.

Numer.Anal., 11 (1974),pp. 321-331.

[10] J. P. HENNART AND R. ENGLAND, A comparison between several piecewise continuous one step

integration techniques, Working Papers for the 1979SIGNUM Meetingon Numerical Ordinary

DifferentialEquations, Champaign,IL,R. D. Skeel, ed.; Report963,Dept.ofComputerScience, Univ.of Illinois atUrbana-Champaign, 1979, pp. 33.1-33.4.

[11] F. R. DE HOOG AND R. M. M. MATTHEIJ, On dichotomy and well-conditioning in BVP, SIAM J. Numer. Anal.,24 (1987),pp. 89-105.

[12] H. B. KELLER, Numerical Solution _ofTwoPoint BVP, CBMS-NSFRegional Conference Series in Applied Mathematics24,Society for Industrial and AppliedMathematics,Philadelphia,PA,1976. [13] B.KREISSANDH.-O.KREISS,Numericalmethodsforsingular perturbationproblems, SIAM J. Numer.

Anal., 18(1981),pp. 262-276.

[14] J. P.LAMBERT,Computational MethodsinOrdinary_DifferentialEquations, JohnWiley,London,1973. 15] M.LENTINI, M.R.OSBORNE,ANDR. D. RUSSELL,Theclose relationships between methodsforsolving

two-pointboundaryvalueproblems, SIAM J. Numer. Anal.,22 (1985),pp. 280-309.

[16] M. LENTINIANDV. PEREYRA,"PASVA4";anODE boundarysolver_forproblemswith discontinuous

interfacesand algebraic parameters, MatematicaApl.Comput.,2 (1983),pp. 103-118.

[17] R. M. M. MATTHEIJ,Estimates_fortheerrors inthe solution_oflinearboundaryvalueproblems,dueto

perturbation, Computing, 27(1981),pp. 299-318.

[18] The conditioning_oflinear boundary valueproblems, SIAM J. Numer. Anal., 19 (1982), pp. 963-978.

[19] R. M. M. MATTHEIJANDG. W. M. STAARINK, Anefficientalgorithm_forsolvinggenerallineartwo

pointBVP,SIAM J.Sci. Statist. Comput.,5(1984),pp. 745-763.

[20] R. M. M.MATTHEIJ, Decoupling and stability_ofBVPalgorithms,SIAMRev.,27(1985),pp. 1-44. [21] R. E.O’MALLEY,JR.ANDJ. E.FLAHERTY,Analyticaland numerical methodsfornonlinear singular

singularlyperturbedinitialvalueproblems, SIAM J. Appl. Math.,38 (1980),pp. 225-248. [22] R. D. RUSSELL, Acomparisonofcollocationand_{finite differencesfor}two-pointBVP,SIAM J. Numer.

Anal., 14(1977),pp. 19-39.

[23] R. J. SACKERANDG.R. SELL,Singular perturbations and conditionalstability, J.Math. Anal.Appl., 76(1980),pp. 406-431.

[24] L. F.SHAMPINEANDG.W.GEAR,A user’sview_ofsolvingstiffordinarydifferentialequations,SIAM

Rev.,21 (1979),pp. 1-17.

[25] M.VANVELDHUIZEN,Some resultson symmetricdiscretizations_forastiffmodelproblem,Report208, VrijeUniv.,Amsterdam,1982.

[26] G.WANNER,E.HAIRER,ANDS. P. NRSETT,Orderstarsandstabilitytheorems,BIT,18(1978),pp. 475-489.

[27] R.WEISS, The application_ofimplicit Runge-Kuttaand collocation methodstoboundaryvalueproblems, Math.Comp.,28(1974),pp. 449-464.

[28]

,

Ananalysisofthe box and trapezoidal schemesforlinear singularlyperturbed boundaryvalue problems,Math.Comp.,42 (1984),pp. 41.-67.

[29] O. B.WIDLUND,A noteonunconditionally stable linear multistepmethods, BIT,7 (1967),pp. 65-70.