Travelling waves in a reaction-diffusion model for electrodeposition

Available online at www.sciencedirect.com

Mathematics and Computers in Simulation 81 (2011) 1027–1044

Travelling waves in a reaction-diffusion model for electrodeposition Benedetto Bozzini a , Deborah Lacitignola b,c , Ivonne Sgura b,∗ a

Dipartimento di Ingegneria dell’Innovazione, Universitá del Salento – Lecce, Via per Monteroni, I-73100 Lecce, Italy Dipartimento di Matematica “Ennio De Giorgi”, Universitá del Salento – Lecce, Via per Arnesano, I-73100 Lecce, Italy Dipartimento di Scienze Motorie e della Salute, Università di Cassino, Campus Folcara, Loc. S. Angelo, I-03043 Cassino, Italy b

c

Available online 21 October 2010

Abstract In this paper we consider an analytical and numerical study of a reaction-diffusion system for describing the formation of transition front waves in some electrodeposition (ECD) experiments. Towards this aim, a model accounting for the coupling between morphology and composition of one chemical species adsorbed at the surface of the growing cathode is addressed. Through a phasespace analysis we prove the existence of travelling waves, moving with specific wave speed. The numerical approximation of the PDE system is performed by the Method of Lines (MOL) based on high order space semi-discretization by means of the Extended Central Difference Formulae (D2ECDF) introduced in [1]. First of all, to show the advantage of the proposed schemes, we solve the well-known Fisher scalar equation, focusing on the accurate approximation of the wave profile and of its speed. Hence, we provide numerical simulations for the electrochemical reaction-diffusion system and we show that the results obtained are qualitatively in good agreement with experiments for the electrodeposition of Au–Cu alloys. © 2010 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: High order finite difference schemes; Method of Lines; Travelling waves; Fisher equation; Electrochemical modeling

1. Introduction and the model Electrodeposition of Au–Cu alloys from cyanocomplex baths containing free cyanide has been shown to exhibit electrokinetic instabilities that can lead to compositional heterogeneity in the electrodeposit bulk [4]. Such instabilities derive from a hysteretic current–voltage characteristic related to the buildup of CN− concentration in the catholyte and attending variations the Cu(I)- cyanocomplex nobility [6]. The technological drawbacks deriving from such instabilities can be efficiently removed by changing the solution chemistry: (i) removing the free-CN− and (ii) changing the Cu(I)cyanocomplex for a Cu(II)-EDTA compound [5]. Notwithstanding the technological improvement, a rich dynamic scenario is still attached to the presence of adsorbed CN− at the growing metal interface [9]. In this paper we consider a general model coupling the surface morphology and CN− concentration, in order to rationalise the formation of the morphological patterns sometimes developing in electrodeposition (ECD). Some dynamic details will be discussed with model parameter values typically corresponding to the more industrially widespread free-CN− solutions. This phenomenon is a special case, original within the realm of metal electrochemistry, of a more general, well-known type of chemical and electrochemical dynamics (electrocatalysis [18], corrosion of Cu [15], Fe-group metals [18] and Ag [8]).

∗

Corresponding author. Tel.: +39 0832 297591; fax: +39 0832 297594. E-mail addresses: [email protected] (B. Bozzini), [email protected] (D. Lacitignola), [email protected] (I. Sgura).

0378-4754/$36.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2010.10.008

1028

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

The proposed model reveals a surprisingly rich phenomenology: in [9] theoretical and numerical investigations were performed, showing the initiation of spatial patterns induced by diffusion, i.e. related to the mechanism introduced by Turing in [23]. In this paper, we show the existence of travelling wave solutions, corresponding to transition front waves, sometimes occurring in the experiments. We present an analytical study of the model, its numerical approximation by accurate space discretization and comparisons between simulations and electrochemical experiments. To the aim of completeness, we now briefly recall the model derivation proposed in [7,9]. First, we stress the following basic assumptions: (i) the electrodeposit surface can be regarded as isotropic in the substrate plane, so that we will consider a 1D model in space (ii) one chemical species adsorbed at the surface of the growing cathode is addressed, so that we deal with a system of two reaction-diffusion equations, one for the morphology and one for the chemistry. The equation for the morphological dynamics is: dη d2η = Ds∗ 2 + S ∗ , dτ dξ

(1)

with S ∗ = α1

η2 − β1 ηθ. 1+η

Here η(ξ, τ) is the dimensionless electrode shape (i.e. the intersection of the electrodeposit surface with a plane normal to the substrate), ξ is the dimensionless space coordinate, τ is the dimensionless time and Ds∗ is the dimensionless surface diffusion coefficient of adatoms. In the source term S∗ , θ(ξ, τ) is the surface coverage with the adsorbed chemical species, the parameters α1 and β1 are strictly positive and weight the two terms in S∗ accounting for: (i) localization of the ECD process and (ii) effects on the ECD rate of surface chemistry, i.e. the presence of adsorbates at the growing cathode, respectively. The presence of adsorbable species in the ECD bath gives rise to the fact that θ(ξ, τ) develops at a growing electrochemical interface as a function of space and time. θ(ξ, τ) is controlled by the nature of the adsorbable species and of the surface active sites. The surface coverage dynamics can be described, as customary in chemical kinetics, in terms of a material balance with a source term containing positive and negative contribution related to adsorption and desorption. More precisely, the dimensionless form for the surface chemical dynamics can be expressed as: 2 dθ ∗ d θ = Dsc + Sc∗ , dτ dξ 2

(2)

∗ is the surface diffusion coefficient of CN− and S ∗ is the chemical source term, given by where Dsc c ∗ ∗ Sc∗ = KADS (η, θ)(1 − θ) − KDES (η, θ)θ ∗ ∗ and KDES which represent the adsorption and desorption rate constants, respectively. and defined in terms of KADS By coupling Eqs. (1) and (2), one obtains the following general model in dimensionless form: ⎧ ∂η ∂2 η ⎪ ⎪ = + ρf (η, θ), ⎨ ∂τ ∗ ∂ξ 2 (3) 2 ⎪ ⎪ ⎩ ∂θ = d ∂ θ + σg(η, θ), ∂τ ∗ ∂ξ 2

with f (η, θ) =

η2 − ηθ, 1+η

∗ ∗ ∗ (η, θ) − [KADS (η, θ) + KDES (η, θ)]θ, g(η, θ) = KADS

(4)

which is defined for (ξ, τ ∗ ) ∈ [0, L] × [0, T ], with L a characteristic length of the electrode, T a characteristic time of the electrodeposition process, and where τ ∗ = Ds∗ τ

d=

∗ Dsc , Ds∗

=

α1 , β1

σ=

1 , Ds∗

ρ=

β1 . Ds∗

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

1029

We also require (3) to be equipped with zero-flux boundary conditions and the following initial conditions: η(ξ, 0) = η0 ,

θ(ξ, 0) = θ0 ,

ξ ∈ [0, L] .

