- Email: [email protected]
Transition from pushed-to-pulled fronts in piecewise linear reaction–diffusion systems
ARTICLE IN PRESS
Physica A 375 (2007) 51–64 www.elsevier.com/locate/physa
Transition from pushed-to-pulled fronts in piecewise linear reaction–diffusion systems V. Me´ndez, V. Ortega-Cejas, E.P. Zemskov, J. Casas-Va´zquez Grup de Fı´sica Estadı´stica, Departament de Fı´sica, Facultat de Cie`ncies, Edifici Cc,Universitat Auto`noma de Barcelona, E-08193 Bellaterra (Cerdanyola), Spain Received 14 December 2005; received in revised form 3 July 2006 Available online 4 October 2006
Abstract The front dynamics in reaction–diffusion equations with a piecewise linear reaction term is studied. A transition from pushed-to-pulled fronts when they propagate into the unstable state is found using a variational principle. This transition occurs for a critical value of the discontinuity position in the reaction function. In particular, we study how the transition depends on the properties of the reaction term and on the delay time. Our results are in good agreement with the numerical solutions of the model. r 2006 Elsevier B.V. All rights reserved. Keywords: Pulled fronts; Pushed fronts; Reaction–diffusion
1. Introduction Piecewise linear models for reaction kinetics are used to get exact analytic solutions for the speed and wave profiles in reaction–diffusion equations [1]. These are just an approach of a more realistic but analytically intractable nonlinear reaction–diffusion equations. The piecewise linear emulation consists in considering the reaction term as two or more linear segments reaching the value 0 at each equilibrium state. In the last decade these linear models have been widely employed, for instance, to find a speed selection mechanism for both singlecomponent [2] and two-components [3] parabolic reaction–diffusion (PRD) equations, to study Allee effect in the spread of invading organisms [4] and the effect of transport memory in reaction–diffusion fronts [5] and to analyze the wave propagation in discrete [6] and excitable [7] media and under forcing in bistable systems [8]. In Refs. [1], the piecewise linear emulation is applied to bistable kinetics by considering f H u þ yðu aÞ, a ¼ const. This allows to find analytic solutions for traveling waves like kinks, wave trains, solitary waves and two-dimensional spiral waves described by activator-inhibitor systems. We deal, however, with a more fundamental problem which appears only in fronts propagating into unstable states: the pushed–pulled transition. Although this is an already known problem for nonlinear reaction term with cubic nonlinearities [9] we are interested here to study this problem when one considers the piecewise linear emulation of these nonlinear Corresponding author.
E-mail address: [email protected] (V. Me´ndez). 0378-4371/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2006.09.007
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
52
reaction terms. In Refs. [2,3,5], the authors deal also with fronts propagating into unstable states but no attention is paid into the pushed–pulled transition because these works only deal with analytic solutions. The main advantage of dealing with piecewise linear reaction terms is to get analytical solutions for the front profiles and their speeds but, as we show in this work, these solutions only hold when the values of the parameters lie in the pushed region. However, from the variational principle and numerical solutions developed in this paper, we are able to find the front speed in both regions (pushed and pulled). We consider a function [10] ( nu; 0puoa; f ðuÞ ¼ (1) bð1 uÞ; aoup1; which is an approximation of a general reaction term f ðuÞ with two equilibrium states: u ¼ 0 and u ¼ 1 for the Fisher–Kolmogorov reaction–diffusion equation qt u ¼ qxx u þ f ðuÞ.
(2)
In this work, we deal with a front propagating into the unstable (n40) or metastable (no0) state u ¼ 0 connected to the stable state u ¼ 1 by means of Eq. (2), and by its generalization to the case with a time delay, with f ðuÞ given by Eq. (1). We calculate the front profiles and the front speed equation, once Eq. (2) has been transformed into an ordinary differential equation by using the traveling wave coordinate z ¼ x ct, being c the asymptotic front speed. We also apply a variational principle to the problems (1) and (2) to calculate lower and upper bounds for the front speed and perform the numerical integrations employing an explicit finite difference method to obtain the front profiles and c. Both techniques allow us to find a pushed-to-pulled transition when the front propagates into the unstable state. We say the front is pulled when it moves into an unstable state with constant speed (after a transient) given by the linear analysis around the unstable state. We refer to this speed as linear speed. On the other hand, we say the front is pushed when its speed after the transient is higher than the linear one [9]. We end this work by illustrating how to select the parameters n, b, and a in (1) to emulate some specific examples of nonlinear reaction terms. We show that the piecewise model captures the essential qualitative dynamical properties of the system under emulation whenever it preserves the relevant properties of the real nonlinear function. 2. Variational principle Variational principles have been used to find bounds for the speed of reaction–diffusion fronts. They are very useful to obtain a bound for the speed when the linear marginal stability fails and the speed selected is strongly affected for the nonlinear terms. In our model, the monotonically decreasing fronts u ¼ uðz ¼ x ctÞ propagating with constant speed c join u ¼ 1 with u ¼ 0: The front satisfies u00 þ cu0 þ f ðuÞ ¼ 0; where the prime symbol stands for derivative respect to z and fills the boundary conditions uðz ! 1Þ ¼ 0 and uðz ! 1Þ ¼ 1. From a variational principle, the asymptotic front speed is given by [11] " R1 # 2 0 fg du c ¼ max 2 R 1 , (3) g 2 =hÞ du ðg 0 where the maximum is taken over all positive decreasing auxiliary functions gðuÞ in ð0; 1Þ (h ¼ dg=du40) for which the integrals involved in Eq. (3) exist. It is often impossible to know the form of the function g ¼ gðuÞ which maximizes the relationship in brackets in Eq. (3) and one can find only lower bounds R1 fg du 2 cL ¼ 2 R 1 0 . (4) 2 0 ðg =hÞ du ItR is easy toR get an upper R 1 bound for the speed independent of g [12]. From R 1 the Schwarz R 1 inequality 1 1 ð 0 ug duÞ2 p 0 ðg2 =hÞ du 0 u2 h du and assuming gð1Þ ¼ 0 in the integration by parts 0 u2 h du ¼ 2 0 ug du, then " R1 # " R1 # fg du fg du f ðuÞ 2 p4 max 2 R 01 p4 sup c ¼ max 2 R 1 0 g g u u2ð0;1Þ ðg2 =hÞ du ug du 0
0
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
and the upper bound is sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi f ðuÞ cU ¼ 2 sup u u2ð0;1Þ
53
(5)
in agreement with the result of Aronson and Weinberger [13]. In this work, we are going to deal with the family of functions gðuÞ ¼ ua1 1 which are decreasing functions if ao1 and fulfil gð1Þ ¼ 0: To get the nearest lower bound to true speed value we will take the maximum value over all the values for the parameter a such that the integrals in (4) are convergent. So, calculating the lower bound from Eq. (4) for f ðuÞ in form (1) one gets c2L ¼ max1oao1 Gða; a; n; bÞ where aþ1 2ð1 þ aÞð3 aÞ a a2 1 b 1 a ðn þ bÞ a þ Gða; a; n; bÞ ¼ . (6) þb a 1a 2 a 1þa aþ1 2
2.1. Propagation into unstable state (n40) ðLÞ For n40, there exists a critical value aðLÞ c , such that when aoac the maximum of G, i.e., of cL ; is attained at ðLÞ a ¼ a with 1oa o0. However, when aXac the maximum is reached always at a ¼ 1 and just at this point one has Gð1; a; n; bÞ ¼ 4n: To compute the critical position of the discontinuity from the lower bound aðLÞ c ; we take dG=da ¼ 0 at a ¼ 1 and find ðLÞ 2 ðLÞ ðLÞ 2 ðLÞ 2aðLÞ c ðn þ bÞ½ðac Þ 2 ln ac 4b½1 þ ðac Þ ac ðn 6bÞ ¼ 0.
To sum up, we find from the variational principle that the lower bound for the speed is given by ( pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Gða ; a; n; bÞ; aoaðLÞ c ðn; bÞ; pffiffiffi cL ¼ 2 n; a4aðLÞ c ðn; bÞ;
(7)
(8)
where Gða ; a; n; bÞ44n; a is such that dG=da ¼ 0 at a ¼ a and aðLÞ c ðn; bÞ is solution of Eq. (7). Let us now to evaluate the upper bound for the front speed. From Eqs. (1) and (5) we get 8 pffiffiffi b > > ; a4acðUÞ ¼ > < 2 n; nþb r ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi (9) cU ¼ > 1a > > b ; otherwise; 2 : a ðUÞ with aðLÞ c oac . This result shows, together with (8), that the speed of the front is linearly selected for a sufficiently large values of a and it is nonlinearly selected for sufficiently low values of a. Therefore, there exists a transition from pushed-to-pulled front when a goes from 0 to 1. The transition is found at a critical value of the position of the discontinuity. From the variational principle we have found two possible values: aðLÞ c from the lower bound (Eq. (8)) and aðUÞ from the upper bound (Eq. (9)). c
2.2. Propagation into metastable state (no0) For no0, the variational principle shows that there is no linearly selected speed and neither a transition pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi from pushed-to-pulled fronts. The speed is always nonlinearly selected and it is given by c ’ max1oao1 G : For large enough values of a the maximum is attained at a ¼ 1: In this region there exists a maximum value of ^ such that for a ¼ a^ one has Gð1; a; ^ n; bÞ ¼ 0 and as a result c ¼ 0: From Eq. (6) it is obtained that a^ a, namely a, satisfies the equation 2½2b aðn þ bÞa ln a þ a2 ðn þ bÞ þ bð3 4aÞ ¼ 0: For a4a^ one finds Go0 and the variational principle is not able to supply lower bounds for the front speed, because the square root of the negative function G appears. In consequence, our variational principle only can predict positives values for c.
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
54
3. Analytic solution In this section we find the analytic solution of the equation u00 þ cu0 þ f ðuÞ ¼ 0; where the reaction term is given by Eq. (1). As the differential equation for the monotonic asymptotic fronts is linear, the solution involves two exponential functions for each piece. In consequence, four integration constants have to be found from boundary and matching conditions. In the first piece f ðuÞ ¼ nu; where 0puoa and we take u ¼ u1 ðzÞ with z0 pzo1: In the second piece f ðuÞ ¼ bð1 uÞ with aoup1 and we have u ¼ u2 ðzÞ with 1ozpz0 : The boundary conditions for each piece are then u1 ðz ! 1Þ ¼ 0 and u1 ðz0 Þ ¼ a for the first piece and u2 ðz ! 1Þ ¼ 1 and u2 ðz0 Þ ¼ a for the second piece. Therefore, the solutions for the front profiles read pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi c c2 4n lþ z l z u1 ðzÞ ¼ A1 e þ A2 e with l ¼ and c42 n, 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c c2 þ 4b þ . ð10Þ u2 ðzÞ ¼ 1 þ B1 em z þ B2 em z with m ¼ 2 Consider first the front propagation into the unstable state, i.e., the case with n40: As l o0 the boundary condition u1 ðz ! 1Þ ¼ 0 does not allow us to find any integration constant. For the second piece the boundary condition u2 ðz ! 1Þ ¼ 1 requires B2 ¼ 0 because mþ 40 and m o0: One can find two of the three integration constants from two boundary conditions u1;2 ðz0 Þ ¼ a and one matching condition for the derivative, u01 ðz0 Þ ¼ u02 ðz0 Þ; one can obtain the third one. All integration constants may be determined but there is no any additional condition to determine the equation for the front speed; c is then another free parameter and no unique front exists. Alternatively, when the front propagates into the metastable state (no0), then l o0 but lþ 40 and A1 ¼ 0 according to the boundary condition u1 ðz ! 1Þ ¼ 0: Thus, in this case, the front speed may be determined from the matching condition for the continuity of the derivative at z0 and an unique front exists. In order to determine the front speed for the propagation into the unstable state in the context of the piecewise linear model, we assume initially the boundary conditions u1 ðz ¼ LÞ ¼ 0 and u2 ðz ¼ LÞ ¼ 1; instead of u1 ðz ! 1Þ ¼ 0 and u2 ðz ! 1Þ ¼ 1; and finally we will take L ! 1 [2]. The expressions for the integration constants read
þ
el L el L ; A ¼ a , 2 þ þ el ðz0 LÞ el ðz0 LÞ el ðz0 LÞ el ðz0 LÞ þ em L em L B1 ¼ ð1 aÞ m ðz þLÞ ; B ¼ ð1 aÞ . ð11Þ 2 þ e 0 em ðz0 þLÞ em ðz0 þLÞ emþ ðz0 þLÞ From the matching condition for the continuity of the derivative in z0 ; one finds, after some algebra an equation for the family of possible front speeds pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 1 c2 þ 4bðz0 þ LÞ . c ¼ a c2 4n coth (12) c 4nðz0 LÞ þ ð1 aÞ c2 þ 4b coth 2 2 A1 ¼ a
Taking the limit at L ! 1 one gets the selected speed which is independent on z0 , according to the translational invariance of the model equation [2]. Thus, the speed equation is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ a c2 4n þ ð1 aÞ c2 þ 4b. (13) Taking L ! 1; the integration constants become A1 ¼ B2 ¼ 0; A2 ¼ a expðl z0 Þ and B1 ¼ ð1 aÞ expðmþ z0 Þ and the exact solution for the front profiles are finally ( ael ðz0 zÞ ; z 2 ½z0 ; 1Þ; uðzÞ ¼ (14) mþ ðz0 zÞ 1 ð1 aÞe ; z 2 ð1; z0 : pffiffiffi It is worthy to stress that if we assume co2 n in Eq. (10) the solution to u1 ðzÞ involves sines and cosines while u2 ðzÞ remains the same. By imposing the matching conditions for u1 ðzÞ and u2 ðzÞ at z ¼ z0 and z ¼ L one gets the corresponding integration constants. By requiring the matching condition for the derivatives at z ¼ z0 one obtains the speed equation. However, when one takes the limit L ! 1 it is not possible to get an equation equivalent to (13) to get the selected speed.
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
55
3.1. Propagation into unstable state (n40) Eq. (13) is an implicit relationship for the selected front speed. It is important to note that this equation has pffiffiffi real solutions only if c42 n, as one can easily check from the first term of the right-hand side of Eq. (13). This means that the selected front speed from (13), once it is explicitly expressed in terms of n, b and a, has to pfound ffiffiffi be higher than the linear speed 2 n. Solving Eq. (13) for c one finds the selected front speed a2 ðn þ bÞ þ bð1 2aÞ c ¼ pffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . a a2 ðn þ bÞ ðn þ 2bÞa þ b
(15)
pffiffiffi Introducing this result in the restriction c42 n, this inequality turns into rffiffiffiffiffiffiffiffiffiffiffi n ðeÞ . aoac ¼ 1 nþb
(16)
This means that the exact solution (15) is restricted to values of a such that a 2 ð0; aðeÞ For c : p ffiffiffi this range of values of a the front is pushed because the speed given in (15) is higher than the linear speed 2 n as we have required. From the variational principle applied to (1) and (2) we have shown in Eqs. (8) and (9) that a monotonic front and its speed exist for any value of a in (0,1) exhibiting a pushed–pulled transition. This result has been also supported by numerical solutions. The exact solution is not able to describe this transition because (15) does not hold for a 2 ðaðeÞ c ; 1Þ and nothing can be known only from the exact solution. To illustrate this fact, we plot in Fig. 1 the front speed obtained from the exact solution (15) and the lower and upper bounds found from the variational principles (8) and (9), respectively, versus the position a of the discontinuity in the reaction term. We depict with symbols the numerical solution for the front speed computed directly from Eqs. (1) and (2) with an initial condition in the form of the Heaviside function between 0 and 1. In dashed lines we depict the exact solution, which is in very good agreement with the numerical calculations but only for ðeÞ a 2 ð0; aðeÞ c , because this curve stops at the point a ¼ ac and cannot go beyond this point. The bounds are
6
ν = 0.5
c
5
cU
4 3
(U)
cL
2 (a)
ac
6
(L)
ac
(e)
ac
ν=1
c
5
cU
4
(U)
ac
3 2
cL 0.0
(b)
0.2
(L)
ac
(e)
ac
0.4
0.6
0.8
1.0
a
Fig. 1. Dimensionless selected front speed from the exact solution (15) (dashed lines), the lower (cL ) and upper (cU ) bounds from the variational principle (VP) and the numerical solutions (symbols) versus the position a of the discontinuity in the reaction term. In (a) n ¼ 12 and b ¼ 1. In (b) n ¼ b ¼ 1.
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
56
depicted in solid lines. The lower bound shows also a good agreement with the numerical solution. Both bounds join within ½acðUÞ ; 1 describing the pushed–pulled transition because from Eqs. (8) and (9) pffiffiffi the range ðUÞ c ¼ cL ¼ cU ¼ 2 n for a 2 ½ac ; 1. While the lower bound describes the transition at a ¼ aðLÞ c ; the upper ðLÞ bound describes the transition at a ¼ aðUÞ 4a . In Table 1, we compare both values, for the critical position c c of the discontinuity, with the obtained from the analytic solution (Eq. (15)) and from the numerical solution. ðeÞ ðnumÞ It is observed a good agreement between aðLÞ . In Fig. 2, we plot the c and ac with the numerical solution ac critical value of a in the pushed-to-pulled transition obtained from the lower bound (aðLÞ c from Eq. (7)) and from the exact solution (aðeÞ from Eq. (15)) together with numerical computations. It is shown that both exact c and variational results yield very similar results and in agreement with numerical solutions. It is also interesting to note that the critical position moves towards 0 as the slope of the first segment n increases keeping fixed the slope of the second segment b. Analogously, for fixed n when b increases, it moves towards 1. In Fig. 3, we represent two front profiles using the exact solution (14) (solid lines) together with the profiles obtained from numerical integrations (symbols) exhibiting an excellent agreement.
Table 1 Table of the results for the critical position of the discontinuity for different values of the slope of the first segment n for the case b ¼ 1 n
acðLÞ
aðUÞ c
aðeÞ c
aðnumÞ c
0.1 0.3 0.5 0.7 0.9 1
0.67 0.51 0.42 0.35 0.30 0.28
0.91 0.77 0.67 0.59 0.53 0.50
0.70 0.52 0.42 0.36 0.31 0.29
0.75 0.55 0.45 0.35 0.35 0.25
ðUÞ is obtained from the upper bound (Eq. (9)), aðeÞ aðLÞ c is the critical value obtained from the lower bound (Eq. (7)), ac c is obtained from the ðnumÞ analytic solution (Eq. (15)) and ac is computed numerically.
1.0
0.8 β=2
ac
0.6 β=1 0.4
β = 0.5
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
ν Fig. 2. Dimensionless critical discontinuity position versus the slope n of the first piece for three different values of the slope b of the ðLÞ second piece. The dashed curves correspond to aðeÞ c ; Eq. (15) and the solid curves correspond to ac ; Eq. (7). The symbols correspond to numerical solutions (circles for b ¼ 0:5; squares for b ¼ 1 and triangles for b ¼ 2).
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
57
1.0 ν = 0.5, β = 1
u(z)
0.8 0.6 0.4 0.2 0.0
(a) 1.0 ν=β= 1
u(z)
0.8 0.6 0.4 0.2 0.0 -14
-12
-10
(b)
-8
-6
-4
-2
0
2
4
z
Fig. 3. Front profiles versus z ¼ x ct z0 : With solid lines we depict the results obtained from the exact solution (Eq. (14)) and with symbols we depict the numerical solutions. For (a) we take n ¼ 0:5, a ¼ 0:1 and for these values the speed is c ¼ 2:91. In this case we take a front profile at t ¼ 120 with z0 ¼ 50:35: For (b) we take n ¼ 1, a ¼ 0:1 and for these values the speed is c ¼ 3:02 and we take a front profile at t ¼ 120 with z0 ¼ 50:75:
3.2. Propagation into metastable state (no0) When the piecewise linear model describes a front propagating into a metastable state, no critical value for a is found, neither from variational principles nor from numerical solutions and the analytic solution for the speed (15) holds for any a: In this situation, the front never selects the linear speed. These results are in agreement with that known for fronts propagating into metastable states described by reaction terms with cubic nonlinearities [9]. Therefore, there is no transition from pushed-to-pulled front. The front speed from the variational principle (lower and upper bound found from Eqs. (8) and (9)) and from the analytic solution (found from Eq. (15)) are plotted in Fig. 4 together with numerical calculations. In this case, the front speed computed from the exact solution is in a better agreement with numerical integrations than the bounds obtained from variational principle. The variational principle for no0 can only predict positive bounds for the front speed, as we have shown in Section 2.2. Therefore, it is unable to account, contrarily to the analytic result (15), for the change of sign in c, characteristic of some bistable systems. 4. Effect of delay In this section, we extend our analysis to study how depends ac on a time delay w. For simplicity we take b ¼ 1 in Eq. (1). To do this, we consider the hyperbolic reaction–diffusion equation [14] wqtt u þ qt u ¼ qxx u þ f ðuÞ þ wqt f ðuÞ,
(17)
where the delay time is included in w. In fact, w is the quotient between the delay time and the characteristic reaction time. As above, let us to begin with the variational principle. The front satisfies ð1 wc2 Þu00 þ c½1 wf 0 ðuÞu0 þ f ðuÞ ¼ 0 and the boundary conditions uðz ! 1Þ ¼ 0 and uðz ! 1Þ ¼ 1: Hence, from (17) one finds that the monotonic fronts are solutions of ð1 wc2 Þp
dp c½1 wf 0 ðuÞp þ f ðuÞ ¼ 0, du
(18)
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
58
6 ν = -1
5 4 c
cU
analytic
3 2
cL
1 0 (a)
0.0
0.1
0.2
0.3
0.4
6 ν = -0.5
5
c
4
cU
analytic
3 2
cL
1 0 0.0
0.1
0.2
0.3
0.4
a
(b)
Fig. 4. Selected front speed versus a. In solid lines, we plot the speed obtained from the exact solution (15) and with dashed curves we plot lower and upper bounds obtained from the variational principle. Symbols corresponds to numerical solutions: (a) corresponds to the case n ¼ 1 and b ¼ 1; and (b) corresponds to n ¼ 0:5 and b ¼ 1:
with pð0Þ ¼ pð1Þ ¼ 0 in (0,1). Multiplying Eq. (18) by gðuÞ and integrating between u ¼ 0 and u ¼ 1 we obtain, Z 1 Z 1 1 gf du ¼ fðpÞ du with fðpÞ cgp½1 wf 0 ðuÞ ð1 wc2 Þhp2 (19) 2 0 0 after integrating by parts. As c; p; g; and h are positive functions for any fixed u in (0,1), fðpÞ has a maximum at pmax ¼
cg½1 wf 0 ðuÞ 40 hð1 wc2 Þ
(20)
if 1 co pffiffiffi w
and
1 wf 0 ðuÞ40
(21)
or analogously if woM 1
with M ¼ max f 0 ðuÞ. u2½0;1
In consequence, fðpÞp
c2 g2 ½1 wf 0 ðuÞ2 2hð1 wc2 Þ
for any fixed u in (0,1). Finally, it follows R1 gf du c2 X2 R 1 g2 0 1 wc2 ½1 wf 0 ðuÞ2 du 0 h
(22)
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
and the lower bound is given by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R1 gf du A with A ¼ 2 R 1 g2 0 . cL ¼ 1 þ wA ½1 wf 0 ðuÞ2 du
59
(23)
0 h
To prove that this is a variational principle we must to ensure that there exists a function g for which the equality in (22) holds. From Eq. (20) we see that the case of equality is attained when g^ fulfills the following differential equation: pð1 wc2 Þ g^ . 0 ¼ c½1 wf 0 ðuÞ g^ ^ obtained by integrating this equation is The maximizing g; Z u c 1 wf 0 ðuÞ ^ ¼ exp du gðuÞ 1 wc2 u0 p
(24)
with 0ou0 o1: We must ensure now that the integrals in Eq. (22) are convergent. The problems may appear when one computes the integrals near u ¼ 0. In this limit plu with pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c½1 wf 0 ð0Þ þ c2 ½1 þ wf 0 ð0Þ2 4f 0 ð0Þ l¼ 2ð1 wc2 Þ pffiffiffiffiffiffiffiffiffi ffi whenever c42 f 0 ð0Þ=½1 þ wf 0 ð0Þ to ensure that l 2 R. Then, from (24) we get k ^ gðuÞu
with k ¼
c½1 wf 0 ð0Þ lð1 wc2 Þ
(25)
near u ¼ 0 and gf ðg2 =hÞ½1 wf 0 ðuÞ2 u1k . Finally, the integrals in Eq. (22) exist if ko2 which is fulfilled pffiffiffiffiffiffiffiffiffiffi provided that c42 f 0 ð0Þ=½1 þ wf 0 ð0Þ. In consequence one can conclude that the asymptotic front speed is given by ( ) R1 gf du c2 ¼ max 2 R 1 g2 0 . (26) g 1 wc2 ½1 wf 0 ðuÞ2 du 0 h
It is not difficult to obtain an upper bound for the speed independent of g. From the inequality Z 1 2 Z 1 2 g g 0 0 2 2 ½1 wf ðuÞ duX inf ½1 wf ðuÞ du u2ð0;1Þ h 0 0 h R1 R1 R1 R1 and Schwarz inequality ð 0 ug duÞ2 p 0 ðg2 =hÞ du 0 u2 h du and assuming gð1Þ ¼ 0 one obtains 0 ug dup R 1 the 2 0 ðg2 =hÞ du, then ( ) R1 c2 0 gf du ¼ max 2 R 1 g2 0 2 g 1 wc2 0 h ½1 wf ðuÞ du " R1 # 2 0 fg du p max R 1 2 inf u2ð0;1Þ ½1 wf 0 ðuÞ2 g 0 ðg =hÞ du f ðuÞ supu2ð0;1Þ u p4 inf u2ð0;1Þ ½1 wf 0 ðuÞ2 and the upper bound is f ðuÞ sffiffiffiffiffiffiffiffiffiffiffiffiffiffi supu2ð0;1Þ B u with B ¼ 4 . cU ¼ 1 þ wB inf u2ð0;1Þ ½1 wf 0 ðuÞ2
(27)
ARTICLE IN PRESS 60
V. Me´ndez et al. / Physica A 375 (2007) 51–64
In this work, we are going to deal with the family of functions gðuÞ ¼ ua1 1 which are decreasing functions if ao1 and fulfil gð1Þ ¼ 0: To get a lower bound nearest to true speed value we will take the maximum value over all the values for the parameter a such that the integrals in (23) are convergent. So, calculating the lower pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bound from Eq. (23) one gets cL ¼ max1oao1 A=ð1 þ wAÞ where aþ1 a a2 1 1 aa þa ðn þ 1Þ þ að1 þ aÞ 2 aþ1 2 a . (28) Aða; w; a; nÞ ¼ ð1 aÞ aþ1 3a aþ1 3a a a 1 a 1 a 2 2 2 þ þ 1þa ð1 wnÞ þ ð1 þ wÞ aþ1 3a aþ1 3a ðLÞ When n40, there exists a critical value aðLÞ c , such that when aoac the maximum of A, i.e., of cL ; is attained at a ¼ a with 1oa o0 and the front speed is higher than that linearly selected (pushed front). However, when aXacðLÞ the maximum is reached always at a ¼ 1 and just at this point Að1; w; a; nÞ ¼ 4n=ð1 wnÞ and the pffiffiffi front speed is that linearly selected (pulled front) c ¼ 2 nð1 þ wnÞ1 . Therefore, this shows that the lower bound for the front speed describes a pushed-to-pulled transition. However, for no0 the lower bound does not describe transition as already occurs for w ¼ 0. In the following we will take n ¼ 1 without loosing of generality. The critical position of the discontinuity from the lower bound aðLÞ c is such that dA=da ¼ 0 at a ¼ 1, so that it follows the transcendent equation 2 ðLÞ ðLÞ ðLÞ ðLÞ 2 ðLÞ 5 ðLÞ 3 ðLÞ 2 ðLÞ ½4 þ 4ðaðLÞ c Þ ð1 ac Þ 5ac þ 8ac ln ac ð1 þ w Þ þ w½4ðac Þ 8ðac Þ 8ðac Þ þ 22ac 8 ¼ 0. (29)
Note that for w ¼ 0 in (29), one recovers (7) taking n ¼ 1: To sum up, we have found that the lower bound for the speed is given by 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Aða ; w; a; 1Þ > ðLÞ > > < 1 þ wAða ; w; a; 1Þ; aoac ðwÞ; (30) cL ¼ 2 > ðLÞ > > ; a4a ðwÞ; :1 þ w c where Aða ; w; a; 1Þ44=ð1 wÞ; a is such that dA=da ¼ 0 at a ¼ a and aðLÞ c ðwÞ is solution to Eq. (29). Let us now to evaluate the upper bound for the front speed. From Eq. (27) we get 8 2 1 > ðUÞ > > rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; apac ¼ ; > 2 a > < ð1 wÞ2 þ 4w 1 a (31) cU ¼ > > 2 1 > > ðUÞ > aXac ¼ ; : 1 þ w; 2 This result shows, together with (30), that the speed of the front is exactly the linearly selected for a sufficiently large values of a and it is nonlinearly selected for sufficiently low values of a. Therefore, there exists a transition from pushed-to-pulled front when a goes from 0 to 1. From the variational principle we have found ðUÞ two possible values: aðLÞ from the upper bound (Eq. (31)). c from the lower bound (Eq. (30)) and ac The HRD equation (17) with (1) can be analytically solved as we did in Section 3 with the parabolic case w ¼ 0. The solution is formally the same as for the parabolic equation under an appropriate redefinition of the parameters. So, by imposing matching conditions for u and its derivative at the discontinuity one gets finally the following selected speed equation: " rffiffiffiffiffiffiffiffiffiffiffiffiffi# 1 2a þ 2a2 1 1 þ w2 ðeÞ c ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi if aoac ¼ 1 . (32) 1þw 2 1 2a að1 aÞð1 w2 Þ þ wð1 2aÞ
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
61
It is worthy to stress, as occurs for the parabolic case, that the exact speed given in Eq. (32) only holds for aoaðeÞ c and the front is always pushed. However, from the variational principle it is shown in Eqs. (30) and (31) that the speed exists for any value of a in (0,1) exhibiting a transition from pushed-to-pulled front. As already occurs in (15) and (16) the exact solution (32) does not inform about what happens beyond aðeÞ c . Anyway, it seems clear the inclusion of the time delay does not affect the existence of the transition at all. In Fig. 5, we plot the front speed obtained from the exact solution (32) together with the lower and upper bounds found from the variational principle (30) and (31), respectively, versus the position of the discontinuity in the reaction term a. Three cases have been plotted: w ¼ 0 (parabolic), w ¼ 0:1 and w ¼ 0:5. We depict with symbols the numerical solutions for the front speed computed from an initial condition in the form of the Heaviside function between 0 and 1. In dashed lines we depict the exact solution, which is in very good agreement with the numerical calculations. The bounds (solid lines) exhibiting the lower bound show also a good agreement with the numerical solution but both bounds have the feature of describing the pushed–pulled transition. First, contrarily to the parabolic case w ¼ 0; there is no divergence for the speed, in the limit a ! 0; pffiffiffi because the speed is always bounded above by 1= w, a characteristic property of the time delayed hyperbolic reaction–diffusion equation [14]. This fact may be also checked analytically from Eq. (32) taking the limit a ! 0: Second, it is also worthy of mention how delay time affects the pushed-to-pulled transition. We have ðUÞ just seen the lower bound describes the transition at a ¼ aðLÞ and the exact c ; the upper bound at a ¼ ac ðeÞ solution at a ¼ ac . We show in Fig. 5 the position of these critical values and how they move to the left as parameter w increases. So, we can talk about the enlargement of the pulled region as delay time increases, meanwhile the pushed region shrinks. In Fig. 6 we plot the front speed versus w in the pulled region (Fig. 6(a)) and in the pushed region (Fig. 6(b)). In the pulled region, it is shown that the speed is independent of the discontinuity position a and the numerical solutions (symbols) are in a very good agreement with the theoretical result c ¼ 2=ð1 þ wÞ found from the variational analysis (solid line). In the pushed region, the speed is different for different values of a. In solid lines, we plot our theoretical results obtained for the lower bound (30) and in dashed lines the result obtained from the exact solution (32). A good agreement is also found between numerical and theoretical results.
4.0
VP Exact
3.5
=0
3.0
2.5
(L)
(e)
ac ac 2.0
= 0.1 (L)
1.5
(e)
ac
ac
(L)
ac
(e)
ac
= 0.5 1.0 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
Fig. 5. Dimensionless selected front speed from the exact solution (32) (dashed lines), the lower (cL ) bound (solid lines) (30) and the numerical solutions (symbols) versus the position a of the discontinuity in the reaction term for w ¼ 0; 0:1 and 0:5 with n ¼ b ¼ 1: It is also depicted the position of the critical a for each case.
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
62
2.0
: a = 0.3 : a = 0.5 : VP
1.8 1.6 1.4
6a
1.2 0.0
0.1
0.2
0.3
0.4
0.5
4.50
: a = 0.1 : a = 0.05 : VP : Exact
3.75 3.00 2.25
0.6
6b
1.50 0.00
0.05
0.10
0.15
0.20
0.25
0.30
Fig. 6. Dimensionless selected front speed versus w. In (a) it is plotted the lower bound in the pulled region (solid line) together with the numerical results (symbols). It is observed that in this region the front speed does not depend on a. In (b) it is plotted the lower bound (30) (solid lines) and the exact result (15) for two different values of a. In both cases n ¼ b ¼ 1:
5. Piecewise linear emulation examples In this section, we illustrate how to emulate a nonlinear reaction term f ðuÞ with two linear segments like in (1). We consider, for simplicity the parabolic case. First of all, we have to establish how to choose the parameters n; b and a. As we deal with fronts connecting the states u ¼ 0 and 1 we are interested in nonlinear reaction terms such that f ð0Þ ¼ f ð1Þ ¼ 0: Near these states, the reaction–diffusion equation is approximately linear and f ðuÞ ’ f 0 ð0Þu near u ¼ 0 and f ðuÞ ’ f 0 ð1Þð1 uÞ near u ¼ 1. So, it seems convenient to take n f 0 ð0Þ and b f 0 ð1Þ. Moreover, the piecewise emulation must resemble as closely as possible the nonlinear function in order to keep up its main properties such symmetries, position of extremum, etc. As an example, we consider first the generalized logistic reaction term f ðuÞ ¼ uð1 uk Þ with kX1: It is known that Eq. (2) with this reaction term has only pulled fronts propagating into the unstable state u ¼ 0 with speed c ¼ 2. Let us now to check the results of our piecewise emulation. According to previous considerations, now we have n ¼ 1 and b ¼ k: In this instance, we choose the discontinuity point a as the point where f ðuÞ reaches its maximum value, i.e., a ¼ ð1 þ kÞ1=k . For the logistic case (k ¼ 1), this choice preserves 1=2 the mirror symmetry about u ¼ 12: From the analytic result (15) one gets aðeÞ and for any c ¼ 1 ð1 þ kÞ ðeÞ kX1 one has a4ac and the front is always pulled, as is indeed the case. From (7) may be checked that a4aðLÞ c for any kX1 and from (8) the speed predicted by the lower bound is 2 as for the numerical solution on (2) with (1). However, the upper bound (9) grows with k. This means that the analytic solution and the lower bound of the variational principle of the piecewise emulation yield the correct result. Let ffiffiffiuÞðu þ bÞ with b40 where it is known that the front speed is pffiffiffi us now pffiffiffi to consider the case f ðuÞ ¼ uð1p 1= 2 þ b 2 for 0obo12 (pushed region) or 2 b for bX12 (pulled region). Our piecewise linear emulation takes pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n ¼ b, b ¼ 1 þ b and a ¼ ð1 b þ b2 þ b þ 1Þ=3. The analytic result predicts a pushed front for 0obo0:173: From the variational principle, we obtain a lower bound predicting a pushed front for 0obo0:14 and a pulled one otherwise. The upper bound, however, predicts a pushed front for any b. These results are depicted in Fig. 7. It is necessary to comment that we have illustrated how to emulate nonlinear reaction terms using reaction term composed by two-segment. Our results are reasonably good compared to the exact known results. However, if one desires a higher accuracy in the piecewise-linear emulation one should consider more segments. This was already done in Ref. [15] where piecewise linear reaction terms of
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
63
2.2 2.0
cU
1.8 1.6 1.4
c
1.2 1.0 0.8
: analytic
cL
0.6 0.4 0.2
:
1 +b 2 2
:
2 b
0.0 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
b Fig. 7. Plot of the selected front speed obtained from the piecewise linear emulation for the nonlinear reaction term f ðuÞ ¼ uð1 uÞðu þ bÞ: The results obtained from the variational principle are depicted with dashed lines (upper and lower bounds), the analytic result is depicted by solid line and the speed for the front under emulation is depicted with symbols.
five segments were used to emulate reaction terms with cubic and quintic nonlinearities showing an excellent agreement with the exact results. 6. Conclusions One-component reaction–diffusion model with two-piece linear reaction functions have been analyzed. The development of variational principles, analytical solution and numerical calculations jointly let us show the existence of a pushed-to-pulled transition as the discontinuity in the reaction term a goes from 0 to 1 when fronts propagate into the unstable state. On the other hand, when the front propagates into metastable state, there is not such transition and fronts are always pushed. Specifically, for unstable state u ¼ 0 (n40) the variational principle points out that for lower values of discontinuity position a than a critical one, the front speed is higher than the linear speed and the front is then pushed. Otherwise, the linear speed is selected (because upper and lower bounds from the variational principle are equal) and the front is pulled. In consequence, the variational principle rigorously prove the existence of a transition from pushed-to-pulled front propagating into the unstable state when a goes from 0 to 1. This result is in agreement with the numerical solutions. Furthermore, piecewise linear reactions allow to work out the analytic solution for the front speed and we have seen how this solution holds only for discontinuity position ranging from 0 to a critical value. We have compared the different critical values of a from the variational principle and the analytic solution with numerical computations exhibiting a relatively good agreement. Nevertheless, when the front propagates into metastable state (no0), there is no linear speed selection and neither pushed-to-pulled transition. The analytic solution for the speed is, in this case, in agreement with the numerical solutions for any a, and may predict, contrarily to the variational principle, negative values for the speed. When a time delay is considered in our reaction–diffusion system, all the above arguments are still valid but it plays an important role on the width of the pushed and pulled regions: as delay increases, the pushed region shrinks while the pulled region expands. This fact has an intuitive explanation: time delay slows the front down and the system stays for more time near unstable state. Therefore, the linear part around u ¼ 0 becomes more important as the delay increase, hence the pulled region enlarges.
ARTICLE IN PRESS 64
V. Me´ndez et al. / Physica A 375 (2007) 51–64
We have finally analyzed how to emulate nonlinear reaction terms by two-piece linear reaction terms for fronts propagating into an unstable state. We have considered, in particular, the emulation of the reaction terms f ðuÞ ¼ uð1 uk Þ with kX1 and f ðuÞ ¼ uð1 uÞðu þ bÞ showing in both cases how the two-piece emulation captures the known qualitative behaviour of the front speed when it resembles close enough the nonlinear reaction. Acknowledgment This work has been supported by the MCYT under Grant FIS2006-12296-C02-01 and by Generalitat de Catalunya under Grant 2005 SGR-00087. References [1] H.P. McKean, Adv. Math. 4 (1970) 209; J. Rinzel, J.B. Keller, Biophys. J. 13 (1973) 1313; Y.B. Zeldovich, A.P. Aldushin, S.I. Kudayev, in: J.P. Ostriker (Ed.), Selected Works of Y. B. Zeldovich, vol. I, Princeton University Press, Princeton, 1992, p. 320; A. Ito, T. Ohta, Phys. Rev. A 45 (1992) 8374; V.S. Zykov, Simulations of Wave Processes in Excitable Media, Manchester University Press, New York, 1987. [2] S. Theodorakis, E. Leontidis, Phys. Rev. E 62 (2000) 7802. [3] S. Theodorakis, E. Leontidis, Phys. Rev. E 65 (2002) 026122. [4] M. Wang, M. Kot, M.G. Neubert, J. Math. Biol. 44 (2002) 150. [5] V. Me´ndez, A. Compte, Physica A 260 (1998) 90; K.K. Manne, A.J. Hurd, V.M. Kenkre, Phys. Rev. E 61 (2000) 4177; G. Abramson, A.R. Bishop, V.M. Kenkre, Phys. Rev. E 64 (2001) 066615. [6] G. Fath, Physica D 116 (1998) 176; A. Lahiri, P. Majumdar, M.S. Roy, Phys. Rev. E 65 (2002) 026106; A. Tonnelier, Phys. Rev. E 67 (2003) 036105. [7] A.M. Pertsov, M. Wellner, J. Jalife, Phys. Rev. Lett. 78 (1997) 2656; Y.B. Chernyak, J.M. Starobin, R.J. Cohen, Phys. Rev. Lett. 80 (1998) 5675; Y.B. Chernyak, J.M. Starobin, R.J. Cohen, Phys. Rev. E 58 (1998) R4108. [8] E.P. Zemskov, K. Kassner, S.C. Mu¨ller, Eur. Phys. J. B 34 (2003) 285; E.P. Zemskov, K. Kassner, S.C. Mu¨ller, Phys. Rev. E 70 (2004) 056208; R. Bakanas, Phys. Rev. E 69 (2004) 016103; R. Bakanas, Phys. Rev. E 71 (2005) 026201. [9] W. van Saarloos, Phys. Rep. 386 (2003) 29. [10] P. Clavin, A. Li na´n, in: M.G. Velarde (Ed.), Nonequilibrium Cooperative Phenomena in Physics and Related Fields, NATO Advanced Study Institute Series B: Physics, vol. 116, Plenum, New York, 1984, p. 291. [11] R.D. Benguria, M.C. Depassier, Phys. Rev. Lett. 77 (1996) 1171. [12] R.D. Benguria, M.C. Depassier, Phys. Rev. E 66 (2002) 026607. [13] D.G. Aronson, H.F. Weinberger, Adv. Math. 30 (1978) 33. [14] J. Fort, V. Me´ndez, Rep. Prog. Phys. 65 (2002) 895. [15] S. Theodorakis, E. Svoukis, Phys. Rev. E 68 (2003) 027201.
Physica A 375 (2007) 51–64 www.elsevier.com/locate/physa
Transition from pushed-to-pulled fronts in piecewise linear reaction–diffusion systems V. Me´ndez, V. Ortega-Cejas, E.P. Zemskov, J. Casas-Va´zquez Grup de Fı´sica Estadı´stica, Departament de Fı´sica, Facultat de Cie`ncies, Edifici Cc,Universitat Auto`noma de Barcelona, E-08193 Bellaterra (Cerdanyola), Spain Received 14 December 2005; received in revised form 3 July 2006 Available online 4 October 2006
Abstract The front dynamics in reaction–diffusion equations with a piecewise linear reaction term is studied. A transition from pushed-to-pulled fronts when they propagate into the unstable state is found using a variational principle. This transition occurs for a critical value of the discontinuity position in the reaction function. In particular, we study how the transition depends on the properties of the reaction term and on the delay time. Our results are in good agreement with the numerical solutions of the model. r 2006 Elsevier B.V. All rights reserved. Keywords: Pulled fronts; Pushed fronts; Reaction–diffusion
1. Introduction Piecewise linear models for reaction kinetics are used to get exact analytic solutions for the speed and wave profiles in reaction–diffusion equations [1]. These are just an approach of a more realistic but analytically intractable nonlinear reaction–diffusion equations. The piecewise linear emulation consists in considering the reaction term as two or more linear segments reaching the value 0 at each equilibrium state. In the last decade these linear models have been widely employed, for instance, to find a speed selection mechanism for both singlecomponent [2] and two-components [3] parabolic reaction–diffusion (PRD) equations, to study Allee effect in the spread of invading organisms [4] and the effect of transport memory in reaction–diffusion fronts [5] and to analyze the wave propagation in discrete [6] and excitable [7] media and under forcing in bistable systems [8]. In Refs. [1], the piecewise linear emulation is applied to bistable kinetics by considering f H u þ yðu aÞ, a ¼ const. This allows to find analytic solutions for traveling waves like kinks, wave trains, solitary waves and two-dimensional spiral waves described by activator-inhibitor systems. We deal, however, with a more fundamental problem which appears only in fronts propagating into unstable states: the pushed–pulled transition. Although this is an already known problem for nonlinear reaction term with cubic nonlinearities [9] we are interested here to study this problem when one considers the piecewise linear emulation of these nonlinear Corresponding author.
E-mail address: [email protected] (V. Me´ndez). 0378-4371/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2006.09.007
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
52
reaction terms. In Refs. [2,3,5], the authors deal also with fronts propagating into unstable states but no attention is paid into the pushed–pulled transition because these works only deal with analytic solutions. The main advantage of dealing with piecewise linear reaction terms is to get analytical solutions for the front profiles and their speeds but, as we show in this work, these solutions only hold when the values of the parameters lie in the pushed region. However, from the variational principle and numerical solutions developed in this paper, we are able to find the front speed in both regions (pushed and pulled). We consider a function [10] ( nu; 0puoa; f ðuÞ ¼ (1) bð1 uÞ; aoup1; which is an approximation of a general reaction term f ðuÞ with two equilibrium states: u ¼ 0 and u ¼ 1 for the Fisher–Kolmogorov reaction–diffusion equation qt u ¼ qxx u þ f ðuÞ.
(2)
In this work, we deal with a front propagating into the unstable (n40) or metastable (no0) state u ¼ 0 connected to the stable state u ¼ 1 by means of Eq. (2), and by its generalization to the case with a time delay, with f ðuÞ given by Eq. (1). We calculate the front profiles and the front speed equation, once Eq. (2) has been transformed into an ordinary differential equation by using the traveling wave coordinate z ¼ x ct, being c the asymptotic front speed. We also apply a variational principle to the problems (1) and (2) to calculate lower and upper bounds for the front speed and perform the numerical integrations employing an explicit finite difference method to obtain the front profiles and c. Both techniques allow us to find a pushed-to-pulled transition when the front propagates into the unstable state. We say the front is pulled when it moves into an unstable state with constant speed (after a transient) given by the linear analysis around the unstable state. We refer to this speed as linear speed. On the other hand, we say the front is pushed when its speed after the transient is higher than the linear one [9]. We end this work by illustrating how to select the parameters n, b, and a in (1) to emulate some specific examples of nonlinear reaction terms. We show that the piecewise model captures the essential qualitative dynamical properties of the system under emulation whenever it preserves the relevant properties of the real nonlinear function. 2. Variational principle Variational principles have been used to find bounds for the speed of reaction–diffusion fronts. They are very useful to obtain a bound for the speed when the linear marginal stability fails and the speed selected is strongly affected for the nonlinear terms. In our model, the monotonically decreasing fronts u ¼ uðz ¼ x ctÞ propagating with constant speed c join u ¼ 1 with u ¼ 0: The front satisfies u00 þ cu0 þ f ðuÞ ¼ 0; where the prime symbol stands for derivative respect to z and fills the boundary conditions uðz ! 1Þ ¼ 0 and uðz ! 1Þ ¼ 1. From a variational principle, the asymptotic front speed is given by [11] " R1 # 2 0 fg du c ¼ max 2 R 1 , (3) g 2 =hÞ du ðg 0 where the maximum is taken over all positive decreasing auxiliary functions gðuÞ in ð0; 1Þ (h ¼ dg=du40) for which the integrals involved in Eq. (3) exist. It is often impossible to know the form of the function g ¼ gðuÞ which maximizes the relationship in brackets in Eq. (3) and one can find only lower bounds R1 fg du 2 cL ¼ 2 R 1 0 . (4) 2 0 ðg =hÞ du ItR is easy toR get an upper R 1 bound for the speed independent of g [12]. From R 1 the Schwarz R 1 inequality 1 1 ð 0 ug duÞ2 p 0 ðg2 =hÞ du 0 u2 h du and assuming gð1Þ ¼ 0 in the integration by parts 0 u2 h du ¼ 2 0 ug du, then " R1 # " R1 # fg du fg du f ðuÞ 2 p4 max 2 R 01 p4 sup c ¼ max 2 R 1 0 g g u u2ð0;1Þ ðg2 =hÞ du ug du 0
0
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
and the upper bound is sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi f ðuÞ cU ¼ 2 sup u u2ð0;1Þ
53
(5)
in agreement with the result of Aronson and Weinberger [13]. In this work, we are going to deal with the family of functions gðuÞ ¼ ua1 1 which are decreasing functions if ao1 and fulfil gð1Þ ¼ 0: To get the nearest lower bound to true speed value we will take the maximum value over all the values for the parameter a such that the integrals in (4) are convergent. So, calculating the lower bound from Eq. (4) for f ðuÞ in form (1) one gets c2L ¼ max1oao1 Gða; a; n; bÞ where aþ1 2ð1 þ aÞð3 aÞ a a2 1 b 1 a ðn þ bÞ a þ Gða; a; n; bÞ ¼ . (6) þb a 1a 2 a 1þa aþ1 2
2.1. Propagation into unstable state (n40) ðLÞ For n40, there exists a critical value aðLÞ c , such that when aoac the maximum of G, i.e., of cL ; is attained at ðLÞ a ¼ a with 1oa o0. However, when aXac the maximum is reached always at a ¼ 1 and just at this point one has Gð1; a; n; bÞ ¼ 4n: To compute the critical position of the discontinuity from the lower bound aðLÞ c ; we take dG=da ¼ 0 at a ¼ 1 and find ðLÞ 2 ðLÞ ðLÞ 2 ðLÞ 2aðLÞ c ðn þ bÞ½ðac Þ 2 ln ac 4b½1 þ ðac Þ ac ðn 6bÞ ¼ 0.
To sum up, we find from the variational principle that the lower bound for the speed is given by ( pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Gða ; a; n; bÞ; aoaðLÞ c ðn; bÞ; pffiffiffi cL ¼ 2 n; a4aðLÞ c ðn; bÞ;
(7)
(8)
where Gða ; a; n; bÞ44n; a is such that dG=da ¼ 0 at a ¼ a and aðLÞ c ðn; bÞ is solution of Eq. (7). Let us now to evaluate the upper bound for the front speed. From Eqs. (1) and (5) we get 8 pffiffiffi b > > ; a4acðUÞ ¼ > < 2 n; nþb r ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi (9) cU ¼ > 1a > > b ; otherwise; 2 : a ðUÞ with aðLÞ c oac . This result shows, together with (8), that the speed of the front is linearly selected for a sufficiently large values of a and it is nonlinearly selected for sufficiently low values of a. Therefore, there exists a transition from pushed-to-pulled front when a goes from 0 to 1. The transition is found at a critical value of the position of the discontinuity. From the variational principle we have found two possible values: aðLÞ c from the lower bound (Eq. (8)) and aðUÞ from the upper bound (Eq. (9)). c
2.2. Propagation into metastable state (no0) For no0, the variational principle shows that there is no linearly selected speed and neither a transition pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi from pushed-to-pulled fronts. The speed is always nonlinearly selected and it is given by c ’ max1oao1 G : For large enough values of a the maximum is attained at a ¼ 1: In this region there exists a maximum value of ^ such that for a ¼ a^ one has Gð1; a; ^ n; bÞ ¼ 0 and as a result c ¼ 0: From Eq. (6) it is obtained that a^ a, namely a, satisfies the equation 2½2b aðn þ bÞa ln a þ a2 ðn þ bÞ þ bð3 4aÞ ¼ 0: For a4a^ one finds Go0 and the variational principle is not able to supply lower bounds for the front speed, because the square root of the negative function G appears. In consequence, our variational principle only can predict positives values for c.
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
54
3. Analytic solution In this section we find the analytic solution of the equation u00 þ cu0 þ f ðuÞ ¼ 0; where the reaction term is given by Eq. (1). As the differential equation for the monotonic asymptotic fronts is linear, the solution involves two exponential functions for each piece. In consequence, four integration constants have to be found from boundary and matching conditions. In the first piece f ðuÞ ¼ nu; where 0puoa and we take u ¼ u1 ðzÞ with z0 pzo1: In the second piece f ðuÞ ¼ bð1 uÞ with aoup1 and we have u ¼ u2 ðzÞ with 1ozpz0 : The boundary conditions for each piece are then u1 ðz ! 1Þ ¼ 0 and u1 ðz0 Þ ¼ a for the first piece and u2 ðz ! 1Þ ¼ 1 and u2 ðz0 Þ ¼ a for the second piece. Therefore, the solutions for the front profiles read pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi c c2 4n lþ z l z u1 ðzÞ ¼ A1 e þ A2 e with l ¼ and c42 n, 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c c2 þ 4b þ . ð10Þ u2 ðzÞ ¼ 1 þ B1 em z þ B2 em z with m ¼ 2 Consider first the front propagation into the unstable state, i.e., the case with n40: As l o0 the boundary condition u1 ðz ! 1Þ ¼ 0 does not allow us to find any integration constant. For the second piece the boundary condition u2 ðz ! 1Þ ¼ 1 requires B2 ¼ 0 because mþ 40 and m o0: One can find two of the three integration constants from two boundary conditions u1;2 ðz0 Þ ¼ a and one matching condition for the derivative, u01 ðz0 Þ ¼ u02 ðz0 Þ; one can obtain the third one. All integration constants may be determined but there is no any additional condition to determine the equation for the front speed; c is then another free parameter and no unique front exists. Alternatively, when the front propagates into the metastable state (no0), then l o0 but lþ 40 and A1 ¼ 0 according to the boundary condition u1 ðz ! 1Þ ¼ 0: Thus, in this case, the front speed may be determined from the matching condition for the continuity of the derivative at z0 and an unique front exists. In order to determine the front speed for the propagation into the unstable state in the context of the piecewise linear model, we assume initially the boundary conditions u1 ðz ¼ LÞ ¼ 0 and u2 ðz ¼ LÞ ¼ 1; instead of u1 ðz ! 1Þ ¼ 0 and u2 ðz ! 1Þ ¼ 1; and finally we will take L ! 1 [2]. The expressions for the integration constants read
þ
el L el L ; A ¼ a , 2 þ þ el ðz0 LÞ el ðz0 LÞ el ðz0 LÞ el ðz0 LÞ þ em L em L B1 ¼ ð1 aÞ m ðz þLÞ ; B ¼ ð1 aÞ . ð11Þ 2 þ e 0 em ðz0 þLÞ em ðz0 þLÞ emþ ðz0 þLÞ From the matching condition for the continuity of the derivative in z0 ; one finds, after some algebra an equation for the family of possible front speeds pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 1 c2 þ 4bðz0 þ LÞ . c ¼ a c2 4n coth (12) c 4nðz0 LÞ þ ð1 aÞ c2 þ 4b coth 2 2 A1 ¼ a
Taking the limit at L ! 1 one gets the selected speed which is independent on z0 , according to the translational invariance of the model equation [2]. Thus, the speed equation is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ a c2 4n þ ð1 aÞ c2 þ 4b. (13) Taking L ! 1; the integration constants become A1 ¼ B2 ¼ 0; A2 ¼ a expðl z0 Þ and B1 ¼ ð1 aÞ expðmþ z0 Þ and the exact solution for the front profiles are finally ( ael ðz0 zÞ ; z 2 ½z0 ; 1Þ; uðzÞ ¼ (14) mþ ðz0 zÞ 1 ð1 aÞe ; z 2 ð1; z0 : pffiffiffi It is worthy to stress that if we assume co2 n in Eq. (10) the solution to u1 ðzÞ involves sines and cosines while u2 ðzÞ remains the same. By imposing the matching conditions for u1 ðzÞ and u2 ðzÞ at z ¼ z0 and z ¼ L one gets the corresponding integration constants. By requiring the matching condition for the derivatives at z ¼ z0 one obtains the speed equation. However, when one takes the limit L ! 1 it is not possible to get an equation equivalent to (13) to get the selected speed.
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
55
3.1. Propagation into unstable state (n40) Eq. (13) is an implicit relationship for the selected front speed. It is important to note that this equation has pffiffiffi real solutions only if c42 n, as one can easily check from the first term of the right-hand side of Eq. (13). This means that the selected front speed from (13), once it is explicitly expressed in terms of n, b and a, has to pfound ffiffiffi be higher than the linear speed 2 n. Solving Eq. (13) for c one finds the selected front speed a2 ðn þ bÞ þ bð1 2aÞ c ¼ pffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . a a2 ðn þ bÞ ðn þ 2bÞa þ b
(15)
pffiffiffi Introducing this result in the restriction c42 n, this inequality turns into rffiffiffiffiffiffiffiffiffiffiffi n ðeÞ . aoac ¼ 1 nþb
(16)
This means that the exact solution (15) is restricted to values of a such that a 2 ð0; aðeÞ For c : p ffiffiffi this range of values of a the front is pushed because the speed given in (15) is higher than the linear speed 2 n as we have required. From the variational principle applied to (1) and (2) we have shown in Eqs. (8) and (9) that a monotonic front and its speed exist for any value of a in (0,1) exhibiting a pushed–pulled transition. This result has been also supported by numerical solutions. The exact solution is not able to describe this transition because (15) does not hold for a 2 ðaðeÞ c ; 1Þ and nothing can be known only from the exact solution. To illustrate this fact, we plot in Fig. 1 the front speed obtained from the exact solution (15) and the lower and upper bounds found from the variational principles (8) and (9), respectively, versus the position a of the discontinuity in the reaction term. We depict with symbols the numerical solution for the front speed computed directly from Eqs. (1) and (2) with an initial condition in the form of the Heaviside function between 0 and 1. In dashed lines we depict the exact solution, which is in very good agreement with the numerical calculations but only for ðeÞ a 2 ð0; aðeÞ c , because this curve stops at the point a ¼ ac and cannot go beyond this point. The bounds are
6
ν = 0.5
c
5
cU
4 3
(U)
cL
2 (a)
ac
6
(L)
ac
(e)
ac
ν=1
c
5
cU
4
(U)
ac
3 2
cL 0.0
(b)
0.2
(L)
ac
(e)
ac
0.4
0.6
0.8
1.0
a
Fig. 1. Dimensionless selected front speed from the exact solution (15) (dashed lines), the lower (cL ) and upper (cU ) bounds from the variational principle (VP) and the numerical solutions (symbols) versus the position a of the discontinuity in the reaction term. In (a) n ¼ 12 and b ¼ 1. In (b) n ¼ b ¼ 1.
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
56
depicted in solid lines. The lower bound shows also a good agreement with the numerical solution. Both bounds join within ½acðUÞ ; 1 describing the pushed–pulled transition because from Eqs. (8) and (9) pffiffiffi the range ðUÞ c ¼ cL ¼ cU ¼ 2 n for a 2 ½ac ; 1. While the lower bound describes the transition at a ¼ aðLÞ c ; the upper ðLÞ bound describes the transition at a ¼ aðUÞ 4a . In Table 1, we compare both values, for the critical position c c of the discontinuity, with the obtained from the analytic solution (Eq. (15)) and from the numerical solution. ðeÞ ðnumÞ It is observed a good agreement between aðLÞ . In Fig. 2, we plot the c and ac with the numerical solution ac critical value of a in the pushed-to-pulled transition obtained from the lower bound (aðLÞ c from Eq. (7)) and from the exact solution (aðeÞ from Eq. (15)) together with numerical computations. It is shown that both exact c and variational results yield very similar results and in agreement with numerical solutions. It is also interesting to note that the critical position moves towards 0 as the slope of the first segment n increases keeping fixed the slope of the second segment b. Analogously, for fixed n when b increases, it moves towards 1. In Fig. 3, we represent two front profiles using the exact solution (14) (solid lines) together with the profiles obtained from numerical integrations (symbols) exhibiting an excellent agreement.
Table 1 Table of the results for the critical position of the discontinuity for different values of the slope of the first segment n for the case b ¼ 1 n
acðLÞ
aðUÞ c
aðeÞ c
aðnumÞ c
0.1 0.3 0.5 0.7 0.9 1
0.67 0.51 0.42 0.35 0.30 0.28
0.91 0.77 0.67 0.59 0.53 0.50
0.70 0.52 0.42 0.36 0.31 0.29
0.75 0.55 0.45 0.35 0.35 0.25
ðUÞ is obtained from the upper bound (Eq. (9)), aðeÞ aðLÞ c is the critical value obtained from the lower bound (Eq. (7)), ac c is obtained from the ðnumÞ analytic solution (Eq. (15)) and ac is computed numerically.
1.0
0.8 β=2
ac
0.6 β=1 0.4
β = 0.5
0.2
0.0 0.0
0.2
0.4
0.6
0.8
1.0
ν Fig. 2. Dimensionless critical discontinuity position versus the slope n of the first piece for three different values of the slope b of the ðLÞ second piece. The dashed curves correspond to aðeÞ c ; Eq. (15) and the solid curves correspond to ac ; Eq. (7). The symbols correspond to numerical solutions (circles for b ¼ 0:5; squares for b ¼ 1 and triangles for b ¼ 2).
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
57
1.0 ν = 0.5, β = 1
u(z)
0.8 0.6 0.4 0.2 0.0
(a) 1.0 ν=β= 1
u(z)
0.8 0.6 0.4 0.2 0.0 -14
-12
-10
(b)
-8
-6
-4
-2
0
2
4
z
Fig. 3. Front profiles versus z ¼ x ct z0 : With solid lines we depict the results obtained from the exact solution (Eq. (14)) and with symbols we depict the numerical solutions. For (a) we take n ¼ 0:5, a ¼ 0:1 and for these values the speed is c ¼ 2:91. In this case we take a front profile at t ¼ 120 with z0 ¼ 50:35: For (b) we take n ¼ 1, a ¼ 0:1 and for these values the speed is c ¼ 3:02 and we take a front profile at t ¼ 120 with z0 ¼ 50:75:
3.2. Propagation into metastable state (no0) When the piecewise linear model describes a front propagating into a metastable state, no critical value for a is found, neither from variational principles nor from numerical solutions and the analytic solution for the speed (15) holds for any a: In this situation, the front never selects the linear speed. These results are in agreement with that known for fronts propagating into metastable states described by reaction terms with cubic nonlinearities [9]. Therefore, there is no transition from pushed-to-pulled front. The front speed from the variational principle (lower and upper bound found from Eqs. (8) and (9)) and from the analytic solution (found from Eq. (15)) are plotted in Fig. 4 together with numerical calculations. In this case, the front speed computed from the exact solution is in a better agreement with numerical integrations than the bounds obtained from variational principle. The variational principle for no0 can only predict positive bounds for the front speed, as we have shown in Section 2.2. Therefore, it is unable to account, contrarily to the analytic result (15), for the change of sign in c, characteristic of some bistable systems. 4. Effect of delay In this section, we extend our analysis to study how depends ac on a time delay w. For simplicity we take b ¼ 1 in Eq. (1). To do this, we consider the hyperbolic reaction–diffusion equation [14] wqtt u þ qt u ¼ qxx u þ f ðuÞ þ wqt f ðuÞ,
(17)
where the delay time is included in w. In fact, w is the quotient between the delay time and the characteristic reaction time. As above, let us to begin with the variational principle. The front satisfies ð1 wc2 Þu00 þ c½1 wf 0 ðuÞu0 þ f ðuÞ ¼ 0 and the boundary conditions uðz ! 1Þ ¼ 0 and uðz ! 1Þ ¼ 1: Hence, from (17) one finds that the monotonic fronts are solutions of ð1 wc2 Þp
dp c½1 wf 0 ðuÞp þ f ðuÞ ¼ 0, du
(18)
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
58
6 ν = -1
5 4 c
cU
analytic
3 2
cL
1 0 (a)
0.0
0.1
0.2
0.3
0.4
6 ν = -0.5
5
c
4
cU
analytic
3 2
cL
1 0 0.0
0.1
0.2
0.3
0.4
a
(b)
Fig. 4. Selected front speed versus a. In solid lines, we plot the speed obtained from the exact solution (15) and with dashed curves we plot lower and upper bounds obtained from the variational principle. Symbols corresponds to numerical solutions: (a) corresponds to the case n ¼ 1 and b ¼ 1; and (b) corresponds to n ¼ 0:5 and b ¼ 1:
with pð0Þ ¼ pð1Þ ¼ 0 in (0,1). Multiplying Eq. (18) by gðuÞ and integrating between u ¼ 0 and u ¼ 1 we obtain, Z 1 Z 1 1 gf du ¼ fðpÞ du with fðpÞ cgp½1 wf 0 ðuÞ ð1 wc2 Þhp2 (19) 2 0 0 after integrating by parts. As c; p; g; and h are positive functions for any fixed u in (0,1), fðpÞ has a maximum at pmax ¼
cg½1 wf 0 ðuÞ 40 hð1 wc2 Þ
(20)
if 1 co pffiffiffi w
and
1 wf 0 ðuÞ40
(21)
or analogously if woM 1
with M ¼ max f 0 ðuÞ. u2½0;1
In consequence, fðpÞp
c2 g2 ½1 wf 0 ðuÞ2 2hð1 wc2 Þ
for any fixed u in (0,1). Finally, it follows R1 gf du c2 X2 R 1 g2 0 1 wc2 ½1 wf 0 ðuÞ2 du 0 h
(22)
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
and the lower bound is given by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R1 gf du A with A ¼ 2 R 1 g2 0 . cL ¼ 1 þ wA ½1 wf 0 ðuÞ2 du
59
(23)
0 h
To prove that this is a variational principle we must to ensure that there exists a function g for which the equality in (22) holds. From Eq. (20) we see that the case of equality is attained when g^ fulfills the following differential equation: pð1 wc2 Þ g^ . 0 ¼ c½1 wf 0 ðuÞ g^ ^ obtained by integrating this equation is The maximizing g; Z u c 1 wf 0 ðuÞ ^ ¼ exp du gðuÞ 1 wc2 u0 p
(24)
with 0ou0 o1: We must ensure now that the integrals in Eq. (22) are convergent. The problems may appear when one computes the integrals near u ¼ 0. In this limit plu with pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c½1 wf 0 ð0Þ þ c2 ½1 þ wf 0 ð0Þ2 4f 0 ð0Þ l¼ 2ð1 wc2 Þ pffiffiffiffiffiffiffiffiffi ffi whenever c42 f 0 ð0Þ=½1 þ wf 0 ð0Þ to ensure that l 2 R. Then, from (24) we get k ^ gðuÞu
with k ¼
c½1 wf 0 ð0Þ lð1 wc2 Þ
(25)
near u ¼ 0 and gf ðg2 =hÞ½1 wf 0 ðuÞ2 u1k . Finally, the integrals in Eq. (22) exist if ko2 which is fulfilled pffiffiffiffiffiffiffiffiffiffi provided that c42 f 0 ð0Þ=½1 þ wf 0 ð0Þ. In consequence one can conclude that the asymptotic front speed is given by ( ) R1 gf du c2 ¼ max 2 R 1 g2 0 . (26) g 1 wc2 ½1 wf 0 ðuÞ2 du 0 h
It is not difficult to obtain an upper bound for the speed independent of g. From the inequality Z 1 2 Z 1 2 g g 0 0 2 2 ½1 wf ðuÞ duX inf ½1 wf ðuÞ du u2ð0;1Þ h 0 0 h R1 R1 R1 R1 and Schwarz inequality ð 0 ug duÞ2 p 0 ðg2 =hÞ du 0 u2 h du and assuming gð1Þ ¼ 0 one obtains 0 ug dup R 1 the 2 0 ðg2 =hÞ du, then ( ) R1 c2 0 gf du ¼ max 2 R 1 g2 0 2 g 1 wc2 0 h ½1 wf ðuÞ du " R1 # 2 0 fg du p max R 1 2 inf u2ð0;1Þ ½1 wf 0 ðuÞ2 g 0 ðg =hÞ du f ðuÞ supu2ð0;1Þ u p4 inf u2ð0;1Þ ½1 wf 0 ðuÞ2 and the upper bound is f ðuÞ sffiffiffiffiffiffiffiffiffiffiffiffiffiffi supu2ð0;1Þ B u with B ¼ 4 . cU ¼ 1 þ wB inf u2ð0;1Þ ½1 wf 0 ðuÞ2
(27)
ARTICLE IN PRESS 60
V. Me´ndez et al. / Physica A 375 (2007) 51–64
In this work, we are going to deal with the family of functions gðuÞ ¼ ua1 1 which are decreasing functions if ao1 and fulfil gð1Þ ¼ 0: To get a lower bound nearest to true speed value we will take the maximum value over all the values for the parameter a such that the integrals in (23) are convergent. So, calculating the lower pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bound from Eq. (23) one gets cL ¼ max1oao1 A=ð1 þ wAÞ where aþ1 a a2 1 1 aa þa ðn þ 1Þ þ að1 þ aÞ 2 aþ1 2 a . (28) Aða; w; a; nÞ ¼ ð1 aÞ aþ1 3a aþ1 3a a a 1 a 1 a 2 2 2 þ þ 1þa ð1 wnÞ þ ð1 þ wÞ aþ1 3a aþ1 3a ðLÞ When n40, there exists a critical value aðLÞ c , such that when aoac the maximum of A, i.e., of cL ; is attained at a ¼ a with 1oa o0 and the front speed is higher than that linearly selected (pushed front). However, when aXacðLÞ the maximum is reached always at a ¼ 1 and just at this point Að1; w; a; nÞ ¼ 4n=ð1 wnÞ and the pffiffiffi front speed is that linearly selected (pulled front) c ¼ 2 nð1 þ wnÞ1 . Therefore, this shows that the lower bound for the front speed describes a pushed-to-pulled transition. However, for no0 the lower bound does not describe transition as already occurs for w ¼ 0. In the following we will take n ¼ 1 without loosing of generality. The critical position of the discontinuity from the lower bound aðLÞ c is such that dA=da ¼ 0 at a ¼ 1, so that it follows the transcendent equation 2 ðLÞ ðLÞ ðLÞ ðLÞ 2 ðLÞ 5 ðLÞ 3 ðLÞ 2 ðLÞ ½4 þ 4ðaðLÞ c Þ ð1 ac Þ 5ac þ 8ac ln ac ð1 þ w Þ þ w½4ðac Þ 8ðac Þ 8ðac Þ þ 22ac 8 ¼ 0. (29)
Note that for w ¼ 0 in (29), one recovers (7) taking n ¼ 1: To sum up, we have found that the lower bound for the speed is given by 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Aða ; w; a; 1Þ > ðLÞ > > < 1 þ wAða ; w; a; 1Þ; aoac ðwÞ; (30) cL ¼ 2 > ðLÞ > > ; a4a ðwÞ; :1 þ w c where Aða ; w; a; 1Þ44=ð1 wÞ; a is such that dA=da ¼ 0 at a ¼ a and aðLÞ c ðwÞ is solution to Eq. (29). Let us now to evaluate the upper bound for the front speed. From Eq. (27) we get 8 2 1 > ðUÞ > > rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; apac ¼ ; > 2 a > < ð1 wÞ2 þ 4w 1 a (31) cU ¼ > > 2 1 > > ðUÞ > aXac ¼ ; : 1 þ w; 2 This result shows, together with (30), that the speed of the front is exactly the linearly selected for a sufficiently large values of a and it is nonlinearly selected for sufficiently low values of a. Therefore, there exists a transition from pushed-to-pulled front when a goes from 0 to 1. From the variational principle we have found ðUÞ two possible values: aðLÞ from the upper bound (Eq. (31)). c from the lower bound (Eq. (30)) and ac The HRD equation (17) with (1) can be analytically solved as we did in Section 3 with the parabolic case w ¼ 0. The solution is formally the same as for the parabolic equation under an appropriate redefinition of the parameters. So, by imposing matching conditions for u and its derivative at the discontinuity one gets finally the following selected speed equation: " rffiffiffiffiffiffiffiffiffiffiffiffiffi# 1 2a þ 2a2 1 1 þ w2 ðeÞ c ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi if aoac ¼ 1 . (32) 1þw 2 1 2a að1 aÞð1 w2 Þ þ wð1 2aÞ
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
61
It is worthy to stress, as occurs for the parabolic case, that the exact speed given in Eq. (32) only holds for aoaðeÞ c and the front is always pushed. However, from the variational principle it is shown in Eqs. (30) and (31) that the speed exists for any value of a in (0,1) exhibiting a transition from pushed-to-pulled front. As already occurs in (15) and (16) the exact solution (32) does not inform about what happens beyond aðeÞ c . Anyway, it seems clear the inclusion of the time delay does not affect the existence of the transition at all. In Fig. 5, we plot the front speed obtained from the exact solution (32) together with the lower and upper bounds found from the variational principle (30) and (31), respectively, versus the position of the discontinuity in the reaction term a. Three cases have been plotted: w ¼ 0 (parabolic), w ¼ 0:1 and w ¼ 0:5. We depict with symbols the numerical solutions for the front speed computed from an initial condition in the form of the Heaviside function between 0 and 1. In dashed lines we depict the exact solution, which is in very good agreement with the numerical calculations. The bounds (solid lines) exhibiting the lower bound show also a good agreement with the numerical solution but both bounds have the feature of describing the pushed–pulled transition. First, contrarily to the parabolic case w ¼ 0; there is no divergence for the speed, in the limit a ! 0; pffiffiffi because the speed is always bounded above by 1= w, a characteristic property of the time delayed hyperbolic reaction–diffusion equation [14]. This fact may be also checked analytically from Eq. (32) taking the limit a ! 0: Second, it is also worthy of mention how delay time affects the pushed-to-pulled transition. We have ðUÞ just seen the lower bound describes the transition at a ¼ aðLÞ and the exact c ; the upper bound at a ¼ ac ðeÞ solution at a ¼ ac . We show in Fig. 5 the position of these critical values and how they move to the left as parameter w increases. So, we can talk about the enlargement of the pulled region as delay time increases, meanwhile the pushed region shrinks. In Fig. 6 we plot the front speed versus w in the pulled region (Fig. 6(a)) and in the pushed region (Fig. 6(b)). In the pulled region, it is shown that the speed is independent of the discontinuity position a and the numerical solutions (symbols) are in a very good agreement with the theoretical result c ¼ 2=ð1 þ wÞ found from the variational analysis (solid line). In the pushed region, the speed is different for different values of a. In solid lines, we plot our theoretical results obtained for the lower bound (30) and in dashed lines the result obtained from the exact solution (32). A good agreement is also found between numerical and theoretical results.
4.0
VP Exact
3.5
=0
3.0
2.5
(L)
(e)
ac ac 2.0
= 0.1 (L)
1.5
(e)
ac
ac
(L)
ac
(e)
ac
= 0.5 1.0 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
Fig. 5. Dimensionless selected front speed from the exact solution (32) (dashed lines), the lower (cL ) bound (solid lines) (30) and the numerical solutions (symbols) versus the position a of the discontinuity in the reaction term for w ¼ 0; 0:1 and 0:5 with n ¼ b ¼ 1: It is also depicted the position of the critical a for each case.
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
62
2.0
: a = 0.3 : a = 0.5 : VP
1.8 1.6 1.4
6a
1.2 0.0
0.1
0.2
0.3
0.4
0.5
4.50
: a = 0.1 : a = 0.05 : VP : Exact
3.75 3.00 2.25
0.6
6b
1.50 0.00
0.05
0.10
0.15
0.20
0.25
0.30
Fig. 6. Dimensionless selected front speed versus w. In (a) it is plotted the lower bound in the pulled region (solid line) together with the numerical results (symbols). It is observed that in this region the front speed does not depend on a. In (b) it is plotted the lower bound (30) (solid lines) and the exact result (15) for two different values of a. In both cases n ¼ b ¼ 1:
5. Piecewise linear emulation examples In this section, we illustrate how to emulate a nonlinear reaction term f ðuÞ with two linear segments like in (1). We consider, for simplicity the parabolic case. First of all, we have to establish how to choose the parameters n; b and a. As we deal with fronts connecting the states u ¼ 0 and 1 we are interested in nonlinear reaction terms such that f ð0Þ ¼ f ð1Þ ¼ 0: Near these states, the reaction–diffusion equation is approximately linear and f ðuÞ ’ f 0 ð0Þu near u ¼ 0 and f ðuÞ ’ f 0 ð1Þð1 uÞ near u ¼ 1. So, it seems convenient to take n f 0 ð0Þ and b f 0 ð1Þ. Moreover, the piecewise emulation must resemble as closely as possible the nonlinear function in order to keep up its main properties such symmetries, position of extremum, etc. As an example, we consider first the generalized logistic reaction term f ðuÞ ¼ uð1 uk Þ with kX1: It is known that Eq. (2) with this reaction term has only pulled fronts propagating into the unstable state u ¼ 0 with speed c ¼ 2. Let us now to check the results of our piecewise emulation. According to previous considerations, now we have n ¼ 1 and b ¼ k: In this instance, we choose the discontinuity point a as the point where f ðuÞ reaches its maximum value, i.e., a ¼ ð1 þ kÞ1=k . For the logistic case (k ¼ 1), this choice preserves 1=2 the mirror symmetry about u ¼ 12: From the analytic result (15) one gets aðeÞ and for any c ¼ 1 ð1 þ kÞ ðeÞ kX1 one has a4ac and the front is always pulled, as is indeed the case. From (7) may be checked that a4aðLÞ c for any kX1 and from (8) the speed predicted by the lower bound is 2 as for the numerical solution on (2) with (1). However, the upper bound (9) grows with k. This means that the analytic solution and the lower bound of the variational principle of the piecewise emulation yield the correct result. Let ffiffiffiuÞðu þ bÞ with b40 where it is known that the front speed is pffiffiffi us now pffiffiffi to consider the case f ðuÞ ¼ uð1p 1= 2 þ b 2 for 0obo12 (pushed region) or 2 b for bX12 (pulled region). Our piecewise linear emulation takes pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n ¼ b, b ¼ 1 þ b and a ¼ ð1 b þ b2 þ b þ 1Þ=3. The analytic result predicts a pushed front for 0obo0:173: From the variational principle, we obtain a lower bound predicting a pushed front for 0obo0:14 and a pulled one otherwise. The upper bound, however, predicts a pushed front for any b. These results are depicted in Fig. 7. It is necessary to comment that we have illustrated how to emulate nonlinear reaction terms using reaction term composed by two-segment. Our results are reasonably good compared to the exact known results. However, if one desires a higher accuracy in the piecewise-linear emulation one should consider more segments. This was already done in Ref. [15] where piecewise linear reaction terms of
ARTICLE IN PRESS V. Me´ndez et al. / Physica A 375 (2007) 51–64
63
2.2 2.0
cU
1.8 1.6 1.4
c
1.2 1.0 0.8
: analytic
cL
0.6 0.4 0.2
:
1 +b 2 2
:
2 b
0.0 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
b Fig. 7. Plot of the selected front speed obtained from the piecewise linear emulation for the nonlinear reaction term f ðuÞ ¼ uð1 uÞðu þ bÞ: The results obtained from the variational principle are depicted with dashed lines (upper and lower bounds), the analytic result is depicted by solid line and the speed for the front under emulation is depicted with symbols.
five segments were used to emulate reaction terms with cubic and quintic nonlinearities showing an excellent agreement with the exact results. 6. Conclusions One-component reaction–diffusion model with two-piece linear reaction functions have been analyzed. The development of variational principles, analytical solution and numerical calculations jointly let us show the existence of a pushed-to-pulled transition as the discontinuity in the reaction term a goes from 0 to 1 when fronts propagate into the unstable state. On the other hand, when the front propagates into metastable state, there is not such transition and fronts are always pushed. Specifically, for unstable state u ¼ 0 (n40) the variational principle points out that for lower values of discontinuity position a than a critical one, the front speed is higher than the linear speed and the front is then pushed. Otherwise, the linear speed is selected (because upper and lower bounds from the variational principle are equal) and the front is pulled. In consequence, the variational principle rigorously prove the existence of a transition from pushed-to-pulled front propagating into the unstable state when a goes from 0 to 1. This result is in agreement with the numerical solutions. Furthermore, piecewise linear reactions allow to work out the analytic solution for the front speed and we have seen how this solution holds only for discontinuity position ranging from 0 to a critical value. We have compared the different critical values of a from the variational principle and the analytic solution with numerical computations exhibiting a relatively good agreement. Nevertheless, when the front propagates into metastable state (no0), there is no linear speed selection and neither pushed-to-pulled transition. The analytic solution for the speed is, in this case, in agreement with the numerical solutions for any a, and may predict, contrarily to the variational principle, negative values for the speed. When a time delay is considered in our reaction–diffusion system, all the above arguments are still valid but it plays an important role on the width of the pushed and pulled regions: as delay increases, the pushed region shrinks while the pulled region expands. This fact has an intuitive explanation: time delay slows the front down and the system stays for more time near unstable state. Therefore, the linear part around u ¼ 0 becomes more important as the delay increase, hence the pulled region enlarges.
ARTICLE IN PRESS 64
V. Me´ndez et al. / Physica A 375 (2007) 51–64
We have finally analyzed how to emulate nonlinear reaction terms by two-piece linear reaction terms for fronts propagating into an unstable state. We have considered, in particular, the emulation of the reaction terms f ðuÞ ¼ uð1 uk Þ with kX1 and f ðuÞ ¼ uð1 uÞðu þ bÞ showing in both cases how the two-piece emulation captures the known qualitative behaviour of the front speed when it resembles close enough the nonlinear reaction. Acknowledgment This work has been supported by the MCYT under Grant FIS2006-12296-C02-01 and by Generalitat de Catalunya under Grant 2005 SGR-00087. References [1] H.P. McKean, Adv. Math. 4 (1970) 209; J. Rinzel, J.B. Keller, Biophys. J. 13 (1973) 1313; Y.B. Zeldovich, A.P. Aldushin, S.I. Kudayev, in: J.P. Ostriker (Ed.), Selected Works of Y. B. Zeldovich, vol. I, Princeton University Press, Princeton, 1992, p. 320; A. Ito, T. Ohta, Phys. Rev. A 45 (1992) 8374; V.S. Zykov, Simulations of Wave Processes in Excitable Media, Manchester University Press, New York, 1987. [2] S. Theodorakis, E. Leontidis, Phys. Rev. E 62 (2000) 7802. [3] S. Theodorakis, E. Leontidis, Phys. Rev. E 65 (2002) 026122. [4] M. Wang, M. Kot, M.G. Neubert, J. Math. Biol. 44 (2002) 150. [5] V. Me´ndez, A. Compte, Physica A 260 (1998) 90; K.K. Manne, A.J. Hurd, V.M. Kenkre, Phys. Rev. E 61 (2000) 4177; G. Abramson, A.R. Bishop, V.M. Kenkre, Phys. Rev. E 64 (2001) 066615. [6] G. Fath, Physica D 116 (1998) 176; A. Lahiri, P. Majumdar, M.S. Roy, Phys. Rev. E 65 (2002) 026106; A. Tonnelier, Phys. Rev. E 67 (2003) 036105. [7] A.M. Pertsov, M. Wellner, J. Jalife, Phys. Rev. Lett. 78 (1997) 2656; Y.B. Chernyak, J.M. Starobin, R.J. Cohen, Phys. Rev. Lett. 80 (1998) 5675; Y.B. Chernyak, J.M. Starobin, R.J. Cohen, Phys. Rev. E 58 (1998) R4108. [8] E.P. Zemskov, K. Kassner, S.C. Mu¨ller, Eur. Phys. J. B 34 (2003) 285; E.P. Zemskov, K. Kassner, S.C. Mu¨ller, Phys. Rev. E 70 (2004) 056208; R. Bakanas, Phys. Rev. E 69 (2004) 016103; R. Bakanas, Phys. Rev. E 71 (2005) 026201. [9] W. van Saarloos, Phys. Rep. 386 (2003) 29. [10] P. Clavin, A. Li na´n, in: M.G. Velarde (Ed.), Nonequilibrium Cooperative Phenomena in Physics and Related Fields, NATO Advanced Study Institute Series B: Physics, vol. 116, Plenum, New York, 1984, p. 291. [11] R.D. Benguria, M.C. Depassier, Phys. Rev. Lett. 77 (1996) 1171. [12] R.D. Benguria, M.C. Depassier, Phys. Rev. E 66 (2002) 026607. [13] D.G. Aronson, H.F. Weinberger, Adv. Math. 30 (1978) 33. [14] J. Fort, V. Me´ndez, Rep. Prog. Phys. 65 (2002) 895. [15] S. Theodorakis, E. Svoukis, Phys. Rev. E 68 (2003) 027201.