Reactive-Diffusive-Advective Traveling Waves in a Family of Degenerate Nonlinear Equations.

This paper deals with the analysis of existence of traveling wave solutions (TWS) for a diffusion-degenerate (at D(0) = 0) and advection-degenerate (at h'(0) = 0) reaction-diffusion-advection (RDA) equation. Diffusion is a strictly increasing function and the reaction term generalizes the kinetic part of the Fisher-KPP equation. We consider different forms of the convection term h(u): (1)  h'(u) is constant k, (2)  h'(u) = ku with k > 0, and (3) it is a quite general form which guarantees the degeneracy in the advective term. In Case 1, we prove that the task can be reduced to that for the corresponding equation, where k = 0, and then previous results reported from the authors can be extended. For the other two cases, we use both analytical and numerical tools. The analysis we carried out is based on the restatement of searching TWS for the full RDA equation into a two-dimensional dynamical problem. This consists of searching for the conditions on the parameter values for which there exist heteroclinic trajectories of the ordinary differential equations (ODE) system in the traveling wave coordinates. Throughout the paper we obtain the dynamics by using tools coming from qualitative theory of ODE.

1. Introduction

The strong effect produced by the addition of the nonlinear convective term kuu on the solutions behavior of the classical Fisher-KPP equation u = u + u(1 − u) is well documented in the literature (see [1-3]). We mean those for the nonlinear reaction-diffusion-advection equation u = u − kuu + u(1 − u). This is particularly remarkable, when the diffusion is negligible compared to the convective effects. In such a case, the solutions can exhibit shock-like behavior (see [4-6]). In the above equation the term ku is called the advection “speed.” It has also been proved (see [3]) that the previous equation has monotonic decreasing traveling wave solutions (TWS) u(x, t) = ϕ(x − ct) ≡ ϕ(ξ), where c is the speed of the wave, satisfying the boundary conditions lim⁡ϕ(ξ) = 1, lim⁡ϕ(ξ) = 0 with 0 < ϕ(ξ) < 1, ∀ξ ∈ (−∞, +∞), if and only if c ≥ c(k), where

In [7], the author carries out a Painlevé analysis to get some approximate solutions for the above-mentioned equation.

In a series of papers authored by Malaguti et al. [8-11], the existence of TWS of the monostable (the reaction term has two equilibria, one is asymptotically stable and the other is unstable) reaction-diffusion-convection equation,was investigated. The constant diffusion case was studied in [12] and the nondegenerate case (D(u) > 0  ∀u ∈ (0,1)) in [11]; they proved that (2) admits decreasing TWS.

In [8] the authors looked at the case, where D(u) is such that D(0) = 0, D(1) > 0 (simply degeneracy) and D(1) = 0 (double degeneracy) and D′(0) = 0 and D′(1) = 0. Although they take h(u) to be nonlinear, specific properties of h(u) are not explicitly stated. In particular, their equation is not necessarily a convection-degenerate one. Note also that even though an application of their results to the evolution of a bacterial colony is presented, this equation does not contain convection term, which makes no real application of the convection-diffusion problem. In [10] continuous dependence of the threshold wave speed and of the traveling wave profiles is studied with respect to the diffusion, reaction, and convection terms. In [12] degenerate convection was considered. More recently the authors considered aggregation (i.e., the diffusion term changes sign) [9], D(u) > 0 for u ∈ (0, α), D(u) < 0, for u ∈ (α, 1), and a bistable term [13]. See [14] for a review.

Gilding and Kersner [15] and separately Mansour [16] looked at the particular case u + bu u = (au u ) + cu(1 − u ) arising in the study of pattern formation by bacterial colonies. Here a ≥ 0, b, c ≥ 0, and k > 0 are constants. Kamin and Rosenau [17] also looked at a similar special case, with D(u) = h(u) = u and g(u) = u(1 − u). In both papers, the set of wave speeds from which the equation admits a wavefront are studied.

The incorporation of a more general nonlinear advection term in the above equation also has been discussed in [3]. In such case, that equation takes the form u = u − [h(u)] + u(1 − u), where the convection “speed” is h′(u).

The aim of this paper is the investigation of the existence of TWS for the one-dimensional nonlinear degenerate RDA equationwhere the functions D, h, and g are defined on the interval [0,1] and there, they satisfy the following conditions:The degeneracy (at u = 0) could have two sources: the diffusion term D(0) = 0 and the case of item (2)(b), for the advective term, where h′(0) = 0.

D ∈ C [0,1] 2 with D(0) = 0, D(u) > 0  ∀u ∈ (0,1]; D′(u) > 0  ∀u ∈ [0,1] and D′′(u) ≠ 0  ∀u ∈ [0,1].

h ∈ C [0,1] 1 with h′(u) > 0  ∀u ∈ (0,1]. Two cases will be considered for h′(0). Namely,

h′(0) > 0,

h′(0) = 0.

g ∈ C [0,1] 2 with g(0) = g(1) = 0, g(u) > 0  ∀u ∈ (0,1); g′(0) > 0 and g′(1) < 0.

The inclusion of the first-order spatial derivative term h′(u)u in (3) transforms the parabolic degenerate nature of (3) with h′(u) ≡ 0 into a hyperbolic-like type. In fact, given that the partial derivative ∂F/∂u = D(u) of the nonlinear operatorvanishes at u = 0, F is not elliptic precisely at u = 0. Because of that, the nonlinear operator,is not parabolic at u = 0. See [18].

The degeneracy of the equation involves two important features of its solutions. One is the finite speed of propagation throughout the space. The other is that, for general rule, we do not expect that all the initial and boundary conditions problem associated with (3) possesses a classical solution, that is, smooth enough solution.

The TWS analysis we carried out through this paper uses a dynamical systems approach, which is different to that used by other authors [12, 15] and focuses on the qualitative behavior of the trajectories of a phase portrait as the involved parameters change. Additionally, in order to show the TWS whose existence we prove, we numerically solve the initial and boundary value problems associated with the full RDA in each considered case.

In ecological terms, (3) could describe the space-temporal dynamics of one species living in a one-dimensional habitat subject to the following factors: a density-dependent diffusion term D(u) which produces a pressure on the individuals of the population to migrate from crowded areas to sparse ones (for more details on this interpretation of D(u), see [19] and references therein), a nonlinear advective term h′(u) which “pushes” the population towards the direction −u (a sort of “directed wind”; see [3]), and a density-dependent growth rate g(u) which, by its qualitative features given in item (3), gives the dynamics of a habitat with limited resources (logistic growth). The carrying capacity of the habitat, in nondimensionalized form, is one.

The derivation of (3) can be done by using the microscopic individual behavior (random walks approach) which can be seen somewhere else (see [20] or [21]). Here we are omitting the details.

Note that taking h′(u) = V(u) + uV′(u), where A is an arbitrary constant, we obtain the equation mentioned and studied in [22].

Among the possible space-time patterns which could be described by (3), are those of traveling wave type, that is, solutions of the form u(x, t) = ϕ(x − ct), where c is the wave speed. These can be interpreted as waves of invasion of the population into the habitat.

Our analysis is based on the assumption that to look for TWS in a functional space is equivalent to search the set of parameters (in which the speed c is included) for which a two-dimensional system of ODE possesses heteroclinic trajectories. This system comes from the restatement of the original problem into the appropriate traveling wave variable.

The TWS analysis for (3) we present in this paper will be carried out in stages corresponding to the different levels of complexity the function h′ exhibits. These are the cases we consider.

Case 1 .

h′(u) = k.

Case 2 .

h′(u) = ku.

Case 3 .

No specific form for h′(u). This function must satisfy the quite basic requirements as stated in item (2).

These three cases are studied separately. Thus, the analysis of each one is the contents of the following three sections of this paper.

2. The TWS Analysis for h′(u) = k

With the choice h′(u) = k, the RDA equation (3) becomeswhere the density-dependent diffusion coefficient D and the kinetic part g satisfy the conditions listed in the previous section.

Note that because of the qualitative features of D and g on the interval [0,1], the pair of functions u 0(x, t) ≡ 0 and u 1(x, t) ≡ 1 are homogeneous and stationary solutions of (6) for all (x, t) ∈ ℝ × ℝ +. Their role in choosing the boundary condition for the TWS of (6) will be clear later.

2.1. A Quick Review on the k = 0 Case

For k = 0, (6) becomes

The traveling wave dynamics of (7) has been already studied (see [19]). In this reference, the authors proved that for D and g satisfying the conditions stated before, there exists a unique value, c > 0, of c, such that (7) has

no traveling wave solutions for 0 < c < c ,

a unique traveling wave solution of sharp type for c = c ,

a traveling wave solution of smooth decreasing front type satisfying the boundary conditions lim⁡ϕ(ξ) = 1 and lim⁡ϕ(ξ) = 0 with 0 < ϕ(ξ) < 1  ∀ξ ∈ (−∞, +∞), for each c > c .

The profile of the TWS mentioned in items (2) and (3) can be found in [19].

2.2. The Analysis for k ≠ 0

In this subsection, we investigate the existence of decreasing TWS for (6). The specific TWS satisfy the boundary conditions lim⁡ϕ(ξ) = 1 and lim⁡ϕ(ξ) = 0 with 0 < ϕ(ξ) < 1  ∀ξ ∈ (−∞, +∞).

Our analysis starts with the physical interpretation of the convective term −ku in (6). The diffusive substance is pushed out with speed k towards the direction of −u .

Suppose that is a traveling wave solution of (6) satisfying the appropriate boundary conditions. We then set , by substituting into (6); we obtain the second-order ODE:where the symbol  ′ on ϕ means derivative respect ξ 1 and on D means derivative with respect to ϕ. Note that, denoting by , actually (8) is of the formThis second-order ODE, except that the derivative  ′ is with respect to ξ 1, has exactly the same form as the corresponding second-order ODE in the traveling wave coordinate for the TWS for (7) (see [19]). The following argument justifies this simplification. Let T : ℝ 2 → ℝ 2 be a linear transformation such that for (x, t) ∈ ℝ × ℝ + The following proposition holds. Let w(x′, t) be such that u(x, t) = w(x′, t).

Proposition 1 .

If u(x, t) is solution of (6), then w(x′, t) satisfies

Proof

This follows by using the chain rule; we obtain u = −kw + w and u = w . Then we substitute u and u into (6) to arrive to the above equation.

Remark 2 .

Formally, the meaning of Proposition 1 is as follows: the convective effect in (6) can be suppressed by simply traveling in a space moving system of coordinates, which moves parallel to the x-axis with the convective speed k.

Hence, in the light of the previous reasoning the existence of TWS analysis for (6) is essentially reduced to the methodology developed in [19]. Therefore, by adapting and reinterpretation of the results given there, the following proposition holds in the present case.

Proposition 3 .

If the functions D and g satisfy the conditions (1) and (3) and ϕ(x, t) is solution of (7) on −∞ < x < +∞, t > 0, then for each k ≠ 0, ψ(x, t) ≡ ϕ(x − kt, t) is solution of (6).

This follows by using the hypothesis and the chain rule. In fact, by hypothesis for each (x, t) such that −∞ < x < +∞, t > 0Computing the relevant derivativesand using (8) and (9) appropriately, we obtain the equation for which ψ is solution, as stated in the proposition.

In Figure 1 we show the corresponding traveling wave profiles for two values of k. These were obtained by numerically solving (6) with the numerical routine NAG P03PCF (NAG Matlab Toolbox). As initial conditions we took 0.5(1.0 + tanh⁡((0.4 − x)/0.1)).

Figure 1

Front traveling waves solutions of (6) at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = 1. The values for k are indicated in each figure. The initial condition used was u(x, 0) = 0.5(1.0 + tanh⁡((0.4 − x)/0.1)), over the spatial domain is 0 < x < 140.

For k > 0, the speed of the TWS of (6) is faster than the TWS when k = 0. For negative k, the speed is slower than TWS when k = 0. In particular, this is true for the sharp type solution. See in Figure 2 an illustration of how the sharp type traveling wave evolves. The sharp traveling wave shown here (and in all remaining cases through the paper) was computed using the PDE numerical solver of Matlab (pdepe) using compact support initial conditions (in a larger domain −50 < x < 40, u(x, 0) = 1 for x < −20 and u(x, 0) = 0 for x ≥ 20). Note that the issue of accurately numerically computing a degenerate traveling wave represents on its own a research topic in numerical analysis; this is beyond the scope of this paper and we refer the interested reader to [23].

Figure 2

Evolution of a compact support function as initial condition under the reaction-diffusion convection process described by (6). The numerics show that as the time grows, the approximative solution tends to the sharp traveling wave solution. Note the agreement between the calculated speed of the sharp solution with that for which we have the saddle-saddle heteroclinic trajectory; see Figure 3. The profiles are plotted at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = k = 1, compact support initial conditions (in a larger domain −50 < x < 40, u(x, 0) = 1 for x < −20, and u(x, 0) = 0 for x ≥ 20).

We also computed numerically the speed of the sharp wave (see Figure 3). To compute this we follow how the wave's location changes with each time step; that is for a fixed value u ≡ 0.5 we compute x 1, x 2,… so that x is x, where u(x, t ) = 0.5. We plot this curve versus time and compute the slope to approximate the wave speed; we have obtained 1.8 as a result. We note that the corresponding sharp wave for k = 0 has a speed of c ≈ 0.8; that is, when k = 1, c = c + 1 = 1.8, making an illustrative example of Proposition 3. Furthermore, for these parameter values (β = k = 1), the phase portrait numerical results show that the value of c for which we have a unique heteroclinic connexion between P 1 and P (corresponding to the sharp traveling wave) is c ≈ 1.82.

Figure 3

Approximation of the wave speed of the profiles evolving to a sharp traveling wave; to compute this we follow how the wave's location changes with each time step. For a fixed value u ≡ 0.5 we compute x 1, x 2,… so that x is the x, where u(x, t ) = 0.5. We plot this curve versus time and compute the slope to approximate the wave speed and obtain 1.8 as a result. For these parameter values (β = k = 1), the phase portrait numerical results show that the value of c for which we have a unique heteroclinic connexion between P 1 and P is for c ≈ 1.82.

3. Traveling Wave Solutions Analysis with h′(u) = ku

Here, (3) takes the formwhere the functions D and g satisfy conditions (1) and (3) of Section 1. As in the equation analyzed in Section 2, the functions u 0(x, t) ≡ 0 and u 1(x, t) ≡ 1 are homogeneous and stationary solutions of (14) for all (x, t) ∈ ℝ × ℝ +.

We will impose the same conditions as in Section 2 for the possible TWS of (14). The TWS analysis we are going to do starts by assuming the existence of positive c and ϕ such that u(x, t) = ϕ(x − ct) ≡ ϕ(ξ) is solution of (14). Hence, by substituting u(x, t) = ϕ(x − ct) into (14), we get the nonlinear second-order ODE:which, by setting v = ϕ′, can be written as the singular (at ϕ = 0) two-dimensional ODE system:The singularity can be removed by using a standard reparametrization of (16) (see [19]). Thus, by introducing the new parameter τ such that for ϕ > 0,we obtain the nonsingular system:where the dot on ϕ and v denotes the derivative of these variables with respect to τ. This new system is topologically equivalent to (16) in the region

In this region, for positive c system (18) has three equilibria: P 0 = (0,0), P 1 = (1,0), and P = (0, −c/D′(0)). Hence, according to the conditions of interest, the problem of showing the existence of the TWS satisfying such conditions in the full nonlinear PDE transforms into a dynamical systems problem. This is searching for the existence of the parameter values for which there exist heteroclinic trajectories of (18) connecting P 1 with P 0 or with P . The analysis is conducted by stages.

3.1. Local Dynamics

First, let us determine the local behavior of the trajectories of (18) in a neighborhood of each equilibria. For this aim, we obtain the Jacobian matrix of the vector field defined in (18) at any point (ϕ, v). This isThe evaluation of (20) at P 0 gives usfrom which we have tr⁡J[f 1, f 2](0,0) = −c < 0 for all positive values of c and det⁡J[f 1, f 2](0,0) = 0. Given that the eigenvalues of (21) are λ 1 = 0 and λ 2 = −c, then P 0 is a nonhyperbolic point of codimension one (see [24]). The corresponding eigenvectors are v 1 = (−c/g′(0), 1) and v 2 = (0,1).

Because the Hartman-Grobman Theorem is not applicable here, the local dynamics of (18) around P 0 does not follow from the corresponding linear approximation. In such a case, we must use the higher order terms in the Taylor series of the vector field (f 1, f 2) around P 0. In fact, we should obtain the normal form of (18) and then use the Center Manifold Theorem (see [24]). This tells us the local dynamics of (18) around P 0 can be essentially reduced to that around its center manifold. By proceeding as we already mentioned it, we conclude that P 0 is a saddle-node point (see [19] for a similar analysis).

Evaluating (20) at P 1 we obtain the Jacobian matrixFrom here, tr⁡J[f 1, f 2](1,0) = −(c − k) and det⁡J[f 1, f 2](1,0) = g′(1)D(1) < 0; therefore P 1 is a hyperbolic saddle point for all positive values of c and k. The corresponding eigenvalues and the eigenvectors arerespectively.

At P , (20) reduces toThen, it follows tr⁡J[f 1, f 2](0,− = 0 and det⁡J[f 1, f 2](0,− = −c 2 < 0; therefore P is a hyperbolic saddle point for all c ≠ 0 and k. The eigenvalues and eigenvectors are λ 1 = −c, λ 2 = c; v 1 = (2cD′(0)2/(D′′(0)c 2 + kD′(0)c + D′(0)2 g′(0)), 1) and v 2 = (0,1), respectively.

3.2. The Nullclines: Towards the Global Dynamics

In order to study the global dynamics associated with system (18), we should understand how its nullclines behave as the involved parameters change. The horizontal nullclines of system (18) are the horizontal and the vertical axis of the ϕv plane. The vertical nullcline has these two branchesFrom its respective expression, it follows Note that for all positive c, V 2(0) < 0; meanwhile, given the positiveness of D′(1), the sign of V 2(1) changes according with the sign of the term −(c − k).

3.3. Dynamics for Extreme Values of c

Here we are going to analyze the dynamics of (18) by considering extreme (including c = 0) values of c. This is done by considering two separate cases.

3.3.1. For c = 0

For c = 0 system (18) becomeswhose equilibrium points (in the region of interest) are P 0 = (0,0) and P 1 = (1,0). Here P 0 comes from the collapse of P into the origin and this point becomes a nonhyperbolic equilibrium of codimension two; meanwhile P 1 stills as a hyperbolic saddle point.

The vertical nullcline branches of system (27) are Figure 4 contains a panel of figures illustrating the behavior of (28) in a representative case, where D(u) = u + βu 2 and g(u) = u(1 − u) for fixed positive β and different values of k.

Figure 4

The branches of the vertical nullcline of system (27) for c = 0 and several values of k. Here β = 2.

As it can be seen in Figure 4, for small values of k, V 1 and V 2 (given by (28)) are not defined on the whole interval [0,1]. They are, however, well defined on this interval for big enough values of k and fixed positive β. Moreover, the behavior of (28)—according to Figure 4—can be classified in three main categories. These are illustrated in each row of the mentioned figures.

For each positive k, let us introduce the following notation:The following proposition holds.

Proposition 4 .

For each positive k system (27) does not have closed trajectories on the following sets:

{(ϕ, v)∣0 ≤ ϕ ≤ 1, −∞ < v < 0},

{(ϕ, v)∣0 ≤ ϕ ≤ 1, v < v (ϕ)},

{(ϕ, v)∣0 ≤ ϕ ≤ 1, v > v (ϕ)}.

This follows by a straightforward application of Bendixson's Negative Criterion. In fact, the divergence, div, of the vector field which defines system (27) isGiven D is a strictly increasing function on [0,1], for item (1)   is positive, the same sign for item (2); meanwhile in the third case. Then the proof follows.

As a consequence of Proposition 4 and the Poincaré-Bendixson Theorem on each set this proposition states the ω and α limit sets of the trajectories are equilibrium points.

In what follows we are going to use the behavior of both branches of the vertical nullcline for the determination of the phase portrait of system (27).

Proposition 5 .

For small enough positive values—including c = 0—of c system (18) does not have nontrivial heteroclinic trajectory in the strip {(ϕ, v)∣0 ≤ ϕ ≤ 1, −∞ < v < +∞}.

By trivial heteroclinic trajectories, we mean those of (18) for which ϕ = 0  ∀τ ∈ (−∞, ∞), as it is the case for each trajectory running on the negative vertical axis of the phase portrait, connecting P with P 0 which exists for all c > 0. Because of the physical interpretation of ϕ, we are not interested in those.

Phase portraits of system (27) for c = 0 and for the same positive values of k as those in Figure 4 can be found in Figure 5.

Figure 5

Phase portraits of system (27) for c = 0 and for the same positive values of k as those in Figure 4. Red points are equilibrium points. Green curves are nullclines and blue curves are trajectories.

3.3.2. For c > 0

Let us introduce the following notation:where the maximum is taken on the closed interval [0,1]. Through this subsection we are going to distinguish two main cases:

0 < c < M ,

c ≥ M .

For values of c satisfying c ≥ M , V 1 and V 2 are defined1 for all ϕ ≥ 0; in particular they do so on the interval [0,1].

In another side, for values of c such that 0 < c < M these branches of the vertical nullcline are not defined on the whole interval [0,1] but they do so on the union of subintervals contained within it. Figure 6 illustrates the behavior of (25) in the same representative2 case as in previous subsection for fixed positive β and k. Here, the positive c varies.

Figure 6

Behaviour of the two branches of the vertical nullclines of system (18) for fixed k = 2 and c > 0: (a)–(g) for 0 < c < M , (h) for c = M , and (i) for c > M . Here β = 2.

Now, we use this behavior to determine the phase portrait of (18). We proceed by considering extreme values of c. We can prove the following proposition.

Proposition 6 .

For each value of c such that c ≥ M system (18) has a heteroclinic trajectory (ϕ (τ), v (τ)), connecting the equilibria P 1 and P 0, that is, satisfying

0 < ϕ (τ) < 1 and v (τ) < 0 for all τ ∈ (−∞, +∞),

lim⁡(ϕ (τ)  , v (τ)) = (1,0) and lim⁡(ϕ (τ), v (τ)) = (0,0).

For c ≥ M the vertical null-clines look like in Figures 6(h) or 6(i). On each one of these branches the vector field, being horizontal, points out towards the left; meanwhile on the region {(ϕ, v)∣0 < ϕ < 1, V 2(ϕ) < v < V 1(ϕ)}, the vector field points left up. This behavior allows us to construct a positive invariant region for such a vector field in a similar fashion as that carried out in [25]. In fact, we can select a function f : [0,1] → ℝ belonging to the set C [0,1] 1 satisfying (a) f(0) = v 0 < 0, f(1) = 0; (b) f′(ϕ) > 0  ∀ϕ ∈ (0,1); and V 2(ϕ) < f(ϕ) < V 1(ϕ) in such a way the restriction of the vector field (18) on the graph of f points inwards the region {(ϕ, v)∣0 < ϕ < 1, f(ϕ) < v < 0}. Then, by a straightforward application of the Poincaré-Bendixson Theorem we have that any trajectory of (18)—in particular that leaving P 1 following the left branch of the unstable manifold at P 1—entering this region, must end at one equilibrium point. Given that there is not any other possibility such trajectory must end at P 0. Hence the proof follows.

The phase portrait of (18) for fixed k = 2 and c > 0 can be seen in Figure 7.

Figure 7

Phase portrait of (18) for fixed k = 2 and c > 0. In (a)–(c) there are no heteroclinic trajectories. (d) shows the existence of a saddle-saddle heteroclinic trajectory for c ≈ 1.63. (e)–(g), for the indicated value of c, show a saddle (P 1) saddle-node (P 0) connection. Here β = 2. Red points are equilibrium points, green curves are nullclines, and blue curves are trajectories.

3.4. A Monotonicity Property

Let (ϕ, v) be any fixed (but arbitrary) point belonging to the regionFor each pair (c, k), let us denote by θ(ϕ, v; c, k) the angle formed by the vector field (18) with the positive semihorizontal ϕ-axis. Then

Proposition 7 .

The angle θ(ϕ, v; c, k) is a monotone decreasing function of the parameter c and a monotone increasing function with respect to the parameter k.

Calculating the partial derivative with respect to c in the equalitywe obtainmeanwhile the corresponding partial derivative with respect to k isThen the proof follows.

Both monotonicity properties contained in the above proposition have important implications in refining the analysis for searching heteroclinic trajectories of system (18) connecting P 1 with P or with P 0. In particular, for fixed k > 0, if we continuously decrease the parameter c starting from M , the left branch of the unstable manifold at P 1, W (P 1), will move continuously downwards within the region ℱ 1; meanwhile, the right branch of the stable manifold at P , W (P ), moves continuously upwards as c decreases in the same region. By continuity of the vector field with respect to the parameter c and using shooting arguments, there exists a unique, c > 0, value of c for which both manifolds touch each other resulting in a saddle-saddle heteroclinic trajectory connecting the equilibrium points P 1 with P . This reasoning constitutes the proof of the following lemma.

Lemma 8 .

For the functions D(u) = u + βu 2 and g(u) = u(1 − u) with β > 0 and for h′(u) = ku with k > 0 there exists a unique c —depending on k—positive value of c for which system (18) has a unique heteroclinic trajectory connecting the equilibria P 1 and P . Moreover

By increasing any of the parameters c or k, such a trajectory is destroyed and for each (c, k) with either c > c or k > k a heteroclinic trajectory connecting P 1 with P 0 emerges.

On the contrary, by decreasing any of these parameters there are not heteroclinic trajectories for system (18) at all.

This follows from Propositions 4–7.

Note that there are two reasons why system (18) is not structurally stable (see Peixoto's theorem in [24]). These are the existence of a nonhyperbolic point and a saddle-saddle heteroclinic trajectory. In particular, any small perturbation of such system, for example, by varying the parameter c in a small neighborhood, V (c ), of the critical value c , involves strong dynamical changes including the destruction of the saddle-saddle trajectory and the emergence of a saddle (P 1) saddle-node (P 0) connexion or the disappearance of heteroclinic trajectories at all.

Theorem 9 .

For the functions D(u) = u + βu 2 and g(u) = u(1 − u) with β > 0 and for h′(u) = ku with k > 0, given k there exists a unique critical value, c (depending on k), of c such that (14) has

no traveling wave solutions for c such that 0 < c < c ,

a unique traveling wave solution of sharp type for c = c ,

a monotonic decreasing traveling wave solution for each c such that c > c .

The previous analysis demonstrates the existence of the associated heteroclinic trajectories (see Lemma 8). Searching for travelling wave solutions of the PDE (3), given that such solutions have the particular form u(x, t) = ϕ(x − ct) ≡ ϕ(ξ), is equivalent to showing the existence of heteroclinic trajectories of the associated ODE system (18). Therefore the theorem follows.

In Figure 8 we show two front traveling wave profiles for two values of k. In Figure 9 the sharp traveling wave can be seen.

Figure 8

Front traveling waves solutions of (14) at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = 1. The values for k are indicated in each figure.

Figure 9

Evolution of a compact support function as initial condition under the reaction-diffusion convection process described by (14). The numerics show that as the time grows, the approximative solution tends to the sharp traveling wave solution. Profiles are plotted at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = 1 and h′(u) = ku, with k = 1.

We close this section by numerically exploring the influence of changes in k on the critical value, c , of c for which the r-d-a equation has a sharp type solution. This is the content of the next subsection.

3.5. The Speed c Depending on k in a Particular Case

As we already mentioned in Section 1, the equation u = u − kuu + u(1 − u) has a monotone decreasing TWS connecting the states u 1(x, t) ≡ 1 and u 0(x, t) ≡ 0 for each c ≥ c(k), where the explicit form of c(k) is given by (1). This result by Murray was our motivation for seeking the corresponding relationship between the speed for which our reaction-diffusion-convection equation (14) has TWS and the parameter k, in particular for those of sharp type.

In this subsection we illustrate this relationship through a particular case. To this aim, we choose D(u) = u + βu 2 and g(u) = u(1 − u) with β > 0. Through this subsection our approach is from a numerical point of view. We carried numerical simulations of the phase portrait of the corresponding nonsingular ODE system in the traveling wave coordinate. The goal of this was to illustrate, for different values of k, the corresponding critical values, c , of c for which a unique saddle-saddle heteroclinic trajectory exists (one for each k). As we already know, associated with this trajectory, there exists a unique TWS of sharp type for (6). In Figure 10(a) we present a numerical approximation of how c depends on k through the corresponding phase portrait, with D(u) and g(u) as before.

Figure 10

“Numerical experimental data” for c for different values of k with h′(u) = ku. The data as they were obtained (a); the fitting (b).

This information tells us that, on this range of the “numerical experiments,” the speed c is a growing function of k. Moreover, we can distinguish two qualitative parts: one exponential for k ≤ 8 and another linear for k ≥ 8. We carried out the corresponding fittings. These are our results for each phase: respectively. See Figure 10(b).

In Figure 11 we show how c changes as function of the parameter k for h′(u) = k and h′(u) = ku for comparison.

Figure 11

Plot c as function of the parameter k: for h′(u) = k and h′(u) = ku.

As result of the numerical experiments, we can see that for k < 0 the critical value c is small but positive, but once k increases, the values of c increase faster with k. The physical interpretation of this is as follows: for k < 0 the “wind” h′(u) = k acts in the opposite direction to which the wave travels; meanwhile for k > 0, the advective term pushes in the same direction as the traveling wave goes. As a result of this, the critical values of the speed, for which there exists the sharp wave, increase as k increases.

4. Traveling Wave Solutions Analysis in the General RDA Equation

The particular cases discussed in previous sections give us some insights in order to carry out the TWS analysis for the general RDA equation:where D, g, and h are real functions defined on the interval [0,1]. The first two functions satisfy conditions (1) and (3) stated at the beginning of Section 1; meanwhile h, in addition to satisfying the conditions in item (2), is such that h′(0) might have different signs. In particular when h′(0) = 0 (39) is degenerate in both the diffusion and the advective terms.

For the analysis of TWS we proceed in a standard way: let us assume u(x, t) = ϕ(x − ct) ≡ ϕ(ξ) with c > 0 is solution of the RDA equation (39). Thus, by substituting ϕ in (39) we obtain a nonlinear second-order ODE equation for ϕ which, by introducing v = ϕ′(ξ), can written as a singular (at ϕ = 0) nonlinear ODE system. The singularity can be removed by introducing the parameter τ in a similar fashion as we did in previous sections (see [19, 26]). The result is the following nonsingular and nonlinear ODE system: where the dot on ϕ and v denotes the derivative with respect to τ. This system and the singular system are topologically equivalent on the stripe

The analysis of the system ((40a) and (40b)) starts by obtaining its nullclines. The horizontal nullcline is the coordinate axis of the ϕv plane. The vertical nullcline has the following two branches:Given the conditions on D, h, and g, the functions V 1 and V 2 are defined on the whole interval [0,1] whenever the inequality,holds which, in turn, gives us a bound for c for which the two branches (42) are defined on [0,1]. This isWhen the above inequality does not hold, V 1 and V 2 are not defined on the whole interval [0,1]. In fact, they are defined on disjoint intervals belonging to the interval [0,1].

From the explicit form of V 1 and V 2 it follows V 1(0) = 0, V 1(1) = 0, andSince we assumed both D′(0) and D′(1) are positive, depending on c compared with h′(0) (or with h′(1)), the following cases might occur:

V 2(0) < 0 for c > h′(0), V 2(1) < 0 for c > h′(1),

V 2(0) = 0 for c = h′(0), V 2(1) = 0 for c = h′(1),

V 2(0) > 0 for c < h′(0), V 2(1) > 0 for c < h′(1).

The equilibrium points of (40a) and (40b) areGiven the positiveness of D′(0), depending on the sign of (h′(0) − c) the third equilibrium is located on the

positive vertical v-semiaxis, for c < h′(0),

the origin for c = h′(0),

negative vertical v-semiaxis, for c > h′(0).

4.1. Local Dynamics

The linear local analysis of the system ((40a) and (40b)) is based on the Jacobian matrix, J[f 1, f 2], of the vector field (f 1, f 2) evaluated at the equilibria. The Jacobian matrix at any point (ϕ, v) is

At P 0, det⁡J[f 1, f 2](0,0) = 0 and tr⁡J[f 1, f 2](0,0) = −(c − h′(0)). Hence, whenever c ≠ h′(0), P 0 is a nonhyperbolic point of codimension one. Evaluating J[f 1, f 2]( at P 1 we obtain det⁡J[f 1, f 2](1,0) = g′(1)D(1) < 0 and tr⁡J[f 1, f 2](1,0) = −(c − h′(1)); hence P 1 is a hyperbolic saddle point. Finally, the evaluation of J[f 1, f 2] at P ( gives us tr⁡J[f 1, f 2] = 0 and det⁡J[f 1, f 2] = −(c − h′(0))2 from which it follows that for c ≠ h′(0), det⁡J[f 1, f 2] < 0. Thus, whenever c ≠ h′(0), P ( is a hyperbolic saddle. On the contrary, for c = h′(0) the equilibrium P ( is a nonhyperbolic point of codimension two.

The Jacobian at P 0 iswith eigenvalues λ 1(P 0) = 0, λ 2(P 0) = h′(0) − c and eigenvectors v 1 = (−(c − h′(0))/g′(0), 1), v 2 = (0,1), respectively.

The Jacobian at P 1 iswith eigenvalues ,   and eigenvectors v 1 and v 2 given byrespectively.

The Jacobian at P ( iswith eigenvalues λ 1(P () = h′(0) − c, λ 2(P () = c − h′(0) and eigenvectors respectively.

4.2. Degeneracy Just in the Diffusion Term at u = 0

Here D(0) = 0 and h′(0) ≠ 0. Thus, for a fixed k and positive c—as we assumed—depending on the sign of h′(0) its effect on P ( is as follows:

For h′(0) > 0 we have two subcases:

0 < h′(0) < c; here P ( is closer to P 0 than P does; this is on the negative vertical v-axis.

0 < c < h′(0); P ( switches to the positive vertical axis of the ϕv plane.

For h′(0) < 0 the equilibrium runs away on the vertical negative axis. Here P (, for all c > 0, is more far away from P 0 on the vertical negative axis.

4.2.1. h′(0) > 0

We will take h′(u) = α + γu (α > 0) as an example of the first case. The behavior of the phase portrait depends on how h′(0) compares with c. In this example this means how α compares to c. For our example h′(0) = α = 1, we show in Figure 12 the phase portrait for several values of c. In particular, as stated before, we have two major subcases: 0 < c = 1 < h′(0) and 0 < h′(0) = 1 < c. For completeness, we start by showing the phase portrait for c = 0. In Figures 12(a) and 12(b) there are three equilibrium points and a heteroclinic connexion from P 0 to P 1, in Figure 12(c) the heteroclinic connexion is from P to P 1, in Figure 12(d) there are only trivial heteroclinic connexions, in Figure 12(e) there are only two equilibria, in Figure 12(f) there are only trivial connexions, Figure 12(g) corresponds to c , that is, one connexion from P 1 to P , and in Figure 12(h) for each c > c there is connexion from P 1 to P 0.

Figure 12

Phase portrait of (40a) and (40b) for eight values of c; (a)–(h) contemplate the first subcase (1 = h′(0) > c); (e) shows c = 1 (only two equilibria); and (f)–(h) contain the second subcase (c > h′(0) = 1). Here k = γ = α = 1, β = 2. Red points are equilibrium points. Green curves are nullclines and blue curves are trajectories.

In Figure 13 we show the corresponding front traveling wave profiles for two values of k. These were obtained by numerically solving (39) for h′(u) = α + γu with the numerical routine NAG P03PCF (NAG Matlab Toolbox). Figure 14 shows the sharp traveling wave. We also plot how c changes with k in this case (see Figure 16). Although we are concentrating on heteroclinic connections leaving P 1, we display the phase portrait behavior for a broad range of values for k to show the richness of the dynamics. To this end, we would also like to show how the heteroclinic connexion from P 0 to P 1 from Figure 12(b) translates into a front traveling from right to left (the initial condition in that case is marked with dots) connecting 0 with 1; see Figure 15.

Figure 13

Front traveling waves solutions of (39) for h′(u) = α + γu at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = 1. The values for k are indicated in each figure.

Figure 14

Evolution of a compact support function as initial condition under the reaction-diffusion convection process described by (39). The numerics show that as the time grows, the approximative solution tends to the sharp traveling wave solution. Profiles are plotted for h′(u) = α + γu at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = k = 1.

Figure 16

Plot of how c changes with k for h′(u) = α + γu. For a comparison, refer to Figure 10, where h′(u) = ku.

Figure 15

Front traveling waves solutions of (39) for h′(u) = α + γu at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = k = 1.

4.2.2. h′(0) < 0

This case is basically the same as before. We have two main subcases according to the sign of c − h′(0). Figure 17 shows how c changes with k for h′(u) = α + γu. We do not plot the corresponding phase portraits or traveling waves as they are qualitatively similar to the ones in the previous subsection. We do, however, show how the advection parameter changes with k in Figure 17; note how here also there is linear dependency to later become an exponential one. For a comparison, refer to Figure 10.

Figure 17

Plot of how c changes with k for h′(u) = α + γu, with α = −1 < 0. For a comparison, refer to Figure 10, where h′(u) = ku.

4.3. Degeneracy in Both Diffusion and Advection Terms at u = 0

Here D(0) = h′(0) = 0. Then the equilibrium P ( coincides with the equilibrium, P , which is allocated on the seminegative vertical v-axis for the corresponding system when no advection term is included. Here we will take h′(u) = u + ku 2 as an example.

Mutatis mutandis one of the monotonicity properties contained in Proposition 7 holds for the general case. In fact, by using the same notation as that used there, the following proposition holds.

Proposition 10 .

For fixed (but arbitrary) (ϕ, v) ∈ ℱ 1 the angle, θ(ϕ, v; c), formed by the positive ϕ semihorizontal axis and the vector field which defines the system ((40a) and (40b)), is a monotone decreasing function of the parameter c.

We haveand then, by simply calculating the derivative with respect to c in the above equality, we obtainfor all (ϕ, v) ∈ ℱ 1; then the proof follows.

Of course, all the implications this proposition has are assumed by the corresponding ODE system which in turn implies the existence of sharp and monotone traveling wave solutions.

In Figure 18 we show the corresponding front traveling wave profiles for two values of k for h′(u) = u + ku 2 and in Figure 19 the sharp type.

Figure 18

Front traveling waves solutions of (39) for h′(u) = u + ku 2 at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = 1. The values for k are indicated in each figure.

Figure 19

Evolution of a compact support function as initial condition under the reaction-diffusion convection process described by (39). The numerics show that as the time grows, the approximative solution tends to the sharp traveling wave solution. Profiles are for h′(u) = u + ku 2 at equal time intervals t = 1,2, 3,… with g(u) = u(1 − u) and D(u) = u(1.0 + βu) for β = 1.

Figure 20 shows how c changes with k for h′(u) = u + ku 2.

Figure 20

Plot of how c changes with k for h′(u) = u + ku 2. For a comparison, refer to Figure 10, where h′(u) = ku.

5. Discussion

We have presented a study on the effect of the incorporation of a nonlinear advection term in the degenerate reaction-diffusion equation (3). When investigating the traveling wave behavior, we found that the “advection speed” influences the type and the speed of the possible traveling waves. The aim of this paper was the investigation of the existence of TWS for the one-dimensional nonlinear degenerate RDA equation (3). The degeneracy of the equation causes its solution to possess finite speed of propagation throughout the space.

The TWS analysis for (3) was carried out for three cases: (1)  h′(u) = k, (2)  h′(u) = ku, and (3) no specific form for h′(u), as long as it satisfies the quite basic requirements as stated in item (2). In all these cases, (3) has

no traveling wave solutions for 0 < c < c ,

a unique traveling wave solution of sharp type for c = c ,

a traveling wave solution of smooth decreasing front type satisfying the boundary conditions lim⁡ϕ(ξ) = 1 and lim⁡ϕ(ξ) = 0 with 0 < ϕ(ξ) < 1  ∀ξ ∈ (−∞, +∞), for each c > c .

We also numerically solved the initial and boundary value problems associated with the full RDA in each considered case. Moreover, we show how the advection speed impacts the type and speed of the traveling waves. In particular, the unique c = c for which there are sharp traveling waves depends initially linearly on the advection speed to later increase exponentially.