Motivated by experimental evidence, it is possible to explicitly choose the adsorption and desorption rates, in the following form: ∗ KADS = A exp(aη + bθ),

∗ KDES = A1 exp(a1 η + b1 θ),

(5)

where A, A1 , a, a1 , b, b1 are positive parameters. More precisely, we consider the case, A = / A1 , a = / a1 , b = b1 , and choose a1 = a +  ln (A/(A1 )) and A/A1 = /2. In the next, we deal with the reaction-diffusion model (3), whose reaction terms f and g are specialized as: f (η, θ) =

η2 − ηθ, 1+η

g(η, θ) =

[A exp(aη + bθ)]( − θ − θ21−η ) . 

(6)

The rest of the paper is organized as follows: in Section 2, we will perform the analytical investigations on the PDE system (3)–(6), to derive general conditions for the existence of travelling wave solutions. Hence, in Section 3 we deal with numerical approximation of the electrodeposition model by the Method of Lines (MOL) based on high order finite difference schemes in space, which are extension of the well-known central difference formula for second order derivatives. Moreover, in this section we present the advantage of using these schemes by solving the well-known scalar Fisher equation, focusing on the accurate approximation of the wave front and its speed. In Section 4 we present some numerical experiments for the PDE reaction-diffusion system for significant parameter values that yield transition front waves of invasion and pursuit and evasion type, as announced by the theoretical analysis. In the last part of the section, we outline that these results are qualitatively in good agreement with electrodeposition experiments of gold–copper (Au–Cu) alloys. 2. Travelling wave solutions Systems of reaction-diffusion equations from applied sciences usually support travelling wave solutions, see e.g. [20] and references therein. In this section, we want to show that for the model (3) travelling waves should be expected for a continuous range of velocities. The waves are of transition front type, as the ones discussed by Fisher, Kolmogorov et al. [13,17] for a scalar reaction diffusion equation. In some cases, they may be not necessarily monotone [12,20]. In the present case, a distinguished value of the wavespeed c, say c∗ = c∗ (d, ρ, ), is found such that for c ≥ c∗ the travelling waves are of one type, i.e. approaching the steady state solution monotonically, whereas for c < c∗ they approach the steady state solution in a damped oscillatory manner. The former type of waves, representing invasion processes, are called waves of invasion [22]; the latter type of waves are referred as waves of pursuit and evasion (see [20] for more details). We introduce the new variable z = ξ − cτ ∗ , where c is a positive constant and −∞< z < ∞. To search for profiles representing an invasion process, we have to look for positive, monotonically, increasing functions satisfying:  lim η = 0; lim θ = ; z→−∞ z→−∞ +2 (7) 1  lim η = ; lim θ = . z→+∞ z→+∞  +1 Differently, profiles representing waves of pursuit and evasion, satisfy (7) but approach the steady state through damped oscillations. By considering, η(ξ, τ ∗ ) = η(ξ − cτ ∗ ) = η(z),

θ(ξ, τ ∗ ) = θ(ξ − cτ ∗ ) = θ(z),

the PDE system (3) is converted into the following system of ordinary differential equations:   η + c η + ρf (η(z), θ(z)) = 0, 

dθ + c θ  + σg(η(z), θ(z)) = 0,

(8)

1030

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

with boundary values η(−∞) = 0, η(+∞) =

1 , 

 , +2  θ(+∞) = , +1

θ(−∞) =



where represents the first derivative with respect to the z variable and η( − ∞)= 0 means lim z→−∞ η(z) = 0 and so on. The study of the travelling waves for the electrodeposition model (3), involves a four-dimensional phase-space. The boundary value problem (8) is in fact equivalent to the following system of four differential equations, ⎧  η = u, ⎪ ⎪ ⎪ ⎪ ⎨ θ  = v, (9) u = −cu − ρf (η, θ), ⎪ ⎪ ⎪ c σ ⎪ ⎩ v = − v − g(η, θ), d d with the equivalent boundary conditions:  η(−∞) = 0, θ(−∞) = , u(−∞) = 0, v(−∞) = 0, +2 1  η(+∞) = , θ(+∞) = , u(+∞) = 0; v(+∞) = 0.  +1 Travelling wave solutions for system (3), can thus be interpreted as heteroclinic trajectories of the dynamical system (9) in the four-dimensional phase space (η, θ, u, v), connecting the steady state E1 = (0, (/ + 2), 0, 0) to the steady state E2 = ((1/), (/ + 1), 0, 0). A necessary condition for the existence of such trajectories is that the steady state E1 must exhibit an unstable manifold whereas the steady state E2 must have a stable manifold. The Jacobian matrix for system (9) is: ⎡ ⎤ 0 0 1 0 ⎢ 0 0 0 1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ k ρ η −c 0 ⎥ ⎣ ⎦ c ψp ψq 0 − d where 2 η  +  η2 − θ − 2 η θ − θ η2 σ Aeaη+bθ , , ψ = d (1 + η)2 p = −a + aθ  + aθ 21−η  − θ 21−η   ln 2, q = −b + bθ  + bθ 21−η  +  + 21−η  .

k=−

By linearizing about the steady state E1 , one obtains the characteristic polynomial, ( dλ2 +  cλ − σ Aeb/(+2)  − 2 σ Aeb/(+2) )( λ2 + 2 λ2 +  cλ + 2 cλ − ρ ) = 0, which roots are:

1 λ1,2 = (− c ∓ 2 c2 + 4 2 dσ Aeb/(+2) + 8  dσ Aeb/(+2) ), 2d and

−c( + 2) ± c2 ( + 2)2 + 4 ρ ( + 2) λ3,4 = . 2( + 2) These expressions confirm that the steady state E1 is hyperbolic and exhibits a 2-dimensional unstable manifold, spanned by those eigenvectors corresponding to the real and positive eigenvalues [14]. Besides, as the eigenvalues are real, it is not possible to find oscillatory behavior in the neighboring of E1 .

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

1031

By linearizing about the steady state E2 , the following characteristic polynomial is obtained: p2 (λ) = λ4 + B3 λ3 + B2 λ2 + B1 λ + B0 ,

(10)

where the coefficient Bi are given by: (d + 1) c c2 A( + 1)σeβ 2 ρ − , B2 = + , d d d ( + 1)2 c2 ρ A( + 1)cσeβ Aσρeβ , B − = (ln 2 − ) B1 = , 0 d ( + 1)d d( + 1)2 B3 =

with β = 

a+a+b2 (+1) .

 +1

3

If we choose the parameter values such that:

e−β =

Aσ , ρ

(11)

