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.


Introduction
The strong effect produced by the addition of the nonlinear convective term on the solutions behavior of the classical Fisher-KPP equation = + (1 − ) is well documented in the literature (see [1][2][3]). We mean those for the nonlinear reaction-diffusion-advection equation = − + (1 − ). 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][5][6]). In the above equation the term is called the advection "speed." It has also been proved (see [3]) that the previous equation has monotonic decreasing  (1) 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][9][10][11], the existence of TWS of the monostable (the reaction term 2 The Scientific World Journal 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 ( ( ) > 0 ∀ ∈ (0, 1)) in [11]; they proved that (2) admits decreasing TWS.
In [8] the authors looked at the case, where ( ) is such that (0) = 0, (1) > 0 (simply degeneracy) and (1) = 0 (double degeneracy) and (0) = 0 and (1) = 0. Although they take ℎ( ) to be nonlinear, specific properties of ℎ( ) 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], ( ) > 0 for ∈ (0, ), ( ) < 0, for ∈ ( , 1), and a bistable term [13]. See [14] for a review.
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 = − [ℎ( )] + (1 − ), where the convection "speed" is ℎ ( ).
The inclusion of the first-order spatial derivative term ℎ ( ) in (3) transforms the parabolic degenerate nature of (3) with ℎ ( ) ≡ 0 into a hyperbolic-like type. In fact, given that the partial derivative / = ( ) of the nonlinear operator ( , , ) ≡ ( ) + ( ) 2 − ℎ ( ) vanishes at = 0, is not elliptic precisely at = 0. Because of that, the nonlinear operator, is not parabolic at = 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 ( ) 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 ( ), see [19] and references therein), a nonlinear advective term ℎ ( ) which "pushes" the population towards the direction − (a sort of "directed wind"; see [3]), and a density-dependent growth rate ( ) 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.
Among the possible space-time patterns which could be described by (3), are those of traveling wave type, that is, solutions of the form ( , ) = ( − ), where is the wave speed. These can be interpreted as waves of invasion of the population into the habitat.
The Scientific World Journal 3 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 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 ℎ exhibits. These are the cases we consider.
Case 3. No specific form for ℎ ( ). 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.

The TWS Analysis for ℎ ( )=
With the choice ℎ ( ) = , the RDA equation (3) where the density-dependent diffusion coefficient and the kinetic part satisfy the conditions listed in the previous section.
Note that because of the qualitative features of and on the interval [0, 1], the pair of functions 0 ( , ) ≡ 0 and 1 ( , ) ≡ 1 are homogeneous and stationary solutions of (6) for all ( , ) ∈ R × R + . Their role in choosing the boundary condition for the TWS of (6) will be clear later.
The profile of the TWS mentioned in items (2) and (3) can be found in [19].
Our analysis starts with the physical interpretation of the convective term − in (6). The diffusive substance is pushed out with speed towards the direction of − .
Proof. This follows by using the chain rule; we obtain = − + and = . Then we substitute and 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 -axis with the convective speed .
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.   Proof. This follows by using the hypothesis and the chain rule. In fact, by hypothesis for each ( , ) such that −∞ < < +∞, > 0 Computing the relevant derivatives and 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 . 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 − )/0.1)).
For > 0, the speed of the TWS of (6) is faster than the TWS when = 0. For negative , the speed is slower than TWS when = 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 < < 40, ( , 0) = 1 for < −20 and ( , 0) = 0 for ≥ 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].
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 : 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 value ≡ 0.5 we compute 1 , 2 , . . . so that is , where ( , ) = 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 = 0 has a speed of * =0 ≈ 0.8; that is, when = 1, * =1 = * =0 +1 = 1.8, making an illustrative example of Proposition 3. which we have a unique heteroclinic connexion between 1 and (corresponding to the sharp traveling wave) is * ≈ 1.82.

