In this paper, we study the stability of traveling waves of a certain chemotaxis model. In a chemotaxis model, formation of traveling fronts is a central theme of the study: let a long tube be filled with a nutrient that is to be consumed by a certain micro-organism, such as bacteria. Let a certain amount of bacteria population be placed at the leftmost entrance of the tube. Conversion of nutrient to a bacteria population at some rate would take place and eventually results in the exhastion of nutrient at the place. If the bacteria is capable of moving from one place to another, then those that move right will find a new opportunity of further growth with a high nutrient density at the new place. As a proof that this does take place, traveling wave solutions where the nutrient front retreats, and that of the bacteria population propagates, comes with an important physical relevance.

Kim et al. [8] introduced a following chemotaxis model:

Here Ω⊂ℝn is a simply connected bounded open set with smooth boundary, u is the population density, n is the nutrient density, r0>0 is the growth rate, and ϵ>0 is a small diffusion constant for the nutrient. Initial and boundary data are assigned as above, where zero flux boundary conditions are imposed for both u and n. ν is the outward normal vector. For the nonlinearty γ:[0,∞)→ℝ, we assume

This chemotaxis model has been suggested from the consideration of non-fickian flux law for diffusion, which is an independent and substantial subject of the random walk in a heterogeneous medium. See the series of studies [7, 5, 6]. The key difference from the classical Keller-Segel model [4] is to have a flux function as in the first equation of (1.1)

As a consequence, migration of population is from the place where p is high to the place p is low. The bacteria are supposed to sense this pressure p by sensing the values of u and n at the place. The value of p is decided by the relative values of n and u. For instance, if p is decreasing in n and increasing in u then migration occurs from the place where the population is crowded, relative to the nutrient, to the place that is less so. This type of chemotaxis model has attracted considerable attention in recent years for a variety of applications (see [7, 1, 14]). Existence theorems on related models also have been actively studied by those authors in [12, 13, 3, 2, 11, 9, 10].

The measurements of this pressure p is done by measurements of n and u at the place by the bacteria and the value of p is decided by the relative values of n and u.

Consider one space dimensional Cauchy problem

Kim et al. [8] proved the existence of traveling wave solutions of (1.3) for the case ϵ=0 that are of the form

Plugging in the ansatz (1.4) to (1.3), one obtains the system of ordinary differential equations

The existence of traveling waves are proved by phase space analysis of (1.6).

Because those traveling waves are constructed by solving the system (1.6), not the system (1.3), once existence is established, the stability against a small perturbation as a solution of (1.3) is subjected to a further study. Any of traveling waves that is not robust against a small perturbation would have little physical significance.

In general, nonlinear stability is a difficult problem, in particular around this non-steady state solution. The objective of this paper is thus to study numerical nonlinear stability of traveling waves: we numerically construct the traveling wave solutions of (1.6) and the computed data are prepared as an initial data of the system (1.3). In this preparation, a small perturbation is added. The system (1.3) with the perturbed initial data is then numerically solved until a finite time T>0. The precise test procedure will be specified in Section 4.

The remainder of this paper is organised as follows. In Section 2 we summarize the existence results of traveling waves presented in [8] for completeness of our discussion. In Section 3 we present numerical integration results done by python SciPy Library for a few relevant instances of traveling waves. In Section 4 we present the numerical stability results, and finally we give conclusions in Section 5.

In this section, the system of odes (2.1) is further discussed, which is a summary of what are presented in [8]. For notational simplicity overlines on variables are suppressed.

Since nu=cn′ for the case ϵ=0, the first equation of (2.1) is integrated once to have

for some integration constant P₀. By imposing the condition (u,n)→(0,n*) as ξ→∞ and using new variables (p,n)=(γ(n)u,n), the system can be arranged in the form

Here, we have used the assumption (1.2) that γ(n)≥γ0>0 for every n.

2.1. Equilibrium points and linearization

It is straightforward to verify that (2.1) has only two equilibrium points that are

This implies that u_* in (1.5) is not an independent parameter but is decided as

Linearization around each equilibrium point is summarized below.

1. (p0,n0) is a saddle point:

2. If c>c*=2r0n*γ(n*) then (p1,n1) is a stable node with two real and negative eigenvalues:

Phase space analysis [8] based on the discussion so far gives the following result:

Theorem 2.1. Fix n*>0 and c>c*. There exists a heteroclinic orbit that is an intersection of the unstable manifold of (p0,n0) and the stable manifold of (p1,n1). The heteroclinic orbit entirely lies in the first quadrant of pn-plane.

As stated in Theorem 2.1, the traveling wave exists as a heterocilinic orbit of (2.1) that is a saddle-to-stable node connection. If (p^,n^) is any point on the stable manifold of (p0,n0), integrating numerically (2.1) from (p^,n^) forward in time, capturing the saddle-to-stable node connection, is not a stiff problem.

Based on that, in this section we numerically compute the saddle-to-stable node connection. We specify the procedure below.

1. We first let δ>0 sufficiently small and set

where X₂ is the eigenvector in Section 2.1.

2. Let M>0 be sufficiently large. We use SciPy Library to integrate (2.1) on interval [0,M] with initial data