as a consequence, the coefficient B1 vanishes and we have   1 c2 3 ρ 2 2 ρ . B2 = 1 − = (ln 2 − ) + , B 0 d d ( + 1)2 ( + 1)4 d To investigate on the eigenvalues of the characteristic polynomial (10), we preliminary observe that B3 is positive since all the involved parameters are positive constants; B2 is positive if c2 > 2 ρ/( + 1)2 (1 − d) and negative otherwise. Moreover, B0 is positive if  < ln 2 holds. In the following, we consider a region in the parameter space such that d and  satisfy d < 1 and  < ln 2. This is perfectly consistent with the present experimental set-up. In fact, the physical meaning of d is the ratio of DCN , the diffusion coefficient of CN− on the Au surface, to DAu , that of Au adatoms on the Au surface itself: DCN /DAu . Within an Arrhenius framework, the diffusion coefficient DX of species X can be regarded to depend exponentially on the site-hopping energy ES : DX ≈ exp(−ES /RT), where R is the gas constant and T the absolute temperature. For our present purposes, ES can be identified with the adsorption energy: since, in the case of free-CN− baths, with a high concentration of CN− at the equilibrium, we have ECN > EAu (more specifically, as a rule of thumb: ECN ≈ 2EAu , see e.g. [16,24]), then it follows that d < 1. For this choice of the parameters d and , one finds through the Descartes’ rule of sign, that the characteristic polynomial (10) has no real positive roots and two or zero real negative roots. Since we are interested in the stable manifold of the steady state E2 , we have to look at the negative roots of the characteristic polynomial (10). To this aim, it is useful to consider (10), as a c-family of λ polynomials, i.e.   2 (1 + d) c 3 (ln (2) − ) 3 ρ2 c (12) p2 (λ, c) = λ4 + λ + + 2 ρ (1 − d −1 )( + 1)−2 λ2 + d d ( + 1)4 d A qualitative picture of (12) for different values of c (with the other parameters fixed), is represented in Fig. 1. Here it is shown that a critical value of the parameter c, say c∗ = c∗ (, ρ, d), exists such that for c < c∗ this polynomial has two complex conjugate eigenvalues with negative real part (and two complex conjugate eigenvalues with positive real part) whereas for c > c∗ it has two real negative eigenvalues (and two complex conjugate eigenvalues with positive real part). In any case, the steady state E2 exhibits a 2-dimensional stable manifold, spanned by those eigenvectors corresponding to the eigenvalues with negative real part. Complex conjugate roots with negative real parts ensure a solution for the boundary value problem (9) characterized by the fact that the resulting trajectory will approach the steady state E2 with damped oscillations, thus representing waves of pursuit and evasion type. Oscillatory behaviors around E2 will be excluded when the characteristic polynomial (12) has two real and negative roots: this circumstance will instead produce waves of invasion type. Hence, we have to find the critical value c = c∗ , discriminating these two types of waves. By considering ∂p2 (λ, c) =0 ∂λ

1032

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044 100

80

c < c* c = c* c > c*

60

40

20

0

−20 −5

−4

−3

−2

−1

0

1

2

3

4

5

λ Fig. 1. Qualitative graphs of the characteristic polynomial (12) for different values of c. The other parameters have been fixed as d = 0.4,  = 0.5, ρ = 40. For these parameter values c∗ = 0.3040 holds.

and solving for λ, one gets the turning points of (12):  3 1 9 B32 − 32 B2 . λ = 0, λ∓ (c) = − B3 ∓ 8 8 Now, substituting λ− (c) in (12), c∗ (, ρ, d) is obtained as solution of the following algebraic equation: √ √ (ln 2 − ) 3 ρ2 9 (d + 1)3 c3 δ 9 δ1 (d + 1)2 c2 δ1 (d + 1) c δ δ1 2 27 (d + 1)4 c4 =0 − + + − + − 512 d4 512 d3 32 d2 16d 4 ( + 1)4 d (13) where δ=9

(d + 1)2 c2 c2 32 2 ρ (1 − d −1 ) − 32 , − d2 d ( + 1)2

δ1 =

c2 (1 − d −1 ) . + 2 ρ d ( + 1)2

Note that, for the previous choices of  and d (i.e.  < ln 2 and d < 1), δ is a positive quantity, so that (13) is algebraically well defined. For c ≥ c∗ , one obtains the two real and negative eigenvalues for (12) so that the oscillatory behavior about E2 will not occur. Therefore, c ≥ c∗ is a necessary condition for the existence of invasion travelling wave solutions. Fig. 2 gives an indication of the curve c = c∗ (d), when the parameters  = 0.5 and ρ = 40 have been fixed. In the region A, the inequality c ≥ c∗ (d) holds and so, for this continuous range of travelling wave speeds, we might expect the existence of invasion travelling wave solutions. In the region B, the inequality c < c∗ (d) holds and so, in this continuous range of travelling wave speeds, we might expect the existence of pursuit and evasion travelling wave solutions. The existence of a continuous range of acceptable travelling wave speeds, for any fixed value of the diffusion parameter d, implies the following questions: given an initial profile, which of these travelling wave solutions, in practice, may be obtained? And which is the related speed? From a theoretical point of view, stability arguments for travelling wave solutions might provide the answer [19–21]. However, because of the peculiar form of the reaction terms in system (3) and because of the many parameters involved, a rigorous mathematical approach to such a problem is very complex to be faced and falls out of the aim of this paper. For this reason, in Section 4 we will show that numerical simulations provide qualitative and quantitative answers to these questions. In fact, for a fixed d, the numerical approximation of the reaction-diffusion system (3), for quite general initial conditions, identifies numerical profiles moving with an approximated wave speed cnum (d) that can be: (i) invasion travelling waves, if the velocity is c = cnum ≥ c∗ or (ii) pursuit and evasion travelling waves, if the velocity is c = cnum < c∗ . More precisely, by fixing the parameter values of  and ρ as in Fig. 2, we find that the type of waves which is observable

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

1033

4 3.5 3 2.5

c

A

2 1.5 B

1 0.5 0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

d Fig. 2. Case  = 0.5, ρ = 40. The curve c = c∗ (d) separates the two regions A and B, in the (d,c) parameter space. We denote with A, the region in the (d,c) parameter space where c ≥ c∗ (d) holds and with B, the region in the (d,c) parameter space where c < c∗ (d) holds.