Traveling Wave Solutions Analysis with ℎ ( )=
Here, (3) takes the form where the functions and satisfy conditions (1) and (3) of Section 1. As in the equation analyzed in Section 2, the functions 0 ( , ) ≡ 0 and 1 ( , ) ≡ 1 are homogeneous and stationary solutions of (14) for all ( , ) ∈ R × R + . 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 and such that ( , ) = ( − ) ≡ ( ) is solution of (14). Hence, by substituting ( , ) = ( − ) 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 system (18) has three equilibria: 0 = (0, 0), 1 = (1, 0), and = (0, − / (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 1 with 0 or with . The analysis is conducted by stages.
Because the Hartman-Grobman Theorem is not applicable here, the local dynamics of (18) around 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 6 The Scientific World Journal the vector field ( 1 , 2 ) around 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 0 can be essentially reduced to that around its center manifold. By proceeding as we already mentioned it, we conclude that 0 is a saddle-node point (see [19] for a similar analysis).

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 branches From its respective expression, it follows Note that for all positive , 2 (0) < 0; meanwhile, given the positiveness of (1), the sign of 2 (1) changes according with the sign of the term −( − ).

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

For
whose equilibrium points (in the region of interest) are 0 = (0, 0) and 1 = (1, 0). Here 0 comes from the collapse of into the origin and this point becomes a nonhyperbolic equilibrium of codimension two; meanwhile 1 stills as a hyperbolic saddle point.
The vertical nullcline branches of system (27) are  For each positive , let us introduce the following notation: The following proposition holds.
Proof. This follows by a straightforward application of Bendixson's Negative Criterion. In fact, the divergence, div, of the vector field which defines system (27) is Given is a strictly increasing function on [0, 1], for item (1) div ⃗ ( , V) is positive, the same sign for item (2); meanwhile div ⃗ ( , V) < 0 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).
Proof. 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 with 0 which exists for all > 0. Because of the physical interpretation of , we are not interested in those. Phase portraits of system (27) for = 0 and for the same positive values of as those in Figure 4 can be found in Figure 5.

3.3.2.
For > 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: (1) 0 < < , For values of satisfying ≥ , 1 and 2 are defined 1 for all ≥ 0; in particular they do so on the interval [0, 1].
In another side, for values of such that 0 < < 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 representative 2 case as in previous subsection for fixed positive and . Here, the positive varies. Now, we use this behavior to determine the phase portrait of (18). We proceed by considering extreme values of . We can prove the following proposition. Proof. For ≥ 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, 2 ( ) < 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 : [0, 1] → R belonging to the set 1 [0,1] satisfying (a) (0) = V 0 < 0, (1) = 0; (b) ( ) > 0 ∀ ∈ (0, 1); and 2 ( ) < ( ) < 1 ( ) in such a way the restriction of the vector field (18) on the graph of points inwards the region {( , V) | 0 < < 1, ( ) < V < 0}. Then, by a straightforward application of the Poincaré-Bendixson Theorem we have that any trajectory of (18)-in particular that leaving 1 following the left branch of the unstable manifold at 1 -entering this region, must end at one equilibrium point. Given that there is not any other possibility such trajectory must end at 0 . Hence the proof follows.
The phase portrait of (18) for fixed = 2 and > 0 can be seen in Figure 7.

A Monotonicity
Property. Let ( , V) be any fixed (but arbitrary) point belonging to the region Proof. Calculating the partial derivative with respect to in the equality we obtain meanwhile the corresponding partial derivative with respect to is Then 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 1 with or with 0 . In particular, for fixed > 0, if we continuously decrease the parameter starting from , the left branch of the unstable manifold at 1 , ( 1 ), will move continuously downwards within the region F 1 ; meanwhile, the right branch of the stable manifold at , ( ), moves continuously upwards as decreases in the same region. By continuity of the vector field with respect to the parameter and using shooting arguments, there exists a unique, * > 0, value of for which both manifolds touch each other resulting in a saddle-saddle heteroclinic trajectory connecting the equilibrium points 1 with * . This reasoning constitutes the proof of the following lemma.   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 in a small neighborhood, ( * ), of the critical value * , involves strong dynamical changes including the destruction of the saddle-saddle trajectory and the emergence of a saddle ( 1 ) saddle-node ( 0 ) connexion or the disappearance of heteroclinic trajectories at all. Proof. 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 ( , ) = ( − ) ≡ ( ), 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 . In Figure 9 the sharp traveling wave can be seen.
We close this section by numerically exploring the influence of changes in on the critical value, * , of for which the r-d-a equation has a sharp type solution. This is the content of the next subsection.

The Speed * Depending on
in a Particular Case. As we already mentioned in Section 1, the equation = − + (1 − ) has a monotone decreasing TWS connecting the states 1 ( , ) ≡ 1 and 0 ( , ) ≡ 0 for each ≥ ( ), where the explicit form of ( ) 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 , in particular for those of sharp type. In this subsection we illustrate this relationship through a particular case. To this aim, we choose ( ) = + 2 and ( ) = (1 − ) 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 , the corresponding critical values, * , of for which a unique saddle-saddle heteroclinic trajectory exists (one for each ). 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 * depends on through the corresponding phase portrait, with ( ) and ( ) as before.
This information tells us that, on this range of the "numerical experiments," the speed * is a growing function of . Moreover, we can distinguish two qualitative parts: one exponential for ≤ 8 and another linear for ≥ 8. We carried out the corresponding fittings. These are our results for each phase: * ( ) = 1.007151763 0.1662829302⋅ , for ≤ 8, (37) * ( ) = 0.5010815735 ⋅ + 0.2125604852, respectively. See Figure 10(b). In Figure 11 we show how * changes as function of the parameter for ℎ ( ) = and ℎ ( ) = for comparison. As result of the numerical experiments, we can see that for < 0 the critical value * is small but positive, but once increases, the values of * increase faster with . The physical interpretation of this is as follows: for < 0 the "wind" ℎ ( ) = acts in the opposite direction to which the wave 14 The Scientific World Journal travels; meanwhile for > 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 increases.

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 , , and ℎ 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 ℎ, in addition to satisfying the conditions in item (2), is such that ℎ (0) might have different signs. In particular when ℎ (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 ( , ) = ( − ) ≡ ( ) with > 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 , ℎ, and , the functions 1 and 2 are defined on the whole interval [0, 1] whenever the inequality, holds which, in turn, gives us a bound for for which the two branches (42) are defined on [0, 1]. This is When the above inequality does not hold, 1  (1) . (45) Since we assumed both (0) and (1) are positive, depending on compared with ℎ (0) (or with ℎ (1)), the following cases might occur: ( ,ℎ (0)) = 0 and det [ 1 , 2 ] ( ,ℎ (0)) = −( − ℎ (0)) 2 from which it follows that for ( ,ℎ (0)) < 0. Thus, whenever ̸ = ℎ (0), ( ,ℎ (0)) is a hyperbolic saddle. On the contrary, for = ℎ (0) the equilibrium ( ,ℎ (0)) is a nonhyperbolic point of codimension two.

16
The Scientific World Journal  For completeness, we start by showing the phase portrait for = 0. In Figures 12(a) and 12(b) there are three equilibrium points and a heteroclinic connexion from 0 to 1 , in Figure 12(c) the heteroclinic connexion is from to 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 * , that is, one connexion from 1 to , and in Figure 12(h) for each > * there is connexion from 1 to 0 .
In Figure 13 we show the corresponding front traveling wave profiles for two values of . These were obtained by numerically solving (39) for ℎ ( ) = + with the numerical routine NAG P03PCF (NAG Matlab Toolbox). Figure 14 shows the sharp traveling wave. We also plot how * changes with in this case (see Figure 16). Although we are concentrating on heteroclinic connections leaving 1 , we display the phase portrait behavior for a broad range of values for to show the richness of the dynamics. To this end, we would also like to show how the heteroclinic connexion from 0 to 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.

4.2.
2. ℎ (0) < 0. This case is basically the same as before. We have two main subcases according to the sign of − ℎ (0). Figure 17 shows how * changes with for ℎ ( ) = + . 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 in Figure 17; note how here also there is linear dependency to later become an exponential one. For a comparison, refer to Figure 10.

Degeneracy in Both Diffusion and Advection Terms at
= 0. Here (0) = ℎ (0) = 0. Then the equilibrium ( ,0) coincides with the equilibrium, , which is allocated on the seminegative vertical V-axis for the corresponding system when no advection term is included. Here we will take ℎ ( ) = + 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) ∈ F 1 the angle, ( , V; ), 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 .
Proof. We have and then, by simply calculating the derivative with respect to in the above equality, we obtain for all ( , V) ∈ F 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 for ℎ ( ) = + 2 and in Figure 19 the sharp type. Figure 20 shows how * changes with for ℎ ( ) = + 2 .

Discussion
We have presented a study on the effect of the incorporation of a nonlinear advection term in the degenerate The Scientific World Journal 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) ℎ ( ) = , (2) ℎ ( ) = , and (3) no specific form for ℎ ( ), as long as it satisfies the quite basic requirements as stated in item (2). In all these cases, (3) has (1) no traveling wave solutions for 0 < < * , (2) a unique traveling wave solution of sharp type for = * , 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 = * for which there are sharp traveling waves depends initially linearly on the advection speed to later increase exponentially.  and̃= −1/ . Both of them are outside of the region of interest. Moreover, given that the vertical axis of the V plane is an invariant set of the system, the dynamics around them does not influence the dynamics on the region of interest.