On the Solution of Convection-Diffusion Boundary Value Problems Using Equidistributed Grids.
The effect of using grid adaptation on the numerical solution of model convection-diffusion equations with a conservation form is studied. The grid adaptation technique studied is based on moving a fixed number of mesh points to equidistribute a generalization of the arc-length of the solution. In particular, a parameter-dependent monitor function is introduced which incorporates fixed meshes, approximate arc-length equidistribution, and equidistribution of the absolute value of the solution, in a single framework. Thus the resulting numerical method is a coupled nonlinear system of equations for the mesh spacings and the nodal values. A class of singularly perturbed problems, including Burgers's equation in the limit of small viscosity, is studied. Singular perturbation and bifurcation techniques are used to analyze the solution of the discretized equations, and numerical results are compared with the results from the analysis. Computation of the bifurcation diagram of the system is performed numerically using a continuation method and the results are used to illustrate the theory.  It is shown that equidistribution does not remove spurious solutions present on a fixed mesh and that, furthermore, the spurious solutions can be stable for an appropriate moving mesh method.
Introduction
The solution of boundary value problems (BVPs) using finite difference
methods consists of discretizing the differential equations on a finite mesh
and solving the system of equations for the values of the functions at these
mesh locations. Adaptive meshes are widely utilized for problems containing
strong gradients and for problems where important features of the solution
develope on small length scales. References for applications in computational
fluid dynamics are given in [23] and [16], [17],[18], [25] give more
general references for the application and analysis of mesh apaption. Adaptive
meshes are essential if accurate and reliable estimates of the solution
structure are to be made and much effort is put into designing suitable
meshes for complex geometries. However, the analysis of the reliability of
adaptive meshing strategies with regard to giving accurate solutions and rejecting
spurious ones is still in its infancy. The fundamental question we ask
is this: how well do the qualitative properties of the discrete system with
mesh adaption reflect the qualitative properties of the continuous system it
is approximating.
Our studies of this question will focus on adaptive numerical discretizations
of reaction-convection-diffusion BVPs of conservation form:
where
f(u;
We look at both steady states and time evolving solutions. The steady-states
satisfy the identity
where k is a suitable constant. We note that even if f is a linear function
of u, the addition of an adaptive strategy to solve (1.1) typically leads to a
nonlinear discrete system.
An important example of such an equation is the viscous Burgers' equation
where
f(u; (1:1),(1:2). This is an important test problem for many methods in computational
fluid dynamics and has a unique, monotone decreasing and symmetric
steady solution given by
where m is a constant given by the unique solution of
This steady solution is a global attractor for the time dependent problem.
As F the steady state of (1.1) is invariant under the action of
the map
All continuous steady solutions of (1.1) under (1.2), are necessarily symmetric
and satisfy
so that
If " is small then there is a sharp transition layer in the solution where u(x)
changes sign. This layer is centred at 1and is of O(") in width. Outside
the transition layer juj is close to A.
In contrast, when the resulting inviscid equation has many steady
weak solutions which satisfy u are given by
For these solutions the transition layer is at x 0 where x 0 is arbitrary. Only
one of these is the weak limit of a viscous solution. However, exponentially
small (in ") changes in the symmetry of the boundary conditions of the
viscous problem can lead to solutions with significantly displaced transition
layers. This phenomenon is referred to as super sensitivity in the papers of
[12], [13], [14] and [24]. For the time dependent problem, if " is small, a
solution will rapidly evolve to one with u x close to zero almost everywhere.
This solution again has a thin transition layer of width O("), where u changes
sign, which is centred at a point determined by the initial conditions. Over
much longer timescales this transition layer moves slowly towards the point
and so the solution converges to the symmetric steady state. This
phenonenon is often called meta-stability [3].
When a method using a fixed mesh at points x
adapted to the solution) is used to find the steady states of (1:1) it is well
known [7] here that this may give rise to spurious solutions. Such solutions
frequently place the transition layers in an incorrect position. Typically,
however, when working with problems that have sharp transition layers, an
adaptive strategy is used which places mesh points x j into regions where the
computed solution has a high gradient. Such strategies are hoped to produce
reliable answers and to preserve qualitative properties of the solution such
as symmetry and monotonicity. In this paper we consider strategies which
move a fixed number of mesh points so as to equidistribute a positive monitor
function of the solution. A natural question which we ask is: can such an
adaptive stategy for finding the steady solutions of (1.1) place the transition
layer in the correct position and give monotone solutions? Briefly the answer
to this question is "no, in general". Indeed, we show that even when the
transition layer is well resolved by the mesh, (so that many points are moved
into the region where u changes sign), it can still appear in entirely the wrong
position.
To demonstrate this we solve for the steady state of (1.1) using an adaptive
strategy in which a fixed number of J mesh points x j are placed in such a
manner that a generalisation of the the arc-length of the solution is equidistributed
over each interval [x sufficiently large compared with
1=J then the numerical method has a unique, monotone, symmetric solution
with a transition layer at reduced then asymmetric
solutions bifurcate from this branch when . Close to
" SB then both the true and spurious solutions are monotone decreasing in
x. These asymmetric solutions persist as " is reduced. Further asymmetric
solutions arise at subsequent symmetry breaking bifurcations from the
main branch. For small " these asymmetric solutions are perturbations of
the many weak solutions arising in the continuous problem when
and they exhibit transition layers placed almost arbitrarily in the interval
[0; 1]. Similar behaviour (including the existence of a related value " SB ) is
observed for the case of a uniform grid. Thus a unique numerical solution
which correctly places the transition layer exists only . We make
some estimates of " SB and show that its value is not significantly reduced
by using the arc-length based adaptive strategy. However, we show that
if a strategy is used which ensures that the estimate ju
equidistributed by the mesh, then the resulting solution is both monotone
decreasing and symmetric for all values of ": However such an equidistribution
method is typically not robust and will not perform well in some
computations due to the non-smoothness of the monitor function.
As a second investigation we look at the time dependent problem (1.1)
using an adaptive method which moves the mesh points x j (t) as the solution
evolves in time to again equidistribute its arc-length. The discrete solutions
of the steady state of (1:1) using the adaptive mesh are steady states of this
method. We show through some numerical experiments that the discrete
solution with the transition layer placed at
and that for a range of initial data, the adaptive solution converges to a
problem with an asymmetrically placed transition layer. Thus the spurious
solutions have a marked effect on the dynamics of the computed profile, with
the meta-stability of the asymmetric states of the continuous problem being
replaced by the stability of these states in the discrete problem.
The basic reason for this rather peculiar behaviour, is that in the continuous
problems super sensitivity implies that the location of the transition layer is
uniquely determined by exponentially small perturbations to the boundary
conditions. For small " the effect of these perturbations is easily overcome
by the errors introduced by the discretisation process. An adaptive method
which places mesh points in regions of high gradient only does not seem to
affect this process. An interesting question is whether the correct location
of the transition layer could be obtained by using a numerical method which
places points close to the boundary, and we leave this for future investigation.
We expect that the ideas in this paper are also relevent to problems which
are only algebraically sensitive to the small diffusion parameter, but that
then the effects will be less dramatic.
The remainder of this paper is organised as follows. In Section 2 we derive
a suitable numerical method for the steady state of (1:1) which incorporates
a mesh-adaption procedure based upon a generalisation of the well-known
arc-length method. In Section 3 we analyse this method when applied to
the examples given. In particular we determine the nature of the (discrete)
solution branches as parametrised by " and estimate the location of bifurcation
points where the non-spurious computed solution loses stability. We
also study the form of the solution for " small. In Section 4 we describe some
numerical calculations of the steady state branches and compare them with
the analytic estimates in Section 3. Finally in Section 5 we make some time-dependent
calculations of the solutions of (1:1), based upon the moving-mesh
strategy of [18], and show how the spurious steady-state solutions affect the
underlying dynamics of the problem.
Adaptive Mesh Discretizations of the Steady BVP
Our derivation of discretization of the steady state of the equations is based
upon the conservation form
where the constant k is to be determined. In Section 5 when we look at the
time dependent formulation.
Let U be the discretized solution variable u and H the vector of mesh-
spacings so that
Thus the h j determine the (non-uniform) mesh distribution and J is the
number of grid cells. Hence x
our approximations
to u(x j ). Discretizing the conservation equation over the interval [x
and taking the mean of F (u) over this interval we obtain,
"
Although we will mainly work with this discretisation, we observe that subtracting
two succesive such equations gives the quivalent formulation
'
=h =2
If the mesh is specified (2.2) represents (J equations for (J unknowns
U and k. Our interest, however, lies in the case where the mesh
is adaptive so that the J mesh parameters forming the vector H are calculated
as part of the solution. Thus we must specify further equations for the
mesh. A common approach to mesh adaption is to use an equi-distribution
principle see [11], [16], [17], [18], [25]. In these methods, there is a function
M(u; commonly called a monitor function such that, if x i is the
location of the i th mesh-point (with x
M(u;
In practice, of course, the integrals are approximated on the mesh.
Different monitor functions lead to different meshes with varying properties
and the 'correct' choice of monitor function is not always clear. In [4] a
monitor function for a certain two-point boundary value problem is shown
to give 'optimal' accuracy in a suitable norm. Alternatively, it is proposed
in [1], that for certain partial differential equations invariant under various
rescaling groups, the most appropriate monitor function to give the correct
dynamics of the solution is one which generates meshes which scale in a
similar manner to the exact solution. These various considerations notwith-
standing, the most commonly used monitor function is arc-length for which
M(u;
see [25], [23]. Here we consider a monitor function which is a generalisation
of this one so that
M(u;
This formulation reduces to arc-length equidistribution if
uniform mesh if equi-distribution with respect to ju x j if
varying ff we may relate the properties of various mesh-adaptation
principles, including the limiting case of a uniform mesh.
To equidistribute the mesh we thus aim to choose H such that
x
x dx:
Approximating the arc length to be constant within the interval yields,
A natural discretization of this equation gives
equations as j ranges from need
the additional constraint that
Thus (2:4), (2:5) give a further J equations for the unknowns Coupled
with (2:2) this gives us a system of 2J equations and unknowns. Since
(2:4) is independent of j, the quantity on each side is a constant s 2 . Thus,
Thus (2:2), (2:6) form 2J equations in the 2J
and s. This set of equations can be written as
By a solution we will always mean a real vector
In section 5 we generalise these ideas to include both the time dependence
of u and a dynamic strategy for moving the mesh points.
3 An Analysis of the Discretized Model Equations
To analyse the solutions of (2:2),(2:6) considered as functions of " and ff
we apply techniques from bifurcation theory, determining the existence of a
symmetric solution branch and then looking for asymmetric solutions bifurcating
from it as " varies. We can then compare the qualitative properties
of the discrete solution with those of the continuous solution. To start our
investigation we take " to be large and establish the existence of a unique a
monotone decreasing and symmetric solution of the discrete equations. This
we refer to as the main solution branch and we establish some general properties
of the symmetric and monotone decreasing solutions of (2:2),(2:6). As "
is reduced the main solution branch remains symmetric but loses monotonic-
ity. However, further asymmetric solutions of (2:2),(2:6) arise. In particular,
we show that there are symmetry breaking bifurcation points on the main
branch leading to such asymmetric solutions. It is interesting that (as " is
reduced) these occur before the main branch loses monotonicity. If
however, then the smmetry breaking points do not exist. Finally, we apply
the Implicit Function Theorem to show that provided " is sufficiently small
then we may construct many solutions of the finite difference scheme, some
of which may have arisen at the symmetry breaking points.
3.1 The solution for " AE 1
If we consider solutions which are bounded in the limit of "
then in this limit equation (2:2)
reduces to
Using the boundary conditions, u
Therefore, to satisfy the last equation in (2:2) and (2:5) we require that
Thus, to leading order in 1="
To obtain the location of the mesh points we substitute these identities into
the equation (2:6) to give
Thus h j is a constant which must therefore equal 1=J . Hence we have
Hence, the unique bounded asymptotic solution of (2:2) for " !1 is given
by:
This is , in fact, the exact solution of the steady solutions of (1.1) at
sampled on a uniform grid.
3.2 Properties of the Main Solution Branch
The solution constructed above for " AE 1, has a uniform mesh and is both
monotonic so that u j ? u j+1 and symmetric so that
Both the symmetry and monotonicity of the discrete solution are in qualitative
agreement with the true solution of the continuous problem. As
F (u) satisfies the symmetry condition F reduced
the equations (2:2),(2:6) will continue to admit a symmetric solution satisfying
(3.2). For large enough " this solution will also be monotonic with
As " is reduced then, in general, there will be a first value of
" at which monotonicity is lost and a possible further value of " at which
asymmetric solutions (not satisfying (3.2)) bifurcate from it.
We now study how the monotone and symmetric solutions vary as " re-
duces. In particular we examine the way in which the mesh changes from
being uniform and obtain some bounds for the values of h j ; these are useful
subsequently in finding bifurcation points. For this we assume that the
function F (u) satisfies the following conditions
where G(u) is bounded for finite u. This includes Burgers' equation for
which
Lemma 3.1 Let
be a a symmetric, monotone-decreasing
solution of (2:2),(2:6) satisfying (3.2) with
Suppose further that F (u) satisfies (3.3),(3.4) and ff ? 0. Then, for
we have
and
Corollary 3.2 Under the same conditions as in Lemma 3.1, it follows that,
and q
Proof of Lemma 3.1
It follows from (2:2) that for a monotone decreasing solution, the function
a
is negative and from (3.3) it is monotone decreasing if j ! J=2. Consequently
a 2
j is monotone increasing if j ! J=2. Now combining (2:2),(2:6) we have
ff
a 2
It follows immediately if ff ? 0, that as s is independent of j then as a 2
j is
increasing then h j+1 is decreasing in j. We immediately deduce from (2.6)
that increasing in j. If
and h j is constant, we deduce that as a 2
j is increasing then so is \Delta j . This
proves the lemma. 2.
Proof of Corollary 3.2 Inequality (3:7) follows simply from the observations
that h j is symmetric in j, is decreasing for J=2 and that
Inequality (3:8) is similar and uses the elemetary identity that
To obtain the estimate (3:9) we make use of the geometric result that
the least arc-length of the solution occurs when the points u j lie on a
straight line of uniform gradient joining the two end-points, which has length
Similarly, the greatest arc-length occurs when the solution
is two horizontal line-segments through joined
by a vertical line-segment through x 2 (0; 1): Such a curve has arc-length
In general the solution does not stay monotonic as " is reduced. As a final
result on the monotonic solutions we give an estimate for the value of " at
which they lose monotonicity
Lemma 3.3 If the symmetric solution first loses monotonicity at
" decreases, then at
Consequently, there exists a u   2 [u 2 ; A] with
Proof The solution first loses monotonicity when As
is monotone increasing by Lemma 3.1, this must occur first when
that A. It then follows from taking in (2:2) that (A). The
first result then follows by considering (2:2) when and the second by
dividing both sides of the identity by u and applying the Mean Value
Theorem. 2
3.3 Symmetry Breaking Bifurcations
We now show that (2:2),(2:6) has further asymmetric solutions which bifurcate
from the symmetric branch described above. Hence the uniqueness
and symmetry of the solution is lost as " reduces. In particular we show
that there are symmetry breaking points at which asymmetric solutions (not
satisfying (3.2)) bifurcate from the symmetric branch at a non zero value of
". To establish this we investigate the linearisation of (2:2),(2:6) about the
symmetric solution and establish the existence of an asymmetric null eigenvector
at a specific value of " SB . This is the condition for a point on the
symmetric branch to be an infinitesimal symmetry breaking bifurcation.
Such a point is actually a symmetry breaking bifurcation point provided
that certain non-degeneracy conditions are satisfied. Generically this is so,
and we observe (numerically) that additional spurious solutions do indeed
bifurcate at these points. Furthermore, these spurious solutions appear to
persist for all " ! 0. Owing to the action of the symmetry group, all such
bifurcations are necessarily pitchfork bifurcations and lead to two solution
branches in a neighbourhood of " SB such that
Numerical calculations demonstrate that u \Gamma and u only exist for
The value of " SB is determined in the following theorem.
Theorem 3.4. If F (u) satisifies (3:3),(3:4) and ff ! 1, there is a symmetry
breaking bifurcation on the symmetric solution branch of (2:2),(2:6) at a
SB satisfies the equation
Equivalently
Further infinitesimal symmetry breaking bifurcations from the main symmetric
branch will occur for (smaller) values of " (j)
SB which satisfy
Lemma 3.5 If F (u) satisfies (3.3),(3.4) and ff ! 1 then there must be a
value of " satisfying (3:11). Such a value satisfies
is the greatest value of " at which the main solution branch loses
monotonicity.
Corollary 3.6. Consider a solution of the discrete Burgers' equation with
(i) If
Furthermore2J
(ii) For all ff 2 [0; 1), " (1)
SB satisfies the following inequalities
Note The particular case corresponds to a solution of (1.4) using
a non-adaptive mesh and (3:13) was first derived in [20] using a different
method from the one we shall now give. The estimates (3:14),(3:15) are
new. In Section 5 we compare these bounds with numerically computed
values of the symmetry breaking point.
Corollary 3.7. For the resulting solution of the discrete Burgers' equation,
if
then an infinitesimal symmetry breaking bifurcation has occured on the main
branch for a value " SB ? ":
We can draw some useful conclusions from the previous lemmas. If we first
consider a mesh without adaptivity, numerical evidence indicates that if " ?
1=2J then the computed solution is non spurious and if " ! (1\Gamma2=J)=2J then
spurious solutions exist in agreement with (3:14). Both of these estimates
are reasonable. The natural length scale of the transition layer is O(") and
to obtain an adequate resolution of the transition we must have a mesh such
that 1=J is of the same magnitude as ". Thus spurious solutions are likely
to exist if this is not the case and we have to take a minimum value of J
to limit the effects of spuriosity. For the value of " in (3:16) it is likely
that there are additional non-symmetric solutions to the original equations.
This result (perhaps rather surprisingly) indicates that using an adaptive
mesh does not significantly improve on this result and a comparable value
of J may have to be taken to avoid spuriosity. This is true even though the
adaptive method may place many more points in the actual transition layer
than the non-adaptive one. Note, however, that the case
such spurious bifurcations. This is examined further in Section 3.5.
Proof of Theorem 3.4
To prove the theorem we consider a linearisation of (2:2),(2:6) about a symmetric
solution. For the original system to have a bifurcation point it is
necessary that the linearisation should have a non-zero null-eigenvector ffi
with
corresponds to the linearization of u j , g j to h j , t to s and l to
respectively. A linearization of (2:2),(2:6) implies that these components
satisfy the following system of linear equations:
st
together with the further boundary condition that
As
that if l admits an anti-symmetric solution which
breaks the symmetry condition (3:2) and satisfies
By [10] it follows that if such a solution is non-zero then " SB is an in-
finitessimal symmetry breaking bifurcation point, and asymmetric solutions
of (2:2),(2:6) bifurcate from the main branch at this point. Thus we must
show that there exists such an " where a solution of the form (3:21) is indeed
non-zero.
Setting in (3:17),(3:18) and taking
and
The reduced system (3.22),(3.23) has a non-zero solution for which w 1 6= 0
and
det
0: (3.24)
Calculating the determinant in (3.24) and simplifying the resulting expression
using the identity
gives the expression (3.10). On substituting the definition of s we have the
equivalent condition
We can now construct the resulting null-eigenvector satisfying (3:21). In
this case, the equation for w 2 and h 2 is given by
!/
This problem will have a solution (w that the matrix in (3:26)
has a non-zero determinant that is
which, on rearranging is equivalent to the condition
Here we note that the equivalent term D 1 is identically zero at the bifur-
cation. We shall show presently that a bifurcation first occurs when u
is monotone decreasing and h i decreases monotonically if i ! J=2 thus if
immediately follows that as G(u) is monotone that D 2 6= 0.
Hence a solution w of (3:26) exists. By a similar procedure we may
calculate the subsequent terms (w
We note that if D 1 6= 0 but D then a further bifurcation will occur with
a null-eigenvector for which w Thus a series of
symmetry-breaking bifurcations occurs when each of the terms D i is equal
to zero. This completes the proof. 2
We now show that such a value of " must exist.
Proof of Lemma 3.5
We look at the quantity
By the estimates of the last section, and 1=s are bounded above as
" !1 and hence D 1 is negative for large ". On the contrary, as " decreases,
then if the solution remains monotone, these values are also bounded below
(and are positive) so that D 1 is positive for small ". Thus there must exist
a value of " at which D Now suppose that the main solution loses
. At such a point U and hence from (2:6)
we have that Thus at this point
Now, from Lemma 3.3 we have
where u   ! A. Thus, as G(u) is monotone increasing in u and h
have that D 1 is positive if . Thus D 1 which first becomes zero at a
where the solution is still monotone. 2
Proof of Corollaries 3.6 and 3.7.
The results in Corollaries 3.6 and 3.7 follow immediately from an application
of the estimates of u Corollary 3.2 to the identity (3:11). 2
3.4 Asymptotic Solutions When "  1
The symmetry breaking bifurcation generates solutions additional to the
main branch for small values of ". In fact, for very small values of " there
are a very large number of solutions of (2:2),(2:6) As stated in the introduc-
tion, these are numerical perturbations of the many weak solutions which
exist when closely related to the super-sensitivity of the underlying
problem. To determine these we firstly recast (2:2) in the form
We show that the combined system (3:29),(2:6) has solutions when
and use the Implicit Function Theorem to deduce the existence of solutions
perturbations of the
Setting
together with the mesh defining equations (2:6) and the boundary condition
We see immediately that this system has a solution if either
or
The values of j for which h correspond to points in the solution where
there is a transition layer. We shall assume that the transition layer occurs
for of the mesh. In the transition
layer, it follows from (2:6) that if h decreasing solution
satisfies
For values of j for which j and are not in the transition layer we have
from applying (3.32) twice that
and hence from the monotonicity of F (u) we deduce that
Thus there is a constant ffi for which
A solution of this form is illustrated in Figure 1.
From this result and from (3.32) we may deduce that, since u
These results imply that h j is constant over all intervals outside the transition
layer and takes the value
so that
By (3:33) we have
and this gives another relationship between s and ff since uM and u L are
determined by the parity of L and M:
The solution of the system when characterised completely by
the values of H ,ffi and s each of which is given uniquely by the values of
and M and (hence of J T ). Indeed, the solution is independent of the
function F (u) We observe that if L is even then u if L is odd
then similar results for uM . This gives us three cases of
solutions to consider.
Case 1: L even & J-M even.
In this case it follows from (3.35) that u \GammaA. Thus summing
the expression (3.33) over the transition layer we deduce that
Substituting this expression into (3.37) gives the following expression for
s
The expression (3.37) is valid provided that
0: (3.41)
For example, when implies that
J:
We see that for this example the values of ffi; s and H are (up to a sign of
uniquely by the value of J T and do not depend upon the precise
values of L and M . Thus the solution derived can have a transition layer
anywhere in the interval, provided that L is even.
Case 2: L odd & J-M odd. We have u so that
ff:
Substituting into ??3.45 implies that ffi satisfies the quadratic equation
In general this equation has two solutions, and either may be observed. these
in Section 4.
Case 3: All combinations not included above. The analysis for this case
follows similarly and we have that ffi satisfies the quadratic equation
We now turn our attention to the perturbations of the above solutions when
" is small. To do this we define the solution when to be U
etc. so that H in the transition layer. As u j is everywhere bounded, it
is clear that the full problem described by (3:29),(2:6) is an O(") perturbation
of that of (3:30),(3:30). Hence an application of the Implicit Function
Theorem implies that a solution of (3:29),(2:6) will take the general form
outside the transition layer (3.45)
and
inside the transition layer. (3.46)
We state the full result as follows.
Theorem 3.8 Let ff 2 [0; 1). In each of the cases 3, the solution of
perturbs smoothly in the manner defined by (3.44),(3.45),(3.46)
provided that the corresponding variable \Gamma i defined in (3:41) for case 1, (with
similar expressions for the other two cases) satisfies the condition
Furthermore, if (3.47) is satisfied then to leading order in ", the mesh in the
transition layer takes the form
where take the values calculated in the three cases above.
Proof
To prove that we can apply the Implicit Function Theorem and deduce this
result we must show that the operator derived by linearising (3:30),(2:6)
about the solution (U equivalently, that it has
no (non-zero) null eigenvector. The proof of this is straightforward but
technical and details of this calculation are given in are given in [2]
3.5 Existence and Uniqueness of Solutions for
When we have from (2:6) that ju Consequently, if the
solution is monotone, we have that the value of
Lemma 3.1 is constant. As the solution can only lose monotonicity when
we deduce that a loss of monotonicity can only occur if
for all j. Thus u j is constant for all j, violating the boundary conditions.
We deduce that the main solution branch remains monotone for all ". As
" tends to zero the only possible monotone solution of (3:30) is one for
which otherwise and there is a transition layer
with 1: Similarly, from (3:25) we see that a symmetry
breaking bifurcation can only occur when
the analysis in the previous section does not rule out the existence of non-monotone
spurious asymmetric solutions when small values of "
we can deduce that such solutions can not bifurcate from the main branch.
More significantly they do not alter the stability of the main branch. Thus
taking reduces the impact of the spurious solutions. We shall see
numerical evidence for this in the next section.
4 Numerical Results for the Examples
This section presents the numerical solution of (2:2),(2:6) obtained through
continuation in both " and ff using the package AUTO [6]. In all cases J
was kept fixed at
To compute a solution branch parameterised by " we take as a starting
solution the monotone symmetric solution computed in Section 3.1 with
10. The value of " is then reduced and all bifurcations and bifurcating
branches computed. At the fixed value of \Gamma3 we also look at the
solutions as functions of ff. Further, disconnected, branches can be found
for very small values of " by using the asymptotic solution (when
a starting solution when " is small.
4.1 Continuation in "
We take ff to be fixed at the two values of 0; 1=2 and compute the bifurcation
diagrams in the two cases. As a suitable measure of the solution we take
the value of u 5 as this is effective in distinguishing between solutions with
different locations of the transition layer.

Figures

2 and 3 then show the corresponding bifurcation diagrams of the
solution of (2:2),(2:6). In Figure 4 we also present the value of h 1 as a
function of " in the case of ff = 1=2: In this figure we see that h 1 tends to
one of the three limits of 1=7; 1=8 and 0 as " tends to zero.
It is clear from Figures 2 and 3 that when there are symmetry
breaking bifurcations from the main branch, as predicted by Theorem
3.4. leading to solution branches which are asymmetric. It is also clear that
these solutions tend toward an asymptotic limit as " tends to zero. Using
the notation of Section 3.4 we may thus characterise these solutions by the
values of M and L such that transition layer occurs for j between L+ 1 and
M . These values are indicated on Figure 2. It is interesting to note that the
solutions which have the transition layer at and at
are the continuations of the two solutions which bifurcate
from the main solution branch at the greatest value of ". Subsequent
solution branches which have transition layers closer to
from the main branch at smaller values of ": In Figures 5 and 6 we present
the complete set of solutions for
These clearly show the oscillatory structure predicted in Section 3.4.
If then the first bifurcation occurs when
bifurcations occur when
. Thus the value of " for the first of these bifurcation points
actually increases when the adaptive procedure is used. Presumably, this
is because fewer points are now placed close to the boundary, reducing its
effect on determining the location of the transition layer. If then the
location of the first of these bifurcation points can be estimated from (3:14)
which predicts that 0:004787
2:45. This is consistent with the
above computations. In Figure 7 we give a graph comparing the estimates
of (3:14) with the actual location of the bifurcation point. It is interesting
to note from this graph that " SB has a maximum value at about
before decreasing to zero as ff tends towards one.
If we now take fixed for all ". However, the value of h 1
increases towards 1=2 as " tends towards zero. This behaviour can be seen
in

Figure

8 which plots h 1 as a function of " in this case. As predicted in
Section 3.5 we see that there are no bifurcations from the main solution
branch. The resulting profile when \Gamma3 is given in Figure 9. From
this figure we can see that the transition layer occurs at the right location
but the grid points are concentrated in a small region.
4.2 Continuation in ff
For continuation in the parameter ff we have taken one of the solutions at
\Gamma3 from Figure 3. as a starting solution. The bifurcation
diagram with ff as the parameter is shown in Figure 10 for which u 2 is
a convenient value to plot. This again shows a complex form with many
solutions, although the structure of the bifurcations in this case is less clear
than the continuation with ":
4.3 Asymptotic Properties of the Solutions
As remarked above, it is clear that as " tends to zero that the computed
solutions are converging to the asymptotic solution predicted in Section 3.4.
We now make this comparison more precise. The solutions at
from two different branches of the bifurcation diagram are shown in Figure
11. In the figure, solution 1 has and represents
the symmetric solution on the main branch. In contrast, solution 2 has
14 and has a maximal degree of asymmetry. Since
solution 1 corresponds to the case 1 and Solution 2
corresponds to the case 3 described in Section 3.4.
To verify the asymptotic analysis for Solution 1, we determine j j u1 \Gammau 0from the numerical computations and compare it with the value of ffi given
by (3:40). This analysis implies that if
so that in general
\Gamma3 and at
\Gamma5 . These values show that which is consistent with
the asymptotic value.
Solution 2 is an example of case 3 of Section 3.4. For this solution the
equation (3:43) yields two values for ffi namely
Clearly from the figure it is the latter value which is seen as a solution
bifurcating from the main branch. Presumably the other solution appears
at a fold bifurcation.
5 Time Like Iterations to Steady State
We now consider the implications of our steady state analysis for the solution
of the full time dependent problem (1:1), in particular determining which
of the solutions constructed is stable. To do this we augment the partial
differential equation with additional equations allowing the mesh points to
evolve and for which (2:6) represents a steady state. Thus we can study
the dynamical effect of the steady solution computed in Section 3. We
consider meshes which evolve to equidistribute the quantity M(u; x) and
for simplicity in this section we only consider discretisations of Burgers'
equation with A = 1.
To apply such a moving mesh technique we follow [11] and introduce a new
computational coordinate  2 [0; 1] and a mesh function x(; t) so that the
mesh point is at x(i=J; t). We then introduce a second partial differential
equation to describe the evolution of x. As the value of x is now a function
of time, we must recast Burgers' equation in a moving coordinate system,
so that
To evolve the mesh we then couple this with the following moving mesh
equation given in [11]:
is a relaxation parameter, M(x(; t); t) is the quantity to
be equidistributed and y() determines the initial mesh distribution. In this
section we take
From (5:2) it can be seen
that the solution equidistributes the quantity M if it reaches steady state.
The time evolution is introduced as a fast relaxation to equidistribution
avoiding the solution of differential algebraic equations for the evolution
problem. Our purpose in this section is to study the stability and basins of
attraction of equilibria under this natural evolution process. For the stability
properties of the equation (5:2) see [17]. Other time-like iteration schemes
for the moving mesh equations are studied in [22].
5.1 Discretization of the Equation
Following [11], equations (5:1) and (5:2) are discretized using the method
of lines. The spatial derivatives are, as before, approximated using central
differences and the resulting stiff system of ODEs, for the solution variables
and mesh locations, is solved using a stiff solver.
We set
h l (t)
and define
The equation (5:1) is discretized on a nonuniform mesh in x. to obtain
dt
Here, the operator P is as given in (2:3) so that
and we take a central difference discretisation for the convective operator to
give
The equations (5:2) are then discretized on a uniform mesh in , to give
where
s
This gives
dX
dt
where A is the inverse of a tridiagonal matrix with diagonal elements \Gamma2
and off-diagonal elements 1, and R is a th element
This can also be written as
dt
A jk R k
Using this in (5:4) we get
dt
A jk R k (U;
Then the equations (5:6) and (5:7) represent a system of ODEs,
dZ
dt
the system for the evolution in time of u and x at uniform grid locations
in . At steady state, these equations reduce to the mesh adaptation
introduced earlier, in the case
5.2 Numerical Results of the Unsteady Computations
The solution branches in the bifurcation diagram of Figure 2 represent the
equilibrium points of the system (5:8). Thus the stability of these solution
branches can be analysed using a linearized stability analysis of the system
(5:8). This analysis can be performed by AUTO and the resulting stability
properties of the branches when are sketched in Figure 12 where the
solid lines represent the stable branches and the broken lines represent the
unstable ones. The main branch undergoes supercritical and subcritical bifurcations
alternately at the sucessive symmetry breaking bifurcation points
In particular it loses stability as " is reduced below "
Thus the effect of the spurious solutions is to destabilise the main symmetric
solution branch. The secondary branch bifurcating from the third bifurcation
point at " loses stability through a Hopf bifurcation at
\Gamma3 . Since the underlying partial differential equation cannot
undergo a Hopf bifurcation, the oscillations resulting from the numerically
observed Hopf bifurcation arise directly from the adaptive mesh strategy.
To study the stability characteristics of the branches under large perturbations
we now consider evolution of the system (5:8) from specific initial
data using the Backward Euler method. We have computed the solution
at . For this value of " there are seven equilibrium solu-
tions, one on the main branch and two on each of the secondary branches
bifurcating from " . In this case the branch and the secondary
branch bifurcating from " 2 are unstable, whereas the branches emanating
from bifurcations at " 1 and " 3 are stable. Each of the stable branches has
a non-zero basin of attraction and it is of interest to ask which initial data
is attracted to which branch. To gain insight into theis question we chose
as initial data for the unsteady computation, small perturbations of the
unstable equilibrium solutions. We now define
where M j+ 1is defined as in (5:5). Thus
of the deviation of the solution from equidistribution. For all our computations
we find that increases initially as the solution moves away
from the unstable branch. However, once the solution is close to a stable
branch the deviation from equidistribution decays exponentially. The convergence
history of the unsteady computation with starting solutions close
to the unstable branches at \Gamma2 is shown in Figure 13. To obtain
the results summarised in Curve 1, initial data for both U i and for the mesh
is taken which is a perturbation of a symmetric solution on the main branch.
To obtain the results in Curves 2,3, two nearby sets of initial data are chosen
which are both perturbations of the asymmetric unsymmetric unstable
solution bifurcating from the main branch at " 2 but which lead to different
time histories. In these time histories the solutions represented by Curves
1 and 2 both coverge slowly toward the stable secondary branch emanating
from " 3 . The oscillations in these solutions, and the slow convergence, is due
to the existence of complex eigenvalues of the linearized system which have
real part close to zero. This follows from the fact that " is close to the Hopf
bifurcation point at " H . In contrast, Curve 3 exhibits rapid convergence to
the secondary branch emanating from " 1 . Thus the secondary curve bifurcating
from " 2 in some sense divides the initial data which evolves toward
the two equilibrium solutions.
6 Conclusions
The results of our analysis and of our computations show that the use of
an adaptive strategy does not prevent the existence of spurious asymmetric
solutions of the discretisation of the partial differential equation (1.1).
Furthermore, these solutions can be stable whilst the main symmetric solution
branch (which is the approximation to the true solution) is unstable.
Hence, they have a profound and misleading effect upon the dynamics of
the solution. If we measure the 'goodness' of an adaptive strategy by the
smallness of the value " 1 of " at which the first spurious branch of solutions
bifurcates from the main branch, then it would seem that using a uniform
mesh (with do rather better than taking an adaptive mesh with
almost any value of ff ! 1. Presumably this is because the adaptive methods
reduce the density of mesh points close to the boundaries where the effects of
super-sensitivity indicate that important features of the solution need to be
resolved. We conclude, that an adaptive strategy for problems of the form
which concentrates on significant features such as the transition layer
but ignores subtle effects at the boundaries can give rise to most misleading
results.



--R

Moving mesh methods for problems with blow-up
On the solution of convection-diffusion boundary value problems by grid adaption
Invariant manifolds for metastable patterns in u t
SIAM J.
A suggestion for the numerical solution of the steady Navier-Stokes equations
AUTO: Software for continuation and bifurcation problems in ordinary differential equations.
The global dynamics of discrete semilinear parabolic equations.
Singularities and groups in bifurcation theory.
Moving mesh partial differential equations (MMPDEs) based on the equidistribution principle.
Singularities and groups in bifurcation theory.
Moving mesh partial differential equations (MMPDEs) based on the equidistribution principle.
Viscous shock motion for advection diffusion equations.
Shock layer movement for Burgers equation.

Convergent and spurious solutions of non-linear elliptic equations
Mesh selection for discrete solution of boundary value problems in ordinary differential equations.
Moving mesh techniques based upon equidis- tribution
Adaptive mesh selection strategies for solving boundary value problems.
Steady shock tracking and Newton's method applied to one-dimensional duct flow
Multiple solutions and bifurcation of finite difference approximations to some steady problems of fluid dy- namics
An analysis of local error control for dynamical systems.
On the dynamics of some grid-adaptation schemes
A survey of dynamically-adaptive grids in the numerical solution of partial differential equations
small eigenvalues and the sensitivity of metastable motion.

Dynamical approach study of spurious steady-state numerical solutions for nonlinear differential equations
Dynamical approach study of spurious steady-state numerical solutions for nonlinear differential equations
--TR

--CTR
Q. Sheng , A. Q. M. Khaliq, Modified arc-length adaptive algorithms for degenerate reaction-diffusion equations, Applied Mathematics and Computation, v.126 n.2-3, p.279-297, 10 March 2002
continuation;convection-diffusion equations;mesh adaptation;equidistribution