(and thus stable), depends on the chosen value of the parameter d. This seems to stress the very important role of this parameter for the expected phenomenology. In the next section, we present in details the numerical approximation of the reaction-diffusion system (3). The question of the numerical estimation of the travelling wave speed is also widely discussed. 3. Numerical approximation by high order finite difference schemes To solve the PDE reaction-diffusion system given in (3), we apply the classical Method of Lines (MOL), hence we transform the problem in an initial value ODE problem by considering a finite difference (FD) approximations of high order for the spatial derivatives. Towards this aim we apply the FD methods introduced recently in [1,2]. To deal with the approximation of second order derivatives we use the generalization of central differences, called D2ECDF. Moreover, to deal with the zero Neumann boundary conditions, instead of the usual ghost point technique, we apply on the boundary mesh points the analogous generalization of central formulae devised for the approximation of first order derivatives. First of all, in the following subsection, we describe the D2ECDF schemes and their application for a space semidiscretization in both cases of Dirichlet boundary conditions and (zero) Neumann boundary conditions. Some emphasis is stressed on the approximation of the boundary conditions with the same high order of the schemes used in the interior domain, so that the global order is preserved along all the given mesh. In fact, FD schemes of high order on a long stencil require additional values near the boundaries that usually are calculated by lower order methods, and this implies a reduction in the approximation order. In the second subsection, to show the advantage of using high order approximation in space, we solve a simple scalar reaction diffusion equation, known as Fisher equation, which presents a kink travelling wave solution with theoretical minimum velocity c∗ = 2. Indeed, recent studies [3,10] prove that in presence of a cut-off on the solution, say tol, the velocity of the front is indeed shifted and is given by c¯ = c¯ (tol) = 2 − (π/ ln(tol))2 < c∗ . We present the numerical approximation of the Fisher equation by the MOL technique based on the D2ECDF methods of order p = 2, 4, 6 and some comparisons with results obtained by the Matlab function pdepe based on the finite-element method and more accurate in time (ode15s) than space (p = 2). For aim of comparisons, we have solved the resulting IVP-ODE system in time by means of the Matlab solver ode15s. We show that, in dependence of the stepsize h = hx , the numerical fronts of the kink are indeed shifted and for h → 0 all schemes tend to produce a kink with same asymptotic speed cnum ≈ c¯ (tol), where c¯ (tol) is obtained in correspondence of the value of tol related to the absolute error in the approximation of the zero part of the kink.

1034

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

Fig. 3. Representation of all possible stencils on k + 1 points for the approximation of the derivatives y(ν) (xi ).

In particular, we show that the D2ECDF of order six approximates carefully this velocity also for large stepsize h. Moreover, we have observed that the numerical approximation of the kink solution and the delayed behavior observed on the velocity reduction is not affected by the order of approximation in time. Thank to these results, in the Section 4, we solve the PDE reaction-diffusion system in (3) for some significant choices of the parameters by means of all D2ECDF schemes proposed and again we show some comparisons with the solutions obtained by the pdepe Matlab solver. Hence, we discuss the numerical features of the approximated kinks in terms of the diffusion coefficient d.

3.1. Extended Central Difference Formulae (D2ECFD) Let consider on the (space) interval [x0 , xf ] a meshgrid on n interior points with uniform stepsize h = (xf − x0 )/(n + 1). (ν) Let be k > 2. For the approximations yi ≈ y(ν) (xi ), ν = 1, 2, i = 1, . . ., n, we consider a non-compact (long) stencil of k + 1 points with s (initial) values on the left and k − s (final) values on the right of each discretization point xi . In Fig. 3 we show a sketch of this idea. Note that, in the present application related to the Method of Lines (MOL), the previous unknowns are indeed time (ν) (ν) dependent functions, that is we have yi ≡ yi (t) ≈ y(ν) (xi , t) = (∂y/∂t)(xi , t) for t ∈ [0, T]. For sake of exposition, we drop the dependence on t, rather than not explicitly reported. (ν) Hence, by construction, there exist k + 1 possible schemes Fs , such that each formula requires s initial values and k − s final values, for s = 0, . . ., k, given by

(2)

Fs

:

y (xi ) :=

k−s ∂2 y 1  (s) (x , t) ≈ αj+s yi+j , i ∂x2 h2

s = 0, . . . , k

j=−s

(1) Fs

:

k−s ∂y 1  (s) y (xi ) := (xi , t) ≈ βj+s yi+j , ∂x h

(14) s = 0, . . . , k.

j=−s

(s)

(s)

The coefficients αj+s and βj+s can be calculated in order to have for all s the same maximum possible order p = k. The main result given in [1] can be described as follows. (ν) For all ν ≤ k and for s = 0, . . ., k, the coefficients of the formula Fs of order k to approximate y(ν) (xi ) using the stencil [xi−s , . . ., xi+k−s ] are the entries of the (s + 1) th column of the matrix C(ν) = V −1 Kν V,

C(ν) ∈ R(k+1)×(k+1) ,

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

1035

where we have the Vandermonde matrix ⎛

1

⎜ ⎜0 ⎜ ⎜ V =⎜ ⎜ .. ⎜. ⎝ 0

1

1

...

1

2

...

.. . 1

2k

...

1

⎞

⎟ k⎟ ⎟ ⎟ .. ⎟ ⎟ . ⎟ ⎠

kk

⎛

⎞

0

⎜ ⎜ ⎜1 ⎜ and K = ⎜ ⎜ ⎜ ⎜ ⎝

..

.

..

.

⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎠

0 k

0

It is easy to show that C(ν) is a centro-symmetric matrix for derivatives of even order and a skew-centro-symmetric matrix for derivatives of odd order. Hence the corresponding property of symmetry and skew-symmetry is inherited by the corresponding schemes. For example for k = 4, we have ⎛ 35 ⎛ 25 11 1 1 11 ⎞ 1 1 1 1 ⎞ − − − − − ⎜ 12 ⎜ 12 12 12 12 12 ⎟ 4 12 12 4 ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ⎜ ⎟ ⎟ 5 4 1 14 ⎟ 5 2 1 4⎟ ⎜ 26 ⎜ ⎜− ⎜ 4 − − ⎟ − − − ⎟ ⎜ 3 ⎜ 3 3 3 3 ⎟ 6 3 2 3⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ⎜ ⎟ ⎟ 3 3 5 1 1 19 ⎟ ⎜ 19 ⎜ ⎟ (2) (1) C =⎜ ⎟ and C = ⎜ −3 0 − 3 ⎟. − ⎜ 2 ⎜ ⎟ ⎟ 2 2 2 2 2 2 ⎜ ⎜ ⎟ ⎟ ⎜ ⎜ ⎟ ⎟ 1 2 5 1 5 26 ⎟ 4 ⎜ 14 ⎜ 4 ⎟ ⎜− ⎜ − −4 ⎟ − − ⎟ ⎜ 3 ⎜ ⎟ ⎟ 2 3 6 3 3 3 3 ⎟ ⎜ ⎜ 3 ⎟ ⎝ ⎝ ⎠ ⎠ 1 1 1 1 25 11 1 1 11 35 − − − − 4 12 12 4 12 12 12 12 12 12 To construct the matrix operator of even order k, for the approximation of the ν-order derivative along the (space) domain, we will use all the schemes generated combined in matrix form as follows. First of all, we need to fix a main scheme with s¯ initial values to be used in all possible interior points of the mesh, then at the beginning and at the end of the interval the remaining schemes (i.e. columns in C(ν) ) will be suitably used. Following these arguments, in the papers [1,2] we have defined the Extended Central Difference Formulae (ECDF) for which s¯ = k/2, that are high order extension of the well known central difference formula of order p = 2. In particular, they have been called D2ECFD when used for the second order derivative approximation. The main scheme corresponds to the (k/2 + 1) th column in C(ν) and then the s¯ − 1 = k/2 − 1 schemes on the left will be used at the beginning of the interval, while the k/2 − 1 schemes on the right will be used at the end of the interval. If Dirichlet boundary conditions are assigned in the differential problem, the first and last columns are never used. As we will explain later, they will be used when Neumann boundary conditions are assigned. A sketch of the application of D2ECDF schemes for p = k = 4 and s¯ = 2, is given in Fig. 4. Hence, for the approximation of the second order derivatives by D2ECDF, the matrix operator is given by