We present three instances of traveling waves. As discussed earlier, γ that is decreasing in n corresponds to a chemotaxis model where migration of population uccurs from a crowded place to uncrowded place. Note that the value x in the eigenvector X₂ will be negative in this case because γ′(0)<0. γ that is increasing in n also is relevant to a certain application. In this case, the sign of x in the eigenvector X₂ is decided by n*γ′(0)γ(0)−1. In case x>0, one observes an overshoot of p profile before it converges to 0 as ξ→∞. See Figure 1 (c).

Figure 1. Traveling waves of cases A, B, and C.

Considering the above, those three instances are specified below, including various wave speeds and parameters.

Numerically captured heteroclinic orbits are presented in Figure 1. In every case, n has a retreating front and p has a propagating front. In wave C, we see the overshoot of p profiles. On the other hands, in Figure 1 (a) and (b) where we have x<0, we observe p decreasing all the time in ξ converging to 0 and n increasing all the time converging to n_*.

Any traveling wave solution that is not robust under a perturbation has little physical significance. In this section, we conduct the numerical nonlinear stability tests on three captured data in the preceding section. The specification of stability test procedure is presented below.

4.1. Procedure of numerical stability test

1. Points on the heteroclinic orbit are collected, for j=0,⋯,N−1 with h=MN, that are

Note that the heteroclinic orbits have been computed on sufficiently large interval [0,M] so that profiles of orbits are sufficiently flat near end points. See Figure 1.

2. Preparation of (u,n) profiles. Since variables we integrate are p and n in heteroclinic orbit, we take for j=0,⋯,N−1

3. Preparation of initial data U₀. The system we will run is the cauchy problem (1.3) in the wholse space ℝ. In its approximation, computational domain has to be large enough so that within a certain time T>0 traveling wave phenomena can be observed. For this reason, we take an interval [0,M¯] with M¯ possibly larger than M as a computational domain for (1.3). The heteroclinic orbit data is translated and extended to U₀ data as follows.

for i₀ an integer shift factor.

4. Preparation of a perturbed initial data U˜0. For a given U0, we consider two types of absolute and relative perturbations: For i=0,⋯,N¯−1

for some E with ‖E‖∞=1. η¯,η>0 are small parameters. η¯ measures the level of absolute noises and η measures the level of relative noises.

5. Finally, we numerically solve (1.3) in [0,T]×[0,M¯]. The time interval is discretized as

For the parabolic solver, we simply used the explicit scheme:

for m=1⋯K¯−1 and j=1,⋯,N¯−2. Assuming M¯ is large enough, we assigned the Dirichilet boundary condition, that is for every m

4.2. Results

In Figure 2 and 3 are the tests of wave A against an oscillatory sinusoidal perturbation, i.e.,

Figure 2. Test of wave A against the sinusoidal absolute additive noise. For the computation, h=0.1, k=0.005, η¯=0.01 have been used.

Figure 3. Test of wave A against the sinusoidal absolute additive noise. For the computation, h=0.1, k=0.005, η¯=0.05 have been used.

with l=10.

We first consider the severer absolute noises. Figure 2 presents snapshots of p=uγ(n) and n profiles taken at time t=0, 0.5, 1.0, and 5.0 with η¯=0.01. The value of η¯ is the 1% of the maximum ‖(p,n)‖∞. Interesting behaviors against the absolute additive noise are exhibited: Oscillations first seem to be stabilized as time proceeds so that by the time around t=1.0 we see almost plateu for both p and n, but this attenuated oscillation is then amplified in later time. This behavior is more drastically observed with absolute noise at η¯=0.05 in Figure 3.

This illustrates that the shape of traveling wave is maintained for a short time in a stable manner, but due to the reaction of the system, traveling wave can be oscillatory in later time.

Note that the two end states at ±∞ have nearly 0 state for one of n and u, or the added error in Figure 2 and 3 in relative sense are huge. Furthermore, in the second equation (1.3) for n, currently we have no diffusion (ϵ=0), which is an extreme case. Considering those, behaviors seems not drastically bad in the sense that the traveling wave does not disappear.

In Figure 4, 5, and 6 are tests against relative additive noises of 10% levels. We conducted computations until t=5.0. What are in the last row of each of Figure 4, 5, and 6 are the snapshots at t=1.0, 2.0, 3.0, 4.0, and 5.0 from the left to right. We considered an ambitious white noise for these tests where E(xi) is a random variable valued in [-1,1] with a uniform probability. To add this white noise implies that the initial data are discontinuous and highly oscillatory.

Figure 4. Snapshots of u and n.

Figure 5. Snapshots of u and n.

Figure 6. Snapshots of u and n.

We see in Figure 4 (a) the farily clean traveling wave fronts of u against the relative additive white noises. White noise in n cannot be flatten out simply because the equation for n in (1.3) is an ode.

In Figure 6 (a) for wave C, we presented the profile of p in place of u, to illustrates the overshoot before the convergence to 0 as x→∞. Figure 6 (a) clearly shows the recovery to the clean profile of p traveling with an overshoot.

In this paper, numerical stability of traveling waves of (1.3) against perturbations has been tested on three instances. The traveling waves themselves were computed by Python SciPy Library, and are assigned as initial data for the system (1.3) after adding a perturbation. We tested both the absolute and relative perturbations. Considering that the end states of traveling waves contain the zero state in one of u and n, the added absolute additive perturbations near the zero state are huge ones if measured in relative sense. Yet we observed fairly stable behaviors against the absolute additive noise in Figure 2 and 3. Against the relative additive noise, traveling waves seem to be highly robust as seen in Figure 4, 5, and 6.