NASA Technical Memorandum 102820

DynamicalApproachStudyof

SpuriousSteady-State

Numerical Solutionsof

Nonlinear Differential Equations

Part I -The ODEConnectionand Its

Implications for Algorithm

Developmentin Computational

Fluid Dynamics

H. C. Yee, Ames Research Center, Moffett Field, California

P. K. Sweby, University of Reading, Whiteknights, Reading, England

D. F. Griffiths, University of Dundee, Dundee, Scotland

April 1990

NationalAeronautics and

Space Administration

Ames Research Center

MoffettReid, California 94035-1000

DYNAMICAL APPROACH STUDY OF SPURIOUS STEADY-STATE

NUMERICAL SOLUTIONS OF NONLINEAR DIFFERENTIAL EQUATIONS

I. The ODE Connection and It_ Implications ]or Algorithm

Development in Computational Fluid Dynamics

H.C. Yee 1

NASA Ames Research Center, Moffett Field, CA 94035, USA

P.K. Sweby 2

University of Reading, Whiteknights, Reading RG6 2AX, England

and

D.F. Griflit hs a

University of Dundee, Dundee, DD1 4HN, Scotland

ABSTRACT

The goal of this paper is to attempt to give some insight and guidelines on the

application of nonlinear dynamic theory to the better understanding of steady-state

numerical solutions and nonlinear instability in algorithm development for nonlinear

differential equations that display genuinely nonlinear behavior in computational sciences and, in particular, computational fluid dynamics (CFD). This stems from the

fact that, although the study of nonlinear dynamics and chaotic dynamics for nonlinear differential equations and for discrete maps have independently flourished rapidly

for the last decade, there are very few investigators addressing the issue on the connection between the nonlinear dynamical behavior of the continuous systems and the

corresponding discrete map resulting from finite-difference discretizations. This issue

is especially vital for computational sciences since nonlinear differential equations in

applied sciences can rarely be solved in closed form and it is often necessary to replace

them by finite dimensional nonlinear discrete maps. In addition, it is also important to

realize that these nonlinear discrete maps can exhibit a much richer range of dynamical

behavior than their continuum counterparts.

Furthermore, it is also very important to identify some of the implications of what

happens when linear stability breaks down for problems with genuinely nonlinear behavior. Studies indicate that for relatively simple nonlinear ordinary differential equations

1 Research Scientist, Computational Fluid Dynamics Branch.

Lecturer, Department of Mathematics.

3Senior Lecturer, Department of Mathematical Sciences.
(ODEs) and well-known time-discretization with modest step-sizes some schemes can

converge to a spurious (false) steady-state solution in a deceptively "smooth" manner.

In some instances, spurious steady states may appear below the linearized stability

limit of the schemes, and consequently computation may lead to erroneous results. Our

preliminary studies on partial differential equations (PDEs) also show that much of

nonlinear dynamic (e.g. chaotic) phenomena have a direct relation for problems contanning nonlinear source terms such as the reaction-diffusion, the reaction-convection or

the reaction-convection-diffusion equations. Here our object is neither to provide theory nor to illustrate with realistic examples the connection of the dynamical behavior

of practical PDEs with their discretized counterparts, but rather to give insight into

the nonlinear features unconventional to this type of study and to concentrate on the

fundamental ideas. Thus, in order to bring out the special properties, the illustrations

center on simple scalar differential equation (DE) examples in which the exact solutions

of the DEs are known.

2
I. INTRODUCTION

While the applied computational fluid dynamicists are busy developing numerical

solver computer codes, grid generation codes, three-dimensional graphical stereo displays, and stretching the limits of the faster supercomputers in the world to numerically

simulate the various 3-D complex aerodynamic configurations [1], there is a group of

applied mathematicians, physicists, chemists, biologists, applied mechanicians, and meteorologists who are involved in a new science called "chaotic dynamics" (or nonlinear

dynamics). The science of chaotic dynamics has cut across many traditional scientific

disciplines for the last decade since chaotic dynamics is a science of the everyday world.

It offers a way of seeing order and pattern where formerly only the random, the erratic,

and the unpredictable were present. It explains much of the genuinely nonhnear phenomena that were once unexplainable. See references [2-10] for an introduction to this

subject.

Nonlinear Dynamics £_4Chaotic Dynamics: Before the birth of chaotic dynamical theory, traditional study of nonlinear dynamics belonged to the applied mechanics disciplines of mechanical engineering. Modern nonlinear dynamics (since the late seventies)

includes chaotic dynamics. Thus, unless otherwise stated, the term nonhnear dynamics

and chaotic dynamics are used interehangeably. That is, nonlinear dynamics includes

chaotic dynamics and vice versa.

Loosely speaking, the study of asymptotic behavior (steady-state solutions) of nonlinear differential equations (DEs) and nonlinear discrete maps (difference equations)

and how the asymptotes change as parameters of the system are varied is most often referred to as nonlinear dynamic analysis and chaotic dynamic theory. Topics in

this area include bifurcation theory, period doubling cascades resulting in chaos, etc.

Stable chaotic solutions (chaotic attractors) may be defined loosely and simply as stable asymptotes that have infinite period and yet are still bounded. It is emphasized

here that unless otherwise stated, all DEs and discrete maps are nonhnear and consist

of system parameters, and the terms discrete maps and difference equations are used

interchangeably.

Types of Dynamical Systems: Consider an ordinary differential equation (ODE) and a

partial differential equation (PDE) of the forms

du

d--i= aS(u), (1.1)

Ou Of(,,)

_- + Ox - eO_ -5 + aS(u), (1.2)

where o and e are parameters and S is a nonlinear function in u and is independent

of a (and e). The function .f(u) can be linear or nonlinear in u. An ODE of this

form in which t does not appear explicitly in S is called an autonomous dynamical

3
system. One can also consider a function S which is nonlinear in _t and depends on

_. ODEs of this type are called nonantonomous dynamical systems and they are more

difficult to analyze; see references [5,8] for a discussion. The analysis would be more

complicated if S = S(u, ct) is nonlinear in both u and a. In this case, the DE is not only

nonlinear in the dependent variable u (and independent variable t), but also nonlinear

in the parameter space a. One can also consider systems that depend on more than one

parameter and/or systems of equations of the above type.

Next consider nonlinear discrete maps (nonlinear difference equations) of the forms

u n+1 = u '_ + D(u n, u n-I , r), (1.3)

and

=y+' = u? + a(uy, uj l,r). (1.4)

Here r is a parameter, and D is nonlinear in u '_ and u n-1 and linear or nonlinear

in the parameter space r. The situation is similar for the function G. One can also

consider discrete systems that depend on more than one parameter. A typical example

is a discrete map arising from a finite-difference approximation of DEs such as (1.1) or

(1.2). For the ODE, the resulting discrete maps might be nonlinear in ct as well as the

time step At, depending on the ODE solvers. For the PDE, again depending on the

differencing scheme, the resulting discretized counterparts can be nonlinear in a, At,

the grid spacing Az and the numerical dissipation parameters even though the DEs

consist of only one parameter or none.

One can also consider discrete maps (scalar or system) of the forms

u = u" + D(u "+k, ..., us, ..., us-z , rl, r2, ...,

where k, l, m are positive integers and rl ,r2, ...,r,,_ are parameters, and

(1.5)

u_÷l • = uj,_ + ..1, tJ[uj+n+k1 ,...,uj±l,...,uj±l, n ,_-I u_+k, ...,uj,...,uj n ,_-t ,rl,r2,...,r,,). (1.6)

Again, (1.6) can depend on more than the three indices j,j -4-1. Systems (1.4) and (1.6)

are sometimes referred to as a partial-difference equation. The dynamical behavior of

(1.4) and (1.6) can be many orders of magnitude more difficult than (1.3) and (1.5).

Any of the systems (1.1)-(1.6) are examples of dynamical systems.

Important Consideration: It is emphasized here that discrete maps, regardless of their

origin, are dynamical systems on their own right. It is also important to distinguish the

following five types of discrete maps:

1. Discrete maps arise naturally in physical sciences. They commonly arise through

the inability to measure populations at all points in space and time [5,10,11] in population dynamics. They can also arise through the study of periodic excitation of dynamical

systems [12,13] in applied mechanics.
2. Discrete maps arisefrom Poincar6sectionsin ODEs [2].

3. Discrete maps arisefrom discrete approximationsof ODEs.

4. Discrete maps (partial-differenceequations)arisefrom temporal and spatial finite

differenceapproximations of PDEs.

5. Discrete models arise from the "Inverse Problems of Nonlinear Dynamics" in time

series analysis of observable data or experiments [9].

Discrete maps of types 1 and 5 sometimes might not have any relationship with a specific continuum DE. As a matter of fact, there might be no concrete associated governing

equations (continuum or otherwise) to start with for type 5 except the surrogated discrete map arising from the time series analysis. Type 2 arises naturally from the study

of dynamical behavior of nonlinear ODEs. However, types 3 and 4 have an intimate

link (but with a different tie than type 2) between the original governing continuum

DEs and their discretized counterparts. Furthermore, it is important to distinguish the

complexity involved in the analysis of types 3 and 4. Type 4 involves spatial as well as

temporal dynamical behavior.

Note that for discrete maps of types 3 and 4, even though the DEs might be linear in

the parameter space, depending on the numerical methods, the discretized counterparts

might be linear or nonlinear in that parameter space. In addition, extra parameters

which may be linear or nonlinear can also be introduced by the scheme as noted in

the paragraph after equation (1.4). An important concept is that even though the

DE does not depend on any parameter, its discretized counterpart does depend on at

least one parameter. As can be seen in the subsequent sections, the nature of the

dynamical behavior of these discrete maps is strongly influenced by properties of the

numerical method and the types and forms of nonlinearity on the DEs. Furthermore,

when dealing with nonhnear conservation laws of PDEs, the dynamical behavior of the

discretized counterparts is also strongly influenced by elements such as conservation and

nonlinearity of the schemes, and treatment of the source terms [14-18]. These issues are

very crucial for the existence of spurious steady-state numerical solutions which will be

explained in a later section. Here the term nonlinear scheme refers to a case where the

resulting discrete maps are nonlinear when apphed to scalar constant coefficient linear

DEs.

Objectives: Our ultimate objective is to conduct long term basic research on the interdisciplinary field of integrating the theory of nonlinear dynamics with computational

sciences and, in particular, with computational fluid dynamics (CFD). This new approach to CFD is extremely difficult and complex to analyze. A summary of the difficulty involved was discussed in Yee [14] and will be elaborated in sections IV and V. Our

immediate goal is to study the behavior of spurious steady-state numerical solutions for

nonhnear DEs and the dynamical behavior of this type of numerical solutions. Even

within this frame work, the subject is still very young, board, difficult and unfamihar

to computational scientists as well as researchers working in nonhnear dynamics and
nonlinear physics.

The intent of this paper is to give a flavor of the subject, to familiarize the reader

in computationl sciences with this new and exciting area, and most of all, to explain

through simple illustrations why it is so important for computational scientists to learn

about the subject. Some challenging topics for future research are also proposed.

Because of the complexity involved, there is a vast difference in the degree of difficulty

on the study of the subject between the discretized counterparts of nonlinear ODEs and

the discretized counterparts of nonlinear PDEs. In order to achieve our final goal of

studying the dynamical behavior of numerical methods for nonlinear PDEs that arise

from, e.g., computational fluid dynamics (CFD), we have to first fully understand this

subject on the time discretization and later link this knowledge to the study of both

the temporal and spatial nonlinear dynamical behavior of finite-difference methods for

nonlinear PDEs of the nonhomogeneous hyperbolic and parabolic types.

Therefore, the content of this paper will concentrate on the dynamical behavior of

time discretization for ODEs or systems of ODEs obtained from time-splitting [19]

or method of lines [20] for PDEs, and emphasis will be placed on its implication for

algorithm development in CFD and computational sciences in general. Hopefully this

will be part I of a series of many future research papers to come under the same topic.

Our companion paper [21] studies the dynamical behavior of the class of explidt RungeKutta methods in detail. The intent of this paper is to not only serve as a study

of the state-of-the-art of nonlinear dynamical behavior of ODE solvers, but also more

importantly to serve as an introduction and to present new results to motivate this vast,

new yet unconventional concept. Thus the mission of this paper is no__ttto provide the

answer or theory or to illustrate the connection of dynamical behavior of practical PDEs

to their discretized counterparts, but rather to gain insight into the nonlinear features

unconventionaJ to this type of study and concentrate on the fundamentals. In order

to bring out the new features, the iUustrations concentrate on the simple scalar DEs

examples in which the exact solutions of the DEs are known.

Outline: The outline of the paper is as follows: First, a brief background, motivation and

basic ideas will be given. Then some typical characteristics of dynamical systems with

genuinely nonlinear behavior will be discussed. Next, the dynamical behavior of discrete

maps arising from time discretization of ODEs will be studied and the main results and

their implications for computational sciences will be described. Studies on discrete maps

arising from finite-difference approximations of PDEs will not be elaborated. Rather,

the level of complexity involved and state-of-the-art study on this subject will be briefly

described. The paper will conclude with a few recommendations. Remarks will be given

on the popular misconception of residual test for convergence in steady-state solution via

the "time-dependent" approach and the popular misconception of the use of the "Inverse

Problems of Nonlinear Dynamics" to analyze the dynamical behavior of time series data

from a computer code in an attempt to learn about the true physical solution behavior

of the governing PDEs. This application of time series analysis can be misleading and

a wrong conclusion can be reached if the practitioner does not know by other means

6
other than the numerical solutions the exact solution behavior of the PDEs .

7
II. MOTIVATION & RELEVANCE

As discussed in the introduction, dynamical systems occur in the form of DEs and

discrete maps. In order to motivate why the study of numerical analysis will not be

complete without the utilization of the nonlinear dynamic approach, and to convey to

practitioners in computational sciences the importance of distinguishing the difference

between weakly nonlinear problems and genuinely nonlinear problems, this section is

devoted to a discussion of the typical behavior of dynamical systems with genuinely nonlinear behavior and the basic characteristic difference in dynamical behavior between

DEs and discrete maps in general. This discussion leads to the key elements of this

paper, namely: (1) to establish the connection between the DEs and their discretized

counterparts and (2) to convey to computational scientists how one should change the

traditional way of thinking and practices when dealing with genuinely nonlinear problems.

2.1. Typical Characteristics of Dynamical Systems

with Genuinely Nonlinear Behavior

The terms "nonlinear behavior" and "genuinely nonlinear behavior" are used quite

often in the literature and there seems to be no unified exact definition or meaning

[9]. Here these terms are used for nonlinear dynamical systems that exhibit mainly the

following characteristics.

(1) The study of nonlinear dynamics most often emphasizes the importance of obtaining a global qualitative understanding of the character of the system's dynamics

since local analysis is not sufficient to give the global behavior of genuinely nonlinear

dynamical systems. As a matter fact, this is one of the major reasons why sometimes

it required orders of magnitude more work than solving their linear counterparts.

(2) Unlike linear or weakly nonlinear problems, the solutions of genuinely nonlinear

DEs and discrete maps are strongly dependent on initial data, boundary conditions and

system parameters.

(3) Only genuinely nonlinear dynamical systems can have chaotic behavior and one

of the striking characteristics of chaotic behavior is sensitivity of the solution to initial

data. This characteristic is independent of whether the dynamical system is a continuum

or a discrete map.

From here on, the terms "dynamical systems with genuinely nonlinear behavior" and

"genuinely nonlinear dynamical systems" are used interchangeably. For convenience,

the word "genuinely" is omitted in most parts of the paper.

8
2.2. Typical Difference in Dynamical Behavior of

ODEs and Discrete Maps

The study of discrete maps is the discrete analog to the study of ODEs, as the

study of recursion formulas is a discrete analog to the study of series expansions of

functions. Much of the theory of ODEs can carry over to discrete maps with some

slight modifications. However, there are new phenomena occurring in discrete maps

which are absent in differential systems [22,23,12,13].

With repect to the topographical behavior, there are new kinds of behavior of trajectories in the neighborhood of equilibrium points of discrete maps. The behavior of

separatrices associated with a saddle type of equilibrium point for a nonlinear difference system is far more complicated than the behavior of separatrices for a differential

system. See Yee, Hsu and Hsu et al. [12,13,24,25] for details and examples.

With respect to similar equation types, the minimum number of first-order nonlinear

autonomous ODEs is three for the existence of chaotic phenomena. However, a simple

scalar first-order difference equation [26-30] like the logistic map

Vl'L

v n+' = ttv'_(1 - -_--), p a parameter, (2.1)

or its piecewise linear approximation [31]

V n+l "- #V n, V n <__ l

= p, 1 -__v '_ < 3 (2.2)

= g(4 - v n) 3 __ v n.

possesses very rich dynamical behavior such as period-doubling cascades resulting in

chaos. Equation (2.2) has the same behavior as (2.1) except that simple closed form

asymptotic solutions of all periods can be obtained. These characteristic trade differences between ODEs and discrete maps are very general. The discrete maps can arise

from any of the five types as discussed in the introduction. It is in this spirit that we

say that discrete maps can exhibit a much richer range of dynamical behavior than

DEs. The next two sections focus on the typical difference and connection between the

dynamical behavior of ODEs and their discretized counterparts.

2.3. Background and Motivation

Spurious asymptotic numerical solutions such as chaos were observed by Ushike [32]

and Brezzi et al. [33] on the leapfrog method for the logistic ODE

dl/

-- = au(1 - u). (2.3)

dt
In reference [34], Schreiber and Keller discussed the existence of spurious asymptotic

numerical solutions for a driven cavity problem. Some related studies are reported in

[35].

Spurious solutions of Burgers' equation and channel flows have been studied and computed in [36-38]. Many other investigators in the computational sciences (e.g. [39-43])

have observed some kind of strange or chaotic behavior introduced by the numerical

methods but were not able to explain systematically the source, the cause of their

results, or most of all the implication and impact in practical applications in computational sciences. Due to the popularity of searching for chaotic phenomena, it is very

trendy to relate inaccuracy in numerical methods with the onset of chaos. It is emphasized here that inaccuracy in long time integration of discrete maps resulting from

finite discretization of nonlinear DEs comes in other forms prior to the onset of chaotic

phenomena. Stable and unstable spurious steady states and spurious periodic numerical solutions set in before chaotic behavior occurs. These spurious asymptotes of finite

period are just as inaccurate as chaotic phenomena as far as numerical integration is

concerned. In other words, the prelude to chaotic behavior is the key element that we

want to stress (i.e., before the the onset of chaos or a divergent solution) since the result

of operating the time step beyond the linearized stability limit is not always a divergent

solution in genuinely nonlinear behavior; spurious steady-state solutions can occur. As

can be seen at a later section, this behavior is more difficult to detect than chaotic

phenomena in practical computations.

Recently, it has been realized by numerical analysts that numerical methods for ODEs

and PDEs can be considered as dynamical systems. Several papers [44,45] on numerical methods as dynamical systems have appeared in recent years. These investigators

studied the dynamical behavior of the different ODE solvers per se without relating its

dose tie with the ODEs themselves. Although the study of chaotic dynamics for nonlinear differential equations and for discrete maps have independently flourished rapidly

for the last decade, there are very few investigators addressing the issue of the connection between the nonlinear dynamical behavior of the continuous systems and the

corresponding discrete map resulting from finite difference discretizations. This issue

is especially vital for computational sciences since nonlinear differential equations in

applied sciences can rarely be solved in closed form and it is often necessary to replace

them by finite dimensional nonlinear discrete maps. Most often, typical apphed scientists rely on numerical methods to give insight into the solution behavior of nonlinear

DEs. It is not always clear how well a numerical solution can mimic the true physics of

problems that possess genuinely nonlinear types of behavior.

Why is there such a need to study the connection between the continuum and its

discretized counterparts for CFD applications? This stems from the fact that current

supercomputer power can perform numerical simulations on virtually any simple 3-D

aerodynamic configuration and, due to the limited available experimental data, the

applied engineers are relying on or trusting the numerical simulations whole heartedly

to help design our next generation aircraft and spacecraft. However, many of these

10
applied scientists are still using hneaz analysis as their guide to study highly nonlinear

equations, and most often they are not aware of the limitations and pitfalls of many of

the numerical procedures. Furthermore, most of the numerical algorithms in use operate

under the accuracy and stability limit guidelines of the linearized model equation. It

is only appropriate to analyze nonlinear problems with the nonlinear approach; i.e., by

the nonlinear dynamic approach.

The unique dynamical property of the separate dependence of solutions on initial

data for the individual nonlinear DE and its discretized counterpart is especially important for employing a "time-dependent" approach to the steady state with given initial

data in hypersonic CFD. In many CFD computations, the steady-state equations are

PDEs of the mixed type and a time-dependent approach to the steady state can avoid

the comphcation of dealing with elliptic-parabolic or elliptic-hyperbolic types of PDEs.

However, this time-dependent approach has created a new dimension of uncertainty.

This uncertainty stems from the fact that in practical computations, the initial data

are not known and a freestream condition or an intelligent guess for the initial conditions is used. In particular the controversy of the "existence of multiple steady-state

solutions" through numerical experiments will not be exactly resolved until there is a

better understanding of the separate dependence on initial data for both the PDEs and

the discretized equations.

2.4. Connection Between the Dynamical Behavior of the

Continuum and Its Discretized Counterpart

Aside from truncation error and machine round-off error, a more fundamental distinction between the continuum and its discretized counterparts is new behavior in the

form of spurious stable and unstable asymptotes created by the numerical methods.

This is due to the fact that nonlinear discrete maps can exhibit a much richer range of

dynamical behavior than their continuum counterparts as discussed in section 2.2. Some

instructive examples will be given in section III. These new phenomena were partially

explored by the University of Dundee group [46-54], Sanz Serna [55], Iserles [56,57] and

Stuart [58-62]. Their main emphasis was on phenomena beyond the linearized stability

bruit. The main contribution of our current study is (1) the occurrence of spurious

steady-state numerical solutions below the linearized stability limit of the scheme for

genuinely nonlinear problems, (2) the strong dependence of numerical solutions on the

initial data, as well as other system parameters of the DEs such as boundary conditions

and numerical dissipations terms, and (3) the implications for practical computations

in hypersonic CFD.

Before discussing the numerical examples, the next two subsections will give an overall

summary of our current findings (integrating with other relevant recent results). The

discussion is divided into steady-state solutions and asymptotes of any period, and

transient solutions.

11
2.4.1. Steady-state Solutions and Asympotes of Period Higher Than One:

Table 2.1 shows the possible stable asymptotic solution behavior between DEs (ODEs

or PDEs) and their discretized counterparts. Some of the phenomena will be supported

by simple examples in section III. The main connection between the DEs and their

discretized counterparts is that steady-state solutions of the continuum are solutions

of the discretized counterparts but not the reverse. Their main difference is that new

phenomena are introduced by the numerical methods in the form of spurious stable

and unstable asymptotic solutions of any period. In the past, the phenomena of spurious asymptotes were observed largely beyond the linearized stabihty of the schemes.

Some numerical analysts and applied computational scientists were not alarmed and

were skeptical about these phenomena since, theoretically, one is always guided by the

hnearized stabihty limit of the scheme. However, this resaoning is only valid if one is

solving a scalar nonlinear ODE and the initial data are known. Another important

concept is that the result of operating with time steps beyond the hnearized stabihty

limit is not always a divergent solution; spurious steady-state solutions and spurious

asymptotes of higher period can occur.

Our current study also indicated that depending on the form of the nonlinear DEs, all

ODE solvers can introduce spurious asymptotic solutions of some period or all periods.

However, the most striking result is that for certain schemes and depending on the form

of the nonlinear DEs, spurious steady states can occur below the linearized stability

limit. See section III and our companion paper [21] for more details.

Another important factor is that associated with the same (common) steady-state

solution, the basin of attraction (domain of attraction) of the continuum might be

vastly different from the discretized counterparts. This is due entirely to the separate

dependence and sensitivity on initial and boundary conditions for the individual system.

The situation is compounded by the existence of spurious steady states and asymptotes

of period higher than one and possibly chaotic attractors.

Here the basin of attraction of a dynamical system is the domain for which the set of

initial conditions time asymptotically approaches a specific asymptote. Figures 2.1 and

2.2 show the basins of attraction of two popular ODE dynamical systems. Figure 2.1

shows the multiple stable steady states and their basins of attraction for the damped

pendulum equation

du

dt - v, (2.4a)

de

-- = -_v - sin(u) (2Ab) de

for _ = 0A. Figure 2.2 shows the multiple steady states and their basins of attraction

for the simple predator-prey equation

du

-- = -3u + 4u 2 - uv/2 - u a (2.5a) de

12
dv

= -2.1v + uv. (2.5b)

d_

where u is the population of the prey and v is the population of the predator. These

figures are taken from Parker and Chua [8] and were generated by the use of a variable

time step Runge-Kutta-Fehlberg method with built in accuracy check (if the numerical

solutions are approximating the true solution of the ODE). See reference [8] for details.

These figures, although generated numerically, with the built in accuracy check the fixed

points and basins of attraction coincide with the ODEs. The stable fixed points of the

damped pendulum equation are 2nr, n = 0, 1, .... The unstable fixed points (saddles)

are (2n + 1)7r. The separatrices of the saddle points divide the phase plane into the

different basins of attraction for the corresponding stable fixed points. The fixed points

of the predator and prey equation are slightly less regular than the damped pendulum

equation. Figure 2.2 shows two saddle points at u = 1, v = 0 and u = 3, v = 0, one

stable focal point at u = 2.1, v = 2 and one stable nodal point at u = 0, v = 0. Again

the separatrices of the saddle points divide the phase plane into the basins of attraction

for the corresponding stable fixed points.

Intuitively, in the presence of spurious asymptotes, the basin of the true steady states

(steady states of the DEs) can be separated by the basins of attraction of the spurious

asymptotes and interwoven by unstable asymptotes, whether due to the physics (i.e.,

present in both the DEs and the dlscretized counterparts) or spurious in nature (i.e.,

introduced by the numerical methods).

For PDEs, another added dimension is that even with the same time discretization

but different spatial discretizations or vice versa, the basins of attraction can also be

extremely different. However, mapping out the basins of attraction for any nonlinear

continuum dynamical system other than the very simple scalar equations relies on numerical methods. The type of nonlinear behavior and the dependence and sensitivity

to initial conditions for both the PDEs and their discretized counterparts make the

understanding of the true physics extremely difficult when numerical methods are the

sole source. Under this situation, how can one delineate the numerical solutions that

approximate the true physics from the numerical solutions that are spurious in nature?

Hopefully, with our simple illustrations in section III, we can demonstrate the importance of the current subject and, most of all, stress the importance of knowing the

general dynamical behavior of asymptotes of the schemes for genuinely nonlinear scalar

DEs before applying these schemes in practical calculations.

2.4.2. Transienl or Time-Accurate Solutions:

It is a common misconception that inaccuracy in long time behavior poses no consequences on transient or time-accuate solutions. This is not the case when one is dealing

with genuinely noIrlinear DEs. For genuinely nonhnear problems, due to the possible

existence of spurious solutions, larger numerical errors can be introduced by the nu13
merical methods than one can expect from local linearized analysis or weakly nonlinear

behavior. The situation will get more intensified if the initial data of the DE is in the

basin of attraction of a chaotic transient [63-65] of the discretized counterpart. This is

due to the fact that existence of spurious asymptotes transact wrong behavior in finite

time. In fact, it is possible the whole solution trajectory is likely to be erroneous.

We'd like to end this section with a direct quote from Sanz-Serna and VadiUo's paper

[55]. This quote indicates the danger of relying on linearized stability and convergence

theory in analyzing nonlinear dynamical problems. Reference [55] is one of the few papers trying to convey to numerical analysts the flavor of the powerful "nonlinear dynamic

approach". Hopefully, with the current discussion, we can convey to computational fluid

dynamicists the flavor of the importance of the "nonlinear dynamic approach" in CFD

analysis.

"Assume that the convergence of a numerical method has been established; it is

still possible that for a given choice of At, or even for any such a choice, the qualitative behaviour of the numerical sequence u °, u I , ..., u", ... be competely different from that of the theoretical sequence u(to), u(tl), ... u(t,_), ... This discrepancy which refers to n tending to 0% At fixed cannot be ruled out by the convergence requirement, as this involves a different limit process (namely At tending to

0).

The fact that analyses based on linearization cannot accurately predict the qualitative behaviour of u" for fixed At. should not be surprising:

there is a host of nonlinear phenomena (chaos, bifurcations, limit cycles ...) which

cannot possibly be mimicked by a linear model."

14
III. THE ODE CONNECTION

In this section, we review some of the fundamentals and available theory and discuss

our major results. The discussion will have some overlap with our companion paper

[21].

8.1. Preliminaries

Consider an autonomous nonhnear ODE of the form

du

d--_= aS(u), (3.1)

where a is a parameter and S(u) is nonlinear in u. For simplicity of discussion, we

consider only autonomous ODEs where (_ is hnear in (3.1); i.e., a does not appear

explicitly in S.

A fixed point u" of an autonomous system (3.1) is a constant solution of (3.1); that

is

S(u') = 0. (3.2)

Note that the terms "equihbrium points", "critical points", "stationary points",

"asymptotic solutions" (exclude periodic solutions for the current definition), "steadystate solutions" and "fixed points" are sometimes used with slightly different meanings

in the literature, e.g., in bifurcation theory. For the current discussion and for the majority of nonlinear dynamic hterature, these terms are used interchangeably. We might

want to mention that certain researchers reserve the term "fixed point" for discrete

maps only.

Consider a nonlinear discrete map from finite discretization of (3.1)

u n+l = u '_ + D(u'_,r), (3.3)

where r = aAt and D(u n, r) is linear or nonlinear in r depending on the ODE solvers.

Here the analysis is similar if D is a nonlinear function of u n+p, p = 0,1, ..., m. Examples

to illustrate the dependence on the numerical schemes for cases where D is hnear or

nonlinear in the parameter space will be given in the subsequent section.

A fixed point u" of (3.3) (or fixed point of period 1) is defined by u '_+3 = u '_, or

u* = u* + D(u',r) (3.4a)

or

D(u',r) =: 0. (3.4b)

15
One canalso defineafixed point of periodp, where p is a positive integer by requiring

that tt n+p = u n or

u" = EP(u*,r) but u" # Ek(u',r) for 0 < k < p. (3.5)

Here, EP(u',r) means that we apply the difference operator E p times, where

E(un, r) = u n + D(un,r). For example, a fixed point of period 2 means u n+2 = u n or

= (3.6)

In this context, when dealing with discrete systems, the term "fixed point" without

indicating the period means "fixed point of period 1" or the steady-state solution of

(3.3).

In order to illustrate the basic idea, the simplest form of the Ricatti ODE, i.e., the

logistic ODE (2.3) with

S(u)=u(l-u) (3.7)

is considered. For this ODE, the exact solution is

It 0

u(t) = u0 + (1 - u°)e -_''' (3.8)

where u ° is the initial condition. The fixed points of the logistic equation are roots of

u'(1 - u °) = 0; it has two fixed points u" = 1 and u ° = 0.

To study the stability of these fixed points, we perturb the fixed point with a disturbance _, and obtain the perturbed equation

/:d

= +

d_

Next, S(u" + _) can be expanded in a Taylor series around u', so that

(3.9)

d-_= a [S(u') + S,,(u')_+ -_s,,,,( ,," )_ +... ] , (3.10)

where S,,(u °) =-_sl=" . Stability can be detected by examining a small neighborhood

of the fixed point provided if for given a, u ° is not a hyperbolic point [3,7,9] (i.e., if the

real part of aS,,(u °) # 0). Under this condition _ can be assumed small, its successive

powers _,_3, ... can normally be neglected and the following hnear perturbed equation

is obtained

d--_= aS=(u')(. (3.11)

The fixed point u ° is asymptotically stable if aS=(u °) < 0 whereas u" is unstable if

aS,,(u') > 0. If aS,,(u') = 0, a higher order perturbation is necessary.

16
If weperturb the logistic equation around the fixed point with a > 0, one can find

that u* = I is stable and u" = 0 is unstable. It is well known that the general asymptotic

solution behavior of the logistic ODE is that for any u ° > 0, the solution will eventually

tend to u" = 1. Figure 3.1 shows the solution behavior of the logistic ODE.

Now, let us look at three of the well known ODE solvers. These are exphcit Euler

(Euler, forward Euler), leapfrog and Adam-Bashforth. For the ODE (3.1) with S(u) =

u(1 -u), the dynamical behavior of their corresponding discrete maps is well estabhshed.

The explicit Euler is given by

un+l =u"+,.S(u=), (3.12)

and it is after a linear transformation, the well known logistic map [26-30]. The leapfrog

scheme can be written as

u,,+l : u--1 + 2rS(un),

and it is a form of the H_non map [32]. The Adam-Bashforth method given by

=-+, = =- + ,[3s(=") _ ,

(3.13)

(3.14)

is again a variant of the H_non map and has been discussed by Priiffer [44] in detail.

We can determine fixed points of the discrete maps (3.12)-(3.14) and their stability

properties in a similar manner as for the ODE. It turns out that all three of the discrete

maps have the same fixed points as the ODE (3.1) -- a desired property which is important for obtaining asymptotes of nonlinear DE numerically. Here we use asymptotes

to mean fixed points of any period.

The corresponding linear perturbed equation for the discrete map (3.3), found by

substituting _'_ = u" + _'_ in (3.3) and ignoring terms higher than _'_ is

+ At)]. (3.15)

Here the parameter a of the ODE has been absorbed in the parameter At due to the

assumption that a does not appear explicitly in S(u). For stability we require

l1 + AtD,,(u',;kt) < 1. (3.16)

Again, for [1 + AtD_(u*, At)[ = 1, higher order perturbation is necessary. For a fixed

point of period p the corresponding hnear perturbed equation and stabihty criterion are

C +p= CE (u',At). (3.17)

and

p

IE,, (u , At) l < 1, (3.18a)

17
with

At)=gd E(,"+,-' , ...

(3.1Sb)

For S(u) = u(1 - u), the stability of the stable fixed points of period 1 and 2 for

discrete maps (3.12)-(3.14) with r = aAt are

Explicit Euler:

Leapfrog:

stable if 0 < r < 2

stable if 2 < r < x/6.

u" = 1 unstable for all r > 0

chaotic solution exist for all r no matter how small

Adam-Bashforth:

u*=l stable if0<r< 1

period 2 stable if 1 < r < v_.

Figure 3.2 shows the stable fixed point diagram of period 1,2,4, 8 by solving numerically the roots of (3.12) for S(u) = u(1 - u). The r axis is divided into 1,000 equal

intervals. The numeric labelling of the branches denotes their period. The subscript E

on the period 1 branch indicates the stable fixed point of the DE.

Two of these three examples serve to illustrate that the result of operating with a

time step beyond the linearized stability limit of the stable fixed points of the nonlinear

ODEs is not always a divergent solution; spurious asymptotes of higher period can

occur. This is in contrast to the ODE solution, where only a single stable asymptotic

value u* = 1 exists for any a > 0 and any initial data u ° > 0. It is emphasized here that

these spurious asymptotes, regardless of the period, stable or unstable, are solutions in

their own right of the discrete maps resulting from a finite discretization of the ODE.

3.2. Spurious Steady-State Numerical Solutions

For the previous three ODE solvers, we purposely picked the type of schemes that do

not exhibit spurious fixed points [56] but allow spurious fixed point of period higher than

1. In this section, we discuss the existence of spurious steady-state numerical solutions.

Again, it is emphasized here that these spurious steady states, stable or unstable, are

solutions in their own right of the resulting discrete maps. Consider two second-order

Runge-Kutta schemes, namely, the modified Euler (R-K 2) and the improved Euler (RK 2), the fourth-order Runge-Kutta method (R-K 4), and the second and third-order

predictor-corrector method [66-68] of the forms

18
Modified Euler (R-K 2) method:

U n+l -----U n JI- ?'S (U n

Improved Euler (R-K 2) method:

R-K 4 method:

_n+l = U n .+

(3.18)

(3.19)

u '_+1 =u '_+_ kl+..k2+2ka+k4

k I _ S n

1 '

k2 = S u n+ :rkl/_

2 /

/

1

(3.21)

Predictor-corrector method of order m:

u (°) = u n + rS _

u4k+l) = u_ + _ S_

u'_+a = u" + _ S"

+ S(k)],

+ S('_-1)].

k = 0, 1, ...,rn- 1

(3.22)

Using the same procedures, one can obtain the fixed points for each of the above

schemes (3.18) - (3.22). Figures 3.3 - 3.7 show the stable fixed point diagrams of period

1,2,4 and 8 for these five schemes for S(u) = u(1 - u). Some of the fixed points of lower

period were obtained by closed form analytic solution and/or by a symbolic manipulator

such as MAPLE [69] to check against the computed fixed point. The majority are

computed numerically [2,8]. The stability of these fixed points was examined by checking

the discretized form of the appropriate stability conditions. Again the axis is divided

into 1,000 equal intervals. The numeric labelling of the branches denotes their period,

although some labels for period 4 and 8 are omitted due to the size of the labelling areas.

The subscript E on the main period one branch indicates the stable fixed point of the

DE while the subscript S indicates the spurious fixed points introduced by the numerical

scheme. Spurious fixed points of period higher than one are obvious and are not labeled

19
except for special cases.Note that these diagrams, which appear in most parts as solid

lines are actually points, which are only apparent in areas with high gradients.

To contrast the results, similar stable fixed point diagrams are also computed for

S(u) = u(1 - u)(b - u), 0 < b < 1. See figures 3.8 - 3.14. The stable fixed point for

the ODE in this case is u" = b and the unstable ones are u ° = 0 and u* = 1. For any

0 < u ° < 1 and any a > 0, the solution will asymptotically approach the only stable

asymptote of the ODE u* = b.

Note that contrary to the DE, the maximum number of stable and unstable fixed

points (real and complex) for each scheme varied between 4 to 16 for S(u) = u(1 - u)

and 9 to 81 for S(u) = u(1 - u)(b - u), depending on the numerical methods and the r

value. The domains of all of the fixed point diagrams are chosen so that they cover the

most interesting part of the scheme and ODE combinations. Notice that asymptotes

might occur in other parts of the domain as well.

Aside from the striking difference in topography in the stable fixed point diagrams

of the various methods and ODE combinations, all of these diagrams have one similar

feature; i.e., they all exhibit spurious stable fixed points, as well as spurious stable

fixed points of period higher than one. Although in the majority of cases, these occur

for values of r above the hnearized stability bruit, this not always the case, as in the

modified Euler scheme applied to the logistic ODE and du/dt = ctu(1 -u)(b- u),

0 < b < .5, and the R-K 4 apphed to the logistic DE. For these two methods and ODE

combinations, stable spurious fixed points occur below the hnearzed stability limit. In

some of the instances, these spurious fixed points are outside the interval of the stable

and unstable fixed points of the ODEs. Others not only he below the hnearized stability

limit but also in the region between the fixed points of the DEs and so could be very

easily achieved in practice.

One might argue that for the ODEs that we are considering, it is trivial to check

whether an asymptote is spurious or not. For example, if _ is a spurious asymptote of

period one, then S(_) # 0. The main purpose of the current illustration is to set the

baseline dynamical behavior of the scheme so that one can use it wisely in other more

complicated settings such as when nonlinear PDEs are encountered in which the exact

solutions are not known. Under this situation, spurious asymptotes could be computed

and mistaken for the correct steady-state solutions.

Note that fgr the modified Euler method, spurious fixed points of higher periods and

chaotic attractors as well as spurious steady states occur below the linearized stability

bruit. Let Q be the basin of attraction of the fixed point of the ODE and let r* be

the corresponding linearized stability limit value of the scheme. Then there exists a

portion of the basin ft denoted by N c in which Qc C Q and an interval of r with

0 < T < r" which actually belongs to the basin of attraction of the chaotic attractor of

the discretized counterparts. There also exist some other tip C Q and an interval of r

with 0 < r < r ° and p _> 1 an integer, which actually belongs to the basin of attraction

of a stable asymptote of period p of the corresponding discrete map. This leads to the

2O
issueof the dependence of solutions on initial data which will be a subject of the next

subsection.

3.3. Strong Dependence of Solutions on Initial Data

For simple nonlinear ODEs that we are considering, the fixed point diagram is extremely useful for the understanding of the dynamics of the DEs and their discretized

counterparts. However, when fixed points of higher periods and/or complex nonlinear

equations are sought, searching for the roots and testing for stability of highly complicated nonlinear algebraic equations can be expensive and might lead to inaccuracy.

Equally useful for understanding the dynamics are the bifurcation diagram and basin

of attraction of fixed points for both the DEs and the difference schemes. The bifurcation

diagram for the one-dimensional discrete maps displays the iterated solution u n vs. r

after iterating the discrete map for a given number of iterations with a chosen initial

condition (or multiple initial conditions) for each of the r parameter values.

Bifurcation is broadly used to describe significant qualitative changes that occur in the

orbit structure of a dynamical system as the system parameters are varied. In general,

bifurcation theory can be divided into two general classes, namely, Iocal and global.

Local bifurcation theory is concerned with the bifurcation of fixed points of nonlinear

equations and discrete maps. Global bifurcation studies phenomena away from the fixed

points. It studies the interaction between different types of fixed points. One might

define a bifurcation point as being any dynamical system which is structually unstable

[3,8,9]. A fixed point is structurally stable if nearby solutions have qualitatively the

same dynamics. The linearized stability limit of a fixed point of a scheme is the same as

the bifurcation point in the corresponding bifurcation diagram of the resulting discrete

map.

For the numerical computations of the bifurcation diagrams with a given interval of r

and a chosen initial condition (or multiple initial conditions), the r axis is divided into

500 equal spaces. In each of the computations, the discrete maps were iterated with

600 preiterations and the next 200 iterations were plotted for each of the 500 r values.

The domains of the r and u '_ axes are chosen to coincide with the stable fixed point

diagrams shown previously. For our current interest, it is not necessary to distinguish

the difference between a stable fixed point of period 200 and a chaotic attractor.

Figure 3.15 shows the bifurcation diagram of the Euler scheme applied to the logistic

DE with an initial condition u ° > 0. It is of interest to know that in this case the

bifurcation diagram looks practically the same for any _0 > 0. This is due to the fact

that no spurious fixed points or spurious asymptotes of low periods exist for r < 2.627.

Comparing the bifurcation diagram with figure 3.2, one can see that if we computed all

of the fixed points of period up to 200 for figure 3.2, the resulting fixed point diagram

would look the same as the corresponding bifurcation diagram (assuming 800 iterations

21
of the logistic map are sufficient to obtain the converged stable asymptotes of period

upto 200 and a proper set of initial data are chosen to cover the basins of all of the

periods in question). The numeric labelling of the branches in the bifurcation diagram

denote their period, with only the essential ones labelled for identification purposes.

In order to interpret the bifurcation diagram for other ODE and scheme combinations, some knowledge of the fixed point diagram is necessary, at least for the lower

order periods. Otherwise, one cannot identify the exact periodicity of the asymptotes

easily. As can be seen later, a "full" bifurcation diagram cannot be obtained efficiently

without the aid of the stable and unstable fixed point diagram for schemes that exhibit

spurious fixed points of any period, especially lower periods. In most cases, the unstable asymptotes divide the domain into the proper basins of attraction for the stable

asymptotes (spurious or otherwise), and at least one initial data point is used from each

of the basins of attraction before a full bifurcation diagram can be obtained.

In all of the fixed point diagrams 3.3 - 3.14, the bifurcation phenomena can be divided

into three kinds. For the first kind, the paths (spurious or otherwise) resemble period

doubling bifurcations (flip bifurcation) [2-5] similar to the logistic map. See figures 3.2,

3.6 and 3.8 for examples. The second kind occurs, most often, at the main branch 1E,

with the spurious paths branching from the correct fixed point as it reaches the linearized

stability bruit, and quite often even bifurcating more than once (pitchfork bifurcation

or supercritical bifurcation [70,7]), as r increases still further before the onset of period

doubling bifurcations. See figures 3.4, 3.7, 3.9 - 3.11 and 3.13 for examples. The

third kind again occurs most often at the main branch 1E. The spurious paths near the

hnearized stabifity limit of IE would experience a transcritical bifurcation [3,7,9,70]. See

figures 3.3, 3.5, 3.7 and 3.14 for examples. Notice that the occurrence of transcritical

and supercritical bifurcations are not limited to the main branch 1E. See figures 3.11

- 3.14 for examples. The other commonly occurring bifurcation phenomenon is the

subcritical bifurcation which was not observed in our two chosen S(u) functions. With

a shght change in the form of our cubic function S(u), a subcritical bifurcation can be

achieved [70,3,7,9]. The consequence of the latter three bifurcation behaviors is that

bifurcation diagrams calculated from a single initial condition tt ° will appear to have

missing sections of spurious branches, or even seem to jump between branches. This is

entirely due to the existence of spurious asymptotes of some period or more than one

period, and its dependence on the initial data. This occurs even for the Euler scheme

as depicted in figure 3.8. See section 3.4 for further discussion of these four types of

bifurcation phenomena.

Figures 3.16- 3.18 show the bifurcation diagram by the modified Euler method for the

logistic ODE with three different starting initial conditions. In contrast to the exphct

Euler method, none of these diagrams look alike. One can see the influence and the

strong dependence of the asymptotic solutions on the initial data. Figure 3.19 shows

the corresponding "full" bifurcation diagram, their earlier stages resembling the fixed

point diagram 3.3. Figures 3.20 - 3.22 illustrate similar bifurcation behavior for the

corresponding R-K 4 method. Figure 3.12 serves as an example to illustrate that the

effect of overplotting a number of initial data, but not the appropriate ones, would not

22
be sufficient to coverall of the essentialspuriousbranches. Figures 3.23- 3.25show a

similar illustration for S(u) = u(1 - u)(b - u), 0 < b < .5 by the improved Euler, R-K

4 and the modified Euler method. The strong dependence of solutions on initial data

is evident from the various examples in which this type of behavior is very common for

genuinely nonlinear problems.

In order to compute a "full" bifurcation diagram, we must overplot a number of

diagrams obtained by the guide of the stable and unstable fixed point diagram as an

appropriate set of starting initial data. In the case where the fixed point diagrams are

extremely difficult to compute, a brute force method of simply dividing the domain of

interest of the _t'_ axis into equal increments and using these u '_ values as initial data

is employed. The "full" bifurcation diagram is obtained by simply overplotting all of

these individual diagrams on one.

For completeness, figures 3.26 - 3.38 show the "full" bifurcation diagrams for the

corresponding fixed point diagrams shown previously. Figures 3.36 and 3.37 show a

blow up section of figures 3.34 and 3.35. Notice that the exact values of the initial

data are immaterial as long as these values cover all of the basins of attraction of the

essential lower order periods (i.e., at least one initial data point is used from each of

the basins). Here, we use the term "full" bifurcation diagram to mean just that. No

attempt has been made to compute the true full bifurcation diagram since this is very

costly and involves a complete picture of the basins of attraction for the domain of

interest in question.

3.4. Classification of ODE solvers

(According to the Ezistence of Spurious Fized Points)

In reference [56], lserles studied the stability of ODE solvers for nonlinear autonomous

ODE via the dynamical approach. He proved that linear multistep methods (LMM) [66-

68] that give bounded values at infinity always produce correct asymptotic behavior, but

it is not the case with Runge-Kutta methods and some predictor-corrector methods. He

demonstrated that the Runge-Kutta and predictor-corrector methods may lead to false

asymptotes. However, he did not discuss the possibihty of these spurious asymptotes

existing below the linearized stability limit.

For implicit LMM, he assumed the resulting nonlinear algebraic equations are solved

exactly. He also showed the influence of nonlinear algebraic solvers on the size of stability

regions for implicit LMM. His conclusion wa,; that the standard nonlinear algebraic

solver -- the modified Newton-Raphson method

Ik+l) Ik) ,_+i - u,_ _ (3.23)

u'_+l = u'_+l - 1 - _S,,(u,_) '

can drastically degrade the region of stability lbnit as compared to the Newton-Raphson

method

23
_ (k)

Un+ 1 -- Un+ 1

1- )

On the other hand, the direct iteration method

(3.24)

"1

u.+, (k+l) = ,,,, + _r S(,,n)+ S(,,.+,)(k) / (3.25)

J

converges only if the step sizeis of the same order of magnitude as that required for an

explicit method. Thus the advantage of using an implicit method to enhance stability

is lost. Here for clarity of notation, when iteration procedures are involved, u,_ is used

in place of u n of the previous section.

The implications of behavior detailed in Iserles' work [56] range far beyond pure ODE.

For most CFD application, the use of implicit time discretization to "time" march the

solution to steady state is very common. The resulting nonlinear algebraic systems

are solved by either noniterative linearization [71,14] or by some kind of iterative or

relaxation procedures. Very often, applied computational fluid dynamicists experience

a non-convergent solution where the residual will decrease only so far before reaching

a plateau with a time step larger than the explicit method. Therefore the behavior

observed in Iserles' work could explain the degradation in the stability of the implicit

scheme in practice. Indeed, even though the mechanisms involved are far more complicated than those studied here, elements such as spatial discretization dynamical behavior

and nonlinear coupling effect for systems, could well be an explanation.

More recently, ]series and Sanz-Serna [57] established conditions for using a variable

step size analysis to avoid spurious fixed points in a class of Runge-Kutta methods.

Looking at the problem from another perspective, it is very useful to find the cause

of the existence of spurious asymptotes by looking at the form and properties of the

resulting discrete maps, regardless of the methods. We have the following two observations.

(1) Assume that the only parameter that was introduced by a numerical method

is At. Then from Iserles results and our current investigation, one obvious necessary

condition for the existence of spurious steady states of ODE solvers for (3.1) is the introduction of nonlinearity in the parameter space At. This is evident from our examples

and general analysis. For example, if At (or r) is linear in (3.3), then (3.3) can be

written as

u n+l = a n + crS(un), c a constant of the scheme. (3.26)

Therefore any fixed point of (3.3) is a fixed point of (3.1). Without lost of generality, a

similar proof applies to the resulting difference operator D from a p time level scheme.

(2) The second observation is that one can classify the types of spurious steady

state in the form of bifurcation theory near a bifurcation point or a bifurcation limit

24
point. Figures 3.39 and 3.40 show the definition of the various types of branching

points and the stability of solution in the neighborhood of branch points. In other

words, the classification is according to the onset of spurious asymptotes of subcritical,

supercritical or transcritical bifurcations. See figure 3.41 for the definition of the three

types of phenomena.

Assume an ODE solver introduces nonlinearity in the parameter space At for (3.1).

Then a necessary and sufficient condition for the occurrence of spurious steady states

below the linearized stabihty limit on the main branch 1E (stable fixed points of the

DE) is that a transcritical or subcritical bifurcation of the types shown in figures 3.42

and 3.43 exist at the bifurcation point or near st bifurcation hmit point. It is emphasized

here that the existence of spurious fixed points of higher period can be independent of

the existence of spurious steady states (fixed points of period 1).

A detailed analytical analysis on the existence of transcritical, subcritical and supercritical bifurcations for the class of Runge-Kutta methods can be found in our companion paper [21]. Figures 3.44 - 3.54 illustrate the onset of different types of spurious

steady states by showing the stable and unstable fixed points of periods 1 and 2, and

the types of bifurcation phenomena for the modified Euler, Improved Euler and RK 4 and the predictor-corrector schemes of order 2 and 3 for S(u) = u(1 - u) and

S(u) = u(1 - u)(b - u), 0 < b < .5. In order to illustrate the different behavior in

an uncluster fashion, not all of the periods 1, 2 and branching points are labeled. It is

interesting to see the manner in which the onset of the different types of bifurcations

occur, in particular, the birth of the different types of bifurcations away from the 1E

branches.

3.5. Basins of Attraction

Due to the separate dependence and sensitivity on initial data for the individual

DEs and the discretized counterparts, in conjunction with the existence of spurious

steady states and asymptotes of higher periods, even associated with the same (common)

steady-state solution, the basin of attraction of the continuum might be vastly different

from the discretized counterparts.

Take for example, S(u) = u(1 - u). The only stable fixed point of the logistic ODE

is = = 1. The entire domain of the real u'_-axis is divided into two basins of attraction

for the ODE independent of any positive a. Now if one numerically integrates the ODE

by the modified Euler method, extra stable and unstable fixed points can be introduced

by the scheme depending on the value of r. That is for certain ranges of the 7" values,

the u'_-axis is divided into four basins of attraction. But of course for other ranges of r,

higher period spurious numerical solutions exist, more basins of attraction are created

within the same u'_-axis range, etc. Stable and unstable fixed point diagrams such as

figures 3.44 - 3.54 are very useful in the division of the u'_-axis into different basins of

lower periods.

25
3.6. Systems of ODEs

As can be seen from the previous sections, the rich and complicated dynamical behavior of discrete maps resulting from finite discretization of simple nonlinear scalar

autonomous ODEs is very enlightening, educational and useful in giving some indications of the strange behavior encountered in practice. One would naturally ask how

highly coupled nonlinear first-order autonomous systems complicate the issue. After

all, these types of systems occur naturally in physical science and engineering fields.

Examples are

(1) second or higher order nonlinear scalar autonomous or nonautonomous ODEs

arising from mechanical systems,

(2) meteorology,

(3) chemical reaction equations arising from chemistry,

(4) system of ODEs arising from the method of fines approach in reaction-diffusion,

reaction-convection and reaction-convection-diffusion equations.

Future work will be directed towards investigation into the nonlinear dynamical effect

of using ODE solvers for nonlinear system of ODEs. Here, we do not attempt to give

a detailed discussion on this subject, but rather indicate some of the implications from

our experience as well as what is availiable in the literature.

First, the coupling of first-order nonlinear systems arising from a higher-order scalar

nonlinear ODE is very different from the truly nonlinear coupling on systems of firstorder ODEs. This difference carries over to their discretized couterparts. Second, due

to the nonlinear coupling effect, whatever is observed in the nonlinear scalar case will

definitely exists in the coupled system case in a more complex manner. Even with the

help of the center manifold theorem [2-5], nonlinear systems of higher than three firstorder ODEs are still extremely difficult to analyze. One major factor in analyzing the

associated discrete maps from finite discretization of the continuum is that when three

or more time levels of ODE solvers are used, even though the continuum is a first-order

scalar autonomous ODE, the resulting discrete maps are (p - 1)th-order, where p is the

time level. One can extrapolate the complexity involved if nonhnear coupled systems of

higher-order ODEs were discretized by p-time levels of ODEs solvers. Some aspects and

implications of numerical integration of second- and third-order ODEs are discussed in

references [39,40,72]. Some of our preliminary numerical experiments agree with the

above general conclusion.

3.7. Suitability to the Type of Computational Environment

The main approach that we use in this paper is to establish the necessary mathematical reasoning and then to support this reasoning with extensive numerical experiments.

26
Our current study on one-dimensional scalar nonlinear dynamical equations which consist of a single parameter indicates that the understanding of the nonlinear effects encountered when applying finite-difference schemes to nonhnear differential equations is

greatly aided by the analysis of bifurcation diagrams which record the values of successive iterations for a range of parameters. Equally useful are diagrams showing the basins

of attraction of equilibria, both those of the differential equations and the spurious attractors generated by the difference scheme. The generation of such diagrams, however,

is computationally expensive, especially for the basins of attraction where each point

on the diagram represents a different choice of parameters for which many iterations of

the scheme must be performed to determine its significance.

In all the bifurcation diagrams, the computations were performed on the VMS VAX

in double precision. Take for example, figure 3.34. Each dot on the plot represents a

solution obtained by integrating the discretized equation 800 times with each of the 20

prescribed initial data and each of the 500 equally space values of r = aAt. In other

words, we are integrating the same equation for 10,00(l different values of r and initial

data combinations, and also iterating the same equation for each of these combinations

with 800 iterations. The task can therefore be greatly enhanced by parallel computation,

since essentially the same process needs to be applied to each point in a fine two- or

three-dimensional array, each element representing a pixel on a high resolution screen or

plotter. It is therefore a task highly suited to machines, such as the Connection machine,

which have large numbers of processors enabling the entire region or subregions of the

problem to be analysed in one pass rather than in a sequential point-by-point approach.

The intensity of (repetitive) computing involw_d is too great to gain major benefit from

machines such as the CRAY.

For multidimensional systems consisting of several parameters, we envision that the

intensity of repetitive computing to obtain a bifurcation diagram or a basin of attraction

cannot be realized if it is not performed on a massively parallel computer such as the

Connection machine.

27
IV. LEVEL OF COMPLEXITY FOR PDEs

In order to systematically approach the subject of studying spurious steady-state

numerical solutions of nonlinear nonhomogeneous hyperbolic and parabolic PDEs via

the nonlinear dynamic approach, we propose to pursue the subject in three stages. First,

we will attempt to obtain a full understanding of the subject for time discretization of

ODEs. The investigation can give insight into numerical methods employing the Strang

type of operator splittings or methods of lines approach for nonhomogeneous hyperbolic

and parabolic PDEs. The second stage will involve the study of the discrete travelling

wave solutions of the reaction-convection and reaction-convection-diffusion equations.

The third stage will involve the study of the complete temporal-spatial discretizations

of the reaction-convection and reaction-convection-diffusion equations. The last stage of

the proposed plan is extremely difficult to analyze. Some aspects of full discretizations

and discrete travelling wave solutions were investigated by [46-54, 58-62,73,74,10].

The question now is in what specific area will this approach advance the state-of-theart in CFD. Our preliminary study indicated that many existing results for nonlinear

dynamical systems such as chaos, bifurcations, and limit cycles (closed periodic orbits [5]) have a direct application to problems containing nonlinear source terms such

as the reaction-diffusion, reaction-convection or the reaction-convection-diffusion equations. Also they have a direct application to most of the nonlinear shock-capturing

methods such as the total variation diminishing (TVD) schemes [14,75-78]. With the

advent of increasing demand for numerical accuracy, stability, efficiency, and uniqueness

of numerical solutions in modeling such equations, an interdisciplinary approach for the

analysis of these systems and schemes is needed. Besides it is a common practice in CFD

to employ a time-dependent approach to achieve steady state. The separate dependence

of solutions on initial data and system parameters for the individual PDE and its finitedifference equations is the crucial element in determining how well a numerical solution

can mimic the true physics of the problem.

The following is an attempt to give a flavor of the subject and at the same time provide

a justification for the importance of this subject area in CFD algorithm development

for our next generation aerodynamics needs.

4.1. Model Equations

One of the recent areas of emphasis in CFD has been the development of appropriate finite-difference methods for nonequilibrium gas dynamics in the hypersonic range

[14,78-81]. A nonlinear scalar reaction-diffusion model equation would be of the form

81/ 821/

8--t = ¢_z 2 + aS(u), ¢,a system parameters, (4.1)

a nonlinear scalar reaction-convection model equation would be of the form

28
0-7+ - (4.2)

and a nonlinear scalar reaction-convection-diffusion model equation would be of the

form

Ou Of(u) 02u

o-7+ - + (4.3)

Here f(u) a linear or nonlinear function of u. The nonlinear source term (or the reaction

term) S(u) can be very stiff. Note that phenomena such as chaos, bifurcations and limit

cycles only relate to source terms S(u) which are nonlinear in u. Equation (4.3) can be

viewed as a model equation in combustion or as one of the species continuity equations in

nonequilibrium flows (except in this case, the source term is coupled with other species

mass fractions).

The above mode] equations are good starting points in the investigation of correlation

between the theory of chaotic dynamical systems and uniqueness, stability, accuracy and

convergence rate of finite-difference methods for CFD.

4.2. Level of Complexity

The main interest is to investigate what types of new phenomena arise from the

numerical methods that are not present in the original nonlinear PDE, as a function

of the stiff coefficient a, the diffusion coefficient e, and the time step At with a fixed

(or variable) grid spacing As. The time step can vary greatly depending on whether

the time discretization is explicit or implicit. More precisely, one wants to weed out all

undesirable phenomena due to the numerical method (e.g., additional equilibrium points

introduced by the time as well as spatial discretizations, degradation of the domain of

attraction, etc.) and to identify whether the numerical method really describes the true

solution of the PDE under prescribed initial and boundary conditions with a, e, the

time step At and the grid spacing Az being parameters. The study can be divided into

steady and unsteady behavior with or without shock waves.

The major stumbling block is that combustion-related and high speed hypersonic flow

problems usually contain multiple equilibrium states and shock waves that are inherent

in the governing equations. Furthermore, spurious equilibrium states can be introduced

by the time differencing and/or the spatial differencing. In many instances the stable

and unstable equilibrium states, whether due to the physics or spurious in nature, are

interwoven over the domain of interest and are usually very sensitive to the initial

conditions and the time steps (even when the chosen time step is within the hnearized

stability limit as indicated in our study) as well as variation of parameters such as angle

of attack, Reynolds number and coefficients of physical and numerical dissipations and

physical and numerical boundary conditions.

29
The sensitivity of numerical solutions to coefficients of physical and numerical dissipations is evident from the study of Mitchell and Bruch on the reaction-diffusion equation.

Their main result is that diffusion, which is usually perceived as having a stabihzing

effect, is able to produce chaotic as well as divergent numerical solutions. Another interesting result due to Mitchell and Bruch was the production of chaos by decreasing

the space increment or increasing the time increment. They showed that the addition of

diffusion poses severe problems unless waves of constant speed c are assumed, in which

case it reverts to an ODE with z + c/, as the independent variable. The sensitivity

of numerical solutions to numerical boundary condition procedures was discussed in

[82,83].

On the subject of sensitivity and dependence of solutions on initial data, the basin

of attraction might be very different between the PDE and the discretized counterpart.

The basin of attraction might contract or be very different from the basin of attraction

for the original PDEs depending on the numerical methods. In many instances, even

with the same spatial discretization but different time discretizations, the basins of

attraction can also be extremely different. One can extrapolate the complexity involved

when the influence of the various temporal as well as spatial discretizations are sought

on the basins of attractivity.

Table 4.1 summarizes the level of complexity for a systematic approach to these types

of PDE. The check mark on each type of PDE and approach indicate the ones where

some work has been done on this subject. The majority are credited to the University

of Dundee group [46-54] and some related theory by A. Stuart [58-82].

4.3. Involvement in the Study of Full Discretization of PDE

Consider a three-level explicit time differencing and a three-point spatial differencing

of the reaction-convection-diffusion equation (4.3) of the form

u" _ u" .-1 .-1 .-1 cAt, Ax), u_ +1 =u_ +t1( j_l,uj, j+l,uj_a,uj ,u_+l,a, (4.4)

where u_ is the numerical solution at / = nat and z = jAz. Then the study of the

asymptotes of (4.4) amounts to the study of fixed point behavior of period p in time and

period q in space, denoted by (p, q), where p and q are integers. Here the fixed point of

the partial-difference equation (4.4) is defined in a slightly more complicated way than

for the ODE.

,,+1 " and a fixed point For example, a fixed point of period (1,1) is defined as uj+ 1 = uj

_ n+2 " However, a fixed point of period (1,2) is defined of period (2,1) is defined as u i+1 = uj.

. n+p n

as uj+_- ,.,+1 = ujn. Thus, in general a fixed point of period (p,q) is defined as =j+q = uj.

One can see that for p, q > 3, solving the resulting nonlinear algebraic equation is very

involved, especially when physical boundary conditions and physical dissipation terms

as well as numerical boundary conditions [82,83,34] and numerical dissipation [47] are

3O
additional dimensionsof consideration. Current available work involved the studies

beyond the linearized stability limit of the schemes,and assumedthe nonexistenceof

spuriousfixed points of period (1,1). Seereferences[46-54]for details.

4.4. Influence in Dynamical Behavior by Property of the PDEs

and Schemes, and Treatment of the Source Terms

Although the general study of the dynamical behavior of partial-difference equations

for the conservation law [84,85] of (4.3) is an enormous task, if we can isolate certain

restricted subsets of the PDEs and schemes in hand which are immune to the type of

phenomena discussed in section III for time discretization as well as spatial discretization, then we can concentrate on the rest of the unknowns.

As can be seen in section III, the nature of the dynamical behavior of the discretized

counterparts is strongly influenced by properties of the numerical method and the types

and form of nonlinear DEs. Here we want to study the influence on the dynamical behavior of elements such as conservation and nonlinearity of the schemes, and treatment

of the source terms [14-17,78-81] when nonlinear conservation laws of PDEs are sought.

First, take the convection equation (4.2) with S(u) = 0 and consider a conservative

explicit scheme [76,14] which is consistent with the conservation law of the form

'_ Ah 1-

u_ +1 = ttj +_ 3-

where A = At/Az and h_.+} are the numerical flux functions. For a two-time level and

five-point spatial scheme, hj_ 1⁄2 = h(u'), uj_l , wj_2).

We also can consider a two-parameter family of scheme

u}'+l -+---1 + _o hi+1⁄2 .7__ = u., 1 +w 0) [h.j+1⁄2 -- hj_} .

rt--1

+1--77 - ). (4.6)

where 0 < 0 < 1. When 0 = 0, the scheme is explicit and when 0 = w + 1/2, the scheme

is temporally secorJd-order accurate. One can obtain (4.5) from (4.6) by setting 0 = 0

and w = 0. The time differencing belongs to the class of LMM. Under the assumption

that this scheme is conservative and consistent with the conservation law, discrete map

(4.6) will have no spurious steady-state numerical solution since consistency means

= I(='). (4.7)

Thus any steady-state solutions of (4.6) are steady-state solutions of the original PDE.

31
Now the situation is different when S(u) # O. Under this situation, even if the

same time and spatial discretization are employed, one still has to evaluate S properly.

Here S is the function S evaluated at some proper average state _ [14-17] for the

full discretization that is consistent with the scheme [18], and achieves conservation at

jumps. For a discussion on this subject, see references [78, 15-17] for details. The other

crucial aspect is that when S(u) # 0, a full investigation on the dynamical behavior

of the temporal and spatial discretization is necessary. The knowledge gained from the

finite-difference methods analysis for S(u) = 0 does not carry over to the S(u) # 0 case.

4.5. Discrete travelling Waves

Analysis of the dynamical behavior of the full discretization of nonlinear nonhomogeneous PDEs of the hyperbolic and parabolic types is very involved. In this section,

we look at a more restricted class of solutions -- the discrete travelling wave solutions.

Consider a reaction-diffusion equation

Ou 02 _t

& Oz 2

+ (4.8)

Solution u(z,t) depend on the space variable z and on the time t. Every zero of S(u)

constitutes an equilibrium of the PDE. Then a travelling wave solution is a profile U(z)

that travels along the z-axis with propagation speed X. Neither the shape of the wave

nor the speed of propagation changes. To find travelling waves, we seek solutions

resulting in an ODE

=(x,t) = u(x - it), (4.9)

u" + Xu' + = 0. (4.1o)

By solving this ODE, one can calculate asymptotic states for the PDE. Let U1 and

U2 be roots of ,_(u) and hence equilibrium solutions for both the PDE and ODE. The

asymptotic behavior of solution U for x _ 4-oo determines the type of travelling wave.

Every solution with

U(c_) = ul (4.11a)

U(-oo) = u2, (4.11b)

with ul # us, is a front wave of the ODE. This corresponds to a heterochnic orbit

[3] of the ODE, connecting the two stationary points ul and u2. Here for a secondorder autonomous ODE (4.10), when distinct saddles are connected, one encounters

a heteroclinic orbit; also a heteroclinic orbit may also join a saddle to a node or vice

versa. Another type of special orbit is a homoclinic orbit. A homoclinic orbit connects a

saddle point to itself and such orbits have an infinite period. Several heteroclinic orbits

32
may form a closedpath called a homoclinic cycle. Both the heteroclinic and homoclinic

orbits are of great interest in applications because they form the profiles of travelling

wave solutions of many reaction-diffusion problems. See references [3,10,73,74] for a

discussion.

Similarly, one can study discrete travelling wave solutions for the finite discretization

of (4.8). See references [73,74] for a discussion. Understanding of the discrete traveling

wave solutions of the corresponding PDEs only gives insight into a very small subset of

the dynamics of the PDEs. In most cases, it provides no information at all for the fully

discretized equation.

33
V. IMPLICATIONS & RECOMMENDATIONS

Due to the complexity of the large increase in system dimension and the involvement

of multiple floating parameters for finite difference methods in PDEs, we are not certain

that a similar systematic general result can be arrived at for more complex nonlinear

systems. The main indication at this point is from our time discretization study.

5.1. Results Drawn from the ODE Connection Study

Our study illustrates a few very important implications which are very fundamental

in explaining what happens when linear stability breaks down for truly nonlinear problems; i.e., equations that display genuinely nonlinear types of behavior. The important

points are as follows:

(1) There is sensitivity to initial data and strong dependence on discretization parameters such as the time step and the grid spacing Az. Dependence of solutions on

initial condition is important for employing a time-dependent approach to the steadystate with a given initial condition and boundary conditions in hypersonic or combustion

flows, especially when initial data of the governing PDE are not known.

(2) Associated with the same (common) steady-state solution the basin of attraction

of the DEs might be vastly different from the discretized counterparts. This is mainly

due to the dependence and sensitivity on initial conditions and boundary conditions

for the individual systems. In the absence of the influence of the initial and boundary

conditions, the difference in the basins of attraction between the continuum and its

discretized counterparts occurs even when an implicit LMM type of method is used

unless the resulting nonlinear algebraic equations are solved exactly.

(3) Nonunique steady-state solutions can be introduced by the spatial discretization

even though the original PDEs might possess only an unique steady-state solution and a

LMM type of time discretization is used so that no spurious steady-state exists in time.

The tie between temporal and spatial dynamical behavior is more severe when one is

dealing with the nonseparable temporal and spatial finite-difference discretization such

as the Lax-Wendroff type, where the time and spatial difference cannot be separated

from each other. The situation would be more complicated if the governing nonlinear

PDE possesses more than one steady-state solution as well as the spurious ones that

are purely due to the numerical method.

(4) For certain time discretizations, spurious steady-state solutions may occur below

the linearized stability limit of the scheme.

(5) The result of operating with a time step beyond the linearized stability limit is

not always a divergent solution; spurious steady-state solutions can occur.

34
(6) There is a nfisconceptionthat computational instability or inaccuracy can often

becured simply by making Ai smaller. Other elementssuchas(1) - (5) aboveaswell as

the variation of the grid spacings,numerical dissipation terms and system parameters

other than the time stepscaninterfere with the dynamical behavior.

(7) When linearized stability limits are usedas a guide for a time step constraint

for highly couplednonlinear system problems,this time step might exceedthe actual

linearized stability limit of the coupled equations. Therefore all of the situations in

(1) - (6) canoccur. In particular, whenone tries to stretch the maximum limit of the

linearized allowabletime stepfor highly coupledsystems,most likely all of the different

type of spuriousbranchesof supercritical, subcritical and trancritical bifurcations can

be achievedin practice dependingon the initial conditions. That is why the occurrence

of spurious steady-state solutions beyond the linearized stability limit is not just secondary but might be as important as the occurrence of spurious steady states below the

hnearized stability limit.

5.2. Recommendations

II is of utmost importance to know the nonlinear dynamical behavior of the various

schemes before their actual use for practical applications. Otherwise, it might be very

difficult to asses the accuracy (spurious or otherwise) of the solution when the numerical

method is the sole source of the understanding of the physical solutions. When in doubt,

it is always safer to use schemes that do not produce spurious steady-state solutions for

the nonlinear scalar case. Some examples of methods of this type in time discretization

can be listed:

(1) LMM [56] ODE solvers such as the exphcit, implicit Euler, three-point backward

differentation, etc. can be used.

(2) One can use the "Regular" Runge-Kutta methods [57].

(3) Solving the nonlinear algebraic systems arising from implicit LMM method exactly would avoid spurious steady state numerical solutions. Otherwise, the type of

iteration method in solving nonlinear algebraic systems can degrade the basin of attractivity of imphcit LMM [57].

The insight gained from time discretization will only give an indication in separable

schemes or method of lines approaches. Also, the commonly used residual test [86-88]

in the time-dependent approach to the steady state might be misleading. This is the

direct consequence of what was indicated in section 5.1. The popular misconception of

using the inverse problem of nonlinear dynamics to analyze a time series data from a

finite difference method computer code in an attempt to learn about the true physical

solution behavior of the continuum governing PDEs without knowing by other means

the exact solution behavior of the PDEs other than the numerical solutions can also be

misleading. These will be discussed in the next two sections.

35
5.3. Residual Test

Consider a quasilinear PDE of the form

Ou

N = (5.1)

where G is nonllnear in u, u_ and uzz and a and _ are system parameters. For simphcity,

consider a two time level and a (p + q) point grid stencil of the form

_+1 = u)n -- H(uj+q,..., n uj,..., _ uj_p, rt a, e, At, Az) (5.2)

for the PDE (5.1). Let U*, a vector representing (u_+q, ..., u_, ..., u__p) be a steady-state

numerical solution of (5.2). It is a common practice in CFD to use a time dependent

approach such as (5.2) to solve the steady-state equation G(u, u,, u,,,a, e) = 0. The

iteration is stopped when the residual H or some L2 norm of the dependent variable u

between two successive iterates is less than a pre-selected level.

Aside from the various standard numerical error such as truncation error, machine

round-off error, etc. [89], there is a more fundamental question on the validity of the

residual test and/or L2 norm test. If the scheme happens to produce spurious steadystate numerical solutions, these spurious solutions would still satisfy the residual and L2

norm tests in a deceptively smooth manner. Moreover, aside from the spurious solutions

issue, depending on the combination of time as well as spatial discretizations, it is not

easy to check whether G(u*,u_,u*_,_,c) _ 0 even though H(U*,a,c, At, Az) ---, O.

This is contrary to the ODE case, where if u* is spurious in (1.1) then S(u*) _ O.

Among other factors, this is one of the contributing factors in the increase in magnitude

of difficulty for analyzing the dynamical behavior of numerical methods for hyperbolic

and parabolic PDEs.

One might argue that one can judge the accuracy of the scheme by comparing the

numerical solutions with more than one numerical methods and by doing a sequence of

grid refinement and time step reductions. The latter approach might not be feasible at

an acceptable cost. The former might not be foolproof if one does not know the dynamical behavior of the finite difference schemes being used. One important contributing

factor on the use of the Lax-Wendroff types of schemes [90,91] is that these schemes

are more accurate and sometimes more stable when operated on or near the linearized

stability limit.

5.4. The Inverse Problems of Nonlinear Dynamics

The use of the inverse problem of nonlinear dynamics to analyze the dynamical behavior of time series data arising from experimental or observable data has received much

attention in nonlinear physics as well as in many of the engineering disciplines. The

approach is very useful for gaining some insights into the nonlinear dynamical behavior

36
in problems where experimental or observabledata are the main source of information. Often the associated governing equations (continuum or otherwise) do not exist

to start with. There has been an explosion of theory, numerical procedures and computer software addressing this rapidly growing direction [92-95]. There also has been

much recent interest in forecasting algorithms that attempt to analyze a time series

by fitting nonlinear models. The attractive feature of this approach is that when used

correctly on the correct problems one can reduce the complexity of the problem from

un-manageable higher dimensions to a very low dimension. It is therefore a natural

tendency for practioners in computational sciences to apply this approach to analyze

the dynamical behavior of time series data from a finite difference method computer

code in an attempt to learn about the true physical solution behavior of the governing

PDEs. This application of time series analysis can be misleading and can lead to a

wrong conclusion if the practitioner does not know by other means the exact solution

behavior of the PDEs other than from the numerical solutions. Examples of the use of

this type of approach in CFD computations have been presented in references [96-98].

It can be seen from our study that the conclusions drawn from this type of time series

analysis provide very little information, but rather can actually mislead one as to the

true physics of the problem.

VI. CONCLUDING REMARKS

Spurious stable as well as unstable steady-state numerical solutions, spurious asymptotic numerical solutions of higher period, and even stable chaotic behavior can occur

when finite-difference methods are used to solve nonlinear DEs numerically. The occurrence of spurious asymptotes is independent of whether the DE posseses a unique

steady state or has additional periodic solutions and/or exhibits chaotic phenomena.

The form of the nonlinear DEs and the type of numerical schemes are the determining

factor. In addition, the occurrence of spurious steady states is not restricted to the time

steps that are beyond the linearized stabihty limit of the scheme. In many instances,

it can occur below the linearized stability lirrdt. Therefore, it is essential for practitioners in computational sciences to be knowledgeable aboul the dynamical behavior of

finite-difference methods for nonlinear scalar DEs before the actual application of these

methods to practical computations. It is also important to change the traditional way

of thinking and practices when dealing with genuinely nonlinear problems.

In the past, spurious asymptotes were obserw_d in numerical computations but tended

to be ignored because they all were assumed to lie beyond the hnearized stabihty hmits of the time step parameter At. As can be seen from our study, bifurcations to and

from spurious asymptotic solutions and transitions to computational instability not only

are highly scheme dependent and problem dependent, but also initial data and boundary condition dependent, and not limited to time steps that are beyond the linearized

stability limit.

The symbiotic relation among all of these various factors makes this topic fascinating

37
and yet extremely complex. The main fundamental conclusion is that, in the absence of

truncation and machine round-off errors, there are qualitative features of the nonlinear

DE which cannot be adequately represented by the finite-difference methods and vice

versa. The major feature is that convergence in practical calculations involved fixed At

as n _ oo rather than At _ 0 as n _ c¢. It should be emphasized that the resulting

discrete maps from finite discretizations can exhibit a much richer range of dynamical behavior than their continuum counterparts. A typical feature is the existence of

spurious numerical asymptotes that can interfere with stability, accuracy and basins of

attraction of the true physics of the continuum.

38
References

I1]

[2]

[3]

[4]

[5]

I6]

[7]

[8]

f9]

[lO]

[11]

[12]

[131

[14]

[15]

[17]

NASA Computational Fluid Dynamics Conference, NASA Conference Publication 10038, Vol. 1 and 2, March 7-9, 1989.

R.L. Devaney, An Introduction to Chaotic Dynamical Systems, Addison Wesley,

New York, 1987.

R. Seydel, From Equilibrium to Chaos, Elsevier, New York, 1988.

J. Guchenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems,

and .Bifurcations of Vector Fields, Springer-Verlag, New York, 1983.

J.M.T. Thompson and H.B. Stewart, Nonlinear Dynamics and Chaos, John Wiley, New York, 1986.

C.S. Hsu, Cell-to-Cell Mapping, Springer-Verlag, New York, 1987.

M. Kub_cek and M. Marek, Computational Methods in Bifurcation Theory and

Dissipative Structures, Springer-Verlag, New York, 1983.

T.S. Parker and L.O. Chua, Pnactical Numerical Algorithms for Chaotic Systems,

Springer-Verlag, New York, 1989.

E. A. Jackson, Perspectives of Nonlinear Dynamics, Cambridge, Cambridge,

1989.

E. Beltrami, Mathematics for Dynamic Modeling, Academic Press, Orlando,

1987.

R.M. May, "Biological Populations Obej'ing Difference Equations: Stable Points,

Periodic Orbits and Chaos," J. Theoret Biol., Vol. 51, 1975, pp. 511-524.

H.C. Yee, "A Study of Two-Dimensional Nonlinear Difference Systems and Their

Applications," Ph.D. Dissertation, University of Calif, Berkeley, Calif., USA,

1975.

C.S. Hsu, "On Nonlinear Parameteric Excitation Problems," Advances in Applied

Mechanics, Vol. 17, 1977, pp. 245-301, Academic Press, New York.

I-I.C. Yee, "A Class of High-Resolution Explicit and Implicit Shock-Capturing

Methods," NASA TM-101088, Feb. 1989.

P.K. Sweby, " "TVD" Schemes for Inhomogeneous Conservation Laws," Proc.

Nonlinear Hyperbolic Equations, Eds. J. Ballmann & R. Jeltsch, Notes on Numerical Fluid Machanics, 24, 1989, pp. 599-607.

P.K. Sweby, "Source Terms and Conservation Laws: A Preliminary Discussion,"

Numerical Analysis Report 6/89, Dep;_rtment of Mathematics, University of

Reading, England.

B. Engquist and B. Sj6green, "Numerical Approximation of Hyperbolic Conservation Laws with Stiff Terms," CAM Report 89-07, Department of Mathematics,

39
UCLA, Calif. USA, March 1989.

[18] K.W. Wallace, "Spurious Period Doubling Bifurcations in a Discretized ReactionAdvection Model," M.S. Thesis, Faculty of Science, University of Dundee, Scotland, Sept. 1989.

[19] G. Strang, "On the Construction and Comparison of Difference Schemes," SIAM

J. Num. Anal. Vol. 5, 1968, pp. 506-517.

[20] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial-Value Problems, Interscience-Wiley, New York, 1967.

[21] P.K. Sweby, H.C. Yee and D.F. Griffiths, "On Spurious Steady-State Solutions

of Explicit Runge-Kutta Schemes," NASA TM, March 1990.

[22] A.M. Panov, "Behavior of the Trajectories of a System of Finite Difference Equations in the Neighbourhood of a Singular Point," Uch. Zap. Ural. Gos. Univ.

vyp, Vol. 19, 1956, pp. 89-99.

[23] O. Perron, "Uber Stabilit£t und Asymtotisches l-.lberhalten de LSsungen eines

Systems endlicher Differenzengleichungen," J. Reine Angew. Math. Vol. 161,

1929, pp. 41-64.

[24] C.S. Hsu, H.C. Yee and W.H. Cheng, "Determination of Global Regions of

Asymptotic Stability for Difference Dynamical Systems," J. Appl. Mech., Vol.

44, pp., 1977, pp. 147-153.

[25] C.S. Hsu, H.C. Yee and W.H. Cheng, "Steady-State Response of a Nonlinear

System Under Impulsive Periodic Parametric Excitation," J. Sound Vib., Vol.

50, 1977, pp. 95-116.

[26] R.M. May, "Simple Mathematical Models with Very Complicated Dynamics,"

Nature, Vol. 261, 1976, pp. 459-467.

[27] R.M. May, "Biological Populations with Nonoverlaping Generations: Stable

Points, Stable Cycles, and Chaos," Science, Vol. 186, No. 15, 1974, pp. 645-647.

[28] T.Y. Li and J.A. Yorke, "Period Three Implies Chaos," Am. Math. Monthly,

Vol. 82, 1975, pp. 985-992.

[29] E.N. Lorenz, "The Problem of Deducing the Climate from the Governing Equations," Tellus, Vol. 16, 1964, pp. 1-11.

[30] M.J. Feigenbaum, "Quantitative Universality for a Class of Nonlinear Transformations," J. Stat. Phys., Vol. 19, 1978, pp. 25-52.

[31] C.S. Hsu and H.C. Yee, "Behavior of Dynamical Systems Governed by a Simple

Nonlinear Difference Equation," J. Appl. Mech., Vol. 44, 1975, pp. 870-876.

[32] S. Ushiki, "Central Difference Scheme and Chaos," Physica 4D, 1982, pp. 407-

424.

[33] F. Brezzi, S. Ushike and H. Fujii, "Real and Ghost Bifurcation Dynamics in Dif40
ferenceSchemesfor Ordinary Differential Equations," Numerical Methods for Bifurcation Problems, T. Kupper, H.D. Mittleman and H. Weber eds., BirkhauserVerlag, Boston, 1984.

[34] R. Schreiber and H.B. Keller, "Spurious Solution in Driven Cavity Calculations,"

J. Comput. Phys., Vol. 49, No. 1, 1983.

[35] W.J. Beyn and E.J. Doedel, "Stability and Multiplicity of Solutions to Discretizations of Nonlinear Ordinary Differential Equations," SIAM J. Sci. Statist.

Comput., Vol. 2, 1981, pp. 107-120.

[36] R.B. Kellogg, G.R. Shubin, and A.B. Stephens, SIAM J. Numer. Anal. Vol. 17,

No. 6, 1980, pp. 733-739.

[37] A.B. Stephens and G.R. Shubin, SIAM J. Sci. Statist Comput., Vol. 2, 1981,

pp. 404-415.

[38] G.R. Shubin, A.B. Stephens and H.M. Glaz, "Steady Shock Tracking and Newton's Method Applied to One-Dimensional Duct Flow," J. Comput. Phys., Vol.

39, 1981, pp. 364-374.

[39] P.G. Reinhall, T.K. Caughey and D.W. Storti, "Order and Chaos in a Discrete

Duffing Oscillator: Implications on Numerical Integration," Trans. of the ASME,

J. Appl. Mech., 89-APM-6, 1989.

[40] E.N. Lorenz, "Computational Chaos - A Prelude to Computational Instability,"

Physica D, Vol. 35, 1989, pp. 299-317.

[41] A.J. Lichtenberg and M.A. Lieberman, Regular and Stochastic Motion, Appl.

Math. Sci. Bd. 38, 1983, Springer-Verlag, New York.

[42] R.H. Miller, "A Horror Story about Integration Methods," to appear.

[43] W.A. Mulder and B. van Leer, "Implicit Upwind Methods for the Euler Equations," AIAA-83-1930, July 1983.

[44] M. Priiffer, "Turbulence in Multistep Methods for Initial Value Problems," SIAM

J. Appl. Math. Vol. 45, 1985, pp. 32-6!3.

[45] W.-J. Beyn, "On the Numerical Approximation of Phase Portraits Near Stationary Points," SIAM J. Numer. Anal., Vol. 24, No. 5, 1987, pp. 1095-1113.

[461 A.R. Mitchell and D.F. Griffiths, "Beyond the Linearized Stability Limit in Non

Linear Problems," Report NA/88 July 1985, Department of Mathematical Sciences, University of Dundee, Scotland U.K.

[47] A.R. Mitchell and J.C. Bruch, Jr., "A Numerical Study of Chaos in a ReactionDiffusion Equation," Numerical Methods for PDEs, Vol. 1, 1985, pp. 13-23.

[48] A.R. Mitchell, P. John-Charles and B.D. Sleeman, "Long Time Calculations and

Non Linear Maps," Numerical Analysis Report 93, May 1986, Department of

Mathematical Sciences, University of Dundee, Scotland.

41
[49]

[5o]

V.S. Manoranjan, A.R. Mitchell, and B.D. Sleeman, "Bifurcation Studies in

Reaction-Diffusion," J. Comput. App. Math., Vol. 11, 1984, pp. 27-37.

B.D. Sleeman, D.F. Griffiths, A.R. Mitchell and P.D. Smith, "Stable Periodic

Solutions in Nonlinear Difference Equations," SIAM 3. Sci. Stat. Comput., Vol.

9, No. 3, May 1988, pp. 543-557.

[51] D.F. Griffiths and A.R. Mitchell, "Stable Periodic Solutions of a Nonlinear Partial

Difference Equation in Reaction Diffusion," Report NA/113, Jan. 1988, Dept.

Math. and Compt. Science, University of Dundee, Scotland.

[52] A.R. Mitchell, G. Stein and M. Maritz, "Periodic Structure Beyond a Hopf Bifurcation," Comm. Appl. Num. Meth., Vol. 4, 1988, pp. 263-272.

[53] D.F. Griffiths and A.R. Mitchell, "Stable Periodic Bifurcations of an Exphcit

Discretization of a Nonlinear Partial Differential Equation in Reaction diffusion,"

Inst. Math. Applics., J. Num. Analy., Vol. 8, 1988, pp. 435-454.

[54] A.R. Mitchell and S.W. Schoombie, "Nonlinear Diffusion and Stable Period 2

Solutions of a Discrete Reaction-Diffusion Model," J. Comp. Appl. Math., Vol

25, 1989, pp. 363-372.

[55]

[56]

[57]

[58]

J.M. Sanz-Serna and F. Vadillo, "Nonlinear Instability, the Dynamic Approach,"

in Proceedings Dundee, 1985, G.A. Watson and D.F. Griffiths, eds., Pitman,

London.

A. Iserles, "Nonlinear Stability and Asymptotics of O.D.E. Solvers," International Conference on Numerical Mathematics, Singapore, R.P. Agarwal, ed.,

Birkhauser, Basel, 1989.

A. Iserles and J.M. Sanz-Serna, "Equilibria of Runge-Kutta Methods," Numerical Analysis Reports, DAMTP 1989/NA4, May 1989, Univeristy of Cambridge,

England.

A.M. Stuart, "Linear Instability Implies Spurious Periodic Solutions," submitted

to IMA J. Num. Anal., 1988.

[59] A. Stuart, "The Global Attractor Under Discretisation," to appear, Proc. NATO

Conference on Continuation & Bifurcation, 1989.

[60] A. Stuart and A. Peplow, "The Dynamics of the Theta Method," to appear,

SIAM J. Sci. Stat. Comput.

[61]

[62]

[63]

[64]

A. Stuart, " Nonlinear Instability in Dissipative Finite Difference Schemes,"

SIAM Review, Vol. 31, No. 2, 1989, pp. 191-220.

A. Stuart, "Linear Instability Implied Spurious Periodic Solutions," to appear,

IMA J. Numer. Anal.

C. Grebogi E. Ott and J. Yorke, "Chaos, Strange Attractors, and Fractal Basin

Boundaries in Nonlinear Dynamics," Science, Vol. 238, 1987, pp. 585-718.

S.W. McDonald, C. Grebogi E. Ott and J. Yorke, "Fractal Basin Boundaries,"

42
Physica 17D, 1985, pp. 125-153.

[65] C. Grebogi E. Ott and J. Yorke, "Crises, Sudden Changes in Chaotic Attractors,

and Transient Chaos," Physics 7D, 1983, pp. 181-200.

[66] J.D. Lambert, Computational Methods in Ordinary Differential Equations, John

Wiley, New York, 1973.

[67] C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations

(Prentice-Hall, 1971 ).

[68] J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations, John

Wiley, New York, 1987.

[69] MAPLE, algebraic manipulation package, University of Waterloo, Canada, 1988

[70] D. Whitley, "Discrete Dynamical Systems in Dimensions One and Two," Bull.

London Math Soc. Vol. 15, 1983, pp. 177-217.

[71] R.M. Beam and R.F. Warming, "Implicit Numerical Methods for the Compressible Navier-Stokes and Euler Equations," Lecture Notes for Computational Fluid

Dynamics, yon Karman Institute for Fhfid Dynamics, March 29 - April 2, 1982,

Rohde-Saint- Gen_se, Belgium.

[72] J. Maynard Smith, Mathematical ]deas in Biology, Cambridge University Press,

Cambridge, 1968.

[73] H. Ikeda, M. Mimura and Y. Nishiura, "Global Bifurcation Phenomena of Travelling Wave Solutions for Some Bistable Reaction-Diffusion Systems," Nonlinear

Analysis, Theory, Methods and Applications, Vol. 13, No. 5, 1989, pp. 507-526.

[74] T. Hagstrom and H.B. Keller, "The Numerical Calculation of travelling Wave

Solutions of Nonlinear Parabolic Equations," SIAM J. Set. Stat. Comput., Vol.

7, No. 3, 1986.

[75] A. Harten and S. Osher, "Uniformly High-Order Accurate Nonoscillatory

Schemes 1," SIAM J. Num. Analy. Vol. 24, No. 2, 1987, pp. 279-309.

[76] A. Harten, "On a Class of High Resolution Total-Variation-Stable FiniteDifference Schemes," SIAM J. Num. Anal., Vol. 21, 1984, pp. 1-23.

[77] P.L. Roe, "Some Contributions to the Modelling of Discontinuous Flows," Lectures in Applied Mathematics, Vol. 22 (Amer. Math. Soc., Providence, R.I.,

1985).

I78] R.J. LeVeque and H.C. Yee, "A study of Numerical Methods for Hyperbolic

Conservation Laws with Stiff Source Terms," NASA TM-100075, March 1988, J.

Comput. Phys., Vol. 86, No. 1, Jan. 1990.

[79] M. Pandolfi, M. Germano and N. Botta, "Non-Equilibrium Reacting Hypersonic

Flow About Blunt Bodies: Numerical Prediction, AIAA-88-0514, Jan. 1988.

[80] P. Colella, A. Majda, and V. Roytburd, "Theoretical and Numerical Structure

for Reacting Shock Waves," SIAM J. Sci. Stat. Comput. Vol. 7, 1986, pp.

43
[81]

[82]

[83]

[84]

[85]

[86]

[87]

[88]

[89]

[90]

[91]

[92]

[93]

[94]

[95]

1059-1080.

T.R. Young and J.P. Boris, "A Numerical Technique for Solving Stiff Ordinary

Differential Equations Associated with the Chemical Kinetics of Reactive-Flow

Problems," J. Phys. Chem. Vol. 81, 1977, pp. 2424-2427.

tt.C. Yee, "Numerical Approximation of Boundary Conditions with Applications

to Inviscid Equations of Gas Dynamics," NASA TM-81265, 1981.

H.C. Yee, R.M. Beam and R.F. Warming, "Boundary Approximations for Implicit Schemes for One-Dimensional Inviscid Equations of Gasdynamics," AIAA

J., Vol. 20, No. 9, 1982, pp. 1203-1211.

P.D. Lax and B. Wendroff, "Systems of Conservation Laws," Commun. Pure

Appl. Math., Vol. 13, 1960, pp. 217-237.

R.J. LeVeque, "Hyperbolic Conservation Laws and Numerical Methods," Lecture Series on Computational Fluid Dynamics, yon Karman Institute for Fluid

Dynamics, Rhode-St-Genbse, Belgium, March 5-9, 1990.

A. Jameson, W. Schmidt and E. Turkel, "Numerical Solutions of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes,"

AIAA-81-1259, 1981.

J.D. Anderson, Jr, "Introduction to Computational Fluid Dynamics," yon Karman Institute for Fluid Dynamics, 1985 lecture series, Rhode-Saint-Gen_se, Belgium.

G.A. Sod, Numerical Methods in Fluid Dynamics, Cambridge University Press,

Cambridge, 1985.

J.H. Ferziger, "Estimation and Reduction of Numerical Error," Forum on Methods of Estimating Uncertainty Limits in Fluid Flow Computations, ASME Winter Annual Meeting, San Francisco, Dec. 1989.

P.D. Lax and B. Wendroff, "Difference Schemes for Hyperbolic Equations with

High Order of Accuracy," Commun. Pure Appl. Math., Vol. 17, 1964, pp.

381-398.

R.W. MacCormaek, "The Effect of Viscosity in Hypervelocity Impact Cratering,"

AIAA-69-354, Cincinnati, Ohio, 1969.

N.I-I. Packard, J.P. Crutchfield, J.D. Farmer and R.S. Shaw, "Geometry from a

Time Series," Physical Review Letters, Vol. 45, No. 9, 1980, pp. 712-716.

J.P. Eckmann and D. Ruelle, "Ergodic Theory of Chaos and Strange Attractors,"

Rev. Mod. Phys., Vol. 57, No. 3, Part I, 1985, pp. 617-656.

H. Froehling, J.P. Crutchfield, D. Farmer, N.H. Packard and R. Shaw, "On Determing the Dimension of Chaotic Flows," Physica 3D, 1981, pp. 605-617.

E.J. Kostelich and J.A. Yorke, "The Analysis of Experimental Data Using TimeDelay Embedding Methods," Institute for Physical Science and Technology Technical Report, University of Maryland, College Park, Maryland, USA.

44
[96] T.H. Pulliam, "Low ReynoldsNumber Numerical Solutions of Chaotic Flows,"

AIAA-89-1023, Jan. 9-12,1989.

[97] T.H. PuUiam,"Numerical Simulation of Chaotic Flows: Measuresof Chaos," Forum on Chaotic Dynamicsin Fluid Dynamics, ASME Fluids Engineering Spring

Conference,La Jolla, Calif. July 1989.

[981 A. Fortin, M. Fortin and J.J. Gervais, "A numerical Simulational Flow," J. Cornput. Phys., 70, 1987, pp. 295-310.

45
Figure Captions

Table 2.1 Possible stable asymptotic solution behavior for DEs and their discretized

counterparts.

Fig. 2.1 Phase portrait and basins of attraction of the damped pendulum equation

(this figure is taken from reference [8]).

Fig. 2.2 Phase portrait and basins of attraction of the predator-prey equation (this

figure is taken from reference [8]).

Fig. 3.1 Asymptotic solution behavior of the logistic ODE du/dt = au(1 - u) for

a>0.

Fig. 3.2 Stable fixed points of periods 1,2,4,8 of the explicit Euler scheme for the

logistic ODE du/dt = au(1 - u).

Fig. 3.3 Stable fixed points of periods 1,2,4,8 of the modified Euler (R-K 2) scheme

for the logistic ODE du/dt = au(1 - u).

Fig. 3.4 Stable fixed points of periods 1,2,4,8 of the improved Euler (R-K 2) scheme

for the logistic ODE du/dt = au(1 - u).

Fig. 3.5 Stable fixed points of periods 1,2,4,8 of the Runge-Kutta 4th-order (R-K 4)

scheme for the logistic ODE du/dt = au(1 - u).

Fig. 3.6 Stable fixed points of periods 1,2,4,8 of the predictor-corrector scheme of

order 2 for the logistic ODE du/dt = au(1 - u).

Fig. 3.7 Stable fixed points of periods 1,2,4,8 of the predictor-corrector scheme of

order 3 for the logistic ODE du/dt = au(1 - u).

Fig. 3.8 Stable fixed points of periods 1,2,4,8 of the explicit Euler scheme for the

ODE du/df = au(1 - u)(0.5- u).

Fig. 3.9 Stable fixed points of periods 1,2,4,8 of the modified Euler (R-K 2) scheme

for the ODE du/dt = au(1 - u)(0.5- u).

Fig. 3.10 Stable fixed points of periods 1,2,4,8 of the improved Euler (R-K 2) scheme

for the ODE du/dt = au(1 - u)(0.5- u).

Fig. 3.11 Stable fixed points of periods 1,2,4,8 of the Runge-Kutta 4th-order (R-K 4)

scheme for the ODE du/dt = au(1 - u)(0.5 - u).

46
Fig. 3.12

Fig. 3.13

Fig. 3.14

Fig. 3.15

Fig. 3.16

Fig. 3.17

Fig. 3.18

Fig. 3.19

Fig. 3.20

Fig. 3.21

Fig. 3.22

Fig. 3.23

Fig. 3.24

Fig. 3.25

Fig. 3.26

Stable fixed points of periods 1,2,,i,8 of the predictor-corrector scheme of

order 2 for the ODE du/df = au(1 - u)(0.5 - u).

Stable fixed points of periods 1,2,,1,8 of the predictor-corrector scheme of

order 3 for the ODE du/dl = au(1 - zt)(0.5 - u).

Stable fixed points of periods 1,2,4.8 of the modified Euler (R-K 2) scheme

for the ()DE du/dt = au(1 - u)(b -u), b = 0.1, 0.2, 0.3, 0.4.

Bifurcation diagram of the explicit Euler scheme for the logistic ODE du/dt_ =

au(l- ,,).

Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic

ODE du/dt = au(1 - u) with u ° = 2.7.

Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic

ODE dzL/dt = au(1 - u) with u ° = 1.5.

Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic

ODE du/dt = au(1 - u) with u ° = 0.25.

"Full" bifurcation diagram

logistic ODE du/dt = au(1

Bifurcation diagram of the

logistic ODE du/dt = au(1

Bifurcation diagram of the

logistic ()DE du/dt = au(1

of the modified Euler (R-K 2) scheme for the

--It).

Runge-Kutta 4th-order (R-K 4) scheme for the

- u) with u ° = 0.5.

Runge-Kutta 4th-order (R-K 4) scheme for the

- u) with multiple initial data.

"Full" bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for

the logislic ODE du/dt = au(1 - u).

Bifurcation diagrams of the improved Euler (R-K 2) scheme for the ODE

du/dt = ,_u(1 - u)(0.5 - u) for four different sets of initial input data.

Bifurcation diagrams of the Runge-Kutta 4th-order (R-K 4) scheme for the

ODE du/dt = au(1 - u)(0.5 - u) for four different sets of initial input data.

Bifurcation diagrams of the modified Euler (R-K 2) scheme for the ODE

du/dÂ = t,u(1 - u)(0.4 - u) for four different sets of initial input data.

"Full" bifurcation diagram of the improved Euler (R-K 2) scheme for the

logistic ODE du/dt = au(1 - u).

47
Fig. 3.27 "Full" bifurcation diagram of the Adam-Bashforth schemefor the logistic

ODE du/dt = au(1 - u).

Fig. 3.28 "Full" bifurcation diagram of the predictor-corrector scheme of order 2 for

the logistic ODE du/dt = au(1 - u).

Fig. 3.29 "Full" bifurcation diagram of the predictor-corrector scheme of order 3 for

the logistic ODE du/dt = au(1 - u).

Fig. 3.30 "Full" bifurcation diagrams of the modified Euler (R-K 2) scheme for the

ODE du/dt = au(1 - u)(b- u), b = 0.1, 0.2, 0.3, 0.4.

Fig. 3.31 "Full" bifurcation diagrams of the exphct Euler scheme for the ODE du/dt =

,_u(1 - u)(0.5 - _,).

Fig. 3.32 "Full" bifurcation diagram of the modified Euler (R-K 2) scheme for the

ODE du/dt = au(1 - u)(0.5 - u).

Fig. 3.33 "Full" bifurcation diagram of the Adam-Bashforth scheme for the ODE

du / dt = au(1 - u)(0.5 - u).

Fig. 3.34 "Full" bifurcation diagram of the improved Euler (R-K 2) scheme for the

ODE du/dt = au(1 - u)(0.5 - u).

Fig. 3.35 "Full" bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for

the ODE du/dt = au(1 - u)(0.5 - u).

Fig. 3.36 "Full" bifurcation diagram of the improved Euler (R-K 2) scheme for the

ODE du/dt = au(1 - u)(0.5- u) (enlarged).

Fig. 3.37 "Full" bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for

the ODE du/dt = au(1 - u)(0.5- u) (enlarged).

Fig. 3.38 "Full" bifurcation diagram of the predictor-corrector scheme of order 2 for

the ODE du/dt = au(1 - u)(0.5 - u).

Fig. 3.39 Types of branching points.

Fig. 3.40 Stabihty of solutions in the neighborhood of branch points, one-dimensional

case. -- stable, - - - unstable a,b,c,d: limit (regular turning) point;

e,f,g,h: bifurcation (double) points; i_j,k,l: bifurcation-limit (singular turning) points; m,n,o,p,q: additional possible cases when the dimension of u is

greater than one (this figure is taken from reference [7]).

48
Fig. 3.41 Stability of steady-statesolutions arising through three types of bifurcation

phenomena(-- stable,- - - unstable).

Fig. 3.42 Spuriousfixed points arising from transcritical bifurcations.

Fig. 3.43 Spuriousfixed points arising from subcritical bifurcation.

Fig. 3.44 Stable and unstable fixed points of periods 1,2 of the modified Euler (R-K

2) schemefor the logistic ODE du/dt = c_u(1 - u).

Fig. 3.45 Stable and unstable fixed points of periods 1,2 of the improved Euler (R-K

2) scheme for the logistic ODE du/dt = au(1 - u).

Fig. 3.46 Stable and unstable fixed points of periods 1,2 of the Runge-Kutta 4th-order

(R-K 4) scheme for the logistic ODE du/dt = au(1 - u).

Fig. 3.47 Stable and unstable fixed points of periods 1,2 of the predictor-corrector

scheme of order 2 for the logistic ODE du/dt = au(1 - u).

Fig. 3.48 Stable and unstable fixed points of periods 1,2 of the predictor-corrector

scheme of order 3 for the logistic ODE du/dt = au(1 - u).

Fig. 3.49 Stable and unstable fixed points of periods 1,2 of the modified Euler (R-K

2) scheme for the ODE du/dt = au(1 - u)(0.5 - u).

Fig. 3.50 Stable and unstable fixed points of periods 1,2 of the improved Euler (R-K

2) scheme for the ODE du/dt = au(1 - u)(0.5 - u).

Fig. 3.51 Stable and unstable fixed points of periods 1,2 of the Runge-Kutta 4th-order

(R-K 4) scheme for the ODE du/dt = au(1 - u)(0.5 - u).

Fig. 3.52 Stable and unstable fixed points of periods 1,2 of the predictor-corrector

scheme of order 2 for the ODE du/d_ = au(1 - u)(0.5- u).

Fig. 3.53 Stable and unstable fixed points of periods 1,2 of the predictor-corrector

scheme of order 3 for the ODE du/dt = au(1 - u)(0.5 - u).

Fig. 3.54 Stable and unstable fixed points of periods 1,2 of the modified Euler (R-K

2) scheme for the ODE du/dt = c_u(1 - u)(0.2 - u).

Table 4.1 Systematic approach - level of complexity.

49
SOLUTION

TYPE

# OF

ASYMPTOTES

OR

STEADY-STATE

SOLUTIONS

PERIODIC

SOLUTIONS

1

CHAOS

ODEsORPDEs

SINGLE

SINGLE

MULTIPLE

MULTIPLE

NO

YES

i

NO

YES

DISCRETIZEDCOUNTERPARTS

SINGLE

MULTIPLE

SAME# OF MULTIPLE

ADDITIONAL# OF MULTIPLE

YES

YES (+ EXTRA)

i i

YES

YES (+ EXTRA)

Table 2.1 Possible stable asymptotic solution behavior for DEs and their discretized

counterparts.

5O
ODE CONNECTION: GAIN INSIGHT INTO TIME DISCRETIZATIONOF PDEs

SCALAR

SYSTEM -- TIME SPLrI'riNGORMETHODOF LINES

I1. DISCRETETRAVELUNGWAVE: au au a2U

,, _- + C _ = _ E_ + .S(u)

iREACTION-DIFFUSION

SCALAR: |REACTION-CONVECTION

[REACTION-CONVECTION-DIFFUSION

III. FULL DISCRETIZATION(TEMPORAL AND SPATIAL):

SCALAR:{ (S_O) _'] "

LINEARSCHEMEFORSPATIALDISCRETIZATION

SCALAR:

IS=O) { NONLINEARSCHEMEFORSPATIALDISCRETIZATION

SCALAR:

(S_O) { NONLINEARSCHEMEFORSPATIALDISCRETIZATION

Table 4.1 Systematic approach - level of complexity.

51
q

q.

0 o

0

0

|

o.

0

I"

o.

o.=

0

0

|

0

0

|

0

°,_

cd

QJ

0

.

0

"_ o

°_._

cd

r,_ °_

52
0

_!-i:!-i-i-i-i-ii:i_sisi:i_ _...'_. i 1

 iiiiiiiiiiiiiiiiiii i!iiii 

iiiiiil '°°°°°°'°°°°°°°°°° !!iiiiii °°°°%°°°°°°° Ii_

_""'"""" '"'"" 's/"d_ r,,.,

iiiiii!!ii_::SiiiS: fi.iiSF° __ _ • _i

I,,-.-.'.-F "FFFFF: •

.°.'"'F: FFFFFF _ _""_. •

":'FF:'.t FF"FF: • _ _ . •

:iiii!i!iii!!i!!iiiii!i::!!::::i:!:!:i:i:!::-:... _ i k _" m

...................:.:.:.Z.;.:.;.:.:.:.:.:.;.:.:.:.:.:." :.:.:.:.:.:.:.,.:. ".'.

0 _ 0 It_ 0

I

°_,ii

I=I

0

0

ca

0

0

°_

0

_ 0

0 _

_.,u=

e,i

_b

°_

53
I!

A

O

^

O

I!

,B

I

N

©

°_

°f,q

O

O

O

O

°_

O

bA

o_

54
Un

.8-

.4-

1E

4 8

f_ 4 8

.I

[

\

\

\

4 \

8

0 I I I I I

1.8 2.0 2.2 2.4 2.6 2.8

r = a,At

I

3.0

Fig. 3.2 Stable fixed points of periods 1,2,4,8 of the exphcit Euler scheme for the

logistic ODE du/dt = au(1 - u).

55
3.0 -

2.O

Un

1.0

\'\is

8

\ 4

\-- 8

",8

1E

4 8

0 I I I I I

1.0 1.5 2.0 2.5 3.0 3.5

r =orAt

I

4.0

Fig. 3.3 Stable fixed points of periods 1,2,4,8 of the modified Euler (R-K 2) scheme

for the logistic ODE du/dt = a_t(1 - _t).

56
Un

1.,2-

.8

1E

°4-

•0 |

/

1.5

1S 8

8

2 4_ 8

4

4

2_\ 8 8_

1S "_// 8 8_ 4

! I

2.0

\

4

.\

I I I

2.5 3.0

r = aAt

I

3.5

Fig. 3.4 Stable fixed points of periods 1,2,4,8 of the improved Euler (R-K 2) scheme

for the logistic ODE dtt/dl = au(1 - 4).

57
1.2-

.8

Un

.4-

0

2.25

1E

1S 2 4

2 4

__4_ 2_4

18 4

8

I I I I I I

2.75 3.25 3.75

r = _At

I

4.25

Fig. 3.5 Stable fixed points of periods 1,2,4,8 of the Runge-Kutta 4th-order (R-K 4)

scheme for the logistic ODE du/dt = ctu(1 - u).

58
1.,?.

.8

O n

.4

0

1E

t r

2..4

4

\\

4

_2

4 8

4

4

4

4

\:

\

8

4

• ;7

2

I I I I I I I I I I

1.5 2.0 2.5 3.0 3.5 4.0

r = aAt

Fig. 3.6 Stable fixed points of periods 1,2,4,8 of the predictor-corrector scheme of

order 2 for the logistic ODE du/dt --: au(1 - u).

59
.8

tl II

.4

0

1E

1S

/

2 4

.

1s 2

\,

\

\ 4_

=/4

\

4",8

2

2

1.5

I I I I I I I I

2.0 2.5 3.0 3.5

r = a,At

I

4.0

Fig. 3.7 Stable fixed points of periods 1,2,4,8 of the predictor-corrector scheme of

order 3 for the logistic ODE du/dt = art(1 - u).

6O
Un

1o0 -

.8-

.6-

.4

.2

1E

.!

J

/

,/

/

\

I I I I I I I

9 11 13 15

r =ttAt

4. 8

4-8

Fig. 3.8 Stable fixed points of periods 1,2,4,8 of the explicit Euler scheme for the

ODE du/dt = au(1 - u)(0.5- u).

61
U n

1,0--

.6-

.4-

1E

1S _'_ 4/,'8

/

2 ?

.8

4_-. 8

0 I I I I I I I

7 8 9 10 11 12 13 14

r = aAt

Fig. 3.9 Stable fixed points of periods 1,2,4,8 of the modified Euler (R-K 2) scheme

for the ODE du/dt = au(1 - u)(0.5- it).

62
1.0 -

U n

.8

.6-

.4-

.2-

0

1E

2 4

.///-_ 4

// 2

" %2

i:

_'",, 4 2/

.o

(

1S •

_,>/.,,4.  ̧

I 1 I I I I I I I I

7 9 11 13 15 17

r = a/_t

Fig. 3.10 Stable fixed points of periods 1,2,4,8 of the improved Euler (R-K 2) scheme

for the ODE du/dt = _u(1 - ,,)(0.5-- u).

63
1.0 -

.8-

.6-

U n

.4

.2-

1E

lS

lS

2S

2S

2 4

2 4

2s

2 -

!

1s 2 .

i. _4

2.._4

/

lS 2

O I 1 I I I I

10 12 14 16

r = aAt

X

Fig. 3.11 Stable fixed points of periods 1,2,4,8 of the Runge-Kutta 4th-order (R-K 4)

scheme for the ODE du/dt = au(1 - u)(0.5 - u).

64
U n

1.0

.8-

m

o6 -

1E

.4-

.2

.0

6

/

/

k

\

2_.

....i.-<4

://

2 S _

/ 2 ,:

1S/

1S 2

2

I I I I I

8 10 12 14 16

r = orAt

I

18

Fig. 3.12 Stable fixed points of periods 1,2;1,8 of the predictor-corrector scheme of

order 2 for the ODE du/dl = au(1 - u)(0.5 - u).

65
1.0 -

U n

.8-

.6-

.4-

.2

0

6

1E

1S 2

_2 • ....

Is _ 1----_ _

_ 1s _1S_....._I_ 1s 1_2 2

/

/

I

/

\-. 1S ./ti..

"\\ l/ 2

• 2.

-...2 2 "

"2

I I I I I

8 10 12 14 16

r = oat

I

18

Fig. 3.13 Stable fixed points of periods 1,2,4,8 of the predictor-corrector scheme of

order 3 for the ODE du/dt = au(1 - u)(0.5 - u).

66

.8

Un

1E

/

4

8

0 3.0

Fig. 3.15 Bifurcation diagram of the explicit Euler scheme for the logistic ODE du/dt =

c_u(1 -u).

fi8
3.0

2.5

2.0-

u n 1.5

1.0

.5-

0

1.0

- lS

m

".... ̧

1E

I I I I I

1.5 2.0 2.5 3.0 3.5

r = aAt

I

4.0

Fig. 3.16 Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic

ODE du/dI au(1 _t) with u ° "

69
3.0 -

2.5-

2.0-

u n 1.5 -

1.0

1E

.5

0

1.0

,2

I I T I _ 1

1.5 2.0 2.5 3.0 3.5 4.0

r = aAt

Fig. 3.17 Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic

ODE du/dl = au(1 - u) with u0 = 1.5.

7O
3.0

2.5

2.0

u n 1.5

1.0

.5

0

1E

2

2

1i

I

I I I I I I

1.0 1.5 2.0 2.5 3.0 3.5 4.0

r = <xht

Fig. 3.18 Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic

ODE d_/dt = au(] - 71) with _t° = 0.25.

71
3.0

2.5

2.0-

un 1.5-

1.0

.5

1E

"FULL" BIFURCATION DIAGRAM

0 I I I I I

1.0 1.5 2.0 2.5 3.0 3.5

r : _At

I

4.0

Fig. 3.19 "Full" bifurcation diagram of the modified Euler (R-K 2) scheme for the

logistic ODE du/dt = au(1 - u).

72
1,2 -

1.0

.8

un .6

.4

.2

0

1E

"",..

..

• _._

;_.[.

• .I. _ _.._,

|

I I I I I I

.,..,

2.25 2.75 3.25 3.75

r = ¢rAt

4.25

Fig. 3.20 Bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for ihe

logistic ODE du/d_ = otu(1 - u) with u0 = 0.5.

73
1.2-

1.0

.8-

Un .6

.4-

2 "

1E

lS

":...

\

2/

0 I I I I

2.25 2.75 3.25

r = crAt

mj 4 " .

2,.. ,_I_

"\ 4 "i

\"\:/"- , .

4

I I I I

3.75 4.25

Fig. 3.21 Bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for the

logistic ODE dtt/dl = art(1 - _t) with multiple initial data.

?4
1.2

u n

1.0

.8-

.4-

.2-

0

1E

"FULL" BIFURCATION DIAGRAM

2

I I I I I I I

2.25 2.75 3.25 3.75

' r =a,lt

4.25

Fig. 3.22 "Full" bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for

the logistic ODE du/dt = otu(1 - u).

75
1.0

.8

Un

tP

.6

1E

._..-,

.4

.2

0

1.0

.8

.6

.4

.2

0

1E

u° " 0.08 _

u° - 0.51

1S

1S

u° " 0.49 i

1E '¢

1S

FULL

1S

1S

7 9 11 13 15 17 7 9 11 13

r - o_t r - o_t

15 17

Fig. 3.23 Bifurcation diagrams of the improved Euler (R-K 2) scheme for the ODE

du/clf = au(1 - _)(0.5 - u) for four different sets of initial input data.

76
1.0 u° = 0.25 :_ u° = 0.35

Un

0

1.0 u =0.49 ! t FULL I

.8 1_

.2

0

10 12 14 16 10 12 14 16

r = o_t r = o_t

Fig. 3.24 Bifurcation diagrams of the Runge-Kutta 4th-order (R-K 4) scheme for the

ODE du/dt = au(1 - u)(0.5 - u) for four different sets of initial input data.

77
.9

.7

U° - 0.19

1E 1E

.7

un J5

.3

.1

u° - 0.95

/

1E

1S 2 .

7 9 11 13 7

• -e_t

1E

FULL

2

9 11 13

r -o_t

Fig. 3.25 Bifurcation diagrams of the modified Euler (R-K 2) scheme for the ODE

du/dt = au(1 - u)(0.4 - _t) for four different sets of initial input data.

78
1.4

1.2

U n

1.0

.8-

.6-

.4

.2-

1E

1S

4

0 I I I I

1.5 2.0 2.5

r = o'At

4

!

3.0 3.5

Fig. 3.26 "Full" bifurcation diagram of the improved Euler (R-K 2) scheme for the

logistic ODE du/dt = au(1 - u).

79
U n

1.4 --

1.,2

1.0

m

.6-

1E

2

Fig. 3.27 "Full" bifurcation diagram of the Adam-Bashforth scheme for the logistic

ODE du/dt = au(1 - u).

8O
Ufl

1,4 --

1.2 --

1.0

.8

.6-

.4-

.2

0

1E

.. ;

\

\'_

4

'\

4

4

2

I I I I I I I I I I

1.5 2.0 2.5 3.0 3.5 4.0

r = oat

Fig. 3.28 "Full" bifurcation diagram of the predictor-corrector scheme of order 2 for

the logistic ODE du/d1 = au(1 -=).

81
1.4 --

1.0

.8-

U n

.6-

.4-

.2

1E

ls ., , ls

1

\

\

1S

\- lS

0 i I I I I I I I

1.5 2.0 2.5 3.0 3.5

r = _At

I

4.0

Fig. 3.29 "Full" bifurcation diagram of the predictor-corrector scheme of order 3 for

the logistic ODE du/dt = c_u(1 - _t).

82
Un

.8

.6

.4

.2

0

6

.9

.7

b=O.1

.8

/

:/ m6

1E _ "2

0

10 14 18 22

b = 0,3 .9

/

.7

b=0.2

6 10

b=0.4

14 18 22

un .5 .5

1E

.3 _ .3

.1 .1

7 9 11 13 15 7 9 11 13

r = odkt r - oc&t

Fig. 3.30 "Fun" bifurcation diagrams of the modified Euler (R-K 2) scheme for the

ODE d,L/dt_ = ctu(1 - _)(b - u), b :: 0.1, 0.2, 0.3, 0.4.

83
1.0 -

U n

.8-

.6-

.4-

0

1E

4

'..•

i I I i I I I I i

7 9 11 13 15

r = nat

Fig. 3.31 "Full" bifurcation diagrams of the explict Euler scheme for the ODE du/dt =

_,_(1- ,,)(0.5- ,,).

84
1.0 -

U n

.8

.6

.4

.2

0

1E

2

4

2

I I I I I I

7 8 9 10 11 12 13

b

r = aht

14

Fig. 3.32 "Full" bifurcation diagram of the modified Euler (R-K 2) scheme for the

ODE d,_/dt = au(1 - u)(0.5- u).

85
un

1.0

.8

.6

.4

.2

0

1E

4

../"

e

/

!

¢

%

4

i

i I I I I I I I

3 4 5 6 7

r = aAt

I

8

Fig. 3.33 "Full" bifurcation diagram of the Adam-Bashforth scheme for the ODE

du/dt = au(1 - u)(0.5 - u).

86
1.0 -

.8

.6

Un

.4-

0

7

1E

4

4 2

I I I I I I I

9 11 13 15

r = aAt

I

17'

Fig. 3.34 "Full" bifurcation diagram of the improved Euler (R-K 2) scheme for the

ODE du/dt = au(1 - u)(0.5- u).

87
Un

1.0

.8

.6

.4

.2

0 I I I ! I

10 12 14

1S 2

r =aA!

2

!

I

16

Fig. 3.35 "Full" bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for

the ODE du/dt = art(1 - u)(0.5 - _).

88
Un

1S

2

4

4

4

2

2

4 2

4 2

4

0 I I I I I I

13.5 14.5 15.5 16.5

r = a,At

Fig. 3.36 "Full" bifurcation diagram of the improved Euler (R-K 2) scheme for the

0 D E du / dt = au( 1 - _)(0.5 - it) (enlarged).

89
U n

1,0 -

IS g

.6

2

2

.4-

0

13.5

2

is

2

2

1

I I I I I _ I

14.5 15.5 16.5

r :_rAt

Fig. 3.37 "Full" bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for

the ODE d_t/dt = au(1 - _)(0.5- it) (enlarged).

9O
U n

.5

.4

.3-

.2

mR --

1E

\

\

\

\

\

2S

2

0 I I I I

6 8 10 12 14 16

r = _At

I

18

Fig. 3.38 "Full" bifurcation diagram of the predictor-corrector scheme of order 2 for

the ODE du/dt = au(1 - u)(0.5 - _).

91
i| _ EI

--I

°l--q

C_

_D

C_

f_o

°pl

92
_, ._ _"_

___'_

._., .__

_ l -_ ._ _

93
U

%

!

r

/

I

I

I

I

!

SUBCRITICAL

I

I

I

I

TRANSCRITICAL

a (SYSTEM PARAMETER)

Fig. 3.41 Stability of steady-state solutions arising through three types of bifurcation

phenomena (-- stable, - - - unstable).

94
o

°_

C_

f_

°l-q

c_ o

-_

c_

°p-qo

cI) 0

C_

°p.q

95
0

f_

_p.q

o_

T_

o_

c_

_o_

p_

c_O

o_

c_

ol.lq

96
q.

_o

.JD

M

Ill

vi11

0

<1

II

&.

"<1

II

s._

°l-.

0

,=

0

0

0

._

_._

o

_b

97
0m

.Dm

o)

qw,

0!!0

Q0

lu

III

f

i

ol w ! /

lu I

-j U3

II!

C

!

II

L._

O

_d

0

0,.._o

_

oo

o_._ o_

_O

.M

98
_o

JD

m4mJ

(n /j

i | I I

fur#

u"p

a,_,.

I",,,,

o

t

I

,.,..,_

(_1

In

(_1

t,.

(_ II

In

r,-

In

¢_1

o©

u_

o_

C

99
,q.

\

/

,,//

,n.

o.

II

O

o

f-i

o

0

o _

o _

o O

o

,_ 0

_,_

o

100
.Q

tG

(/)

i,-

v,-

f--"

! !

I

C

o,

_r

o.

t_

II

o,

0

0

0

o_

0 _

,.,_

° 0

°_

o._

0

,._ 0

cO 0

_b

101
.D

I0

G_ o

.12m

sn

e=

"OC

88

J_

!

o.

!

w

\/

III

e-

\|

II

II

L_

II

L.

0

o

_o m

_b

o_,,i

102
_Q

.D

M

¢h

o

m

c

"o

a

m

_j,\J, ,_J. j

¢I)

I I I !

lu

I i I I

u

!

(D

ul

e,

t I

O!

-\

j_e_ .,°

I I I

0

II

II

L_

o I

o _

_._

0

II

o_

o "_

t:z,_

_O

_b

°_

103
cq /c_ _lc_

I I I I ! I I I

0

.a

W

4_

W

"D

W

0

mm

.D

m

!

o,

_m

N

i

104

UC_

T_

II

O _

o _

" _ II

II ro_ ._..

0

0

_b
o

dD

al

4IJ

w

I

.Q

ca

C

C

cO

_Q

dD

m

P\ . . /a

VI

-. ............ J

I

, , , , _] , , , ,

I I !

tl IV

!

/

/-I

I

O

II

L_

0

_D

II

i._

_1_

_D

Q

0

0

o_._

,.=

"El "-_

H

o

_ o

105
0

m

,0

lu

• • _ _. /" " "_ V \J_ ,_ ,_ ""_"

I I I l "1I I I I

i

f!

i

0

i

! !

e=

I

II

li,..

O

II

l,J

O

,1,J

t.D

O

0

t_

0

t=)

o

osiq

0

0

m

°_,,_

106
o_

JDm

Ino

J:

oBm

"oc

f_ o

.D

_mJ

(n

_n

IIIIIl

uJ

IIII!I

e,

-C_l

oo

ii

im

0

,Ira,

(D

@4

¢0

N"

II

Irm

O

Vm

CD

0

°_-_o

,.c:; o

rJ_ o

U2

"" II

,t,,qo ,._

_0

°_

107

AI/LS/X Natiomd Atrcnaulk_

Sf_cQ Admrr_nmJon

1. Report No. I

NASA TM- 102820 I

4. Title and Subtitle

Report Documentation Page

2. Government Accession No.

Dynamical Approach Study of Spurious Steady-State Numerical

Solutions of Nonlinear Differential Equations

Part 1 - The ODE Connection and Its Implications for Algorithm

Development in Computational Fluid Dynamics

3. Recipient's Catalog No.

5. Report Date

7. Author(s)

-EL C. Y___K. Sweby (University of Reading, Whiteknights,

,and), and D. F. Griffiths (University of Dundee,

,and)

ration Name and Address

;h Center

, CA 94035-1000

cy Name and Address

mautics and Space Administration

DC 20546-0001

April 1990

6. Performing Organization Code

8. Performing Organization Report No.

A-90149

10. Work Unit No.

505-60

11. Contract or Grant No.

13. Type of Report and Period Covered

Technical Memorandum

14. Sponsonng Agency Code

|otes

tact: H.C. Yee, Ames Research Center, MS 202A-1, Moffett Field, CA 94035-1000

(415) 604-4769 or FTS 464-4769

tted at 12th International Conference on Numerical Methods in Fluid Dynamics, Oxford,

y 9-13, 1990.

_table as well as unstable steady-state numerical solutions, spurious asymptotic numerical solutions of higher

en stable chaotic behavior can occur when finite-difference methods are used to solve nonlinear differential

tnumerically. The occurrence of spurious asymptotes is independent of whether the DE possesses a unique steady

ditional periodic solutions and/or exhibits chaotic phenomena. The form of the nonlinear DEs and the type of

_'mes are the determining factor. In addition, the occurrence of spurious steady states is not restricted to the time

_eyond the linearized stability limit of the scheme. In many instances, it can occur below the linearized stability

re, it is essential for practitioners in computational sciences to be knowledgeable about the dynamical behavior

ence methods for nonlinear scalar DEs before the actual application of these methods to practical computations.

rtant to change the traditional way of thinking and practices when dealing with genuinely nonlinear problems.

t, spurious asymptotes were observed in numerical computations but tended to be ignored because they all were

beyond the linearized stability limits of the time step parameter At. As can be seen from our study, bifurcations

_urious asymptotic solutions and transitions to computational instability not only are highly scheme dependent

lependent, but also initial data and boundary condition dependent, and not limited to time steps that are beyond

stability limit

.. Jggosted by Author(s))

Nonlinear dynamics, Chaotic dynamics, Dynamics of

numerics, Computational fluid dynamics, Hypersonic

reacting flows, Nonlinear ordinary differential

equations, Nonlinear dynamical systems

19. Security Classif. (of this report)

Unclassified

20. Security Classif. (of this page)

Unclassified

18. Distribution Statement

Unclassified-Unlimited

Subject Category - 64

21. No. of Pages

110

22. Price

A06

NASA FORM 1626 OCT86

For sale by the National Technical Information Service, Springfield, Virginia 22161

NASA-TM-11Z?09 _:

Comp. Fluid Dyn., t995, Vol. 4, p. 219-283

Reprints available directly from the publisher

Photocopying permitted by license only

(C, 1995 OPA (Overseas Publishers Association)

Amsterdam B.V. Published under license by

Gordon and Breach Science Publishers SA

Printed in Malaysia

DYNAMICAL APPROACH STUDY OF SPURIOUS

STEADY-STATE NUMERICAL SOLUTIONS

OF NONLINEAR DIFFERENTIAL EQUATIONS II.

GLOBAL ASYMPTOTIC BEHAVIOR OF TIME

DISCRETIZA TIONS*

H. C. YEE

Fluid Dynamics Division, NASA Ames Research Center,

Moffett Field, CA, 94035, USA

P. K. SWEBY *

Department of Mathematics, University of Reading,

Whiteknights, Reading RG6 2AX, England

(Received 19 August 1993; in final form tO February 1994)

SUMMARY

The global asymptotic nonlinear behavior of 11 explicit and implicit time discretizations for four 2 × 2

systems of first-order autonomous nonlinear ordinary differential equations (ODEs) is analyzed. The

objectives are to gain a basic understanding of the difference in the dynamics of numerics between the scalars

and systems of nonlinear autonomous ODEs and to set a baseline global asymptotic solution behavior of

these schemes for practical computations in computational fluid dynamics. We show how "numerical" basins

of attraction can complement the bifurcation diagrams in gaining more detailed global asymptotic behavior

of time discretizations for nonlinear differential equations (DEs). We show how in the presence of spurious

asymptotes the basins of the true stable steady states can be segmented by the basins of the spurious stable

and unstable asymptotes. One major consequence of this phenomenon which is not commonly known is that

this spurious behavior can result in a dramatic distortion and, in most cases, a dramatic shrinkage and

segmentation of the basin of attraction of the true solution for finite time steps. Such distortion, shrinkage

and segmentation of the numerical basins of attraction will occur regardless of the stability of the spurious

asymptotes, and will occur for unconditionally stable implicit linear multistep methods. In other words, for

the same (common) steady-state solution the associated basin of attraction of the DE might be very different

from the discretized counterparts and the numerical basin of attraction can be very different from numerical

method to numerical method. The results can be used as an explanation for possible causes of error, and slow

convergence and nonconvergence of steady-state numerical solutions when using the time-dependent

approach for nonlinear hyperbolic or parabolic PDEs.

KEY WORDS: Spurious steady-state numerical solutions, spurious asymptotes, global asymptotic

behavior, nonlinear ODEs, numerical methods, time discretizations.

* Part of the material appeared in the Proceedings of the 9th GAMM Conference on Numerical Methods

in Fluid Mechanics, Lausanne, Switzerland, Sept. 25 27, 1991. Full text appeared as a NAS Applied

Research Technical Report RNR-92-008, March 1992, NASA Ames Research Center.

t Part of this work was performed as a visiting scientist at the NASA Ames Research Center.

219
220 H.C. YEE AND P. K. SWEBY

1. INTRODUCTION

The tool that is utilized for the current study belongs to a multidisciplinary field of

study in numerical analysis, sometimes referred to as "The Dynamics of Numerics _''.

Here the phrase "to study the dynamics of numerics" {dynamical behavior of a numerical scheme) is restricted to the study of local and global asymptotic behavior and

bifurcation phenomena of the nonlinear difference equations resulting from finite

discretizations of a nonlinear differential equation (DE) subject to the variation of

discretized parameters such as the time step, grid spacing, numerical dissipation

coefficient, etc. In this paper, standard terminologies of nonlinear dynamics, chaotic

dynamics (Guckenheimer and Holmes, 1983; Hale and Kocak, 1991) and computational fluid dynamics (CFD) are assumed. For an introduction to the dynamics of

numerics and its implications for algorithm development in CFD, see Yee et al. (1991}

and Yee (1991) and references cited therein.

1.1 Backyround

The phenomenon that a nonlinear DE and its discretized counterpart can have

different dynamical behavior (asymptotic behavior} was not uncovered fully until

recently. Aside from truncation error and machine round-offerror, a more fundamental

distinction between the DE {continuum) and its discretized counterparts for genuinely

nonlinear behavior is extra solutions in the form of spurious stable and unstable

asymptotes that can be created by the numerical method. Here we use the term

"discretized conterparts" to mean the finite difference equations (or discrete maps)

resulting from finite discretizations of the underlying DE. Also we use the term "spurious

asymptotic numerical solutions" to mean asymptotic solutions that satisfy the discretized counterparts but do not satisfy the underlying ordinary differential equations

{ODEs} or partial differential equations (PDEs). Asymptotic solutions here include

steady-state solutions (fixed points of period one for the discretized equations), periodic

solutions, limit cycles, chaos and strange attractors. See Section 1II and Guckenheimer

and Holmes (1983), Hale and Kocak (1991) and Yee et al. (1991) for definitions.

lserles (1988) showed that while linear multistep methods (LM Ms) for solving ODEs

possess only the fixed points (fixed points of period one) of the original DEs, popular

Runge-Kutta methods may exhibit additional, spurious fixed points. It has been

demonstrated by the authors and collaborators (Yee et al., 1991; Yee, t991; Sweby

et al., 1990; Griffiths et al., 1992; Yee and Sweby, 1993a, 1993b) for nonlinear ODEs,

and Lafon and Yee {1991, 1992) for nonlinear reaction-convection model equations that

such spurious fixed points as well as spurious fixed points of higher periods may be

stable below the linearized stability limit of the scheme, depending on the initial data.

lserles et al. (1990), Hairer et al. (1989} and Humphries {1991} further advanced some

theoretical understanding of the dynamics of numerics for initial value problems of

ODEs. lserles et al. and Hairer et al. classified and gave guidelines and theory on the

types of Runge-Kutta methods that do not exhibit spurious period one or period two

fixed points. H umphries (1991) showed that under appropriate assumptions if stable

_Named after the First IMA Conference on Dynamics of Numerics and Numerics of Dynamics,

University of Bristol, England, July 31 August 2, 1990.
DYNAMICAL STUDYOFSPURIOUSSTEADY-STATE NUMERICAL SOLUTIONS221

spurious fixed points exist as the time-step approaches zero, then they must either

approach a true fixed point or become unbounded. However, convergence in practical

calculations involves a finite time step At as the number of integrations n _ oo rather

than At _ 0, as n _ oo. There appear to be missing links between theoretical development and practical scientific computation. Our aim is to provide some of these missing

links that were not addressed in Iserles (1988), Iserles et al. (1990), Hairer et al. (1989),

Humphries (1991) and our earlier work. In particular, we want to show in more detail

the global asymptotic behavior of time discretizations when finite but not extremely

small At is used. Other aspects that were not addressed in Iserles (1988) for different

iteration procedures in solving the resulting nonlinear algebraic equations are reported

in greater depth in our companion papers (Yee and Sweby, 1993a, 1993b).

1.2 Relevance and Motivations

Although the understanding of the dynamics of numerics of systems of nonlinear ODEs

and PDEs is important in its own right and has applications in the various nonlinear

scientific fields, our main emphasis is CFD applications. Time-marching types of

methods (time-dependent approach) are commonly used in CFD because the steady

PDEs of higher than one dimension are usually of the mixed type. When a timedependent approach is used to obtain steady-state numerical solutions of a fluid flow

or a steady PDE, a boundary value problem (BVP) is transformed into an initialboundary value problem (IBVP) with unknown initial data. If the steady PDE is

strongly nonlinear and/or contains stiff nonlinear source terms, phenomena such as

slow convergence, nonconvergence or spurious steady-state numerical solutions and

limit cycles commonly occur even though the time step is well below the linearized

stability limit and the initial data are physically relevant. One of our goals is to search

for logical explanations for these phenomena via the study of the dynamics of numerics.

Here the term "time-dependent approach" is used loosely to include some of the

iteration procedures (due to implicit time discretizations), relaxation procedures, and

preconditioners for convergence acceleration strategies used to numerically solve

steady PDEs. This is due to the fact that most of these procedures can be viewed as

approximations of time-dependent PDEs (but not necessarily the original PDE that

was under consideration). If one is not careful, numerical solutions other than the

desired one of the underlying PDE can be obtained (in addition to spurious asymptotes

due to the numerics).

One consequence of the existence of stable and unstable spurious asymptotes below

or above the linearized stability limit of the numerical schemes is that these spurious

features may greatly affect the dynamical behavior of the numerical solution in practice

due to the use of a finite time step. As discussed in details in later sections and also in

Yee et al. (1991), Yee and Sweby (1993a, b), Lafon and Yee (1992), Sweby and Yee

(1991), Yee et al. (1992), it is possible that for the same steady-state solution, the

associated basin of attraction of the underlying DEs (which initial conditions lead to

which asymptotic states) might be very different from that of the basin of attraction of

the discretized counterparts due to the existence of spurious stable and unstable

asymptotic numerical solutions. In other words, there is a separate dependence on

initial data for the individual DEs and their discretized counterparts. Here the basin of

attraction is a domain of a set of initial conditions whose solution curves (trajectories)
222 H.C. YEE AND P. K. SWEBY

all approach the same asymptotic state. Also we use the term "exact" and "numerical"

basins of attraction to distinguish "basins of attraction of the underlying DEs" and

"basins of attraction of the discretized counterparts".

In view of the spurious dynamics, it is possible that numerical computations may

converge to an incorrect steady state or other asymptote which appears to be physically

reasonable. One major implication is that what is expected to be physical initial data

associated with the underlying steady state of the DE might lead to a wrong steady

state, a spurious asymptote, or a divergence or nonconvergence of the numerical

solution. In addition, the existence of spurious limit cycles may result in the type of

nonconvergence of steady-state numerical solutions observed in time-dependent

approaches to the steady states. It is our belief that the understanding of the symbiotic

relationship between the strong dependence on initial data and permissibility of

spurious stable and unstable asymptotic numerical solutions at the fundamental level

can guide the tuning of the numerical parameters and the proper and/or efficient usage

of numerical algorithms in a more systematic fashion. It can also explain why certain

schemes behave nonlinearly in one way but not another. Here strong dependence on

initial data means that for a finite time step At that is not sufficiently small, the

asymptotic numerical solutions and the associated numerical basins of attraction

depend continuously on the initial data. Unlike nonlinear problems, the associated

numerical basins of attraction of linear problems are independent of At as long as At is

below a certain upper bound.

Nonunique Steady-State Solutions of Nonlinear DEs vs. Spurious Asymptotes: The phenomenon of generating spurious steady-state numerical solutions (or other spurious

asymptotes) by certain numerical schemes is often confused with the nonuniqueness (or

multiple steady states} of the DE. In fact, the existence of nonunique steady-state

solutions of the continuum can complicate the numerics tremendously (e.g., the basins

of attraction) and is independent of the occurrence of spurious asymptotes of the

associated scheme. But, of course, a solid background in the theory of nonlinear ODEs

and PDEs and their dynamical behavior is a prerequisite in the study of the dynamics of

numerics for nonlinear PDEs. See Yee et al., 1991 for a discussion. It is noted that the

approach and primary goal of our work is quite different from the work of e.g., Beam

and Bailey (1988) and Jameson (1991). The main goal of Beam and Bailey (1988) and

Jameson (1991) was to study the non unique steady-state solutions admitted by the PDE

as the physical parameter is varied. Our primary interest is to establish some working

tools and guidelines to help delineate the true physics from numerical artifacts via the

dynamics of numerics approach. The knowledge gained from our series of studies (Yee

et al., !991; Lafon and Yee, 1991; Lafon and Yee, 1992) hopefully can shed some light on

the controversy about the existence of multiple steady-state solutions through numerical experiments for certain flow types of the Euler and/or Navier Stokes equations.

1.3 Objectives and Outline

The primary goal of the series of papers (present and the companion papers Yee et al.,

1991: Yee and Sweby, 1993a, b; Lafon and Yee, 1991; Lafon and Yee, 1992) is to lay the

foundation for the utilization of the dynamics of numerics in algorithm development for

computational sciences in general and CFD in particular. This is part II of this series of

papers on the same topic. Part 1(Yee et al., 1991} concentrated on the dynamical behavior
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 223

of time discretizations for scalar nonlinear ODEs. The intent of part I was to serve

as an introduction to motivate this concept to researchers in the field of CFD and to

present new results for the dynamics of numerics for first-order scalar autonomous

ODEs.

The present paper, the second of this series, is devoted to the study of the dynamics

of numerics for 2 x 2 systems of ODEs. Here we show how "numerical" basins of

attraction can complement the bifurcation diagrams in gaining more detailed global

asymptotic behavior of numerical methods for nonlinear DEs. We show how in the

presence of spurious asymptotes the basins of the true stable steady states can be

segmented by the basin of the spurious stable and unstable asymptotes. One major

consequence of this phenomenon which is not commonly known is that this spurious

behavior can result in a dramatic distortion and, in most cases, a dramatic shrinkage

and segmentation of the basin of attraction of the true solution for finite time steps.

Such distortion, shrinkage and segmentation of the numerical basins of attraction

will occur regardless of the stability of the spurious asymptotes, and will occur for

unconditionally stable implicit linear multistep methods. In other words, for the same

steady-state solution, the associated basin of the DE might be very different from its

discretized counterparts. The basins can also be very different from numerical method

to numerical method. The present study reveals for the first time the detail interlocking relationship of numerical basins of attraction and the causes of error, and slow

convergence and nonconvergence of steady-state numerical solutions when using the

time-dependent approach.

The article of Lafon and Yee (1991), the third of this series, was devoted to the study

of the dynamics of numerics of commonly used numerical schemes in CFD for a model

reaction-convection equation. The article of Lafon and Yee (1992), the fourth of

this series, was devoted to a more detailed study of the effect of numerical treatment

of nonlinear source terms on nonlinear stability of steady-state numerical solution

for the same model nonlinear reaction-convection BVP. In our companion papers

(Sweby et al., 1990; Griffiths et al., 1991a, 1992b), a theoretical bifurcation analysis of

a class of explicit Runge-Kutta methods and spurious discrete travelling wave phenomenon were presented. In yet another companion paper, Yee and Sweby (1993a), the global

asymptotic nonlinear behavior of three standard iterative procedures in solving

nonlinear systems of algebraic equations arising from four implicit LM Ms is analyzed

numerically.

1.4 Outline

The outline of this paper is as follows. Section II discusses the connection of the

dynamics of numerics for systems of ODEs and numerical approximations of timedependent PDEs. Section Ill reviews background material for nonlinear ODEs and

their numerical methods. Section IV describes four 2 × 2 systems of nonlinear

first-order autonomous model ODEs. Section V describes the 11 time discretizations

and the associated bifurcation diagrams for the four model ODEs. Section VI discusses

the combined basins of attraction and bifurcation diagr'fims for the underlying

schemes. Comparison between a linearized implicit Euler and Newton method is

briefly discussed in Section 6.5. The paper ends with some concluding remarks in

Section VII.
224 H.C. YEE AND P. K. SWEBY

2. THE DYNAMICS OF NUMERICS OF SYSTEMS OF ODEs AND

NUMERICAL APPROXIMATIONS OF TIME-DEPENDENT PDEs

For finite discretizations of PDEs, spurious asymptotes and especially spatiallyvarying spurious steady states can be independently introduced by time and spatial

discretizations (Yee et al., 1991; Lafon and Yee, 1991; Lafon and Yee, 1992). The interaction

between temporal and spatial dynamical behavior is more complicated when one is

dealing with the nonseparable temporal and spatial finite-difference discretizations

such as the Lax-Wendrofftype. The analysis and the different features of the numerics

due to temporal and spatial discretizations can become more apparent by separable

temporal and spatial finite difference methods (FDM). A standard method for obtaining such a FDM is the method of lines (MOL) procedure where the time-dependent

PDE is reduced to a system of ODEs (by replacing the spatial derivatives by finite difference approximations). The resulting approximation is called semi-discrete, since the

time variable is left continuous. The semi-discrete system in turn can be solved by the

desired time discretizations. Similar semi-discrete systems can be obtained by finite

element methods except in this case an additional mass matrix is involved. Besides the

MOL approach, coupled nonlinear ODEs can arise in many other ways when analyzing

nonlinear PDEs. See for example Globus et al. (1991), Hung et al. (1991), Foias et al.

(1985), Temam (1989), Kwak (1991), Schecter and Shearer (1990), and Shearer et al.

(1987). Among these possibilities, the idea of inertial manifold tiM) and approximate

inertial manifold (AIM) for incompressible Navier-Stokes (Foias et al., 1985; Temam,

1989; Kwak, 1991), the relationship between shock waves, heteroclinic orbits of systems

of ODEs (Schecter and Shearer, 1990; Shearer et al., 1987), and flow visualization of

numerical data (Globus et al., 1991; Hung, 1991) are touched upon here.

2.1 Asymptotic Analysis of the Method of Lines Approach

When the ODEs are obtained from a semi-discrete approximations of PDEs, the

resulting system of ODEs contains additional system parameters and discretized

parameters as opposed to physical problems governed by ODEs. Depending on the

number of grid points "J" used, the dimensions of the resulting system of semi-discrete

approximations of ODEs can be very large. Also, depending on the differencing scheme

the resulting discretized counterparts ofa PDE can be nonlinear in At, the grid spacing

Ax and the numerical dissipation parameters, even though the DEs consist of only one

parameter or none. One major consideration is that one might be able to choose

a "safe" numerical method to solve the resulting system of ODEs to avoid spurious

stable steady states due to time discretizations. However, spurious steady states and

especially spatially varying steady states introduced by spatial discretizations in

nonlinear hyperbolic and parabolic PDEs for CFD applications appear to be more

difficult to avoid. In the case of the MOL approach, if spurious steady states due to

spatial discretizations exist, the resulting ODE system has already inherited this

spurious feature as part of the exact solutions of the semi-discrete case. We remark that

spurious stable and unstable asymptotes other than the steady states due to time

discretizations are also more difficult to avoid than spurious steady states. See Sections

V and VI for some illustrations. Taking for example the nonlinear ODE models that are
DYNAMICAL STUDYOF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 225

considered, it is relatively easy to avoid spurious steady states due to time discretizations since, ira numerical steady state U* for the ODE dU/dt = S(U) is spurious, then

S(U*) v_O. This is not the case for spurious asymptotes such as limit cycles.

In addition to the aforementioned considerations, it is well known from the theory

of nonlinear dynamics for ODEs that much of tho established theory and known

behavior of nonlinear dynamics are restricted to lower dimensional first-order ODEs

(or for problems that exhibit lower dimensional dynamical behavior). Moreover, if

higher than two-time level numerical methods are used, the dynamics of these

discretized counterparts usually are richer in structure and more complicated to

analyze than their two-time level cousins. Therefore, in order to gain a first hand

understanding of the subject we restrict our study to 2 x 2 systems of first-order

autonomous ODEs and two-time level numerical methods with a fixed time step, even

though the current study is far removed from the realistic setting. Studies of 3 x 3

systems and general J x J systems are in progress.

Due to the complexity of the subject matter, this paper concerns fixed time step (and

fixed grid spacing) time-marching methods only. The fixed or local variable time step

control method study can also shed some light on identifying whether certain flow

patterns are steady or unsteady. See Yee et al. (1990) for some examples. Proper

regulation of a variable time step to prevent the occurrence of spurious steady-state

numerical solutions will be a subject of future research. In order to isolate the different

causes and cures of slow convergence and nonconvergence of time-marching methods,

our study concerns nonlinearity and stiffness that are introduced by DEs containing

smooth solutions. Nonlinearity and stiffness that are introduced by the scheme, the

coupling effect in the presence of a source term (terms) in coupled system of PDEs, the

highly stretched nonuniform structured and unstructured grids, the discontinuities in

grid interfaces and/or the discontinuities inherent in the solutions, and external flows

that need special boundary conditon treatment with a truncated finite computation

domain are added factors and require additional treatment or different analysis. These

are not considered at the moment. Generalization of our study to include grid adaption

as one of the sources of nonlinearity and/or stiffness introduced by the numerics is

reported in Sweby and Yee (1964)and Buddet al. (1994).

2.2 Inertial Manifold ( IM ) and Approximate Inertial Manifold (AIM)

The concept of IMs was introduced by Foias et al. (1985). See Foias et al. (1985),

Temam (1989) and Kwak (1991) for details of the subject. The key idea of IMs and

AIMs is to establish theories to aid in better understanding of nonlinear phenomena

and turbulence via the study of the interaction of short and long wavelengths of

dissipative systems. Basically, an IM is a finite-dimensional submanifold that contains

all the attractors and invariant sets of an infinite-dimensional dynamical system

described by some dissipative PDEs. It establishes the criterion for the reduction of

long-term dynamics of certain infinite-dimensional problems to a finite system of

ODEs. An attractive feature is that the reduction introduces no error in the problem.

That is, the IM contains all pertinent information about the long-term dynamics of the

original system. One of the main objectives of AIMs is to handle cases where the IMs

are not known to exist. AIMs also can help in finding good algorithms for dealing with
226 H.C. YEE AND P. K. SWEBY

the IMs that are known to exist. AIMs may also help reduce finite but extremely large

systems of ODEs to lower-dimensional problems. In a nut shell, the derivation of IMs

and AIMs is based on the decomposition of the unknown function into large-scale and

small scale components. In the case of fluid dynamics, those structures can be identified

as large and small eddies. Thus an IM or AIM corresponds to an exact or approximate

interaction law between the short and long wavelengths. Kwak (1991) showed that the

long-term dynamics of some two-dimensional incompressible Navier-Stokes equations

can be completely described by a finite system of ODEs. Kwak does so by finding

a nonlinear change of variable that embeds the incompressible Navier-Stokes equations in a system of reaction-diffusion equations that possess an IM. All of the theories

oflMs and AIMs are very involved and interested readers are encouraged to read Foias

et al. (1985), Temam (1989) and Kwak (1991) and the references cited therein.

2.3 Relationship Between Shock Waves and Heteroclinic Orbits of Systems of ODEs

Another example of the importance of understanding the "dynamics" and the "dynamics of numerics" of systems of ODEs is related to the study of shocks using

equilibrium bifurcation diagrams of associated vector fields. This was introduced by

Shearer et al. (1987). The authors find of great interest how one can reduce the study of

admissible shock wave solutions ofa 2 × 2 hyperbolic conservation laws to the study of

heteroclinic orbits of a system of nonlinear ODEs. Further development in this area

can help in constructing suitable approximate Riemann solvers in numerical computations. Schecter and Shearer (1990) studied undercompressive shocks for nonstrictly

hyperbolic conservation laws by adding information to the equilibrium bifurcation

diagrams (introduced by Shearer et al.) about heteroclinic orbits of the vector fields.

The augmented equilibrium bifurcation diagrams are then used in the construction of

solutions of Riemann problems.

2.4 Dynamics of Numerics and Flow Visualizations of Numerical Data

The use of flow visualization of numerical data (numerical solutions of finite discretizations of e.g., fluid flow problems) in an attempt to understand the true flow physics has

become increasingly popular in the last decade. See, Globus et al. (1991) and Hung et al.

(1991) and references cited therein. Many of the techniques rely on the extraction of

the boundary surfaces by analyzing a set of appropriate vector fields. Approximations

are then performed based on this set of vector fields. The study of the topological

features of certain flow physics based on the numerical data is then related to the

study of fixed points of the associated systems of ODEs. Fluid problems with known

flow physics can be used to reveal how well the associated vector fields of the

numerical data can mimic the true physics. It can also help to delineate spurious

flow patterns that are solely due to the numerics. At the present time we are entering

into the regime where CFD is extensively used to aid the understanding of complicated

flow physics that is not amenable to analysis otherwise. In the situation where the

numerical data indicate flow structures which are not easily understood, a good

understanding of the spurious dynamics that can be introduced by the numerics is

needed.
DYNAMICAL STUDYOFSPURIOUSSTEADY-STATE NUMERICAL SOLUTIONS227

3.PRELIMINARIES

Considera2x 2systemoffirst-orderautonomousnonlinearODEsof theform

dU

-_T = S(U), (3.1)

where U and S are vector functions of dimension 2, and S(U) is nonlinear in U.

A fixed point UE of an autonomous system (3.1) is a constant solution of (3.1);

that is

S(U_) = 0, (3.2)

where the subscript "E" stands for "exact" and UE denotes the fixed points of the ODE

as opposed to the additional fixed points of the discretized counterparts (spurious fixed

points) due to the numerical methods which we will encounter later.

Let the eigenvalues of J(UF) = (OS/OU)Iv_ (the Jacobian matrix of S(U) evaluated

at Ue) be 21 and 3.2. Here J(UE) is assumed to be nonzero. The fixed point U E is

hyperbolic if Re(2i) 4: 0, i = 1, 2. If both 2i are real, U Eis a saddle if2122 < 0 and a node if

2122 > 0. If exactly one 2i = 0, then U E is semihyperbolic. If the eigenvalues are

complex, then U e is a spiral. The "tightness" of the spiral is governed by the magnitude

of the imaginary part of the eigenvalues. If the eigenvalues both have a zero real part,

then U_ is non-hyperbolic. Such a fixed point is called a center. Under this situation,

more analysis is needed to uncover the real behavior of(3.1) around a non-hyperbolic

fixed point. The fixed point U E is stable if both 21 and 22 have negative real parts. U_ is

unstable if a 2_ has a positive real part. In the non-hyperbolic case the fixed point is

neutral.

If due to a variation of a parameter of the ODE a fixed point becomes unstable, then,

if at the point of instability the eigenvalues are distinct and real, the resulting bifurcation

will be to another fixed point. Such bifurcation is called a steady bifurcation. If, however,

the eigenvalues are complex, then the bifurcation will be of a Hopf type. This is a

slightly simplified classification, since our main concern in this work is not on the

variation of the ODE parameter. Detailed background information can be found in

(Guckenheimer and Holmes, 1983; Hale and Kocak, 1991).

Consider a nonlinear discrete map from a finite discretization of (3.1)

U,+ 1 = U" + D(U", r), (3.3)

where r=At and D(U",r) is linear or nonlinear in r depending on the numerical

method. A fixed point U o of(3.3) is defined by U "+ 1 = U", or

U D = U o + D(UD, r) (3.4)

or D(Uo, r)=O. A fixed point U D of period p>0 of (3.3) is defined by U"+P= U"

with U "+k 4: U" of k < p. In the context of discrete systems, the term "fixed point"

without indicating the period means "fixed point of period 1" or the steady-state

solution of (3.3). Here we use the term asymptote to mean a fixed point of any period,

a limit cycle (in the discrete sense invariant set), chaos, or a strange attractor.
228 H.C. YEE AND P. K. SWEBY

The type of finite discretization of (3.1) represented in (3.3) assumed the use of

two-time level schemes. Otherwise the vector dimension of (3.3) would be 2(k - l)

instead of 2 where k is the number of the time level of the scheme. Here the vector

function D is assumed to be consistent with the ODE (3.1) in the sense that fixed points

of the ODE are fixed points of the scheme; however, the reverse need not hold. It is this

feature accompanied by other added dynamics, that the discretized counterparts of the

underlying ODE possess a much richer dynamical behavior than the original ODE.

Thus the fixed points U D of D(Uo, r) = 0 may be true fixed points Ue of(3. l) or spurious

fixed points U s. The spurious fixed points U s are not roots of S(U)= O. That is

S(U,) _ O. Spurious asymptotes are asymptotic numerical solutions of (3.3) but not

(3.1).

Letting U" = U o + 6_, then a perturbation analysis on (3.3) yields

6"+1= I+ OD(UD'r)') -_- j "+1 6.o

(3.5)

Assuming OD(Uo, r)/dU _ 0, then the fixed point U o is stable if the eigenvalues of

Jo = I + dD(U a, r)/dU lie inside the unit circle. If both eigenvalues are real and both lie

inside (outside) the unit circle, then the fixed point is a stable (unstable) node. If one is

inside the unit circle and the other outside, then the fixed point is a saddle. If both

eigenvalues are complex, then the fixed point is a spiral. If the eigenvalues lie on the unit

circle, then the fixed point of(3.3) is indeterminant and additional analysis is required to

determine the true behavior of(3.3) around this type of fixed point. For a more refined

definition and the difference in fixed point definition between ODEs and discrete maps,

see Panov et al. (1956), Perron (1929) and Hsu (1987) and references cited therein. The

reader is referred to Guckenheimer and Holmes (1983), Hale and Kocak (1991),

Langford and Iooss (1980), and Werner (1980) for full details on the subject of

bifurcation theory.

An important feature which can arise (for both systems of ODEs (3.1) and their

discretizations) as the result of a Hopf bifurcation is a limit cycle where the trajectory

traverses a closed curve in phase space. In all but a few simple cases such limit cycles are

beyond analysis.

4. MODEL 2 x 2 SYSTEMS OF NONLINEAR FIRST-ORDER

AUTONOMOUS ODEs

Four 2 x 2 systems of nonlinear first-order autonomous model ODEs are considered.

The systems considered with U r = (u, v) or z = u + iv are a

1. Dissipative complex model:

dz

_= z(i + e- lzl 2) (4.1)
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 229

2. Damped Pendulum model:

du

dt v (4.2a)

d_

--= -ev-sin(u) (4.2b) dt

3. Predator-Prey model:

du

dt

dv

dt

--- - 3u +4u z-O.5uv- us

(4.3a)

-- - - 2.Iv + uv (4.3b)

4. Perturbed Hamiltonian System model:

du

dt - e(1 - 3u) + 3 [1 - 2u + u2 - 2v(1 - u)] (4.4a)

dv

d_ = e(1 - 3v) - 3 [1 - 2v + v2 - 2v(1 - v)] (4.4b)

Here e is the system parameter for (4.1), (4.2) and (4.4).

The perturbed Hamiltonian model can be related to the numerical solution of the

viscous Burgers' equation with no source term

c3u 1 O(u 2) c_2u

Ot +2 c3x fl_x 2 fl>0. (4.5)

Let uj(t) represent an approximation to u(xj, t) of (4.5) where xj =jAx, j = 1..... J,

with Ax the uniform grid spacing. Consider the three-point central difference in space

with periodic condition u_+j = uj, and assume ZJ= _ uj = constant, which implies that

J

_,_= _duj/dt = 0. If we take J = 3 and Ax = 1/3 then, with _ = 9fl, this system can be

reduced to a 2 × 2 system of first-order nonlinear autonomous ODEs (4.4) with

U T = (u 1,u2) = (u, v). In this case, the nonlinear convection term is contributing to the

nonlinearity of the ODE system (4.4).

These four equations were selected to bring out the dynamics of numerics for four

different types of solution behavior of the ODEs. The dissipative complex system (4.1)

possesses either a unique stable fixed point or limit cycle with an unstable fixed point

depending on the value of_.. This is the rare situation where the analytical expression of

a limit cycle can be found. The purpose of choosing (4.1) is to illustrate the numerical

accuracy of computing a limit cycle and the spurious dynamics associated with this

type of asymptote. The damped pendulum (4.2), arising from modelling of a physical

process, exhibits a periodic structure of an infinite number of fixed points. The

predator-prey model (4.3), arising from modelling of biological process, exhibits

multiple stable fixed points without a periodic pattern as model (4.2). The perturbed

Hamiltonian model (4.4), which arises as a gross simplification of finite discretization of
230 H.C.YEEAND P. K. SWEBY

the viscous Burgers' equation, exhibits an unique stable fixed point. Following the

classification of fixed points of(3.1) in Section III, one can easily obtain the following:

Fixed Point of (4.1): The dissipative complex model has a unique fixed point at

(u, t,) = (0,0) for _,_<0. The fixed point is a stable spiral if_ < 0. It is a center if_ = 0. For

> 0, the fixed point (0, 0) becomes unstable with the birth of a stable limit cycle with

radius equal to _ centered at (0, 0). Figure 4.1 shows the phase portrait (u v plane) of

system (4.1) for e = - 1 and c,= 1 respectively. Here the entire (u, v) plane belongs to the

basins of attraction of the stable fixed point (0,0) if_, < 0. On the other hand, ift > 0, the

entire (u, t,) plane except the unstable fixed point (0, 0) belongs to the basin of attraction

of the stable limit cycle centered at (0,0).

Fixed Points of (4.2): The damped pendulum (4.2) has an infinite number of fixed

points, namely (kn, 0) for integer k. Ilk is odd, the eigenvalues of the Jacobian J(U E)are

of opposite sign and these fixed points are saddles. Ifk is even, however, two cases must

be considered, depending on the value of _:.If _,< 2 and positive, the eigenvalues are

complex with negative real part and the fixed points are stable spirals. If _,/> 2, the

eigenvalues are real and negative and the fixed points are nodes. If e = 0, the spirals

become centers. Figure 4.2 shows the phase portrait and their corresponding basins of

attraction for system (4.2). The different shades of grey regions represent the various

basins of attraction of the respective stable fixed points for _,= 0.5 and c = 2.5.

Fixed Points oJ"(4.3): The fixed points of the predator-prey equation are less regular

than those for the damped pendulum equation. System (4.3) has four fixed points C0,0),

(0, I), (3,0) and (2.1, !.98). By looking at the eigenvalues of the Jacobian of S, one finds

that C0,0) is a stable node, (2.1, 1.98) is a stable spiral, and (1,0) and (3,0) are saddles.

Figure 4.3 shows the phase portrait and their corresponding basins of attraction for

system (4.3). The different shades of grey regions represent the various basins of

attraction of the respective stable fixed points. The white region represents the basin of

3- t" =-1

--2

--3 .... _2 .... I1.... i .... r .... I .... -3 - - O 1 2 5

U

3-

2

> 0

I

-2

E=I

3 .... 12.... _1 .... ! .... I .... I .... I

-3 - - O I 2 3

U

Figure4.1 Phase Portraits and basins of Attraction Dissipative Complex Equation.
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 231

E =0.5

> 0.0 i

-2,5

-5,0

-10.0 -7,5 -5.0 -Z5 0.0 2.5 5.0 7.5 10.0

U

=2.5

2,5

> 0,0

-5.0

-10.0 -7.5 -5,0 -2.5 0.0 23 5.0 7.5 10.0

U

Figure 4.2 Phase Portraits and Basins of Attraction Damped Pendulum Equation.

>

4 5

Figure 4.3 Phase Portraits and Basins of Attraction Predator-Prey F.quatlon.
232 H.C. YEE AND P. K. SWEBY

Figure 4.4

Space).

1.0= 3

t 0 r- £- [-

U

0=3

-[

£

o <

¢-

£-

[-

0 r- £- [-

U

Phase Portraits and Basins of Attraction Viscous Burger's Equation (Central Difference in

divergent solutions. Note that the trajectories near the unstable separatrices actually

do not merge with the unstable branch ofseparatrices, but only appear to merge due to

the thick drawings of the solution trajectories.

Fixed Points of (4.4): The perturbed Hamiltonian has four steady-state solutions of

which three are saddles and one is a stable spiral at (1/3, 1/3) for e ¢ 0. For e, = 0 the

stable spiral becomes a center. Figure 4.4 shows the phase portrait and their corresponding basins of attraction for system (4.4). The shaded region represents the basins

of attraction for the fixed point (1/3, 1/3) for e = 0 and e = 0.01. The white region

represents the basin of divergent solutions. From here on we refer to (4.4) also as

a viscous Burgers' equation with central difference in space.

5. NUMERICAL METHODS AND BIFURCATION DIAGRAMS

This section describes the 11 time discretizations and their corresponding bifurcation

diagrams for the four model ODEs (4.1) (4.4). The 11 numerical methods are listed in

Section 5.1. Section 5.2 discusses the stability of selected fixed points of the discretized

counterparts of the model ODEs as functions of system parameters. Section 5.3

discusses the bifurcation diagrams as a function of the discretized parameter At with the

system parameter held fixed.

5.1 Numerical Methods

The 9 explicit and two implicit methods considered are the explicit Euler, two

second-order Runge-Kutta, namely, the modified Euler (R-K 2) and the improved

Euler (R-K 2), two third-order Runge-Kutta (R-K 3), a fourth-order Runge-Kutta

(R-K 4), the two and three-step predictor-corrector (Lambert, 1973), and noniterative

linearized forms of the implicit Euler and the trapezoidal methods.
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 233

(1) Explicit Euler (lst-order; R-K l):

U"+I=U"+rS"; S"=S(U"), (5.1)

(2) Modified Euler (R-K 2):

U"+1 =U"+rS(U"+2S"), (5.2)

(3) Improved Euler (R-K 2):

r n

U.+ 1 = U" + _ [S + S(U + rS")], (5.3)

(4) Heun (R-K 3):

r

U.+ 1 = U" + _ (k a + 3k3) (5.4)

k 2 =S n

2r

(5) Kutta (R-K 3):

(6) R-K 4:

n r

U "+I=U +g(k 1+4k2+k3)

k I =S n

k3 = S(U _- rk I + 2rk2),

U n + 1 U n r

= +-_(k_ +2k2+2k3+k4)

k I = S"

k2=S U"+_kl

(5.5)

(5.6)
234 H. C. YEE AND P. K. SWEBY

k3=S (r) Un +_k2

k 4 = S(U _ + rk3),

(7, 8) Predictor-.corrector for m = 2, 3 (PC2, PC3):

U_O)= U _ + rS _

U_k+ 1) = U" + 2 [Sn + S_k)]' k = O, 1..... m - 1

U _+ 1 = U _ + 2 [S" + S _m- 1_],

(9) Adam-Bashforth (2nd-order):

r n

U"+1= u"+_E3s(u )- s(u"-l)],

(10) Linearized Implicit Euler:

U "÷l = U _ + r(l-rJ _)-lS _

\OU,/ and det(l - rJ _) # 0,

(11) Linearized Trapezoidal:

(5.7)

(5.8)

(5.9)

U"+ I = Un + r I- J" S _ (5.10)

\8U,] and det l-_J _ =#0,

where the numeric identifier after the "R-K" indicates the order of accuracy of the

scheme and r = At and det( ) means the determinant of the quantity inside the ().

Schemes (10) and (1 l) are unconditionally stable methods. See Beam and Warming

(1976) and Yee (1989) for the versatility of the linearized implicit Euler and linearized

trapezoidal methods in CFD applications. A comparison between Newton method in

solving the steady part of the ODEs and the linearized implicit method (5.9) for model

(4.4) is included in Section 6.5. Studies on Newton method in solving the steady state

part of the PDE and some iteration procedures in solving the nonlinear algebraic
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 235

equation resulting from four implicit LMMs are reported in a separate paper (Yee and

Sweby 1933a). Although the explicit Euler can be considered as an R-K l, it is also

a LMM. All of the R-K methods (higher than first order) and the predictor-corrector

methods are nonlinear in the parameter space r, and all LMMs are linear in r. As

discussed in Yee et al. (1991), a necessary condition for a scheme to produce spurious

fixed points of period one is the introduction of nonlinearity in the parameter space r. It

can be shown later that this property plays a major role on the shapes and sizes of the

associated numerical basins of attraction of the scheme. For simplicity in referencing,

hereafter we use "implicit Euler" and "trapezoidal" to mean the linearized forms (5.9)

and (5.10), respectively, unless otherwise stated.

5.2 Stability of Fixed Points of Numerical Methods as a Function of System

Parameters

In our later study, we assume a fixed system parameter so that only the discretized

parameter comes into play. However, in order to get a feel for the numerical stability of

these schemes around selected stable fixed points Ue as a function of the system parameter e, Figures 5.1-5.3 show the stability regions of the schemes as a function of the

5-

4

3

2"

1"

0

4

R-K 1

5-

4-

3

2

,

il_ 1 0 ,ill

4-3-2-1012 -5_-3_-1012

£ (

R-K 2

il RK3 l RK4

i i i t 0 r i ' I

-5 -4 -3 -2 -1 • 1 2 -S -4 -3 -2 -1 0 ! 2

E £

Figure 5.1 Stability Regions vs. System Parameters Dissipative Complex Equation.
236 H.C. YEE AND P. K. SWEBY

il PC2

-5 .4 -3 -2 -1 0 1 2

(

5-

-5

PC3

-4 -3 -2 -1 0

(

$-

4-

3"

2"

1

0

Adams-Bashforth

-5 -4 -3 -2 -1 0 ! 2

E

Figure 5.1 (Continued)

I , I

1 2

system parameter e around a selected fixed point for each of the models. The linearized

stability regions for the R-K methods of the same order behave in exactly the same

manner, and the linearized stability regions around stable U e of the linearized

implicit methods are not interesting, since they have the same regions of stability as

the ODEs.

The stability diagrams presented were obtained by numerically solving the absolute

stability polynomials for the various methods, in most cases using Newton iteration.

For the Runge-Kutta schemes (of order p _< 4) the stability (Gritiiths et al., 1992;

Lambert, 1973) condition is that

l+,_r+---+ P[ <1,

(5.11)

where g are the eigenvalues of the Jacobian of S(U). For the Predictor-Corrector of

steps p = 2, 3 the stability condition is that

)tPrP I 1 +_r + ... +-_- <1,

(5ll)
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 237

and for the Adams-Bashforth method the roots p of

p2_ (1+--3 r)_t+_=O (5.13)

satisfy ]#[ < 1. Note that all of these expressions only hold for the U Efixed points of the

system.

In all cases the boundary of the stability region is when unit modulus is attained. The

linearized implicit Euler and trapezoidal methods are unconditionally stable for the

stable exact fixed points U r of the ODE systems we are considering.

These stability regions can be used to isolate the key regions of the e parameter to be

considered for the study of dynamics of numerics later. Due to the enormous number of

possibilities, detailed study can only concentrate on one to two representative system

parameters. Even with such a restriction, as can be seen later, computing the corresponding bifurcation diagrams and basins of attraction is very CPU intensive. Fortunately the computation can be made highly parallel. Figures 5.1-5.3 also can serve as

a spot check on the numerical results presented in the next section.

s R-K 1

4

$

0 1 2 3 4 !

(

R-K 2

,;tlllttl;l,rlt;I/ltltl I

1 2 3 4 S

(

il R-K3 s4-

3"

0 iiliiir I 0

O 1 2 3 4 S

£

R-K 4

;_tl ltrl IIltlTI_I I

I 2 3 4 S

(

Figure 5.2 Stability Regions vs. System Parameters Damped Pendulum Equation.
238 H.C. YEE AND P, K. SWEBY

0 ,;,,i,f,pl,l,rl

0 1 2 3 4 5

£

5-

4-

3"

PC3

!_r,_,llrrrl,,rl I

1 2 3 4 S

5

4

3

2

1

0

Adams-Bashlorth

i1, rl,Ir_l,Tt]lr rl rT,_tT I

1 2 ] 4 S

C

Figure 5.2 (Continued)

5.3 Bifurcation Diagrams

In this section, we show the bifurcation diagrams of selected R-K methods. It illustrates some of the many ways in which the dynamics of a numerical discretization

of 2 x 2 first-order autonomous nonlinear system of ODEs can differ from the

system itself. Note that there is no limit cycle or higher dimensional tori counterpart for the scalar first-order autonomous ODEs. Spurious limit cycles and higher

dimensional tori can only be introduced by the numerics when solving nonlinear

ODEs other than scalar first-order autonomous ODEs (if 2-time level schemes are

used) and/or by using a scheme with higher than two-time level for the scalar firstorder autonomous ODEs. In Section VI, we showed how numerical basins of attraction can complement the bifurcation diagrams in gaining more detailed global

asymptotic behavior of numerical schemes. We purposely present our results in

this order (not showing the basins of attraction) in order to bring out the importance

of basins of attraction for the time-dependent approach in obtaining steady-state

numerical solutions.

Even though the analytical solutions of the these models are known, depending on

the scheme, the dynamics of their discretized counterparts might be very difficult to

analyze. In particular, some analytical linearized analysis (without numerical computations) of fixed points of periods one and two is possible for the predator-prey and the

damped pendulum case. However, analytical analysis for the dissipative complex
DYNAMICAL STUDYOF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 239

model and the perturbed Hamiltonian is not practical. For a detailed analysis of these

selected cases, readers are referred to Sweby and Yee (1991). For the majority of the cases

where rigorous analysis is impractical we study the dynamics of numerics using numerical

experiments.

Note that some global solution behavior of fixed points of the nonlinear discretized

equations (5.1)-(5.10) for (4.1)-(4.4) can be obtained by the pseudo arclength continuation method devised by Keller (1977), a standard numerical method for obtaining

bifurcation curves in bifurcation analysis. A major shortcoming of the pseudo arclength

continuation method is that for problems with complicated bifurcation patterns, it

cannot provide the complete bifurcation diagram without known start up solutions

for each of the main bifurcation branches before one can continue the solution along

a specific main branch. For spurious asymptotes it is usually not easy to locate even just

one solution on each of these branches.

The nature of our calculations requires thousands of iterations of the same equation

with different ranges of initial data on a preselected (u, v) domain and range of the

discretized parameter space At. Since the NASA Ames CM-2 allows vast numbers

(typically 65,536) of calculations to be performed in parallel, our problem is perfect for

computation on the CM-2.

5-

4-

3"

2.

1'

0

0

R-K 1

tl,tl,lllllfl,l,,,lllr+, I

1 2 3 4 5

E

_!1 R-K2

i ,

0 1 2 3 4 S

E

:l

3

2

R-K 3 s4

3

1

ITI;I+IIIIITIIlI+I+IIIll I 0

1 2 3 4 5

R-K 4

irl JllllllllllllTI Jllrlll

1 2 3 4 S

£

Figure 5.3 Stability Regions vs. System Parameters Viscous Burgers' Equation (Central Difference in

Space).
240 H.C. YEE AND P. K. SWEBY

0tll_l[rl_'l'lllll_lJll'll I

0 1 2 3 4 5

£

0

0

PC3

I,l,l,q,rl,l,_j,,,,l,_ll I

1 2 3 4 5

E

5

4

3

2

1

0

Adams-Bashforth

''''ll'/llll_ll'lTli,,T, I

1 2 3 4 5

E

Figure 5.3 (Continued)

To obain a "full" bifurcation diagram, the domain of initial data and the range of the

At parameter are typically divided into 512 equal increments. For each initial datum

and At, the discretized equations are preiterated 3,000-5,000 (more or less depending

on the ODE and scheme) before the next 4,000 6,000 iterations are plotted. The

preiterations are necessary in order for the trajectories to settle to their asymptotic

value. The high number of iterations plotted (overlay on the same plot) is to detect

periodic orbits or invariant sets. Since the results are a three dimensional graph

((At, u, _,)), we have taken slices in a given constant v- and u-plane in order to enhance

viewing the decrease CPU computations. Note that with this method of computing the

bifurcation diagrams, only the stable branches are plotted. Some of the bifurcation

diagrams in a t, = constant plane for the four model ODEs and for the modified Euler,

improved Euler, Kutta and R-K 4 methods are shown in Figures 5.4 5.8. Figure 5.4

shows a typical example of spurious stable fixed points (branches 3 and 4 on the

diagram) occurring below the linearized stability by the modified Euler method. It also

shows the existence of spurious asymptotes such as limit cycles, higher order periodic

solutions and possibly numerical chaos (chaos introduced by numerics). See later

sections and subsections for further details. Selected bifurcation diagrams for the rest of

the numerical methods are illustrated in Section IV with basins of attraction

superimposed (see Figures 6.3-6.5,6.13, 6.14, and 6.19-6.20). See also the original

NASA internal report RNR-92-008, March 1992 for additional illustrations. Due to the
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS

Un

800

241

6 O0

4,00

2.00

0.00

-2.00

-400

/

0.20I 0.140 0.60I 0.80I 1 I I I I I I r=aDt .IXI 1.20 1,40 1,60 1,80 2.00

Figure 5.4 Bifurcation Diagram Predator-Prey Equation Modified Euler.

plotting package, the labels (u,, v,) on all of the figures are (u", v"). In the plots, r = At

unless stated.

The term "full bifurcation" as defined in Yee et al. (1991) is used to mean bifurcation

diagrams that cover the essential lower-order periods in such a way as to closely resolve

the "true" bifurcation diagram of the underlying discrete map for a selected range of

initial data domains. This is necessary since solutions with different initial conditions

will converge to different asymptotic limits. All of the computations shown are "full"

bifurcation diagrams.

The following summarizes the spurious dynamical behavior of the 11 numerical

methods based on selected domains of initial data and ranges of the discretized

parameter r. Numerical results agree with the analytical linearized analysis reported in

Sweby and Yee ( 1991 ).

Bifurcation Diaqrams of Numerical Methods for Model (4.1): For _:= 0, (4.1) is nondissipative (or a Hamiltonian system), and all of the I I numerical methods which are

non-simplectic converge quite slowly to the fixed point (0,0). We conjecture that

simplectic schemes (Sanz-Serna, 1990) would be more appropriate for t: = 0. For

sufficiently small negative (positive) t:, all of the studied schemes converge extremely

slowly to the stable spiral (limit cycle). This is a typical example of slow convergence

of the numerical solution due to the stiffness of the system parameter. While the

bifurcation diagrams for _:_<0 for the various numerical methods are not too

interesting, the bifurcation diagrams for ;: > 0 are very instructive. Figure 5.5 shows the

bifurcation diagrams for the four R-K methods for _:= 1.
242 H C. YEE AND P. K. SWEBY

Un

3.0.

2.0

,.o

o.o.

Modified Euler

J t

/

ols ,is

u n

3,0-

2.0

1.0

Improved Euler

, 0.0

2'.5 015 115 2'.5

r=a(

un

Figure 5.5

Kutta

=:s 3._ ,:5 '

u n

4. R-K 4

3-

2-

o

' 2

r=a_

Bifurcation Diagrams Dissipative Complex Equation, _ = 1, v = 0.0.

Note also that R-K 4 method gives the most overall accurate numerical approximations of the true limit cycle with radius _ centered at (0,0). The Adam-Bashforth,

PC2, PC3, implicit Euler and trapezoidal methods give the least accurate numerical

approximation of the limit cycle for r closer to the linearized stability. The R-K 4 and

Heun methods produced spurious higher-order limit cycles (invariant set of multiple

circles on the diagrams). See Section IV and Figures 6.6 and 6.8 for more details. These

diagrams illustrate the unreliability of trying to compute a true limit cycle with any

sizable r.This should not be surprising since the scheme only gives an O(rp)approximation to the solution trajectories. In addition, since the limit cycle is not a fixed point, we

would expect inaccuracies to be introduced. However, inaccuracies are not easy to

detect in practice, especially when a numerical solution produces the qualitative

features expected. See Section VI and Figures 6.3-6.9 for additional details. All of the

studied explicit methods produce spurious asymptotes.

For e > 0, the trapezoidal method produces no spurious steady states. However, the

implicit Euter method in addition to maintaining an unconditionally stable feature of

the exact limit cycle, also turns the unstable fixed point UE = (0, 0) of the ODE (4.1) into

a stable fixed point for r >t 1. See Figures 6.5, 6.8 and 6.9 for additional details.
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 243

Un Modified Euler Un Improved Euler

5"

4

3

2

1

o ,2. 2:o 2:, 2'8 o lls _o 2:,

Un

5

r=a

un

5-

4

3

2

1

,16 21o 214 ' 21a ,.s 21o 215

Kutta R-K 4

r=a£

Figure 5.6 Bifurcation Diagrams Damped Pendulum Equation, e = 1, v = 0.0.

Bifurcation Diagrams of Numerical Methods for Model (4.2): All of the studied explicit

and implicit methods produce spurious asymptotes. In particular, some of the explicit

methods (even explicit Euler) produce spurious limit cycles for certain e values. For

certain ranges of r and e values the implicit Euler and trapezoidal methods turn the

saddle points ot"(4.2)into an unstable fixed point of different type (see Figure 6.12). For

the modified Euler method, spurious steady states occur below the linearized stability

limit of the scheme. See Section VI and Figures 6.10-6.12 for additional details.

Bifurcation Diagrams of Numerical Methods for Model (4.3): Again, all of the

studied explicit and implicit methods generate spurious asymptotes. Also, some of

the explicit methods produce spurious limit cycles. For certain ranges of the r, the

trapezoidal method turns the saddle points (exact fixed points of (4.3)) into unstable

fixed points of different types, and the implicit Euler method turns the saddle points

into stable fixed points of different type. The numerical results coincide with analytical

analysis by examining the eigenvalues of the Jacobian of the resulting discrete map.

Transcritical bifurcations introduced by the R-K 4 method resulted in the production

of spurious steady state below (and very near) the linearized stability limit of the

scheme. See, Section VI and Figures 6.13 6.18 for additional details.
244 H C, YEE AND P. K. SWEBY

u n

Modified Euler

0.2 0.6 1.0

Un

4- Improved Euler

o" M

114 l'e ' -1 0'2 016 110 ' 114 ' 118 ,r = aDt

4

3

2

1

o

I

u n

Kutta

u n

R-K 4

=

o

02 0;6 110 114 110 ' -1 012 01S 1'.0 114 110 ,r = aDt

Figure 5.7 Bifurcation Diagrams Predator-Prey Equation, v = 0.0.

More than one spurious fixed point below the linearized stability of the scheme was

introduced by the modified Euler method (see Fig. 5.4). From the form of the modified

Euler scheme it is easily seen that as well as the exact fixed points UE of the ODEs, any

other value U s satisfying

r

U s +-_ S(Us) = U e (5.14)

will also be a fixed point of the scheme. As mentioned earlier, we refer to these

additional fixed points as spurious fixed points. Note that the U s on the right-hand side

of(5.14) encompasses both stable and unstable fixed points of the ODE and so, for the

predator-prey equations (since S contains cubic terms in U), there are up to twelve (real)

spurious steady states, three for each exact fixed point Us. In fact there are six such

spurious steady states which lie in the v = 0 plane. All of them occur below the

linearized stability limits of the exact fixed points, although not all are stable there.

Four (stable ones) of the six are shown in the bifurcation diagram of Figure 5.4,

numbered 1, 3, 4, 6 of the bifurcation branch. The other two are unstable. Note also that

the branch numbered 6 is in fact not stable but represents the stable eigen-direction

(separatrix) in the v = 0 plane of a saddle point.
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 245

U n

10

(]5

U n

Modified Euler Is Improved Euler

1.0

05

00

-05

-1.00P 5" i -10 _ T r = aDt

1.0 15 20 2r 5 30 0S 1_0 1S 210 2IS 310

f_

Un

1.5

1_o

05

o.o-

-0.5

-1 o

05

Un

15

lO0,5

c

oo

-05

r = aDt

lO 15 21o 25 30 1°15 20 25 3r0 35

Kutta R-K 4

Figure S.$ Bifurcation Diagrams Viscous Burgers' Equation, c = 0.1. l_= 0.333 (Central Difference in

Space).

Bifurcation Diagrams of Numerical Methodsfi_r Model (4.4): For _:= 0, the ODE (4.4) is

nondissipative and thus for small r, slow convergence was experienced. For r beyond

the linearized limit and with _ = 0 all of the explicit methods produce spurious limit

cycles. For r, > 0 (and not too large) all of the studied 11 explicit and implicit methods

produce spurious asymptotes. Also, all of the explicit methods produce spurious limit

cycles. For r,=0.1, the Kutta and Heun methods introduce spurious asymptotes

(higher than period one) that are below the linearized stability limit of the scheme. See

Figures 6.19 6.25 for additional details.

6. BASINS OF ATTRACTION AND BIFURCATION DIAGRAMS

This section illustrates how basins of attraction can complement the bifurcation

diagrams in gaining more detailed global asymptotic behavior of time discretizations
246 H.C. YEE AND P. K. SWEBY

for nonlinear DEs. Two different representations of the numerical basins of attraction were computed on the NASA Ames CM-2. One representation is bifurcation

diagrams as a function of At with numerical basins of attraction superimposed on

a constant v- or u-plane. The other representation is the numerical basins of attraction

with stable asymptotes superimposed on the phase plane (u,v). Before discussing

numerical results for each of the model ODEs, the next subsection gives some

preliminaries on how to compute and on how to interpret the basins of attraction

diagrams for the CM-2.

6.1 Introduction

To obtain a bifurcation diagram with numerical basins of attraction superimposed on

the CM-2, the preselected domain of initial data on a constant v- or u-plane and the

preselected range of the At parameter are divided into 512 equal increments. Again, for

the bifurcation part of the computations, with each initial datum and At, the discretized

equations are preiterated 3,000-5,000 steps before the next 5,000 iterations (more or

less depending on the problem and scheme) are plotted. The bifurcation curves appear

on the figures as white curve, white dot and white dense dots. While computing the

bifurcation diagrams it is possible to overlay basins of attraction for each value of At

used. For the numerical basins of attraction part of the computation with each value of

At used, we keep track of where each initial datum asymptotically approaches and

color code them (appearing as a vertical strip) according to the individual asymptotes.

While efforts were made to match color coding of adjacent strips on the bifurcation

diagram, it was not always practical or possible. Care must therefore be taken when

interpreting these overlays.

For the basins of attraction on the phase plane (u, v) with selected values of At

and the stable asymptotes superimposed, the (u, v) domain is divided into 512 × 512

points of initial datum. With each initial datum and the selected At, we preiterate

the respective discretized equation 3,000 5,000 steps and plot the next 5,000 steps

to produce the asymptotes. Again, for the basins of attraction part of the computations, for each value of At used, we keep track of where each initial datum

asymptotically approaches and color code them according to the individual asymptotes. All of the selected time steps At shown are based on the bifurcation diagram

with the basins of attraction superimposed. The chosen time steps were selected

to illustrate special features of the different bifurcation phenomena on the (u,v)

plane. Details of the techniques used for detection of asymptotes and basins of

attraction are given in the appendix of Sweby and Yee (1991). Note that in all of

the plots, ifcolor printing is not available, the different shades of grey represent different

colors.

As a prelimiary, and before discussing our major results, we discuss the numerical

basins of attraction associated with modified Euler, improved Euler, Kutta and R-K 4

methods for the two scalar first-order autonomous nonlinear ODEs studied in part 1of

our companion paper (Yee et al., 1991_. The two scalar ODEs are:

du

dt = au(1 - u} (6.1)
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 247

and

du

dt = au( I - u)(0.5 - u). (6.2)

The fixed points for (6.1) with a > 0 are u = 0 (unstablet and u = 1 Istable), and no

additional higher order periodic fixed points or asymptotes exist. The basin of attraction

for the stable fixed point u = 1 is the entire positive plane for all values of a > 0.

The fixed points for (6.2) with a > 0 are u = 0 (unstable), u = 1 (unstablel and u = 0.5

tstable) and no additional higher-order periodic solutions or asymptotes exist. The

basin of attraction for the stable fixed point u = 0.5 is 0 < u < 1for all a > 0. The white

curve, white dots and white dense dots of Figures 6.1 and 6.2 show the bifurcation

diagrams for four of the R-K methods for (6.1) and (6.2). For more details of the dynamics

of numerics for systems (6.1) and (6.2), see Yee et al. (1991). Intuitively, in the presence of

spurious asymptotes the basins of the true stable steady states can be separated by the

numerical basins of attraction of the stable and unstable spurious asymptotes.

Take, for example, the ODE (6.1) where the entire domain u is divided into two basins

of attraction for the ODE independent of any real a. Now if one numerically integrates

the ODE, depending on the scheme and r, extra stable and unstable fixed points of any

order can be introduced by the scheme. The bifurcation part of Figures 6.1 and 6.2,

cannot distinguish the types of bifurcation and the periodicity of the spurious fixed

points of any order. With the numerical basins of attraction and their respective

bifurcation diagrams superimposed on the same plot, the type of bifurcation and to

which initial data asymptote to which stable asymptotes become apparent. Note that

for Figures 6.1 and 6.2 r = aAt.

For example, any initial data residing in the green region in Figure 6.1 for the

modified Euler method belong to the numerical basin of attraction of the spurious

(stable) branch emanating from u = 3 and r = 1. Thus, if the initial data is inside the

green region, the solution can never converge to the exact steady state using even

a small fixed but finite At (all below the linearized stability limit of the scheme). Note

that the green region extends upward as r decreases below 1.Thus for certain ranges of

r values, the domain is divided into four basins (instead of two for the ODE). But of

course higher period spurious fixed points exist for other ranges of r and more basins

are created within the same u domain.

A similar situation exists for the R-K 4 method (Fig. 6.1t, except now the numerical

basins of attraction of the spurious fixed points occur very near the linearized stability

limit of the scheme, with a small portion occurring below the linearized stability limit.

In constrast to the improved Euler method (Fig. 6.1), the green region represents the

numerical basins of one of the spurious stable transcritical bifurcation branches of the

fixed point. The bifurcation curve directly below it with the corresponding red portion

is the basin of the other spurious branch. See Yee et al. (1991 ) or Hale and Kocak (1991)

for a discussion of the different types of bifurcations. With this way of color coding the

basins of attraction, one can readily see tfrom the plots) that for ODE (6.1), the modified

Euler, improved Euler and R-K 4 methods, experience one steady bifurcation before

a period doubling bifurcation occurs (Fig. (6.1)j. Using the PC3 method to solve 16.1)

tfigure not shown; see Yee et al., 1991t, more than one consecutive steady bifurcation

occurs before period doubling bifurcation. For ODE (6.2), the improved Euler experi-
248 H.C.YEEANDP.K.SWEBY

Figure6.1 (See ('olor Plate 1al the back of this issue.}

cnces two consecutive steady bifurcations before a period doubling bifurcation occurs

(Fig. (6.2}!. Using the PC3 method to solve _6.2) (figure not shown; see Yee et al. 1991),

four consecutive steady bifurcations occur before period doubling bifurcations. The

modilied Euler and R-K 4 methods, however, experience only one stead',' bifurcation

before period doubling bifurcations occur.
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 249

Figure 6.2 {See Color Plate II at the back of lhis issue.)

The next section presents similar diagrams for the 2 x 2 systems of model nonlinear

ODEs (4.1) (4.4). In this case, only basins of attraction with bifurcation diagrams

superimposed on t_ = constant planes are shown. Selected results for both representations of numerical basins of attraction are shown in Figures 6.3 6.5, 6.8 6.25 for the

numerical methods. Section 6.6 summarizes similar results presented in Yee and
250 H.C. YEE AND P. K. SWEBY

Sweby (1993a) for iterative procedures in solving nonlinear systems of algebraic

equations arising from four implicit LMMs. In the plots r = At. White dots and white

curves on the basins of attraction with bifurcation diagrams superimposed represent

the bifurcation curves. White dots and white closed curves on the basins of attraction

with the numerical asymptotes superimposed represent the stable fixed points, stable

periodic solutions or stable limit cycles. The black regions represent divergent

solutions.

Note that the streaks on some of plots are either due to the non-settling of the

solutions within the prescribed number of preiterations or the existence of small

isolated spurious asymptotes. Due to the high cost of computation, no further

attempts were made to refine their detailed behavior since our purpose was to show

how, in general, the different numerical methods behave in the context of nonlinear

dynamics.

6.2 Numerical results for the Dissipative Complex Equation

Figures 6.3-6.5, 6.8, 6.9 show selected results for the two representations of numerical

basins of attraction for model (4.1) for e = 1. The exact solution for (4.1) with e = 1 is

a stable limit cycle with unit radius centered at (0, 0). The basin of attraction for the limit

cycle is the entire (u, v) plane except the unstable fixed point (0, 0).

Comparing Figures 6.3-6.5 with Figure 5.5, one can appreciate the added information that the basin of attraction diagrams can provide. As At moves closer to the

linearized stability limit of the limit cycle, the size (red) of the numerical basins of

attraction decreases rapidly. This is due to the existence of spurious unstable asymptotes below as well as above the linearized stability limit. The green region, shown in

Figure 6.5 using the implicit Euler method, is the numerical basin of attraction for the

stabilized fixed point (0, 0). Note how the implicit Euler method turns an unstable fixed

point (0, 0) of the ODE system into a stable one for At t> 1.

Figures 6.6 and 6.7 show the phase trajectories and Figures 6.8 and 6.9 show the

same figures with numerical basins of attraction superimposed for four different At by

the R-K 4 and implicit Euler methods, respectively. Note how little information

Figures 6.6 and 6.7 can provide as compared to Figures 6.8 and 6.9. Note also how

rapidly the size of the basin (redl decreases as At increases for the R-K 4 method. This

phenomenon can relate to practical computations where only a fraction of the

allowable linearized stability limit of At is safe to use if the initial data is not known. For

At = 1.75 and 2, spurious limit cycles of higher order period exist. (The multiple white

circles with only one distinct basin of attraction). In this case, the red regions represent

the basins of the spurious numerical solutions.

Figure 6.9 illustrates the situation where unconditionally stable LM M schemes can

converge to a wrong solution if one picks the initial data inside the green region which

are valid physical initial data for the ODE. Thus even though LM M preserved the same

number of fixed points as the underlying ODE, these fixed points can change type

and stability. This phenomenon is related to the "non-robustness" of implicit methods

sometimes experienced in CFD computations. In this type of computation where the

initial data are not known, the highest probability of avoiding spurious asymptotes is

achieved when a fraction of the allowable linearized stability limit of At is employed.
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 251

Figure 6.3 (Scc Color Pl_tlc 111_t thc back of this issue.)

6.3 Nttmerical Restdts /'or the Damped Penduhtm Equation

Selected results for the studied numerical methods for _:= 1 and _:= 1.5 t_re shown in

Figures 6.10 6.12. Here. for each At value the different colors represent diffcrcnt

numerical b_sins of _tttraction of the respective asymptotes. Observe the striking
252 H.C. YEE AND P. K. SWEBY

Figure 6.4 {SeeColor Plate IVat the back of this issue.)

difference in behavior between the explicit and implicit methods. The shapes and sizes

of the numerical basins of attraction by the implicit Euler method for At = 0.1 shown in

Figure 6.12 appear to be similar to the exact basins of attraction of the DE (4.2). From

the different colors of the basins in Figure 6.10 one can readily identify that spurious

higher than period one and spurious limit cycles exist for the different At values by the
DYNAMICAL STUDY OF" SPURIOUS STEADY-STATE NUMERICAl. SOEUTIONS 253

e-,

©
254 H.C. YEE AND P. K. SWEBY

v n

3-

-3

-3

v n

3"

-3

At = 0.5

Vn

3- At = 1.5

, i UR

_, _ _ -3 -3 -, ; 3

At = 1.75

©

v n

3 At= 2

1

@

-1

, r , -3 iun

-1 1 3 -3 -'1 1 3

Figure 6.6 Phase Trajectories Dissipative Complex Equation, e = 1, R-K 4.

explicit Euler method. For At = 1.4, the explicit Euler produces spurious period two

fixed points. Figure6.11 shows the existence of spurious fixed points below the

linearized stability limit by the modified Euler method and spurious fixed points of

period 4 (the four white dots one each basin) above the linearized stability limit by the

R-K 4 method. Figure 6.12 shows the evolution (birth and death) of spurious fixed

points of higher-order period for the implicit Euler method. This figure illustrates

another situation where unconditionally stable schemes can converge to a wrong

solution even though these schemes preserved the same number and type of fixed points

as the underlying ODE. In this case it is the birth of spurious stable and unstable

asymptotes or even numerical chaos that contributes to the size reduction of the true

basins of attraction of the ODE.

6.4 Numerical Results for the Predator-Prey Equation

Selected results for the two representations of numerical basins of attraction are shown

in Figures 6.13 6.18. Comparing Figures 6.13, 6.15, 6.16, with Figures 5.4 and 5.7, one

can again appreciate the added information that the basin of attraction diagrams can
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 255

Vn V n

3" At 3" = 0.5

-1

At = 1.5

-3 -3 Un

Vn v n

3- 3- At= 2 Al=4

-3 -3 , Un

Figure 6.7 Phase Trajectories Dissipative Complex Equation, e = 1 Linearized Implicit Euler.

provide. Here for all of these figures (except Fig. 6.15 for the last four At values), the

green regions represent the numerical basins of attraction for the stable spiral

(2.1, 1.98) and red regions represent the numerical basins of attraction for the stable

node (0, 0).

The numerical basins ofattraction in Figure 6.18 with At = 0.1 appear to be the same

as the exact basins of attraction of the DE (4.3). The numerical basin of attraction by the

implicit Euler for the fixed point (0,0) with At = 0.1 is larger than the corresponding

exact basin of attraction for the DE (4.3). In this case the numerical basin of attraction

for the divergent solution (black region) is smaller than the true one. The dramatic

difference in shapes and sizes of numerical basins of attraction for the different methods

and solution procedure combinations compared with the exact basin of attraction is

even more fascinating than for the previous two models.

Take, for example, one of the most interesting cases, the modified Euler method.

Figure 6.15 shows how spurious stable fixed points can alter the numerical basins of

attraction of the stable node and spiral of the ODE (4.3) for time steps that are below the

linearized stability limit of both of these stable fixed points of the ODE (see Figs. 5.4
256 H.C. YEE AND P. K. SWEBY

Figure 6.8 [Sce Color F'late VI at thc back of this issue.)

and 6.13). For At = 0.8, the stable node bifurcates into a spurious fixed point. Without

performing the bifurcation analysis one would not be able to detect this particular

spurious fixed point, since the value of the spurious one is so close to the exact fixed

point U_c = 10, 0_. For At = 0.9524, there is the birth of a spurious limit cycle (the white

close curvel. For At = 1.2, spurious higher-order periodic solutions exist. Note that for
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 257

Figure 6.9 (See Color Plate VII at the back of this issue.

the first four At values in Figure 6.15, the fixed points and asymptotic values are colored

black instead of white due to the birth of additional numerical basins of attraction that

are colored white.

The implicit methods change the two saddle points into stable or unstable fixed

points of other types as illustrated in Figures 6.14, 6.17 and 6.18. For the implicit Euler,
258 H.C. YEE AND P. K. SWEBY

Figure 6,1(I ISec c.'t_lt_z-Plalc VIII ill the back of Ibis issuc.I

the two tixed points (2.1, 1.98) and (0,0) arc unconditionally stable and the stahilized

fixed points (I, 0) and (3, 0)(saddles for the original ODE) are almost unconditionally

stable except for small At. This is most interesting in the sense that the numerical basins

of attraction for the stable exact fixed points U_:,of the model (4.3) by the implicit Euler

method were permanently altered for At near or larger than 3 as illustrated in
DYNAMICAL STUDYOFSPURIOUS STEADY-STATE NUMERICAL SOLUTIONS259

-i

t_

_x
260 H.C. YEE AND P. K. SWEBY

Figure 6.12 (See Color Plate X at the back of this issue.)

Figures6.14,6.17. It would be easier to interpret the results in Figure6.14 if one

interchanged the yellow and green colors for At >_ 1. Observe how the newly created

numerical basins of attraction by the stabilized fixed points (1,0) and (3, O) resulted in

the segmentation of the numerical basins of attraction of the stable node (0,0) and

stable spiral (2.1, 1.98). Although the trapezoidal method did not turn the two saddle
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 261

Figure 6.13 (See Color Plate XI at the back of this issue.I

points (1,0) and (3, 0) into stable fixed points of different type, they did turn the two

saddle points into unstable fixed points of different type.

The evolution of the numerical basin of attraction as At changes is very traumatic

these implicit LM Ms. The cause of nonconvergence of these implicit LM Ms may due to

the fact that their numerical basins of attraction are fragmented. Take for example the
262 H.C. YEE AND P. K. SWEBY

rX

e_

o

_3
DYNAMICALSTUDYOFSPURIOUSSTEADY-STATE NUMERICAL SOLUTIONS263

trapezoidal method(Fig.6.18)wheretheschemebecomeseffectivelyunstableforlarge

At.Thesizeofthenumerical basinsofattractionforthestablefixedpointsU E shrink to

almost nonexistence. This phenomenon might be one of the contributing factor to the

unpopularity of the trapezoidal method in CFD. The basins are so fragmented and

small for large At that they are beyond the accuracy of the CM2 to resolve and no

Figure 6.15 (See Color Plate XII1 al the back of this issue.}
264 H.C. YEE AND P. K. SWEBY

Figure 6.15 [Continued) {See Color Plate XIII at the back of this issuc.)

further attempt was made. A better approach in computing these types of basins is to

use interval arithmetic or the enclosure type method (Adams, 1990).

6.5 Numerical Results fi_r the Perturbed Hamiltonian Equation

Selected results for the two representations of numerical basins of attraction of the

various numerical methods for _:= 0.1 are shown in Figures 6.19 6.25. Our studies
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 265

Figure 6.16 (See Color Plate XIV at the back of this issue.)

indicate that all of the studied Runge-Kutta methods exhibit spurious limit cycles and

other spurious periodic solutions. For the Kutta and Heun methods, stable spurious

asymptotes can occur below the linearized stability limit of the scheme. The implicit

methods also exhibit spurious asymptotes. In particular, unstable spurious asymptotes

were produced below the linearized stability limit by all of the studied schemes.
266 H.C.YEEANDP.K.SWEBY

Althoughthisexampleconsistsofanartificiallysmallnumberofgridpoints,it can

shedsomelightontheinterplaybetweeninitialdata,spurioustableandunstable

asymptotes, basinsofattractionandthetime-dependent approachto theasymptotic

numerical solutions. Asolidunderstandingofthisconceptatthefundamental levelcan

Figure6.17ISccColorPlateXVatthebackofthisissue.t
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 267

e-,

x

i----

r_
268 H.C. YEE AND P. K. SWEBY

help to determine the reliability of the time-dependent approach to obtaining steadystate numerical solutions.

In all of Figures 6.19 6.25, red regions represent the numerical basins of attraction

for the stable spiral (1/'3, 1/3) when At is below the linearized stability of the scheme.

When At is above the linearized stability, some of the red regions represent the

numerical basin of attraction of the stable spurious asymptotes. The numerical basins

Figure 6.18 <See Color Plate XVI at the back of this issue.)
DYNAMICAL STUDYOF"SPURIOUSSTEADY-STATE NUMERICAL SOLUTIONS269

A

x

o_

c

_o
270 H.C. YEE AND P. K. SWEBY

Figure 6.19 (SeeColor Plate XVII at the back of this issue.)

of attraction in Figures 6.21 with At = 0.1 appear to be the same as the exact basins of

attraction. Note also that the possibility of the numerical basin of attraction being

larger than the exact one does not always occur when the time step is the smallest. The

numerical basin of attraction for (1/3, 1/3) is larger than the corresponding exact basin

of attraction for At = 1 by the improved Euler and Kutta methods and for At = 0.1 by
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 271

3

x

'-6
272 H.C. YEE AND P. K. SWEBY

Figure 6.21 (See Color Plate XIX at the back of this issue.)

the implicit Euler and trapezoidal

discusses results for the improved

trapezoidal methods.

methods. See Figures 6.21 6.24. The following

Euler, the Kutta, the implicit Euler and the

Impro_'ed Euler Method: This example illustrates the existence of spurious limit cycles

and its effect on the numerical basins of attraction for the exact steady state. Figure 6.21
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 273

Figure 6.22 (See Color Plate XX at lhc back of this issue.)

shows the basins of attraction of the improved Euler method for 4 different At = 0.1, 1,

2.25,2.35 with _:= 0.1. By a bifurcation computation shown in Figure 6.19, we found

that the first two time steps are below the linearized stability limit around the exact

stable steady state tl,,"3 1,,"3), and the last two time steps are above the limit.
274 H.C. YEE AND P. K. SWEBY

Above the linearized stability limit spurious limit cycles and higher dimensional

periodic solutions were observed. Further increasing At resulted in numerical chaostype phenomena and eventually divergence (with additional increase in At). For

At = 2.25 and 2.35, the red or multicolor regions are the basins of the spurious limit

cycle (the irregular white closed curve shown on Fig. 6.21) or other type of spurious

Figure 6.23 [See Color Plate XXI at the back of this issue.t
DYNAMICAL STUDYOFSPURIOUS STEADY-STATE NUMERICAL SOLUTIONS275

v¢

>¢

"-6

v

",o
276 H.C. YEE AND P. K. SWEBY

Figure 6.24 (See Color Plate XXII at the back of this issue.J

asymptote (white dots for Fig. 6.21). For these two time steps the numerical basins for

the exact steady state (1/3, 1/3) by the improved Euler method disappeared. However, if

the initial data are in the red or multicolor region, one gets spurious solution instead of

what the linearized stability predicts divergent solution.
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 277

¢-.

e_

,....

x
278 H.C. YEE AND P. K. SWEBY

Kutta Method: To give an example of the existence of spurious stable asymptotes below

the linearized stability limit of the scheme, as well as the existence of spurious limit

cycles above the linearized stability limit, Figure 6.22 shows the basins of attraction

of the Kutta method for 6 different fixed time steps At =0.1, l, 1.826, 1.85,2.75 and

2.785 (the first four below the linearized stability of the scheme) with e = 0.1. For

At = 1.826, the numerical basin for (1/3, 1/3) has become fractal like with the birth of

fragmented, isolated new basins of attraction due to the presence of spurious periodic

solutions (the three white complicated closed curves with the associated purple, green

and blue basins shown in Fig. 6.22). The last two time steps in Figure 6.22 show the

diappearance of the numerical basin of attraction for the exact steady state with the

birth of basins for the spurious limit cycle.

Implicit Euler Method: This is yet another interesting illustration of the use of an

unconditionally stable implicit method where in practical computations, when the

initial data are not known, the scheme has a higher chance of obtaining a physically

correct solution if one uses a At restriction slightly higher than that for the stability limit

of standard explicit methods (but with larger numerical basins of attraction than the

explicit method counterparts). Figures 6.20 and 6.23 show the two representations of

numerical basins of attraction using the implicit Euler method. These figures show the

generation of stable spurious asymptotes for At >/I. As At increases further, the size of

the same numerical basin decreases and becomes fractal like, and new numerical basins

are generated. The behavior is similar to the predator-prey model (4.3) in a sense that

the numerical basin of attraction for (1/3, 1/3) was permanently altered for At near or

larger than 10. Observe the fragmentation of the numerical basin of attraction for

(1/3, 1/3) by the basins of the spurious asymptotes.

Trapezoidal Method: Figures 6.20 and 6.24 show the two representations of numerical

basins of attraction using the trapezoidal method. As in the implicit Euler case, this

scheme has a higher probability of obtaining a physically correct solution if one uses a

At similar to that of standard explicit methods (but with larger numerical basins of

attraction than the explicit method counterparts). In a manner similar to the implicit

Euler, the numerical basins of attraction for (1/3,1/3) are much larger than the

corresponding exact basin of attraction for At _<2. Their sizes are bigger than the ones

generated by the implicit Euler method with the same At values. The scheme becomes

effectively unstable due to the fragmentation of the numerical basins of attraction.

Again due to the high cost of double precision computations, no further attempts were

made for At large. The computation of these basins requires an interval arithmetic or

the enclosure-type (Adams, 1990) of mathematical operation before a more precise

behavior can be revealed.

Straiqht Newton rs. Other Studied Methods: Figure 6.25 shows the basin of attraction

using Newton method in solving the steady part of the ODE (hereafter referred to as

straight Newton) compared with the implicit Euler at At = 1. One can see that straight

Newton method has a smaller attracting basin for the stable spiral (1/3, 1/3) than the

implicit Euler method for small At. In fact its basin is the same as the implicit Euler

using larger At. Figure 6.25 illustrates the situation where quadratic convergence by the

Newton method can be achieved only if the initial data are in the red regions.
DYNAMICALSTUDYOFSPURIOUSSTEADY-STATE NUMERICAL SOLUTIONS279

Figures6.23and6.25alsoillustratesthefactthatusingverylargeAtbythe(linearized)

implicitEulermethodhasthesamechanceofobtainingthecorrectsteadystateasthe

Newtonmethodif theinitialdataarenotknown.ComparisonofNewtonmethodwith

otheriterationproceduresfortheimplicitEulerandtrapezoidal methodsarereported

inourcompanionpaper(YeeandSweby,1993a).

CombiningthecurrentresultwithYeeandSweby(1993a), wecanconcludethat

contrarytopopularbelief,theinitialdatausingthestraightNewtonmethodmaynot

have to be close to the exact solution for convergence. Straight Newton also exhibits

stable and unstable spurious asymptotes. Initial data can be reasonably removed from

the asymptotic values and still be in the basin of attraction. However, the basins can be

fragmented even though the corresponding exact basins of attraction are single closed

domains. The cause of nonconvergence may just as readily be due to the fact that its

numerical basins of attraction are fragmented.

6.6 Global Asymptotic Behavior of lterative Implicit Schemes

The global asymptotic nonlinear behavior of some standard iterative procedures in

solving nonlinear systems of algebraic equations arising from four implicit linear

multistep methods (LM Ms)in discretizing models (4. l), (4.3) and (4.4) is analyzed in our

companion paper (Yee and Sweby, 1993a). The implicit LMMs include implicit Euler,

trapezoidal, mid-point implicit and three-point backward differentiation methods. The

iterative procedures include simple iteration and full and modified Newton iterations.

The results are compared with standard Runge-Kutta explicit methods, a non-iterative

implicit procedure, and straight Newton method. Here we give a summary of Yee and

Sweby (1993a) so that the reader may get a bigger picture of implicit methods other

than the ones studied in this paper.

Studies in Yee and Sweby (1993a) showed that all of the four implicit LMMs

exhibit a drastic distortion but less shrinkage of the basin of attraction of the true

solution than standard explicit methods studied in this paper. In some cases with

smaller At, the implicit LMMs exhibit enlargement of the basins of attraction of

the true solution. Overall, the numerical basins of attraction of a non-iterative implicit procedure mimic more closely the basins of attraction of the continuum than

the studied iterative implicit procedures for the four implicit LMMs. In general the

numerical basins of attraction bear no resemblance to the exact basins of attraction.

The size can increase or decrease depending on the time step. Also the possible existence

of the largest numerical basin of attraction that is larger than the exact one does

not occur when the time step is the smallest. The dynamics of numerics of the implicit

methods differ significantly from each other, and the different methods of solving

the resulting non-linear algebraic equations are very different from each other

since different numerical methods and solution procedures result in entirely different

nonlinear discrete maps. Although unconditionally stable implicit methods allow

a theoretically large time step At, the numerical basins of attraction {allowable initial

data) for large At some-times are so fragmented and/or so small that the safe (or

practical) choice of At is slightly larger or comparable to the stability limit of standard

explicit methods (but with larger numerical basins of attraction than the explicit

method counterparts). In general, if one uses a At that is a fraction of the stability limit,
280 H.C. YEE AND P. K. SWEBY

one has a higher chance of convergence to the correct asymptote than the standard

explicit methods.

Studies in Yee and Sweby (1993a) also showed that the variable time step control

method can occasionally stabilize unstable fixed points, depending on the initial data,

starting time step and the iterative tolerance value. One shortcoming is that the size of

At needed to avoid spurious dynamics is impractical to use, especially for the explicit

method.

7. CONCLUDING REMARKS

The global asymptotic nonlinear behavior and bifurcation phenomena for the explicit

Euler method, five different multistage Runge-Kutta methods (modified Euler, improved

Euler, Heun, Kutta and 4th-order methods), two and three-step predictor-corrector

methods, Adams-Bashforth method, and implicit Euler and trapezoidal method with

linearization are compared for different model nonlinear ODEs. The five multistage

Runge-Kutta methods and the predictor-corrector methods are nonlinear in the

discretized parameter space At and all LMMs are linear in At. With the aid of the

CM-2, the complex behavior and sometimes fractal like structure of the associated

numerical basins of attraction of these time discretizations are compared and revealed

for the first time.

The numerical results indicate that with sufficiently small At and initial data close to

the steady state (usually not known for the time-marching method), one can have the

highest chance of convergence to the correct asymptote. In general, the initial data can

be far removed from the exact steady state by the studied implicit methods provided

that a fraction of the allowable time step restriction is used. Our study also indicates

that bifurcation to a period two or lower order period soluton is readily detectable in

numerical calculations. However, bifurcation to a limit cycle will not be so obvious

(without a phase portrait representation), especially in the vicinity of the bifurcation

point. Indeed the phenomenon of an artificial time iteration to steady-state of a large

system formed by spatial discretization which nears convergence before the residuals

"plateaus out", could actually be the result of a stable spurious limit cycle around the

Hopf bifurcation point. In addition, the bifurcation of spirals to limit cycles might account

in part for the phenomenon of near {but lack of) convergence in large stiff systems.

For a given initial data and two finite but different At's that are below the linearized

stability limit of the scheme, their numerical solutions might converge to two different

solutions even if no spurious stable steady-state numerical solution is introduced by the

scheme and the initial data are physically relevant. The source of the behavior is due to

the existence of unstable spurious asymptotes or stable asymptotes other than steady

states which have the same detrimental (in terms of robustness) effect. However, in the

case of occurrence of stable spurious steady states, they can be mistaken for the true

steady state in practical computations. In other words depending on the initial data, for

a given At below the linearized stability limit, the numerical solution can (a) converge

to the correct steady state, (b) converge to a different steady state, (c) converge to

a spurious periodic solution, (d) yield spurious asymptotes other than (a) ic), or (e)

diverge, even though the initial data are physically relevant.
DYNAMICAL STUDY OF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 281

Another important finding is that unlike the scalar first-order autonomous ODE

discussed in part I (Yee et al., 1991), the fixed points can change types as the time step is

varied even for two-time-level unconditionally stable implicit LMMs. An unstable

fixed point can become a stable fixed point and can e.g., change from a saddle to a stable

or unstable node (for a fixed system parameter _:).Since these implicit methods can

introduce spurious asymptotes as well, thus even though LMMs preserve the same

number but not the same types of fixed points as the underlying DEs, the numerical

basins of attraction of LM Ms (explicit or implicit) do not always coincide with the exact

basins of attraction of the underlying DEs. One major consequence of this behavior is

that the flow pattern can change type as the discretized parameter is varied. Another

consequence of these phenomena is the fragmentation of the numerical basin of attraction. In general, unconditionally stable implicit LMMs exhibit less shrinkage of the

basin of attraction of the true solution than standard explicit methods. Another

interesting result is that contrary to popular belief, the initial data using the straight

Newton method may not have to be close to the exact steady state for convergence.

However, we believe that one cause of nonconvergence in straight Newton or implicit

LM Ms with large time step may due to the fact that the numerical basins of attraction

are fragmented.

In conclusion, the present results can explain some of the roots of why one cannot

achieve the theoretical linearized stability limit of the typical implicit Euler and

trapezoidal time discretization in practice when solving strongly nonlinear DEs, e.g. in

CFD. The results can also shed some light in bridging some of the gaps between

theoretical convergence criterion (At--* 0, as n---, _c,) and practical scientific computation (finite At as n ---, _).

Acknowledqements

The authors wish to thank A. M. Stuart and D. F. Grilfiths for their valuable discussions during the course of

this research and for suggesting the model equations (4. I ) and 14.4). Special thanks to A. Lafon, M.Vinokur,

T. Coakley and C. Y. McNeil for their critical review of the manuscript. Financial support from T. Lasiuski

for the second author as a visiting scientist at NASA Ames is greaffully acknowledged.

Rel'erences

Adams, E. (1990) "'Periodic Solutions: Enclosure, Verification, and Applications", Computer Arithmetic and

Se!{Validatin_l Numerical Methods, Academic Press, 199 245.

Beam, R. M. and Bailey, H. E. (1988) "Direct Solver for Navier-Stokes Equations", Proceedings of International Conference on Computational Engineering Science, Atlanta, GA.

Beam, R. M. and Warming, R. F. (1976)"An Implicit Finite-Difference Algorithm for Hyperbolic Systems in

Conservation Law Form", J. Comput. Phys., 22, 87 110.

Budd, C. J., Stuart, A. M., Koomullil, G. P. and Yee, H. C. (1994) "Numerical Solution Behavior of Model

Convection-Diffusion BVP with Grid Adaptation", in preparation.

Foias, C., Sell, G. and Temam, R. ( 1985)"Varietes lnertielles des Equations Differentielles Dissipatives", C.R.

Aead. Sci. Paris, Ser. I Math., 301, 139 141.

Globus, A., Levit, C. and Lasinski, T. (1991)"'A Tool for Visualizing the Topology of Three-Dimensional

Vector Fields", NAS Applied Research Branch Report RNR-91-017, NASA Ames Research Center.

Griffiths, D. F., Sweby, P. K. and Yee, H.C. (1992a) "On Spurious Asymplotes Numerical Solutions of

Explicit Runge-Kuna Schemes", IMA J. Numer. Anal., 12, 319 338.

Griffiths, D. F., Stuart, A. M. and Yee. H. C. ( 1992b)"Numerical Wave Propagation in Hyperbolic Problems

with Nonlinear Source Terms", SIAM J. _?fNumer. Analy., 29(5), 1244 1260.
282 H.C.YEEANDP.K.SWEBY

Guckenheimer, J.andHolmes, P.(1983) Nonlinear Oscillations. Dynamical Systems. and Bifurcations of

Vector Fields, Springer-Verlag, New York.

Hale, J. and Kocak, H. (1991) Dynamics and Bifurcations, Springer-Verlag, New York.

Hairer, E., lserles, A. and Sanz-Serna, J. M. (1989)"Equilibria of Runge-Kutta Methods", Numer. Math.

Hsu+ C. S. 11987) Cell-to-CeU Mapping, Springer-Verlag, New York.

Hung, C. M., Sung, C, H. and Chen, C. L. 119911 "Computation of Saddle Point of Attachment", AIAA-91-

1713, AIAA 22nd Fluid Dynamics, Plasma Dynamics and Lasers Conference, Honolulu, Hawaii.

Humphries, A. R. (1992) "Spurious Solutions of Numerical Methods for Initial Value problems", IMA J.

Num, Anal.

Iserles, A. (19881 "'Stability and Dynamics of Numerical Methods for Nonlinear Ordinary Differential

Equations", DAMTP NA 1, University of Cambridge, Cambridge, England.

Iserles, A., Peplow, A. T. and Stuart, A. M. (1990) "A Unified Approach to Spurious Solutions Introduced by

Time Discretisation", Part I: Basic Theory, DAMTP 1990/NA4, Numerical Analysis Reports, University of Cambridge.

Jameson, A. 119911 "Airfoils Admitting Nonunique Solutions to the Euler Equations", AIAA-91-1625.

Keller, H. B. (1977) "Numerical Solution of Bifurcation and Nonliner Eigenvalue Problems", Applications qf

Bili_rcation Theory, P. H. Rabinowitz, ed., Academic Press, 359 384.

Kwak, M. 119911 "Finite Dimensional Inertial Forms for the 3D Navier-Stokes Equations," I MA preprint

Series #828.

l+afon, A. and Yec, H.C. (1991)"Dynamical Approach Study of Spurious Steady-State Numerical

Solutions for Nonlinear Differential Equations", Part Ili: The Effects of Nonlinear Source Terms

and Boundary Conditions in Reaction-Convection Equations, NASA TM-103877, to appear in Intern.

d. ('FD.

Lafon, A. and Yee, H. C. 11992! "'Dynamical Approach Study of Spurious Steady-State Numerical Solutions

of Nonlinear Dift'crential Equations", Part IV: Stability vs. Numerical Treatment of Nonlinear Source

Terms, ONERA-CERT Technical Report DERAT 45/5005.38, to appear in Intern. J. CFD.

Lambert, J. D. ( 19731 Computational Methods in Ordinary Differential Equations, John Wiely, New York.

Langford, W. F. and Iooss, G. (19801 "'Interactions of Hopf and Pitchfork Bifurcations", ISNM 54, Bifurcation Problems and Their Numerical Solution, Workshop on Bifurcation Problems and Their Numerical Solution, Dortmund, ed. H. D. Mittelmann and !1. Weber, Birkhauser Verlag, Basel, 103 134.

Panov, A. M. 119561 "'Behavior of thc Trajectories of a System of Finite Difference Equations in the

Neighbourhood of a Singular Point", Ueh. Zap. Ural. Gos. Unw. vyp, 19, 89 99.

Perron, O. (1929)"'Uber Stabilitiit und Asymtotischcs Uberhalten de L6sungen eines Systems endlicher

Differcnzengleichungen", 3'. Reine Anyew. Math., 161, 41 64.

Richtmyer, R. D. and Morton, K. W. (1967) Di_i, renee Methods [;)r Initial-Value Problems, InterscienceWiely. New York.

Sanz-Scrna, J. M. (19901 "Numerical Ordinary Differential Equations vs. Dynamical Systems," Applied

Math. Comput. Report 1990/3, Universidad de Valladolid.

Schecter, S. and Shearer, M. 11990)"'Undercompressive Shocks for Nonstrictly Hyperbolic Conservation

Laws", IMA Preprint Series #619.

Shearer, M,, Schaeffcr, D. G.. Marchesin, D. and Paes-Leme, P. ( 1987)"Solution of the Riemann Problem for

a Prototype 2 × 2 System of Non-Strictly Ityperholic Conservation Laws", Arch. Rat. Mech. Anal., 97,

299 321).

Strang, G. (1968)+'On the Construction and Comparison of Difference Schemes", SIAM ,I. Num. Anal 5,

506 517.

Swcby, P. K, Yee, tt. C. and Grifliths, Dr F. (1990) "On Spurious Steady-State Solutions of Explicit

Rnngc-Kutta SchemcsZ University of Reading, Department of Mathematics, Numerical Analysis

Report 390, also NASA TM 102819, April 1990.

Swcb_, P. K. and Yee, H.C. 119911 "On Spurious Asymptotic Numerical Solutions of 2 × 2 Systems of

ODEs", Numerical analysis Report 7/91, University of Reading, England.

Swcby, P. K. and Yce, H. C. 119941 "On thc Dynamics of Some (;,rid Adaptation Schemes", Proceedings of

the 4th International Conference on Numerical Grid Generation in CFD and Related Fields,"

Univcrsity Collegc of Swansea, UK, also RIACS Technical Report 94.02. Feb. 1994.

Temam, R 119891 "'Do Inertial Manifolds Apply to Turbulence?", Physica D, 37, 146 152.

Warner, B. (1980) "Turning Points of Branches of Positive Solutions", ISNM 54, B![urcation Problems and

Their Numerical Solution, Workshop on Bifurcation Problems and Their Numerical Solution,

Dortmund, ed. H. D. Mittelmann and H. Weber, Birkhauser Verlag, Basel, 211 226.

Yee tt. C. (1989) "'A Class of High-Resolution Explicit and Implicit Shock-Capturing Methods", VKI

Lecture Series 1989 04 March 6 10, 1989, also NASA TM-101088, Feb. 1989.

Yee, tt. C. I1991) "A Nonlinear Dynamical Approach to Algorithm Development in Hypersonic CFD",

Procccdings of the 4th International symposium on Computational Fluid Dynamics, Davis, Calif.
DYNAMICAL STUDYOF SPURIOUS STEADY-STATE NUMERICAL SOLUTIONS 283

Yee, H.C., Klopfer, G.H. and Montagne, J.-L. (1990) "'High-Resolution shock-capturing Schemes for

Inviscid and Viscous Hypersonic Flows", J. Comput. Phys., 88, 31 61.

Yee, H. C+, Sweby, P. K. and Griffiths, D. F. (1991)"Dynamical Approach Study of Spurious Steady-State

Numerical Solutions for Nonlinear Differential Equations", Part I: The Dynamics of Time Discretizations

and Its Implications for algorithm Development in Computational Fluid Dynamics, NASA TM102820, April 1990, also J. Comput. Phys., 97, 249 310.

Yee, H. C., Sweby, P. K. and Lafon, A. (1992)"Basins of Attraction and the Time-Dependent Approach to

Obtaining Steady-State Numerical Solutions", Proceedings of the 1CFD Conference on Numerical

Methods for Fluid Dynamics, Reading England.

Yee, H. C. and Sweby, P. K. (1993a) "'Global Asymptotic Behavior of Iterative Implicit Schemes", R1ACS

Technical Report 93.11, NASA Ames Research Center, to appear in International J. of Bifurcation and

Chaos, Dec., 1994.

Yee, H. C. and Sweby, P. K+ 11993bi "On the Dynamics of Some lterative Implicit schemes", Proceedings of

the Chaotic Numerics Workshop, Deakin University, Geelong, Australia.

Color Plate I Figure 6.1 (See H.C. Yee and P. K. Sweby.)
Color Plate II Figure 6.2 (See H.C. Yee and P. K. Sweby.)
Color Plate III Figure 6.3 (See H.C. Yee and P. K. Sweby.)
ColorPlateIV Figure6.4(See H.C. Yee and P. K. Sweby.)
_J

r_

_J

C_

ru00
Color Plate VI Figure 6.8 (See H.C. Yee and P. K. Sweby.)
Color Plate VII Figure 6.9 (See H.C. Yee and P. K. Sweby.)
Color Plate VIII Figure 6.10 (See H.C. Yee and P. K. Sweby.)
_P

v_

.,..i

_S

rl

0

N

0
Color Plate X Figure 6.12 (See H.C. Yee and E K. Sweby.)
Color Plate XI Figure 6.13 (See H.C. Yee and E K. Sweby.)
c_

Q_

_J

v_

=

L

0
Color Plate XlII Figure 6.15 (See H.C. Yee and R K. Sweby.)
Color Plate XlII Figure 6.15 (Continued) (See H.C. Yee and P. K. Sweby.)
Color Plate XIV Figure 6.16 (See H.C. Yee and P. K. Sweby.)
Color Plate XV Figure 6.17 (See H.C. Yee and R K. Sweby.)
P_

c_

c_

c_

0

N

L

0

m

0
Color Plate XVI Figure 6.18 (See H.C. Yee and P. K. Sweby.)
¢J

._,

0

g

r,

0
Color Plate XVII Figure 6.19 (See H.C. Yee and P. K. Sweby.)
_J

L. 0
Color Plate XIX Figure 6.21 (See H.C. Yee and P. K. Sweby.)
ColorPlateXX Figure6.22(See H.C. Yee and E K. Sweby.)
Color Plate XXI Figure 6.23 (See H.C. Yee and P. K. Sweby.)
"9.

C_

¢;

r_

0

(0

ee_

.._

r,

f..

@
Color Plate XXII Figure 6.24 (See H.C. Yee and R K. Sweby.)
_J

Lt_

c_Q
CON I mill t E v,E i ,_

D. 26_00

n 7OOO0

ii [lo01Dfl

tl :_;'fl(lll

:!:!i

::l;:r,!

_ _!Lq_j_

, , hib_h_i

Ii 'l/llllll

In q|_l J!.ll I

I • tluut_

I • I_?llOI1

I ii.tllnfl

I ilalloll

i i nllilo

I i ;,nNfl

Color Plate XXIV

Figure 10 Computed surface pressure contours with store located in captive position

and moving through 0.6, 1.0, 1.6, and 2.0 store diameters below the pylon.

(See A. Arabshahi and D. L. Whitfield)