(15)

1036

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

Fig. 4. Sketch of ECDF schemes of order p = k = 4 for the approximation of the derivatives y(ν) (xi ).

˜ = [b0 Similar construction for the operator approximating the first order derivative will yield B T T T If Y = [y1 , y2 , . . ., yn ] and Y˜ = [y0 , Y , yf ] , then our FD schemes can be set in matrix form as Y  (x) ≈

1 ˜ ÃY , h2

Y  (x) ≈

B

bk ].

1˜˜ BY , h

˜ are n × (n + 2) matrices. We define the matrix operator Mp(2) := A and Mp(1) := B to yield a space where à and B discretization of global order p = k for the derivatives involved in the problem. If Dirichlet b.c. are given, then y0 and yf are known, and the vector accounting for the boundary conditions is given by bc ≡ bcD = a0 ⊗ y0 + ak ⊗ yf ∈ Rn ,

(16)

where a0 and ak are the first and last column in the matrix (15). If zero Neumann b.c. are given, as in our reaction-diffusion PDE system, we consider the first and last methods in C(1) , since they do not require initial and final values, to obtain for all t, the approximation of y0 ≈ y(x0 , t) and yn+1 ≈ y(xn+1 , t) as follows: ∂y (x0 , t) = y (x0 ) = 0 ∂x

(1)

⇒

(0)

⇒

(0)

βj yj = 0

(1)

Fk :

(k)

βk yn+1 +

−k 

(k)

βk+j yn+1+j = 0.

j=−1

yn+1 = − (0)

k  j=1

∂y (xn+1 , t) = y (xn+1 ) = 0 ∂x Then we obtain 1 y0 = − (0) v0 Y β0

(0)

β 0 y0 +

F0 :

1

(17)

v Y, (k) f

βk

T

(k)

(k)

T

where v0 = [β1 , . . . , βk , 0, . . . , 0] are vf = [0, . . . , 0, β0 , . . . , βk−1 ] ∈ Rn . This implies that the vector accounting for the approximation of (zero) Neumann boundary conditions is given by   1 1 (18) bc = a0 ⊗ y0 + ak ⊗ yf = − a0 ⊗ (0) v0 + ak ⊗ (k) vf Y = bcN Y β0 βk where bcN ∈ Rn×n .

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

1037

In the next subsection we will apply the D2ECDF schemes of order p = 2, 4, 6 for the semidiscretization of the Fisher equation equipped with Dirichlet b.c. For this reason, in the definition of the corresponding IVP-ODE, we have to consider the contribution due to boundary conditions given by the vector bcD in (16). In the last subsection, we will apply these same schemes to solve the electrochemical reaction-diffusion system equipped with the Neumann conditions. In this case, the semi-discretization must account for the boundary conditions contribution expressed in terms of (18). 3.2. Numerical solution of the Fisher equation The Fisher equation is one of the simplest model describing the space and time evolution of a concentration u(x, t) in a reaction-diffusion phenomenon by means of the propagation of travelling waves connecting a stable and an unstable state. The model is given by ⎧ ⎪ ⎨ ut = uxx + u(1 − u) u(x, 0) = u0 (x), x0 < x < xf (19) ⎪ ⎩ u(x0 , t) = 1, u(xf , t) = 0, t > 0. The kink wave solution connects the uniform solutions u = 1 and u = 0, in principle on an infinite interval, that is for x ∈ [x0 , xf ] for x0 → − ∞ , xf → + ∞. The analytic form of the travelling wave solution u(x, t) = F(x − c t) is not known, but it is possible to show that for a steep initial condition u0 (x) the asymptotic solution attained for large t is that with minimal speed c∗ = 2. This is indeed an ideal situation. In fact, in real models the concentration is never zero, but at least larger than a small value, say tol, which is called cut-off in the literature. The main effect of having u ≥ tol > 0 introduces a perturbation in the tail of the front and this changes its speed. The effect of the cut-off on the travelling wave has been studied by many authors (see e.g. [11]) and in particular for the Fisher equation by [3,10]. In presence of this threshold, the reaction term in (19) has been defined as f (u) = 0,

for

0 < u < tol;

f (u) = u(1 − u),

for

tol ≤ u ≤ 1

(20)

and an estimate for the speed of the corresponding pulled front is given in [10] by c(u) ≈ c¯ = c¯ (tol) := 2 −

π2 (ln tol)2

(21)

Moreover, in [3] the authors refine this estimate and give an upper and lower bound cLOW (tol) < c(u) ≤ cUP (tol), with   1 1 ∗ cUP (tol) = 2 sin(ψ ) = c¯ + O , ψ∗ s.t. ψ∗ tan(ψ∗ ) = | ln(tol)|. (22) 3 2 (ln tol) where c¯ is given by (21). From the numerical point of view, the threshold tol can be related at least to round-off errors, but also to the truncation error introduced by the approximation of the second order derivative in (19) and to the error in the approximation of the zero part of the kink solution. A careful analysis would be needed to relate the results reported, for example, in [3] to the effect of cut-off induced by the numerical approximation. In this paper, we have solved the Fisher equation by the solver pdepe in Matlab and by the MOL method based on the D2ECDF of order p = 2, 4, 6. Then the semidiscretized problem for (19) on a finite interval [x0 , xf ] with space stepsize h = (xf − x0 )/(n + 1) on the grid x = (x1 , . . ., xn ) is given by the ODE-IVP ⎧ ⎨ dY = M Y + bc + f (Y ) p D dt ⎩ Y (0) = Y (0) . with Mp ≡ Mp(2) in (15) and bcD given in (16). For all schemes the time integration is performed on [0, T] = [0, 30] by using the solver ode15s with an input meshgrid of time step ht = 0.05. For the following numerical experiments, we have considered [x0 , xf ] = [0, 80] and as initial profile Y(0) the evaluation on the meshgrid x of the step function u0 (x) = 1

for

x ≤ 10,

u0 (x) = 0

for

10 < x ≤ 80.

1038

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044 1.2

time

time

h=0.8

1

1

0.8

0.8

0.6

0.6

u(x,t)

u(x,t)

1.2

0.4

0.4

0.2

0.2

0

0

−0.2

0

20

40

60

x

80

−0.2

0

20

h=0.1

40

60

80

x

