- Email: [email protected]
Evolution of reaction-diffusion patterns in infinite and bounded domains
Physica A 206 (1994) North-Holland SSDZ
380-400
PHYSICA h!
0378-4371(93)E0585-3
Evolution of reaction-diffusion infinite and bounded domains S.A.
Hassan’,
M.N.
Kuperman2,
Centro Atdmico Bariloche (CNEA) Rio Negro, Argentina
Received
5 April
H.S.
patterns
Wio3
and
and Znstituto Balseiro (UNC),
D.H.
in
Zanette3
8400 Bariloche,
1993
We introduce a semi-analytical method to study the evolution of spatial structures in reaction-diffusion systems. It consists in writing an integral equation for the relevant densities, from the propagator of the linear part of the evolution operator. In order to test the method, we perform an exhaustive study of a one-dimensional reaction-diffusion model associated to an electrical device - the ballast resistor. We consider the evolution of step and bubble-shaped initial density profiles in free space as well as in a semi-infinite domain with Dirichlet and Neumann boundary conditions. The piecewise-linear form of the reaction term, which preserves the basic ingredients of more complex nonlinear models, makes it possible to obtain exact wave-front solutions in free space and stationary solutions in the bounded domain. Short and long-time behaviour can also bc analytically studied, whereas the evolution at intermediate times is analyzed by numerical techniques. We paid particular attention to the features introduced in the evolution by boundary conditions.
1. Introduction The last decade has witnessed a growing interest in the study of formation of spatial structures [l]. In this context one of the main issues has been the study of pattern propagation, not only for fronts connecting different states (as, for instance, in bistable systems), but also for other nonhomogeneous spatial structures [2,3]. For the cases of fluid and reaction-diffusion systems, such studies have addressed problems on unbounded domains on one hand, and on finite ones on the other, where the effects of standard boundary conditions on wave length selection have been analyzed [4]. Some reviews dealing with these topics for systems described by the Ginzburg-Landau equation, as well as for reaction-diffusion systems, can be found in the recent literature [5].
1 Scholarship of CNEA, Argentina. ’ Fellowship of CONICET, Argentina. ’ Member of CONICET, Argentina.
0378.4371/94/$07.00
0
1994 - Elsevier
Science
B.V. All rights
reserved
S.A.
Hassan et al. I Reaction-diffusion
381
patterns in domains
However, the problem of structure propagation is far from being settled. It is known that the results obtained so far are mainly based on singular perturbation-like schemes [6] and/or numerical approaches. In this paper we present a novel
scheme
for studying
spatial
structure
propagation
in an unbounded
domain, as well as in cases with one boundary, leading us to analyze the effects of different types of boundary conditions on the propagation. This scheme is based on a functional-integral approach, leading to an integral equation for the relevant density, making it possible to obtain analytical asymptotic results and providing a convenient numerical procedure to analyze the evolution along the whole range of times. Within this scheme, we aim to study nonlinear reactiondiffusion systems, with one or more interacting species. Here, we introduce the method by discussing a simple one-dimensional bistable model with only one component, and study the propagation of fronts and bubble-like structures in unbounded space and in a semi-infinite domain. The boundary conditions we consider are of total absorption and total reflection. This system has been chosen as a test for our approach, in order to gain some insight before discussing more complicated (many component) systems and more general boundary conditions. The paper is organized as follows. In section 2, the model and the mathematical method are presented. In section 3 we discuss the case of evolution in the infinite domain, showing analytical results on front propagation in the limits of short and long times, and numerical calculations for intermediate times. In section 4 we analyze a bounded, semi-infinite system, studying the existence and stability of stationary solutions and obtaining, as before, analytical results for the asymptotic evolution as well as numerical results for the intermediate regime. Finally, in section 5, we discuss the results obtained so far, and make some final remarks.
2. Model and mathematical The
reaction-diffusion
method model
that
we study
here
describes
the
thermal
behaviour of an electrical device that presents an electrothermal instability the ballast resistor [7-91. This device consists of a wire surrounded by a heat bath at temperature T, and carrying a constant current 1. The spatiotemporal evolution of the temperature distribution along the wire, T(x, t), is governed by the equation
c d,T(x, t) = d,[h(x) d,T(x, t)] - q[T(x, t) - TB] + I%(T) where
h(x) is the heat conductivity,
q
is the coefficient
,
of heat transfer
(1) to the
S.A.
382
bath,
Hassan et al. I Reaction-diffusion
and c and R(T)
length,
respectively.
patterns in domains
are the specific heat and the electrical For a superconducting
resistance
wire, we can take R(T)
per unit = R&T
-
T,) where 8 is the Heaviside step function and T, is the critical temperature. Taking c, A and q as constants and TB = 0, and introducing dimensionless variables,
Eq. (1) becomes
%T(Y> 7) = ~;T(Y, T) - T(y, 7) + T&T-
T,) ,
(2)
with y = xa, r = tqlc and T,, = Z’R,lq. This is a reaction-diffusion equation for the temperature field T(y, 7 ) with a piecewise-linear reaction term f(T) = -T + T&T - T,). Coming back to the variables x and t, and calling 4(x, t) the generic field we are interested in, we shall analyze the reaction-diffusion model described by
a,% t>=+b ,
(3)
f(4) = -4 + d4$(4- 4,) >
(4)
with
where the positive constant parameters &, and +c define the reaction process. In connnection with the most commonly studied reaction-diffusion systems, the field 4 will be thought of as a density. Because of their relevance in the dynamics of the system, it is worthwhile to give a summary on the properties of the reaction term as a function of the density I$. For I_$< 4,, f(4) is simply given by a straight line whose slope equals -1, and has a root at $ = 0. At the critical point, $I = 4,, f(4) has a step discontinuity and, for higher values of the field, it is a straight line with the same slope as before. If $+, < &,, the value of the reaction term just at the right of the discontinuity is negative, f(+,‘) < 0. On the other hand, if +,, > 4,, we have f(+,‘) > 0, and there is a new root at 4 = +h. Therefore, according to how &, and += compare, the reaction dynamics is monostable or bistable. Although our analysis is valid for any value of these parameters, we shall concentrate on the more interesting, bistable case, 4,, > 4,. From this description, it is clear that this reaction-diffusion model, Eqs. (3) and (4), is a piecewise linearized mimic of the Schlogl model [lo], and can be thought of as a simplified version of more general bistability problems [ll]. The general mathematical method used in the following to deal with Eq. (3) is described in detail in appendix A. It consists in writing the solution to the reaction-diffusion equation as
S. A. Hassan et al. I Reaction-diffusion
+ c&
I I dt’
patterns in domains
dx’ X(x, t 1x’, t’) 6(+(x’, t’) - 4,) ,
383
(5)
0
where X(X, t (x’, t’) is the propagator associated to the linear part of Eq. (3) and 4,(x) = 4(x, 0) is the initial condition. The integrals on x’ are performed over the spatial domain under consideration. For a system evolving in free space, the propagator reads
X’O’(x, t 1x’,
t’) =
(- y(t-_xt’~)2 $i& exp (t - I’,)
.
(6)
Besides evolution in free space, we shall also study the problem in a semi-infinite domain, x E [0, +m). To obtain the solution in this bounded case it is necessary to fix a boundary condition at x = 0. Here, we consider two kinds of boundaries. For totally absorbing (Dirichlet) conditions, we impose +(O, t) = 0. The solution to the reaction-diffusion equation is analogous to Eq. (5), by replacing the free-space propagator .X(O)by its antisymmetrized version
.eD’(x, t 1x’,
t’) =
Yc’“‘(x, t 1x’,
t’) -
X(O)(-x, t 1x’,
t’) .
(7)
On the other hand, Neumann boundary conditions, 13,+(0, t) = 0, describe total reflexion at x = 0. In this case, the relevant propagator is given by the symmetric form
YC’N’(X,t 1x’,
t’) =
.9c0’(x, t 1x’,
t’) +
X(O)(-x, t 1x’,
t’) .
(8)
As it stands, Eq. (5) is not an explicit form for the solution to Eq. (3); in fact, the solution itself is involved in the right-hand side, through the nonlinear component ~(c#I- $J. H owever, this implicit expression for 4(x, t) suggests a recursive method to find the solution. The time integral in Eq. (5) can be discretized to a sum over N time intervals of length At, so that t = N At. Introducing the initial condition 4(x, 0) in the integral, we calculate, as a first 4(x, 2 At), step, the solution at time At, 4(x, At) and, successively, $(x,3 At), . . . , up to r$(x, t = N At). We shall use this iterative scheme to perform numerical calculations of the solutions in the cases we are interested in. The relatively simple form of the nonlinearity associated to our reactiondiffusion model makes it possible to construct an explicit expression for its
S.A.
384
Hassan et al. I Reaction-diffusion
patterns in domains
solution from its integral version, Eq. (5), at least if some reasonable hypotheses are supposed to hold. In the first place, assume that the solution 4(x, t) adopts the critical value 4, at only one point x = x,(t) for all t, such that 4(x, t) < qb, for x
replaced by form for the
solution. But once the relevant integrations to calculate 4(x, t) have been performed, the solution still depends on the unknown function x,(t). This critical value of the spatial coordinate has then to be calculated consistently with this solution. From Eq. (5) we get for x,(t),
f +
4
I I dt'
dx’ X(x,(t),
t 1x’, t’) 13(x’ - x,(t’))
,
(9)
0
which can be solved analytically in some limits but, in general, has to be treated numerically. This particular class of solutions, for which 4(x, t) = rpc at only one point x,(t), clearly contains the case of a front structure. The restriction of this front being an increasing function of x at the critical point x, can be obviously relaxed, by changing 0(x - xc) + 0(x, - x). The same method can be adapted to study the second kind of structures that we shall consider here, namely, bubble-type spatial distributions. In this case, it is assumed that the solution is greater than the critical value +C only within a finite interval of the spatial coordinate. This implies that 4(x, t) = $c at two critical points xl(t) and q(t). Now, the nonlinearity in Eq. (5) must be replaced
by a combination of two step functions of the of x,(t)) - 0(x - 4)). F rom this point on, the extension straightforward. We only remark that, for bubble structures, the type of Eq. (9) have to be solved as the consistence solution to the reaction-diffusion model. It is also clear that be extended to more complex situations.
3. Reaction-diffusion
coordinate, 0(x the procedure is two equations of condition for the the method could
in free space
In this work, we consider two initial of structures of interest in the dynamics one is a step-shaped density profile,
conditions, corresponding to two types of reaction-diffusion systems. The first
S.A.
i.e., a front centered
Hassan et al. I Reaction-diffusion
at x = ~~(0). The second
patterns in domains
initial
condition
385
is a bubble-type
distribution,
4(X>0) = hl[Nx -x,(O)) - 0(x - -@))I
(11)
7
i.e., a combination of two fronts, centered at x,(O) and x2(O) (x,(O)
3.1.
Short-time evolution
It is intuitive that, for very short times, the evolution of a front structure is independent of whether there are additional fronts at other points or not. In fact, interaction between fronts is expected to be relevant at intermediate times. This circumstance, which is confirmed by analytical calculations, indicates that it is sufficient to study the short-time evolution of the initial step distribution, Eq. (10). Since the detailed shape of the front at the initial stages is not expected to present particularly interesting features, we concentrate our attention on the evolution of the critical coordinate x,(t), at which the density equals 4,. The initial evolution of the critical coordinate is dominated by the first term in Eq. (9), as the second one is given by an integral Introducing the initial condition (10) and the free-space we conclude that Eq. (9) is satisfied for t+ 0 if x,(t) =x,(O) +
avi
(12)
)
where CY satisfies z = -Erf(a/2) parameters 4, and c#+,: z=l-2$,/4,.
in time from t = 0. propagator, Eq. (6),
and
z is a particular
combination
of the
(13)
Observe that the critical coordinate starts moving at infinite velocity, i,(t) = a/V?. The sign of this velocity depends on the sign of (Y, i.e., on that of z. For z > 0 (4, > 24,) the velocity is negative, and the front, which is an increasing
S.A.
386
Hassan et al. I Reaction-diffusion
patterns in domains
function of x, moves forward. On the other hand, if z < 0 (4, < 24,) the front moves backward. the small-time behaviour determines the For the step initial condition, subsequent evolution. As we shall see in the next section, in fact, the front preserves
the
direction
of its motion
and
approaches
a steady
shape
with
constant velocity. In the case of a bubble initial distribution, if z < 0 the short-time motion of the fronts indicates that the bubble begins to shrink. Since as time elapses there is no way to inhibit this tendency, the bubble finally disappears. For z > 0, instead, the initial tendency of the bubble is to grow. The subsequent interaction of fronts, however, can eventually stop and invert this initial growth, as shown by numerical simulations for intermediate times. 3.2.
Long-time evolution. Wave-front steady solution
Stationary - time independent - solutions are the candidates to describe the asymptotic state approached by the density at long times. In the free-space problem, the only stable stationary solutions to the reaction-diffusion equation are the homogeneous distributions &(x) = 0 and 4,(x) = c$,,, for which the density equals one of the roots of the reaction term, Eq. (4). For the step initial condition, Eq. (lo), it is numerically verified that the density approaches a steady-shaped wave front which moves with constant velocity. This distribution dominates then the long-time behaviour for this initial condition. Wave-front solutions can be found by proposing for the density a spatiotemporal dependence of the form 4(x, t) = 4(x - ut) = 4(e) where .$ =x - ut and u is the (constant) front velocity. The reaction-diffusion equation (3) becomes
where
the primes
indicate
differentiation
with respect
to 5, and 5, is the point
at which the wave front equals 4,. Due to translational invariance, the value 5, can be arbitrarily chosen. In order to represent the long-time evolution the initial condition (lo), we require the solution to vanish for e+ --OOand approach c#+, for [+ +w. Imposing continuity conditions on 4 and derivative at 5 = [,, we obtain
with A, = 3(-u + d/u’ + 4). The continuity value of the velocity:
conditions
at 5, determine
of of to its
also the
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
387
V=-&.
(16)
As in the limit of short times, u is positive for negative z, and vice versa. Thus, the ratio of the signs of u and z is the same for short and long times, and the motion
of the front
preserves
its direction.
For positive
moves backwards and the density at each point 4, = 0. On the other hand, for u < 0 the density
velocities,
the front
tends to the stationary approaches &,.
value
Depending on its initial width and on the value of z, the bubble distribution evolves towards the homogeneous trivial state, 4, = 0, when the interaction of the two fronts at intermediate times leads to their mutual annihilation. Otherwise, the fronts move away from each other and, consequently, become more independent. In this case, each of them adopts asymptotically the form of a wave front of the type of Eq. (15), and moves at constant velocity, whose modulus is given precisely by that of Eq. (16). 3.3.
Numerical study for intermediate times
The transient regime from short to long times is studied here by means of the numerical scheme described in section 2. In Fig. 1, we present the results of calculations of x,(t) for the step initial
1.0
0.5
8%
00.
-0.5
-1 .o
time Fig. 1. Free space evolution of the critical coordinate (lo), with x,(O) = 0, for various values of z. Dashed asymptotic behaviour.
for the step-shaped initial lines indicate the analytical
condition, Eq. result for the
S.A.
388
condition, crossover relatively singular
Hassan
et al. I Reaction-diffusion
patterns
in domains
Eq. (lo), with x,(O) = 0 and for various values of z. We see that the from the short-time to the long-time regime is straightforward. In a short initial
x,(t) m t. Dashed
interval, 0
(16) for each value of z. agreement with the analytical The intermediate evolution
such motion,
as calculated
passes uniform
from its motion,
from Eqs. (15) and
These numerical calculations results presented above. for the bubble initial condition,
are
thus
Eq. (ll),
in
full
is more
interesting. Here we take x,(O) = -x2(O) < 0 and, due to the symmetry of the problem, it is sufficient to study the evolution of only one critical coordinate, say the positive one, x*(t). In Fig. 2, we consider the cases x2(O) = 1 and x2(O) = 5, for some values of z. Similar curves in both cases correspond to the same value of 2. The initial behaviour reproduces the analytical result: the bubble grows or shrinks according to the sign of z. Thus, for z < 0 the bubble will disappear. For z ~0, on the other hand, the interaction of the two fronts becomes relevant at intermediate times and determines their subsequent evolution. Compare, in fact, the curves for z = 0.05. In the first case (x2(O) = l), in spite of its initial growth, the bubble finally shrinks and disappears. In the second case, instead, the fronts are farther from each other, so that their interaction is not sufficient to stop the bubble growth. The fronts will then separate and
6
5
4
2 ry
3
z=o.33
_----
2
_
/
-
-
---
l.z.:~+~~~~___,___“~~o~5, ,,,_.....
1
1,
0 0.0
0.5
,,,I
1.0
..,__,,..._...... ..
.....
/ /,,, 1.5
1
2.0
...‘_”
,,,,
2.5
time Fig. 2. As in Fig. 1, for the bubble-shaped initial condition, Eq. (ll), with x,(O) = -x,(O)
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
389
eventually attain the asymptotic regime described by the front wave (15). As shown in Fig. 2, this is also the case for higher values of z.
4. Reaction-diffusion
with a boundary
in Eq.
condition
In this section, we consider the ballast model in a one-dimensional semiinfinite space, x E [0, +m). Boundary conditions at x = 0 are fixed to determine the solution for a given initial condition c$(.x, 0). Dirichlet boundary conditions,
are required to describe total particle hand, Neumann conditions,
absorption
at the boundary;
on the other
which insure that the particle current vanishes at x = 0, describe total reflection at the origin. Such boundary conditions introduce nontrivial features in the behaviour of the system. In the first place, stationary, time independent nonhomogeneous solutions are now possible. Under conditions on their stability, these stationary density profiles are expected to coincide with the long-time asymptotic solutions. Furthermore, the time evolution of the system is also severely influenced by the boundary conditions: the front-like initial structures considagainst the boundary, ered in the previous section are shown to “bounce” drastically modifying the subsequent studied in the following sections. 4.1.
dynamics.
Details
on these
topics
are
Stationary solutions and their stability
In order to study the stationary solutions of the ballast model with the boundary conditions (17) or (18), we come back to the original form of the reaction-diffusion equation, Eq. (3), and impose a, = 0, obtaining
(19) If we associate time with x and coordinate with 4, this has the form Newton equation describing the motion of a particle under the action potential
V($) = -44’ + 4h(4 - 4%)et4 - A).
of a of a
(20)
390
S.A.
This mechanical
Hassan et al. I Reaction-diffusion
analogy
can be fruitfully
patterns in domains
exploited
stationary reaction-diffusion problem [9]. It is clear that a homogeneous stationary
solution
to find the solution
to the
4,(x) = 0 always
exists.
The character of the nontrivial solutions depends on the boundary condition and on the value of the parameter z, defined in Eq. (13). Regarding their stability, here it will be analyzed in the linear approximation, by studying the evolution
of a small perturbation
applied
to the stationary
4.1.1. Dirichlet boundary condition When the density is required to satisfy Eq. (17), besides a nonhomogeneous solution exists for z > 0. It reads
solution.
the trivial
solution,
forx x, ,
(21)
where
(22) is the stationary critical coordinate, at which the stationary density reaches the critical value $=. This solution clearly satisfies the Dirichlet condition, and grows monotonically in space, approaching c#+,as x* +m. A linear stability analysis indicates that the trivial solution is stable for any z, and that the nonhomogeneous solution (21) is also stable when it exists. Therefore, for z > 0 the system is bistable, and the long-time behaviour of the solutions will be determined by the initial condition. 4.1.2. Neumann boundary condition The requirement of a,~$ = 0 at the origin makes it possible to have a homogeneous stationary solution different from the trivial one. It corresponds to the second root of the reaction term in our reaction-diffusion model,
4,(x) = 4%. For .z < 0, a nonhomogeneous
(23) solution
is proven
to exist also. It reads
forx x, ,
where
(24)
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
391
(25) is, as for Dirichlet
conditions,
the
stationary
critical
coordinate.
density has a finite value at x = 0, namely, 4,(O) = 4,/cash x, < 4,. monotonically for increasing x and asymptotically approaches &,.
Now,
the
It grows
This nonhomogeneous stationary solution proves to be linearly unstable. On the other hand, both homogeneous solutions are stable for any value of z. Thus, independently of the parameters 4, and c$,, (4, > +C), Neumann conditions determine the system to be bistable. Again, its long-time evolution will be selected by the initial condition.
4.2.
Short-time and long-time behaviour
As we could expect, the boundary condition at x = 0 does not affect the short-time evolution of any initial density distribution. Front structures are thus shown to reproduce the initial behaviour observed in free space: they go forward if z > 0 and move backward if z < 0. In some particular situations, also the long-time evolution of the density is well described by a free-space solution. We recall that, in free space, the solution for a step-shaped initial condition evolves towards a limiting steady wave front, moving with velocity u = -22/m. In the semi-infinite domain, this behaviour is reproduced for the same initial condition when the density tends to a vanishing stationary state: the front structure moves towards x-+ +m and, as it separates from the boundary, it approaches the shape and the velocity of the free-space wave front solution, Eqs. (15) and (16). On the other hand, when nonvanishing stationary states are approached from bubbletype initial conditions, the rightmost front moves also towards x+ +m and tends to reproduce the free-space behaviour. These limiting features in the bounded solutions, which exist from an analytical study of Eq. (5), are confirmed calculations presented in the following section.
4.3.
can be proven to by the numerical
Numerical study for intermediate times
We present here numerical results for the bounded system in the intermediate regime. For the step initial condition, Eq. (lo), we take x,(O) = 1; for the bubble initial condition, we choose x,(O) = 1 and x2(O) = 2. These cases exhibit all the qualitatively different features that we were able to observe for other values of the initial critical coordinates. Thus, they constitute a proper summary of more general situations.
392
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
1.03
/
0.60
0
2
4
8
6
10
time Fig. 3. Evolution of the critical x,(O) = 1 and Dirichlet boundary
coordinate condition
for the step-shaped at x = 0.
initial
condition,
Eq.
(lo),
with
4.3.1. Dirichlet boundary condition Fig. 3 displays the evolution of the critical coordinate x,(t) for the step initial condition with Dirichlet boundary conditions, for various values of z. As predicted by the analytical results, for z > 0, x,(t) tends to a limiting value, which can be shown to coincide with the critical coordinate X, of the stationary solution, Eq. (21). For negative z, on the other hand, the front separates from the boundary, approaching the free-space wave-front evolution. At intermediate times nontrivial features appear in the evolution. In particular, for sufficiently small, positive values of z the motion of the critical coordinate is nonmonotonic. It begins by approaching the boundary, but at a certain point it “bounces” and, from then on, it moves backward. For small negative values of z (see the curve for z = -0.03) a sign of this “bounce” is preserved: the time derivative of x,(t) evolves nonmonotonically. Results for the bubble initial condition are plotted in Fig. 4. Here, the lower part of each curve corresponds to xl(t), and the upper one, to x2(t). The bounce effect plays now a determinant role in the survival of the bubble. Consider, for example, the curve for z = 0.30. Up to t- 1.5, the upper critical coordinate x2 seems to have escaped from the influence of the other front, and to approach the wave-front behaviour observed in free space. However, at t ~0.3, xl(t) has bounced back, and approaches now the other critical coordinate at a velocity higher than that of x2(t). Finally, the fronts collide and annihilate at time t,.
S.A.
Hassan et al. I Reaction-diffusion
393
patterns in domains
3.0 2.5 2.0
8
1.5
2 1 .o
0.5
z=oso 0.0
j
0.0
”
1 ” 0.5
”
I ““““‘1 1 .o
”
2.0
1.5
”
1 ”
2.5
”
”
3.0
”
3.5
time Fig. 4. As in Fig. 3, for the bubble-shaped initial condition, Eq. (ll), with x,(O) = 1 and x,(O) = 2. The lower part of each curve corresponds to x,(t), and the upper part, to x,(t). The small plot in the upper-right corner shows the front annihilation time t, as a function of z.
This annihilation time depends on the value of z, and has been plotted in the upper-right small figure. We see that ta(z)+ CCfor z = 0.31. For higher values of z the upper front is able to escape from the other and the bubble survives and grows. In this case, x,(t) approaches the limiting value X, predicted for the stationary solution, and x*(f) acquires a constant velocity whose modulus equals that of Eq. (16). Observe that, even in this case, the intermediate evolution can be highly nontrivial: for z = 0.35, for instance, xl(t) exhibits a kind of damped oscillation as it approaches x,.
4.3.2. Neumann boundary condition We recall that, for Neumann conditions, there exists no stable stationary nonhomogeneous solution, so that the solution must approach one of the two homogeneous states, 4 = 0 or 4 = +,,. As we shall see, this fact determines essential differences with the case of Dirichlet conditions. In Fig. 5, we plotted the critical coordinate for the step initial condition. For sufficiently
negative
values
of z the front moves
as if it were in free space,
and
394
S.A.
0.0
Hassan et al. I Reaction-diffusion
0.5
1.0
1.5
patterns in domains
2.0
2.5
3.0
time Fig. 5. As in Fig. 3, for Neumann
boundary
conditions
at the origin.
the solution tends to the vanishing state. Approaching z = 0, however, the boundary condition begins to play a role in the evolution of the front, and although the critical coordinate moves initially with positive velocity, it eventually inverts its motion and reaches x = 0 in a finite time. From then on, the solution is greater than $, for all x, and approaches asymptotically the homogeneous state $,(x) = 4,. This is also the limiting behaviour displayed for z > 0, where a nonmonotonic time dependence is observed in the derivative of x,(t). Note that the critical coordinate reaches the boundary at infinite velocity. The features observed for the step initial condition determine the evolution of the bubble distribution. Results for this situation are presented in Fig. 6. Now, with respect to the Dirichlet case, the fronts interchange in some sense the role they play in the evolution. In fact, the annihilation of the bubble for z >O is here a consequence of the change of the sign of the upper front velocity. It moves initially with positive velocity but, as observed for the step distribution, it eventually stops and moves towards the origin. Then, it encounters the other front, and their interaction leads to mutual annihilation at a time t,(z). This time diverges at z = 0.216. For higher z, the bubble grows, x,(t) reaches the boundary, x2(t) adopts the wave-front motion, and the solution approaches the nonvanishing homogeneous state. For Neumann boundary conditions, both fronts can display a rather complex time dependence, with the damped oscillatory behaviour already observed with a Dirichlet boundary.
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
395
2.5
0.3 -0.1 4.l 0.0 0.1 0.2
2.0
1.5
z
z
1.0
0.5
0.0 0.0
0.3
0.6
0.9
1.2
1.5
1.8
conditions
at the origin.
time Fig. 6. As in Fig. 4, for Neumann
boundary
5. Conclusions
In this paper we have addressed the problem of pattern propagation within reaction-diffusion models. The system we have chosen for our analysis is a simplified version of the ballast resistor, whose piecewise-linear reaction term mimics the behaviour of other more complicated models (i.e., the Schliigl model and other bistable ones). The scheme we have used is a novel one based on a functional integral approach, leading to an integral equation for the relevant density. This scheme allows us to obtain analytic asymptotic results as well as to analyze the evolution along the whole range of times by means of a suitable numerical procedure. Although our interest is the study of structure propagation iri general nonlinear reaction-diffusion systems composed of one or more species, here the method has been introduced by discussing the simple one-component one-dimensional bistable model indicated above. We have considered the case of fronts connecting different states in a bistable model, as well as bubble-like structures. Our study has involved infinite systems on one hand, and bounded systems on the other. In this latter case, we paid particular attention to the effect of totally absorbing (Dirichlet) and totally reflecting (Neumann) boundary conditions on the propagation of structures. For the case of free-space evolution, numerical results agree quite well with
S.A.
396
Hassan et al.
I Reaction-diffusion
patterns in domains
the analytic limits at short and long times, also showing the tendency to reach the expected asymptotic wave-front behaviour. For bubble-like structures, our results coincide with previous ones, obtained through different methods (see, for instance, [ 121). On the other hand, for the semi-infinite system, we have found that the boundary condition influences the time evolution of the structures in some unexpected ways, for example, with the appearance of bounces and damped oscillatory-like effects. Other relevant aspects have arisen when analyzing the evolution of bubble-like structures. We have seen that, even for situations where the initial condition in the free-space case would lead to a growing structure, the presence of a boundary can cause its collapse. As the behaviour of bubble structures has been used to describe nucleation phenomena [13], these results indicate that boundary effects would influence the nucleation dynamics in small size systems. Here, we have discussed the evolution of the ballast reaction-diffusion model by looking only at the behaviour of the critical coordinate x,(t) - which corresponds to the front position - as the detailed form of the front profiles is not expected to show particularly interesting features. Nevertheless, in Fig. 7 we show a comparison of the density profiles at a fixed time for a bubble initial distribution in free space and with the two boundary conditions considered here. The features introduced by the boundary are apparent. In conclusion, by analyzing the problem described above, we have shown 0.6
,
/
I
0.5 0.4 0.3 0.2
0.1
0.0
x
Fig. 7. Density profiles for a bubble initial condition with x,(O) = 1 and x2(O) = 2, at time t = 0.8, evolving in free space (F), with Dirichlet conditions (D), and with Neumann conditions (N). In this plot, &, = 1 and z = 0.23. Compare with Figs. 4 and 6.
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
391
that the present propagation scheme is very suitable for the study of evolution in reaction-diffusion systems. The insight so far gained will lead us to the analysis of more complex (many component) problems, as well as more general boundary conditions. These topics are the subject of work in progress.
Acknowledgement
We thank Prof. Veronica Griinfeld for her careful reading of the manuscript.
Appendix
A. Path-integral
solution
to reaction-diffusion
equations
We consider here a general nonlinear reaction-diffusion
equation of the form
where we assume that f[+(x, t)] h as such mild behaviour that we may proceed without restrictions. For instance, it is shown in [14] that if f[+(x, t)] as well as and --oo
4(x,t)=/9[x(t)]&(x)exp[-
Ida
(+,+)2+f~$‘~)1)].
2
fo
(27)
However, this is not the most adequate expression to perform quantitative (numerical) studies on the time evolution of the fronts or spatial structures which arise in the system described by Eq. (26). In order to overcome this drawback and to have a representation of the solution suitable for numerical analysis, we use a path integral-like approach. To obtain the propagation for a given initial condition, we integrate Eq. (26) over a very short time interval (t, t + T) to the lowest order in r and use a prepoint discretization scheme [15] to find
a, t +7) = (1+ ~~:>d&t>+ Tf[4+, 91 . Through the usual procedure
(28)
of Fourier transforming
4(x, t = t’ + T) = ~f[+(x, t’)] +
I--$&exp( -
[15], we obtain from here
& (x - x’)‘) 4(x’, t’) . (29)
S.A.
398
Considering propagation
Hassan et al. I Reaction-diffusion
patterns in domains
always the lowest order in T and iterating Eq. (29), we get the of the initial distribution 4(x, t = to) = c#I,(x) up to the point xr and
to the time t, (tf - t, = NT) as
x f[4(Xk, to+ kT)l>
(30)
where xj = x(t,) is the coordinate at tj = t, + j7. In the continuous limit (N-+ 00, r+ 0 and NT = t, - to) and in the usual notation we finally obtain
%I
~~(“)(x,, t, I-%,to) 4(-ql, to)
If
I I dt'
+
dx’
X(“)(xf, t, 1x’,
t’) f[+(x’,
t’)] ,
(31)
10
where X”‘(x, t 1x’, t’) is the propagator diffusion operator, given, in the discrete
associated case, by
to the (linear)
free-space
(32)
with the continuous
limit
X(“)(x,,tr/xo,to)=~9[x]exp(-$
[dsx(s)‘).
(33)
10 It is easy to prove that the result indicated in Eq. (31) is a solution to Eq. (26) and satisfies the given initial condition. It is clear that the above indicated procedure will also work if we separate in our equation all the linear contributions (call it 24) from the nonlinear contributions (call it N(4)). Eq. (26) becomes a&(x,
t) = =W + J(4)
>
(34)
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
399
and has a solution similar in form to Eqs. (31), (32), (33), where now XC’) is the propagator of the linear operator LZ. In the present case (the ballast model), the nonlinear reaction function has the form f(4) = -4 + 4@(4 - 4,) .
(35)
Then we choose as our linear operator ~4=&-4%
(36)
leaving as the nonlinear part J(4)
(37)
= WV+ - 4%).
After evaluating X(O) for the operator
X(O)(Xf, t, Ixo,4J =
&q
in Eq. (36) we get
exp - y;t-_x;y f -(ff - to,), i
(38)
0
which, replaced in Eq. (31), gives the final form used in this work.
References PI P.C. Hohenberg and M.C. Cross, An introduction to pattern formation in nonequilibrium systems, in: Fluctuations and Stochastic Phenomena in Condensed Matter, L. Garrido, ed. (Springer, Berlin, 1986); A.C. Newell, The dynamics and analysis of patterns, in: Complex Systems, D. Stein, ed. (Addison-Wesley, New York, 1989). PI G.T. Dee, J. Stat. Phys. 39 (1985) 705; G.T. Dee and W. van Saarloos, Phys. Rev. Lett. 60 (1988) 2641; M. Nickles, M. Lticke and H. Miiller-Krumbhaar, Europhys. Lett. 9 (1989) 237; W. van Saarloos, Phys. Rev. A 37 (1988) 211; 39 (1989) 6367. [31 P. Hagan, SIAM J. Math. Anal. 57 (1982) 717; J. Rinzel and D. Terran, SIAM J. Appl. Math. 42 (1982) 1111; J.D. Dockery and J.P. Keener, SIAM J. Appl. Math. 49 (1989) 539; H. Ikeda and M. Mimura, SIAM J. Appl. Math. 49 (1989) 515; H. Yamada and K. Nozaki, Prog. Theor. Phys. 84 (1990) 801. [41 L. Kramer and P.C. Hohenberg, Physica D 13 (1984) 357; P.C. Hohenberg, L. Kramer and H. Riecke, Physica D 15 (1985) 402; M.C. Cross, Phys. Rev. Lett. 57 (1986) 2935. PI W. van Saarloos and P.C. Hohenberg, Physica D 56 (1992) 303; E. Meron, Phys. Rep. 218 (1992) 1. 161 J.J. Tyson and J.P. Keener, Physica D 32 (1988) 327. [71 B. Ross and J.D. Lister, Phys. Rev. A 15 (1977) 1246; R. Landauer, Phys. Rev. A 15 (1977) 2117. PI D. Bedeaux and P. Mazur, Physica A 105 (1981) 1.
400 [9] [lo] [ll] [12] [13]
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
C. Schat and H.S. Wio, Physica A 180 (1992) 295. F. Schlogl, Z. Phys. 248 (1971) 446; 253 (1972) 147. Ch. Zulike, A.S. Mikhailov and L. Schimansky-Geier, Physica A 163 (1990) 559. A.S. Mikhailov, Foundations of Sinergetics I (Springer, Berlin, 1991). L. Schimanski-Geier and W. Ebeling, Ann. Phys. (Leipzig) 40 (1983) 10; L. Schimansky-Geier et al., Ann. Phys. (Leipzig) 40 (1983) 277. [14] M. Freidlin, Functional Integration and Partial Differential Equations (Princeton Univ. Press, New Jersey, 1985). [15] F. Langouche, D. Roekaerts and E. Tirapegui, Functional Integration and Semiclassical Expansions (Reidel, Dordrecht, 1982).
380-400
PHYSICA h!
0378-4371(93)E0585-3
Evolution of reaction-diffusion infinite and bounded domains S.A.
Hassan’,
M.N.
Kuperman2,
Centro Atdmico Bariloche (CNEA) Rio Negro, Argentina
Received
5 April
H.S.
patterns
Wio3
and
and Znstituto Balseiro (UNC),
D.H.
in
Zanette3
8400 Bariloche,
1993
We introduce a semi-analytical method to study the evolution of spatial structures in reaction-diffusion systems. It consists in writing an integral equation for the relevant densities, from the propagator of the linear part of the evolution operator. In order to test the method, we perform an exhaustive study of a one-dimensional reaction-diffusion model associated to an electrical device - the ballast resistor. We consider the evolution of step and bubble-shaped initial density profiles in free space as well as in a semi-infinite domain with Dirichlet and Neumann boundary conditions. The piecewise-linear form of the reaction term, which preserves the basic ingredients of more complex nonlinear models, makes it possible to obtain exact wave-front solutions in free space and stationary solutions in the bounded domain. Short and long-time behaviour can also bc analytically studied, whereas the evolution at intermediate times is analyzed by numerical techniques. We paid particular attention to the features introduced in the evolution by boundary conditions.
1. Introduction The last decade has witnessed a growing interest in the study of formation of spatial structures [l]. In this context one of the main issues has been the study of pattern propagation, not only for fronts connecting different states (as, for instance, in bistable systems), but also for other nonhomogeneous spatial structures [2,3]. For the cases of fluid and reaction-diffusion systems, such studies have addressed problems on unbounded domains on one hand, and on finite ones on the other, where the effects of standard boundary conditions on wave length selection have been analyzed [4]. Some reviews dealing with these topics for systems described by the Ginzburg-Landau equation, as well as for reaction-diffusion systems, can be found in the recent literature [5].
1 Scholarship of CNEA, Argentina. ’ Fellowship of CONICET, Argentina. ’ Member of CONICET, Argentina.
0378.4371/94/$07.00
0
1994 - Elsevier
Science
B.V. All rights
reserved
S.A.
Hassan et al. I Reaction-diffusion
381
patterns in domains
However, the problem of structure propagation is far from being settled. It is known that the results obtained so far are mainly based on singular perturbation-like schemes [6] and/or numerical approaches. In this paper we present a novel
scheme
for studying
spatial
structure
propagation
in an unbounded
domain, as well as in cases with one boundary, leading us to analyze the effects of different types of boundary conditions on the propagation. This scheme is based on a functional-integral approach, leading to an integral equation for the relevant density, making it possible to obtain analytical asymptotic results and providing a convenient numerical procedure to analyze the evolution along the whole range of times. Within this scheme, we aim to study nonlinear reactiondiffusion systems, with one or more interacting species. Here, we introduce the method by discussing a simple one-dimensional bistable model with only one component, and study the propagation of fronts and bubble-like structures in unbounded space and in a semi-infinite domain. The boundary conditions we consider are of total absorption and total reflection. This system has been chosen as a test for our approach, in order to gain some insight before discussing more complicated (many component) systems and more general boundary conditions. The paper is organized as follows. In section 2, the model and the mathematical method are presented. In section 3 we discuss the case of evolution in the infinite domain, showing analytical results on front propagation in the limits of short and long times, and numerical calculations for intermediate times. In section 4 we analyze a bounded, semi-infinite system, studying the existence and stability of stationary solutions and obtaining, as before, analytical results for the asymptotic evolution as well as numerical results for the intermediate regime. Finally, in section 5, we discuss the results obtained so far, and make some final remarks.
2. Model and mathematical The
reaction-diffusion
method model
that
we study
here
describes
the
thermal
behaviour of an electrical device that presents an electrothermal instability the ballast resistor [7-91. This device consists of a wire surrounded by a heat bath at temperature T, and carrying a constant current 1. The spatiotemporal evolution of the temperature distribution along the wire, T(x, t), is governed by the equation
c d,T(x, t) = d,[h(x) d,T(x, t)] - q[T(x, t) - TB] + I%(T) where
h(x) is the heat conductivity,
q
is the coefficient
,
of heat transfer
(1) to the
S.A.
382
bath,
Hassan et al. I Reaction-diffusion
and c and R(T)
length,
respectively.
patterns in domains
are the specific heat and the electrical For a superconducting
resistance
wire, we can take R(T)
per unit = R&T
-
T,) where 8 is the Heaviside step function and T, is the critical temperature. Taking c, A and q as constants and TB = 0, and introducing dimensionless variables,
Eq. (1) becomes
%T(Y> 7) = ~;T(Y, T) - T(y, 7) + T&T-
T,) ,
(2)
with y = xa, r = tqlc and T,, = Z’R,lq. This is a reaction-diffusion equation for the temperature field T(y, 7 ) with a piecewise-linear reaction term f(T) = -T + T&T - T,). Coming back to the variables x and t, and calling 4(x, t) the generic field we are interested in, we shall analyze the reaction-diffusion model described by
a,% t>=+b
(3)
f(4) = -4 + d4$(4- 4,) >
(4)
with
where the positive constant parameters &, and +c define the reaction process. In connnection with the most commonly studied reaction-diffusion systems, the field 4 will be thought of as a density. Because of their relevance in the dynamics of the system, it is worthwhile to give a summary on the properties of the reaction term as a function of the density I$. For I_$< 4,, f(4) is simply given by a straight line whose slope equals -1, and has a root at $ = 0. At the critical point, $I = 4,, f(4) has a step discontinuity and, for higher values of the field, it is a straight line with the same slope as before. If $+, < &,, the value of the reaction term just at the right of the discontinuity is negative, f(+,‘) < 0. On the other hand, if +,, > 4,, we have f(+,‘) > 0, and there is a new root at 4 = +h. Therefore, according to how &, and += compare, the reaction dynamics is monostable or bistable. Although our analysis is valid for any value of these parameters, we shall concentrate on the more interesting, bistable case, 4,, > 4,. From this description, it is clear that this reaction-diffusion model, Eqs. (3) and (4), is a piecewise linearized mimic of the Schlogl model [lo], and can be thought of as a simplified version of more general bistability problems [ll]. The general mathematical method used in the following to deal with Eq. (3) is described in detail in appendix A. It consists in writing the solution to the reaction-diffusion equation as
S. A. Hassan et al. I Reaction-diffusion
+ c&
I I dt’
patterns in domains
dx’ X(x, t 1x’, t’) 6(+(x’, t’) - 4,) ,
383
(5)
0
where X(X, t (x’, t’) is the propagator associated to the linear part of Eq. (3) and 4,(x) = 4(x, 0) is the initial condition. The integrals on x’ are performed over the spatial domain under consideration. For a system evolving in free space, the propagator reads
X’O’(x, t 1x’,
t’) =
(- y(t-_xt’~)2 $i& exp (t - I’,)
.
(6)
Besides evolution in free space, we shall also study the problem in a semi-infinite domain, x E [0, +m). To obtain the solution in this bounded case it is necessary to fix a boundary condition at x = 0. Here, we consider two kinds of boundaries. For totally absorbing (Dirichlet) conditions, we impose +(O, t) = 0. The solution to the reaction-diffusion equation is analogous to Eq. (5), by replacing the free-space propagator .X(O)by its antisymmetrized version
.eD’(x, t 1x’,
t’) =
Yc’“‘(x, t 1x’,
t’) -
X(O)(-x, t 1x’,
t’) .
(7)
On the other hand, Neumann boundary conditions, 13,+(0, t) = 0, describe total reflexion at x = 0. In this case, the relevant propagator is given by the symmetric form
YC’N’(X,t 1x’,
t’) =
.9c0’(x, t 1x’,
t’) +
X(O)(-x, t 1x’,
t’) .
(8)
As it stands, Eq. (5) is not an explicit form for the solution to Eq. (3); in fact, the solution itself is involved in the right-hand side, through the nonlinear component ~(c#I- $J. H owever, this implicit expression for 4(x, t) suggests a recursive method to find the solution. The time integral in Eq. (5) can be discretized to a sum over N time intervals of length At, so that t = N At. Introducing the initial condition 4(x, 0) in the integral, we calculate, as a first 4(x, 2 At), step, the solution at time At, 4(x, At) and, successively, $(x,3 At), . . . , up to r$(x, t = N At). We shall use this iterative scheme to perform numerical calculations of the solutions in the cases we are interested in. The relatively simple form of the nonlinearity associated to our reactiondiffusion model makes it possible to construct an explicit expression for its
S.A.
384
Hassan et al. I Reaction-diffusion
patterns in domains
solution from its integral version, Eq. (5), at least if some reasonable hypotheses are supposed to hold. In the first place, assume that the solution 4(x, t) adopts the critical value 4, at only one point x = x,(t) for all t, such that 4(x, t) < qb, for x
replaced by form for the
solution. But once the relevant integrations to calculate 4(x, t) have been performed, the solution still depends on the unknown function x,(t). This critical value of the spatial coordinate has then to be calculated consistently with this solution. From Eq. (5) we get for x,(t),
f +
4
I I dt'
dx’ X(x,(t),
t 1x’, t’) 13(x’ - x,(t’))
,
(9)
0
which can be solved analytically in some limits but, in general, has to be treated numerically. This particular class of solutions, for which 4(x, t) = rpc at only one point x,(t), clearly contains the case of a front structure. The restriction of this front being an increasing function of x at the critical point x, can be obviously relaxed, by changing 0(x - xc) + 0(x, - x). The same method can be adapted to study the second kind of structures that we shall consider here, namely, bubble-type spatial distributions. In this case, it is assumed that the solution is greater than the critical value +C only within a finite interval of the spatial coordinate. This implies that 4(x, t) = $c at two critical points xl(t) and q(t). Now, the nonlinearity in Eq. (5) must be replaced
by a combination of two step functions of the of x,(t)) - 0(x - 4)). F rom this point on, the extension straightforward. We only remark that, for bubble structures, the type of Eq. (9) have to be solved as the consistence solution to the reaction-diffusion model. It is also clear that be extended to more complex situations.
3. Reaction-diffusion
coordinate, 0(x the procedure is two equations of condition for the the method could
in free space
In this work, we consider two initial of structures of interest in the dynamics one is a step-shaped density profile,
conditions, corresponding to two types of reaction-diffusion systems. The first
S.A.
i.e., a front centered
Hassan et al. I Reaction-diffusion
at x = ~~(0). The second
patterns in domains
initial
condition
385
is a bubble-type
distribution,
4(X>0) = hl[Nx -x,(O)) - 0(x - -@))I
(11)
7
i.e., a combination of two fronts, centered at x,(O) and x2(O) (x,(O)
3.1.
Short-time evolution
It is intuitive that, for very short times, the evolution of a front structure is independent of whether there are additional fronts at other points or not. In fact, interaction between fronts is expected to be relevant at intermediate times. This circumstance, which is confirmed by analytical calculations, indicates that it is sufficient to study the short-time evolution of the initial step distribution, Eq. (10). Since the detailed shape of the front at the initial stages is not expected to present particularly interesting features, we concentrate our attention on the evolution of the critical coordinate x,(t), at which the density equals 4,. The initial evolution of the critical coordinate is dominated by the first term in Eq. (9), as the second one is given by an integral Introducing the initial condition (10) and the free-space we conclude that Eq. (9) is satisfied for t+ 0 if x,(t) =x,(O) +
avi
(12)
)
where CY satisfies z = -Erf(a/2) parameters 4, and c#+,: z=l-2$,/4,.
in time from t = 0. propagator, Eq. (6),
and
z is a particular
combination
of the
(13)
Observe that the critical coordinate starts moving at infinite velocity, i,(t) = a/V?. The sign of this velocity depends on the sign of (Y, i.e., on that of z. For z > 0 (4, > 24,) the velocity is negative, and the front, which is an increasing
S.A.
386
Hassan et al. I Reaction-diffusion
patterns in domains
function of x, moves forward. On the other hand, if z < 0 (4, < 24,) the front moves backward. the small-time behaviour determines the For the step initial condition, subsequent evolution. As we shall see in the next section, in fact, the front preserves
the
direction
of its motion
and
approaches
a steady
shape
with
constant velocity. In the case of a bubble initial distribution, if z < 0 the short-time motion of the fronts indicates that the bubble begins to shrink. Since as time elapses there is no way to inhibit this tendency, the bubble finally disappears. For z > 0, instead, the initial tendency of the bubble is to grow. The subsequent interaction of fronts, however, can eventually stop and invert this initial growth, as shown by numerical simulations for intermediate times. 3.2.
Long-time evolution. Wave-front steady solution
Stationary - time independent - solutions are the candidates to describe the asymptotic state approached by the density at long times. In the free-space problem, the only stable stationary solutions to the reaction-diffusion equation are the homogeneous distributions &(x) = 0 and 4,(x) = c$,,, for which the density equals one of the roots of the reaction term, Eq. (4). For the step initial condition, Eq. (lo), it is numerically verified that the density approaches a steady-shaped wave front which moves with constant velocity. This distribution dominates then the long-time behaviour for this initial condition. Wave-front solutions can be found by proposing for the density a spatiotemporal dependence of the form 4(x, t) = 4(x - ut) = 4(e) where .$ =x - ut and u is the (constant) front velocity. The reaction-diffusion equation (3) becomes
where
the primes
indicate
differentiation
with respect
to 5, and 5, is the point
at which the wave front equals 4,. Due to translational invariance, the value 5, can be arbitrarily chosen. In order to represent the long-time evolution the initial condition (lo), we require the solution to vanish for e+ --OOand approach c#+, for [+ +w. Imposing continuity conditions on 4 and derivative at 5 = [,, we obtain
with A, = 3(-u + d/u’ + 4). The continuity value of the velocity:
conditions
at 5, determine
of of to its
also the
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
387
V=-&.
(16)
As in the limit of short times, u is positive for negative z, and vice versa. Thus, the ratio of the signs of u and z is the same for short and long times, and the motion
of the front
preserves
its direction.
For positive
moves backwards and the density at each point 4, = 0. On the other hand, for u < 0 the density
velocities,
the front
tends to the stationary approaches &,.
value
Depending on its initial width and on the value of z, the bubble distribution evolves towards the homogeneous trivial state, 4, = 0, when the interaction of the two fronts at intermediate times leads to their mutual annihilation. Otherwise, the fronts move away from each other and, consequently, become more independent. In this case, each of them adopts asymptotically the form of a wave front of the type of Eq. (15), and moves at constant velocity, whose modulus is given precisely by that of Eq. (16). 3.3.
Numerical study for intermediate times
The transient regime from short to long times is studied here by means of the numerical scheme described in section 2. In Fig. 1, we present the results of calculations of x,(t) for the step initial
1.0
0.5
8%
00.
-0.5
-1 .o
time Fig. 1. Free space evolution of the critical coordinate (lo), with x,(O) = 0, for various values of z. Dashed asymptotic behaviour.
for the step-shaped initial lines indicate the analytical
condition, Eq. result for the
S.A.
388
condition, crossover relatively singular
Hassan
et al. I Reaction-diffusion
patterns
in domains
Eq. (lo), with x,(O) = 0 and for various values of z. We see that the from the short-time to the long-time regime is straightforward. In a short initial
x,(t) m t. Dashed
interval, 0
(16) for each value of z. agreement with the analytical The intermediate evolution
such motion,
as calculated
passes uniform
from its motion,
from Eqs. (15) and
These numerical calculations results presented above. for the bubble initial condition,
are
thus
Eq. (ll),
in
full
is more
interesting. Here we take x,(O) = -x2(O) < 0 and, due to the symmetry of the problem, it is sufficient to study the evolution of only one critical coordinate, say the positive one, x*(t). In Fig. 2, we consider the cases x2(O) = 1 and x2(O) = 5, for some values of z. Similar curves in both cases correspond to the same value of 2. The initial behaviour reproduces the analytical result: the bubble grows or shrinks according to the sign of z. Thus, for z < 0 the bubble will disappear. For z ~0, on the other hand, the interaction of the two fronts becomes relevant at intermediate times and determines their subsequent evolution. Compare, in fact, the curves for z = 0.05. In the first case (x2(O) = l), in spite of its initial growth, the bubble finally shrinks and disappears. In the second case, instead, the fronts are farther from each other, so that their interaction is not sufficient to stop the bubble growth. The fronts will then separate and
6
5
4
2 ry
3
z=o.33
_----
2
_
/
-
-
---
l.z.:~+~~~~___,___“~~o~5, ,,,_.....
1
1,
0 0.0
0.5
,,,I
1.0
..,__,,..._...... ..
.....
/ /,,, 1.5
1
2.0
...‘_”
,,,,
2.5
time Fig. 2. As in Fig. 1, for the bubble-shaped initial condition, Eq. (ll), with x,(O) = -x,(O)
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
389
eventually attain the asymptotic regime described by the front wave (15). As shown in Fig. 2, this is also the case for higher values of z.
4. Reaction-diffusion
with a boundary
in Eq.
condition
In this section, we consider the ballast model in a one-dimensional semiinfinite space, x E [0, +m). Boundary conditions at x = 0 are fixed to determine the solution for a given initial condition c$(.x, 0). Dirichlet boundary conditions,
are required to describe total particle hand, Neumann conditions,
absorption
at the boundary;
on the other
which insure that the particle current vanishes at x = 0, describe total reflection at the origin. Such boundary conditions introduce nontrivial features in the behaviour of the system. In the first place, stationary, time independent nonhomogeneous solutions are now possible. Under conditions on their stability, these stationary density profiles are expected to coincide with the long-time asymptotic solutions. Furthermore, the time evolution of the system is also severely influenced by the boundary conditions: the front-like initial structures considagainst the boundary, ered in the previous section are shown to “bounce” drastically modifying the subsequent studied in the following sections. 4.1.
dynamics.
Details
on these
topics
are
Stationary solutions and their stability
In order to study the stationary solutions of the ballast model with the boundary conditions (17) or (18), we come back to the original form of the reaction-diffusion equation, Eq. (3), and impose a, = 0, obtaining
(19) If we associate time with x and coordinate with 4, this has the form Newton equation describing the motion of a particle under the action potential
V($) = -44’ + 4h(4 - 4%)et4 - A).
of a of a
(20)
390
S.A.
This mechanical
Hassan et al. I Reaction-diffusion
analogy
can be fruitfully
patterns in domains
exploited
stationary reaction-diffusion problem [9]. It is clear that a homogeneous stationary
solution
to find the solution
to the
4,(x) = 0 always
exists.
The character of the nontrivial solutions depends on the boundary condition and on the value of the parameter z, defined in Eq. (13). Regarding their stability, here it will be analyzed in the linear approximation, by studying the evolution
of a small perturbation
applied
to the stationary
4.1.1. Dirichlet boundary condition When the density is required to satisfy Eq. (17), besides a nonhomogeneous solution exists for z > 0. It reads
solution.
the trivial
solution,
forx
(21)
where
(22) is the stationary critical coordinate, at which the stationary density reaches the critical value $=. This solution clearly satisfies the Dirichlet condition, and grows monotonically in space, approaching c#+,as x* +m. A linear stability analysis indicates that the trivial solution is stable for any z, and that the nonhomogeneous solution (21) is also stable when it exists. Therefore, for z > 0 the system is bistable, and the long-time behaviour of the solutions will be determined by the initial condition. 4.1.2. Neumann boundary condition The requirement of a,~$ = 0 at the origin makes it possible to have a homogeneous stationary solution different from the trivial one. It corresponds to the second root of the reaction term in our reaction-diffusion model,
4,(x) = 4%. For .z < 0, a nonhomogeneous
(23) solution
is proven
to exist also. It reads
forx
where
(24)
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
391
(25) is, as for Dirichlet
conditions,
the
stationary
critical
coordinate.
density has a finite value at x = 0, namely, 4,(O) = 4,/cash x, < 4,. monotonically for increasing x and asymptotically approaches &,.
Now,
the
It grows
This nonhomogeneous stationary solution proves to be linearly unstable. On the other hand, both homogeneous solutions are stable for any value of z. Thus, independently of the parameters 4, and c$,, (4, > +C), Neumann conditions determine the system to be bistable. Again, its long-time evolution will be selected by the initial condition.
4.2.
Short-time and long-time behaviour
As we could expect, the boundary condition at x = 0 does not affect the short-time evolution of any initial density distribution. Front structures are thus shown to reproduce the initial behaviour observed in free space: they go forward if z > 0 and move backward if z < 0. In some particular situations, also the long-time evolution of the density is well described by a free-space solution. We recall that, in free space, the solution for a step-shaped initial condition evolves towards a limiting steady wave front, moving with velocity u = -22/m. In the semi-infinite domain, this behaviour is reproduced for the same initial condition when the density tends to a vanishing stationary state: the front structure moves towards x-+ +m and, as it separates from the boundary, it approaches the shape and the velocity of the free-space wave front solution, Eqs. (15) and (16). On the other hand, when nonvanishing stationary states are approached from bubbletype initial conditions, the rightmost front moves also towards x+ +m and tends to reproduce the free-space behaviour. These limiting features in the bounded solutions, which exist from an analytical study of Eq. (5), are confirmed calculations presented in the following section.
4.3.
can be proven to by the numerical
Numerical study for intermediate times
We present here numerical results for the bounded system in the intermediate regime. For the step initial condition, Eq. (lo), we take x,(O) = 1; for the bubble initial condition, we choose x,(O) = 1 and x2(O) = 2. These cases exhibit all the qualitatively different features that we were able to observe for other values of the initial critical coordinates. Thus, they constitute a proper summary of more general situations.
392
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
1.03
/
0.60
0
2
4
8
6
10
time Fig. 3. Evolution of the critical x,(O) = 1 and Dirichlet boundary
coordinate condition
for the step-shaped at x = 0.
initial
condition,
Eq.
(lo),
with
4.3.1. Dirichlet boundary condition Fig. 3 displays the evolution of the critical coordinate x,(t) for the step initial condition with Dirichlet boundary conditions, for various values of z. As predicted by the analytical results, for z > 0, x,(t) tends to a limiting value, which can be shown to coincide with the critical coordinate X, of the stationary solution, Eq. (21). For negative z, on the other hand, the front separates from the boundary, approaching the free-space wave-front evolution. At intermediate times nontrivial features appear in the evolution. In particular, for sufficiently small, positive values of z the motion of the critical coordinate is nonmonotonic. It begins by approaching the boundary, but at a certain point it “bounces” and, from then on, it moves backward. For small negative values of z (see the curve for z = -0.03) a sign of this “bounce” is preserved: the time derivative of x,(t) evolves nonmonotonically. Results for the bubble initial condition are plotted in Fig. 4. Here, the lower part of each curve corresponds to xl(t), and the upper one, to x2(t). The bounce effect plays now a determinant role in the survival of the bubble. Consider, for example, the curve for z = 0.30. Up to t- 1.5, the upper critical coordinate x2 seems to have escaped from the influence of the other front, and to approach the wave-front behaviour observed in free space. However, at t ~0.3, xl(t) has bounced back, and approaches now the other critical coordinate at a velocity higher than that of x2(t). Finally, the fronts collide and annihilate at time t,.
S.A.
Hassan et al. I Reaction-diffusion
393
patterns in domains
3.0 2.5 2.0
8
1.5
2 1 .o
0.5
z=oso 0.0
j
0.0
”
1 ” 0.5
”
I ““““‘1 1 .o
”
2.0
1.5
”
1 ”
2.5
”
”
3.0
”
3.5
time Fig. 4. As in Fig. 3, for the bubble-shaped initial condition, Eq. (ll), with x,(O) = 1 and x,(O) = 2. The lower part of each curve corresponds to x,(t), and the upper part, to x,(t). The small plot in the upper-right corner shows the front annihilation time t, as a function of z.
This annihilation time depends on the value of z, and has been plotted in the upper-right small figure. We see that ta(z)+ CCfor z = 0.31. For higher values of z the upper front is able to escape from the other and the bubble survives and grows. In this case, x,(t) approaches the limiting value X, predicted for the stationary solution, and x*(f) acquires a constant velocity whose modulus equals that of Eq. (16). Observe that, even in this case, the intermediate evolution can be highly nontrivial: for z = 0.35, for instance, xl(t) exhibits a kind of damped oscillation as it approaches x,.
4.3.2. Neumann boundary condition We recall that, for Neumann conditions, there exists no stable stationary nonhomogeneous solution, so that the solution must approach one of the two homogeneous states, 4 = 0 or 4 = +,,. As we shall see, this fact determines essential differences with the case of Dirichlet conditions. In Fig. 5, we plotted the critical coordinate for the step initial condition. For sufficiently
negative
values
of z the front moves
as if it were in free space,
and
394
S.A.
0.0
Hassan et al. I Reaction-diffusion
0.5
1.0
1.5
patterns in domains
2.0
2.5
3.0
time Fig. 5. As in Fig. 3, for Neumann
boundary
conditions
at the origin.
the solution tends to the vanishing state. Approaching z = 0, however, the boundary condition begins to play a role in the evolution of the front, and although the critical coordinate moves initially with positive velocity, it eventually inverts its motion and reaches x = 0 in a finite time. From then on, the solution is greater than $, for all x, and approaches asymptotically the homogeneous state $,(x) = 4,. This is also the limiting behaviour displayed for z > 0, where a nonmonotonic time dependence is observed in the derivative of x,(t). Note that the critical coordinate reaches the boundary at infinite velocity. The features observed for the step initial condition determine the evolution of the bubble distribution. Results for this situation are presented in Fig. 6. Now, with respect to the Dirichlet case, the fronts interchange in some sense the role they play in the evolution. In fact, the annihilation of the bubble for z >O is here a consequence of the change of the sign of the upper front velocity. It moves initially with positive velocity but, as observed for the step distribution, it eventually stops and moves towards the origin. Then, it encounters the other front, and their interaction leads to mutual annihilation at a time t,(z). This time diverges at z = 0.216. For higher z, the bubble grows, x,(t) reaches the boundary, x2(t) adopts the wave-front motion, and the solution approaches the nonvanishing homogeneous state. For Neumann boundary conditions, both fronts can display a rather complex time dependence, with the damped oscillatory behaviour already observed with a Dirichlet boundary.
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
395
2.5
0.3 -0.1 4.l 0.0 0.1 0.2
2.0
1.5
z
z
1.0
0.5
0.0 0.0
0.3
0.6
0.9
1.2
1.5
1.8
conditions
at the origin.
time Fig. 6. As in Fig. 4, for Neumann
boundary
5. Conclusions
In this paper we have addressed the problem of pattern propagation within reaction-diffusion models. The system we have chosen for our analysis is a simplified version of the ballast resistor, whose piecewise-linear reaction term mimics the behaviour of other more complicated models (i.e., the Schliigl model and other bistable ones). The scheme we have used is a novel one based on a functional integral approach, leading to an integral equation for the relevant density. This scheme allows us to obtain analytic asymptotic results as well as to analyze the evolution along the whole range of times by means of a suitable numerical procedure. Although our interest is the study of structure propagation iri general nonlinear reaction-diffusion systems composed of one or more species, here the method has been introduced by discussing the simple one-component one-dimensional bistable model indicated above. We have considered the case of fronts connecting different states in a bistable model, as well as bubble-like structures. Our study has involved infinite systems on one hand, and bounded systems on the other. In this latter case, we paid particular attention to the effect of totally absorbing (Dirichlet) and totally reflecting (Neumann) boundary conditions on the propagation of structures. For the case of free-space evolution, numerical results agree quite well with
S.A.
396
Hassan et al.
I Reaction-diffusion
patterns in domains
the analytic limits at short and long times, also showing the tendency to reach the expected asymptotic wave-front behaviour. For bubble-like structures, our results coincide with previous ones, obtained through different methods (see, for instance, [ 121). On the other hand, for the semi-infinite system, we have found that the boundary condition influences the time evolution of the structures in some unexpected ways, for example, with the appearance of bounces and damped oscillatory-like effects. Other relevant aspects have arisen when analyzing the evolution of bubble-like structures. We have seen that, even for situations where the initial condition in the free-space case would lead to a growing structure, the presence of a boundary can cause its collapse. As the behaviour of bubble structures has been used to describe nucleation phenomena [13], these results indicate that boundary effects would influence the nucleation dynamics in small size systems. Here, we have discussed the evolution of the ballast reaction-diffusion model by looking only at the behaviour of the critical coordinate x,(t) - which corresponds to the front position - as the detailed form of the front profiles is not expected to show particularly interesting features. Nevertheless, in Fig. 7 we show a comparison of the density profiles at a fixed time for a bubble initial distribution in free space and with the two boundary conditions considered here. The features introduced by the boundary are apparent. In conclusion, by analyzing the problem described above, we have shown 0.6
,
/
I
0.5 0.4 0.3 0.2
0.1
0.0
x
Fig. 7. Density profiles for a bubble initial condition with x,(O) = 1 and x2(O) = 2, at time t = 0.8, evolving in free space (F), with Dirichlet conditions (D), and with Neumann conditions (N). In this plot, &, = 1 and z = 0.23. Compare with Figs. 4 and 6.
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
391
that the present propagation scheme is very suitable for the study of evolution in reaction-diffusion systems. The insight so far gained will lead us to the analysis of more complex (many component) problems, as well as more general boundary conditions. These topics are the subject of work in progress.
Acknowledgement
We thank Prof. Veronica Griinfeld for her careful reading of the manuscript.
Appendix
A. Path-integral
solution
to reaction-diffusion
equations
We consider here a general nonlinear reaction-diffusion
equation of the form
where we assume that f[+(x, t)] h as such mild behaviour that we may proceed without restrictions. For instance, it is shown in [14] that if f[+(x, t)] as well as and --oo
4(x,t)=/9[x(t)]&(x)exp[-
Ida
(+,+)2+f~$‘~)1)].
2
fo
(27)
However, this is not the most adequate expression to perform quantitative (numerical) studies on the time evolution of the fronts or spatial structures which arise in the system described by Eq. (26). In order to overcome this drawback and to have a representation of the solution suitable for numerical analysis, we use a path integral-like approach. To obtain the propagation for a given initial condition, we integrate Eq. (26) over a very short time interval (t, t + T) to the lowest order in r and use a prepoint discretization scheme [15] to find
a, t +7) = (1+ ~~:>d&t>+ Tf[4+, 91 . Through the usual procedure
(28)
of Fourier transforming
4(x, t = t’ + T) = ~f[+(x, t’)] +
I--$&exp( -
[15], we obtain from here
& (x - x’)‘) 4(x’, t’) . (29)
S.A.
398
Considering propagation
Hassan et al. I Reaction-diffusion
patterns in domains
always the lowest order in T and iterating Eq. (29), we get the of the initial distribution 4(x, t = to) = c#I,(x) up to the point xr and
to the time t, (tf - t, = NT) as
x f[4(Xk, to+ kT)l>
(30)
where xj = x(t,) is the coordinate at tj = t, + j7. In the continuous limit (N-+ 00, r+ 0 and NT = t, - to) and in the usual notation we finally obtain
%I
~~(“)(x,, t, I-%,to) 4(-ql, to)
If
I I dt'
+
dx’
X(“)(xf, t, 1x’,
t’) f[+(x’,
t’)] ,
(31)
10
where X”‘(x, t 1x’, t’) is the propagator diffusion operator, given, in the discrete
associated case, by
to the (linear)
free-space
(32)
with the continuous
limit
X(“)(x,,tr/xo,to)=~9[x]exp(-$
[dsx(s)‘).
(33)
10 It is easy to prove that the result indicated in Eq. (31) is a solution to Eq. (26) and satisfies the given initial condition. It is clear that the above indicated procedure will also work if we separate in our equation all the linear contributions (call it 24) from the nonlinear contributions (call it N(4)). Eq. (26) becomes a&(x,
t) = =W + J(4)
>
(34)
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
399
and has a solution similar in form to Eqs. (31), (32), (33), where now XC’) is the propagator of the linear operator LZ. In the present case (the ballast model), the nonlinear reaction function has the form f(4) = -4 + 4@(4 - 4,) .
(35)
Then we choose as our linear operator ~4=&-4%
(36)
leaving as the nonlinear part J(4)
(37)
= WV+ - 4%).
After evaluating X(O) for the operator
X(O)(Xf, t, Ixo,4J =
&q
in Eq. (36) we get
exp - y;t-_x;y f -(ff - to,), i
(38)
0
which, replaced in Eq. (31), gives the final form used in this work.
References PI P.C. Hohenberg and M.C. Cross, An introduction to pattern formation in nonequilibrium systems, in: Fluctuations and Stochastic Phenomena in Condensed Matter, L. Garrido, ed. (Springer, Berlin, 1986); A.C. Newell, The dynamics and analysis of patterns, in: Complex Systems, D. Stein, ed. (Addison-Wesley, New York, 1989). PI G.T. Dee, J. Stat. Phys. 39 (1985) 705; G.T. Dee and W. van Saarloos, Phys. Rev. Lett. 60 (1988) 2641; M. Nickles, M. Lticke and H. Miiller-Krumbhaar, Europhys. Lett. 9 (1989) 237; W. van Saarloos, Phys. Rev. A 37 (1988) 211; 39 (1989) 6367. [31 P. Hagan, SIAM J. Math. Anal. 57 (1982) 717; J. Rinzel and D. Terran, SIAM J. Appl. Math. 42 (1982) 1111; J.D. Dockery and J.P. Keener, SIAM J. Appl. Math. 49 (1989) 539; H. Ikeda and M. Mimura, SIAM J. Appl. Math. 49 (1989) 515; H. Yamada and K. Nozaki, Prog. Theor. Phys. 84 (1990) 801. [41 L. Kramer and P.C. Hohenberg, Physica D 13 (1984) 357; P.C. Hohenberg, L. Kramer and H. Riecke, Physica D 15 (1985) 402; M.C. Cross, Phys. Rev. Lett. 57 (1986) 2935. PI W. van Saarloos and P.C. Hohenberg, Physica D 56 (1992) 303; E. Meron, Phys. Rep. 218 (1992) 1. 161 J.J. Tyson and J.P. Keener, Physica D 32 (1988) 327. [71 B. Ross and J.D. Lister, Phys. Rev. A 15 (1977) 1246; R. Landauer, Phys. Rev. A 15 (1977) 2117. PI D. Bedeaux and P. Mazur, Physica A 105 (1981) 1.
400 [9] [lo] [ll] [12] [13]
S.A.
Hassan et al. I Reaction-diffusion
patterns in domains
C. Schat and H.S. Wio, Physica A 180 (1992) 295. F. Schlogl, Z. Phys. 248 (1971) 446; 253 (1972) 147. Ch. Zulike, A.S. Mikhailov and L. Schimansky-Geier, Physica A 163 (1990) 559. A.S. Mikhailov, Foundations of Sinergetics I (Springer, Berlin, 1991). L. Schimanski-Geier and W. Ebeling, Ann. Phys. (Leipzig) 40 (1983) 10; L. Schimansky-Geier et al., Ann. Phys. (Leipzig) 40 (1983) 277. [14] M. Freidlin, Functional Integration and Partial Differential Equations (Princeton Univ. Press, New Jersey, 1985). [15] F. Langouche, D. Roekaerts and E. Tirapegui, Functional Integration and Semiclassical Expansions (Reidel, Dordrecht, 1982).