Fig. 5. Numerical fronts for all methods, when h = 0.8 (left plot) and h = 0.1 (right plot): pdepe (dashed line), D2ECDF of order p = 2 (red line), p = 4 (blue line), p = 6 (black line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

We have applied all the schemes for the space stepsizes h = 0.1, 0.2, 0.4, 0.8. In Fig. 5 we show the wave profiles at the last ten time positions obtained for h = 0.8 (left plot) and h = 0.1 (right plot). In the first case for h = 0.8, it is evident that: (i) the two schemes of order p = 2 exhibit a faster front propagation w.r.t. the higher order schemes, and pdepe is faster than D2ECDF; (ii) the D2ECDF schemes of order p = 4 and p = 6 give almost the same approximation with almost the same velocity. In the case h = 0.1, all the schemes seem to propagate with the same velocity. To give a quantitative measure for the numerical propagation speed cnum of the numerical solutions, we devise an algorithm based on a spline technique, described as follows. Let Y(i) be the numerical approximation of the solution at the time step ti along the space meshgrid x. For each i, we calculate the corresponding spline approximation S(i) = S(Y(i) ) for the data points (Y(i) , x). If S˜ (i) = S(Y˜ (i) ) is the spline evaluation on a uniform grid of function values 0 ≤ (Y˜ (i) )j ≤ 1, j = 1, . . ., m we calculate the speed approximation as Ci =

S˜ (i+1) − S˜ (i) , ti+1 − ti

i = 1, . . . , nt − 1.

Note that each S˜ (i) is a vector of the same dimension of Y˜ (i) . This computation is necessary since in general the approximated values of the solution in xj at the time ti and in xj+1 at the time ti+1 , respectively, are different, that is / (Y (i+1) )j+1 , and then the shift of the front cannot be easily computed. (Y (i) )j = In Fig. 6, for h = 0.1 and the D2ECDF of order p = 6, we show the values (x, Y(i) ) (symbol ‘o’) and the corresponding splines (S˜ (i) , Y˜ (i) ) (symbol red ‘.’) obtained in the last ten time steps ti . In Fig. 7 we show the values Ci obtained by all D2ECDF schemes and by pdepe at the same time steps, when h = 0.8 and h = 0.1. From Fig. 7 it is evident that for h → 0 all schemes tend to the have the speed already obtained for the larger stepsize by the method of order p = 6. By calculating, for example, the speed of the fronts C(r) , for r ∈ {2, 4, 6, pdepe}, in Y˜ (i) = 1/2 and in the last time step i = nt , we can give a quantitative estimate of the convergence w.r.t. the theoretical value c¯ given in (21) and a check w.r.t the upper bound cUP in (22). If we consider as cut-off the threshold tol = 1e−6, which is the default absolute tolerance for each component of the solution used in ode15s (true for all methods considered), we obtain by (21) and (22) the values c¯  1.9489103 and cUP  1.9615511, respectively. Note that, in particular, this tol value can be related to the error approximation on the zero part of the kink. For r ∈ {2, 4, 6, pdepe}, in Table 1 we report the front speed approximation values C(r) (h) obtained for different space stepsize h and the corresponding relative error estimates err(r) = |C(r) − c¯ |/|¯c|. We remark that for large stepsize the solver pdepe and the D2ECDF of order p = 2 have numerical speeds greater than the theoretical upper bound cUP and then also a large relative error w.r.t. c¯ . This gap is reduced for decreasing h. The methods of order p = 4, 6 instead give for all h numerical speeds less than cUP and they approximate well the c¯ value also for large stepsize. Note that the best approximation is obtained by the method of order p = 6 for h = 0.4, that is n + 1 = 200 grid points. These results confirm those shown in Fig. 7.

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

1039

Spline approximation

75 70 65 60

x

55 50 45 40 35 30

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

u(x,t ) i

Fig. 6. Spline approximation in the last ten time steps for the D2ECDF of order p = 6 when h = 0.1.

Thanks to these results, to solve our electrodeposition PDE system, we will apply all D2ECDF schemes with moderate space stepsize and in the Section 4 we will compare again the results with those obtained by pdepe. Towards this aim, in the following subsection, we present the semidiscretized problem for the electrochemical reaction-diffusion system.

p=2

2

c(u)

c(u)

1.96

h=0.8

1.98

p=4

1.97

1.96

h=0.1

1.95

h=0.1

1.95

1.94

1.94 1.93

h=0.8

0

0.2

0.4

0.6

0.8

1.93

1

0

0.2

0.4

u p=6

1.97

1

0.8

1

h=0.8

2.1 h=0.8

c(u)

c(u)

0.8

PDEPE

2.15

1.96 1.95

0.6

u

2.05 2

1.94

h=0.1

1.97 1.95 1.93

1.93 0

0.2

0.4

0.6

u

0.8

1

h=0.1 p=6

0

0.2

0.4

0.6

u

Fig. 7. Speeds of the numerical fronts for all methods. The green dotted line represents the solution obtained by D2ECDF of order p = 6 and is reported on all subplots for aim of comparisons. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

1040

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

Table 1 Front speed approximations C(r) (h) and relative errors err(r) w.r.t. c¯ = c¯ (tol)  1.9489103, tol = 1e − 6. The values in bold correspond to the best performance obtained (minimum error), as described in the text. Cnum (h)

h = 0.8 h = 0.4 h = 0.2 h = 0.1

err(h)

pdepe

p=2

p=4

p=6

pdepe

p=2

p=4

p=6

2.1423 1.9999 1.9621 1.9526

1.9960 1.9627 1.9526 1.9509

1.9460 1.9490 1.9492 1.9506

1.9507 1.9487 1.9495 1.9497

9.9561e−2 2.6503e−2 7.0702e−3 2.2122e−3

2.4479e−2 7.3867e−3 2.2167e−3 1.3211e−3

1.1510e−3 3.6486e−4 4.4843e−4 1.1785e−3

1.2593e−3 1.9415e−4 6.0523e−4 7.4414e−4

3.3. PDE system approximation We sketch the MOL method for the reaction-diffusion system with zero-Neumann boundary conditions as follows. Let be h = L/(n + 1), then η˜ i (t) = η(ξi , t), θ˜ i (t) = θ(ξi , t) on the meshgrid ξ i = i h for i = 0, . . ., n + 1. The zero Neumann b.c. imply: ∂˜η0 ∂˜ηn+1 (t) = (t) = 0, ∂ξ ∂ξ

∂θ˜ 0 ∂θ˜ n+1 (t) = (t) = 0, ∀t. ∂ξ ∂ξ T

If η˜ = [˜η1 (t), . . . , η˜ n (t)]T and θ˜ = [θ˜ 1 (t), . . . , θ˜ n (t)] we have to solve the ODE-IVP system of 2n equations given by ⎧ d η˜ ⎪ ˜ ⎪ = (Mp + bcN )˜η + ρf (˜η, θ) ⎪ ⎪ ⎨ dt d θ˜ ˜ = d(Mp + bcN )θ˜ + σg(˜η, θ) ⎪ ⎪ ⎪ dt ⎪ ⎩ ˜ = θ˜ (0) η˜ (0) = η˜ (0) , θ(0)

(23)

For all t, Mp ≡ Mp(2) is again the matrix operator in (15) for the discretization of the ∂2 y (ξ, t), ∂2 ξ for both y = η and y = θ and bcN is the matrix in (18) accounting for the boundary conditions. In the classical approach, the central finite differences of order p = 2 can be applied together with the ghost point technique to deal with the zero Neumann boundary conditions, that is: y (ξi ) ≈

yi−1 − 2yi + yi+1 , h2

y (ξ0 ) ≈

y1 − y−1 = 0, 2h

y (ξn+1 ) ≈

yn+2 − yn = 0. 2h

In this way y0 and yn+1 are included as unknowns and then the vector of unknown functions Y = Y (t) = η˜ (t), or ˜ belong to Rn+2 . The matrix operator is given by M2 = (1/h2 )tridiag([1, − 2, 1])n+2 and in (23) the Y = Y (t) = θ(t) boundary conditions are given by the matrix bcN such that (bcN )1,2 = 1 = (bcN )n+2,n+1 and zero otherwise. In this paper, we apply the D2ECDF schemes together with the matrix bcN in (18) which will yield an approximation on the boundaries of the same order than in the interior domain. Note that, in our approach we solve the discrete problem on Rn and recover the values of y0 and yn+1 from (17). 4. Numerical examples for the electrochemical model Let be d < 1 and  < ln (2) and fix the other parameters such that  = 0.5, ρ = 40. We observe that, in this case, the condition (13) is function only of the parameters c and d. We choose also the values of the remaining parameters, so that (11) is verified. For example, we fix a = 0.25, b = 2.8, A = 2.5, σ = 0.141340.

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

1041

Fig. 8. Numerical solution of the PDE system (3) obtained by D2ECDF of order p = 6. Left plot: case d = 0.4. Travelling waves of invasion type in 3D for η(ξ, τ ∗ ). Right plot: case d = 0.9. Travelling wave of pursuit and evasion type in 3D for η(ξ, τ ∗ ).

We solve the system with zero-flux boundary conditions (Neumann) and the following initial conditions: ⎧  ⎧ ⎪ 0 if ξ ≤ ξ ⎪ ⎪ 0 ⎨ ⎨  + 2 if ξ ≤ ξ0 η0 = 1 θ0 = ⎪ ⎪  ⎩ ⎪ if ξ > ξ0 ⎩ if ξ > ξ0  +1 for 0 ≤ ξ ≤ L with L = 60 and ξ 0 = 10. We solve the system for t ∈ [0, 20]. The time stepsize is the same used for the solution of the Fisher equation and it is given by ht = 0.05. Several numerical investigations reveal that a critical value d∗ of the parameter d exists such that for 0.1 < d < d∗ ≈ 0.6, invasion travelling wave solutions exist, moving with a specific speed cnum (d) > c∗ (d). For d = d∗ they move with a speed cnum (d) ≈ c∗ (d). For d∗ < d < 1, we find travelling wave solutions approaching the steady state E2 with small damped oscillations and moving with a specific speed cnum (d) < c∗ (d). This suggests that, in the first range of the parameter d, 0.1 < d ≤ d∗ ≈ 0.6, travelling waves of invasion are stable and, on the contrary, in the range d∗ < d < 1, the pursuit and evasion travelling waves are the stable ones. Two examples of these two different scenarios are discussed below into details. As first test, we present the case d = 0.4 as representative of the range 0.1 < d < d∗ ≈ 0.6. For this value of d, the relationship (13) yields c∗ = 0.3040. The numerical results for this case are presented in Figs. 8a and 9a. In Fig. 8a we report the 3D numerical solution obtained by the D2ECDF of order p = 6 for h = 0.15. In Fig. 9a, we show the front propagation obtained by all the schemes, shown with same symbol as in the previous section. In the zoomed plot we note that: (i) the front has a smooth transition; (ii) the schemes show slightly different speeds. A further analysis of the speed approximation is given in Fig. 10. The same remarks reported for the approximation of the Fisher equation hold. In Table 2, we report the approximated values of the wave speeds in the middle of the front on the last time step. All schemes for decreasing stepsizes h converge towards cnum ≈ 1.373 > c∗ = 0.3040. Table 2 Front speed approximations C(r) (h) for all methods and for d = 0.4, d = 0.9. d = 0.4

h = 0.3 h = 0.15 h = 0.075

d = 0.9

pdepe

p=2

p=4

p=6

pdepe

p=2

p=4

p=6

1.3687 1.3719 1.3727

1.3600 1.3698 1.3722

1.3723 1.3730 1.3730

1.3729 1.3730 1.3730

1.0887 0.9122 0.9116

1.0245 0.9078 0.9104

0.8754 0.9112 0.9113

0.8896 0.9113 0.9113

1042

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

d=0.4, c ≈ 1.373

2.5

d=0.9, c ≈ 0.9113

2.5 time 0 ≤ t ≤ 20

2

1.5

1.5

η(ξ,t)

η(ξ,t)

time 0 ≤ t ≤ 20

2

1 2.005

1 2.005

0.5

0.5

2

2 1.995

0

0 1.99 38

−0.5 0

10

20

30

40

ξ

40

50

1.995 29

30

31

32

42

−0.5 0

60

10

20

30

40

ξ

50

60

Fig. 9. Time evolution of the travelling wave front η(ξ, τ ∗ ) for time 0 ≤ t ≤ 20 obtained by pdepe (dashed line) and D2ECDF of order p = 2 (red line), p = 4 (blue line), p = 6 (black line). Left plot: case d = 0.4. The wave profile approaches the steady state E2 monotonically with cnum ≈ 1.373 > c∗ = 0.3040. Right plot: case d = 0.9. The wave profile presents small damped oscillations around the steady state E2 . In this case cnum ≈ 0.9113 < c∗ = 3.7048. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

As second test, we consider d = 0.9, as representative of the case d∗ < d < 1. According to (13), c∗ = 3.7048. The numerical results are shown in the Figs. 8b and 9b. Also in this case, Fig. 8b shows the numerical solution in 3D obtained by the D2ECDF of order p = 6 for h = 0.15. In Fig. 9b, we show the front propagation obtained by all the schemes. In the zoomed plot we note that: (i) the front exhibits damped oscillations before attaining the stable equilibrium; (ii) the schemes show slightly different speeds. In Table 2, we report the approximated values of the wave speeds in the middle of the front on the last time step. All schemes for decreasing stepsizes h converge towards cnum ≈ 0.9113 < c∗ = 3.7048. p=4

p=2 1.376

1.375

h=0.15, 0.075.

h=0.075 1.374

1.37 1.372

c(η)

c(η)

h=0.15 1.365

1.37

h=0.3

h=0.3

1.368

1.36

1.366 1.355

0

0.5

1

η

1.5

2

0

0.5

1.376

1.374

1.374

1.372

1.372

1.37 1.368

1.5

2

1.5

2

PDEPE

1.376

c(η)

c(η)

p=6

1

η

h=0.075

h=0.15 1.37 1.368

h=0.3, 0.15, 0.075

1.366

h=0.3

1.366 0

0.5

1

η

1.5

2

0

0.5

1

η

Fig. 10. Case d = 0.4. Estimates of the numerical speeds of the front for η(ξ, τ ∗ ) obtained by all methods applied and for h = 0.3, 0.15, 0.075. Only estimates at the last ten time steps are shown. For decreasing h, the speeds in the middle of the profile converge towards cnum ≈ 1.373 > c∗ = 0.3040.

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

1043

Fig. 11. (a) Micrograph of the electrodeposited alloy exhibiting unstable growth. (b) Numerical solution of the PDE system (3) corresponding to the travelling wave for the electrodeposit profile η(ξ, τ ∗ ) for all time τ ∗ ∈ [0, 30]. At a fixed time ¯t , the solution along the section aa can be compared with the electrodeposits along any horizontal section of the alloy in Fig.11a.

Fig. 12. Macrograph of the electrodeposited alloy exhibiting unstable growth.

In our opinion, the above results are interesting because they are able to qualitatively describe some essential features related to electrodeposition experiments for Au–Cu alloys, as those shown in Figs. 11 and 12. In particular, the comparison reported in Fig. 11b shows that a travelling wave 1D in space, at a fixed time ¯t along the section aa, can be detected on experimental electrodeposits along any horizontal section of Fig. 11a. This can be an example of invasion travelling wave because of the smooth transition. In Fig. 12, a similar comparison may be performed. This can be considered an example of pursuit and evasion travelling wave, because along 1D vertical section in the middle of the picture, at the interface between Au–Cu electrodeposits, the lighter border may be related to a damped transition. Acknowledgments The authors are very grateful to the anonymous referee for his/her comments and suggestions. References [1] P. Amodio, I. Sgura, High order finite difference schemes for the solution of second order BVPs, J. Comput. Appl. Math. 176 (1) (2005) 59–76. [2] P. Amodio, I. Sgura, High order generalized upwind schemes and the numerical solution of singular perturbation problems, Bit Numer. Math. 47 (2007) 241–257. [3] R.D. Benguria, M.C. Depassier, M. Loss, Upper and lower bounds for the speed of pulled fronts with a cut-off, Eur. Phys. J. B 61 (2008) 331–334. [4] B. Bozzini, G. Giovannelli, P.L. Cavallotti, Investigation into the electrodeposition of Au–Cu-matrix particulate composites, J. Appl. Electrochem. 29 (1999) 685–692. [5] B. Bozzini, G. Giovannelli, P.L. Cavallotti, An investigation into the electrodeposition of AuCu-matrix particulate composites. Part II: Baths not containing free cyanide, J. Appl. Electrochem. 30 (2000) 591–594. [6] B. Bozzini, P.L. Cavallotti, G. Giovannelli, Electrokinetic behavior of gold alloy and composite plating baths, Met. Fin. 100 (2002) 50–60.

1044

B. Bozzini et al. / Mathematics and Computers in Simulation 81 (2011) 1027–1044

[7] B. Bozzini, I. Sgura, L. D’Urzo, Modelling of morphological control of electrodeposited Cu by adsorption of aminoacids and oligopeptides, in: A. El Nemr (Ed.), New Developments in Electrodeposition and Pitting Research, Research Signpost Pub., Fort P.O., Trivandrum-695 023, Kerala, India, 2007. [8] B. Bozzini, G. Giovannelli, C. Mele, Electrochemical dynamics and structure of the Ag/AgCl interface in chloride-containing aqueous solutions, Surf. Coat. Technol. 201 (2007) 4619–4627. [9] B. Bozzini, D. Lacitignola, I. Sgura, A reaction-diffusion model of spatial pattern formation in electrodeposition, J. Phys.: Conf. Ser. 96 (2008) 012051. [10] E. Brunet, B. Derrida, Shift in the velocity of a front due to a cutoff, Phys. Rev. E 56 (1997) 2597–2604. [11] W.A. Day, G. Saccomandi, On the propagation of the bulk of a mass subject to periodic convection and diffusion, Quart. Appl. Math. 57 (1999) 561–572. [12] S.R. Dunbar, Travelling wave solutions of diffusive Lotka-Volterra equations, J. Math. Biol. 17 (1983) 11–32. [13] R.A. Fisher, The wave of advance of advantageous genes, Ann. Eugen. 7 (1937) 355–369. [14] J. Guckenheimer, P. Holmes (Eds.), Nonlinear Oscillations Dynamical Systems and Bifurcations of Vector Fields, Springer-Verlag, NY, 1983. [15] J.L. Hudson, T.T. Tsotis, Electrochemical reaction dynamics: a review, Chem. Eng. Sci. 49 (1994) 1493–1572. [16] D.M. Kolb, Physical and electrochemical properties of metal monolayers on metallic substrates, in: H. Gerischer, C.W. Tobias (Eds.), Advances in Electrochemistry and Electrochemical Engineering, vol. 11, Wiley, New York, 1978, p. 125. [17] A. Kolmogoroff, I. Petrovsky, N. Piscounoff, Etude de l’équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique, Bulletin de l’Université d’Etat à Moscou, Serie Internationale, vol. I, 1937. [18] K. Krischer, Principles of temporal and spatial pattern formation in electrochemical systems, in: B.E. Conway, J.O. Bockris, R. White (Eds.), Modern Aspects of Electrochemistry, Plenum, New York, 1999. [19] M. Lucia, C.B. Muratov, M. Novaga, Linear vs. nonlinear selection for the propagation speed of the solutions of scalar reaction-diffusion equations invading an unstable equilibrium, Comm. Pure Appl. Math. 57 (2004) 616–636. [20] J.D. Murray, Mathematical Biology II; Spatial models and biomedical applications, 3rd Ed., Springer, 2002. [21] B. Sandstede, Stability of travelling waves, in: B. Fiedler (Ed.), Handbook of Dynamical Systems II, Elsevier, Amsterdam, 2002, pp. 983–1059. [22] L.T. Takahashi, N.A. Maidana, W. Ferreira, C. Pulino, P.H.M. Yang Jr., Mathematical models for the Aedes Aegypti dispersal dynamics: travelling waves by wing and wind, Bull. Math. Biol. 67 (2005) 509–528. [23] A.M. Turing, The chemical bases of morphogenesis, Phil. Trans. R. Soc. Lond. B 237 (1952) 37–72. [24] H. Wise, J. Oudar, Material Concepts in Surface Reactivity and Catalysis, Dover Publications Inc., Mineola NY, 1990.