Stability and Numerical Analysis of the Hébraud-Lequeux Model for Suspensions

Abstract

We study both analytically and numerically the stability of the solutions of the Hébraud-Lequeux equation. This parabolic equation models the evolution for the probability of finding a stress 𝜎 in a mesoscopic block of a concentrated suspension, a non-Newtonian fluid. We prove a new result concerning the stability of the fixed points of the equation, and pose some conjectures about stability, based on numerical evidence.

1. Introduction

Non-Newtonian (or complex) fluids are very common in nature and industry, appearing for instance, in foods, biofluids, personal care products, pharmacology and bioengineering, electronics and optical materials, and energy and plastic production. In fact, one could say that Newtonian (or simple) fluids are rather an exception (if not an idealization), even though they include such a prominent member as water. Attending to their rheologic properties, complex fluids are classified in different categories, including suspensions, colloids, polymer melts, liquid crystals, gels and foams, among others.

In this paper, we will consider numerical approximations for the Hébraud-Lequeux equation for suspensions [1] and, more specifically, for the following model (proposed in [2]):𝜕_𝑡𝑝=−𝑏(𝑡)𝜕_𝜎𝑝+𝐷𝑝(𝑡)𝜕²_𝜎𝜎1𝑝−𝑇₀𝜒_{ℝ⧵[−𝜎𝑐,𝜎𝑐]}(𝜎)𝐷𝑝(𝑡)𝑝+𝛼𝛿₀(𝜎)𝑝(),𝑝≥0,0,𝜎=𝑝⁰(𝜎), (1.1) with 𝑝=𝑝(𝑡,𝜎),𝑡≥0,𝜎∈ℝ, and 𝐷(𝑝(𝑡))≡𝐷(𝑝(𝑡,⋅))=𝛼/𝑇₀∫_{|𝜎|>𝜎𝑐}𝑝(𝑡,𝜎)𝑑𝜎. Furthermore, 𝜒_𝐼 denotes the characteristic function of the interval 𝐼 and 𝛿₀ is the Dirac delta function supported at the origin.

We point out the important fact that the function 𝑝(𝑡,𝜎) is a probability density function for every fixed 𝑡, so it must satisfy that 𝑝(𝑡,𝜎)≥0 and ∫_ℝ𝑝(𝑡,𝜎)𝑑𝜎=1, for any 𝑡≥0.

Existence of solutions of (1.1) is given in [3]. In [4] a simplified version of this model is studied, in which the term 𝑏(𝑡)𝜕_𝜎𝑝 is neglected, proving that in such case the weak solutions of the equation generate a multivalued dynamical system possessing a global attractor with respect to a suitable phase space.

In [5] we considered a lattice dynamical system of the simplified model, obtained by discretizating the derivative 𝜕²_𝜎𝜎𝑝 and using a Simpson's quadrature for the integral 𝐷(𝑝). Also, the Dirac delta function 𝛿 is approximated by a parabolic-type function. Observe that the variable 𝜎 belongs to the whole real line, so the lattice dynamical system consists of an infinite number of ordinary differential equations. We proved that this system defines a continuous semigroup of operators in an appropriate space and, as well, the existence of a global compact connected attractor.

Moreover, by a suitable truncation, finite-dimensional approximations of the lattice system are given, proving that each of the approximating systems possesses a global compact connected attractor and also the upper semicontinuous convergence of these attractors to the one generated by the lattice system.

Finally, an implicit Euler scheme is implemented in the finite-dimensional approximations. Thus, a numerical algorithm approximating the solutions of the original equation is obtained. In [5] the resulting discrete dynamical system is studied, proving the existence of a global compact connected attractor and the upper semicontinuity of the sequence of attractors with respect to the previous finite-dimensional approximations.

In this paper, using this algorithm, we present several numerical approximations of the solutions of the equation.

At the same time, we prove a new result concerning the stability of the fixed points of the equation and give also some conjectures which are supported by the numerical simulations.

More precisely, let us consider the simplified equation (i.e., 𝑏(𝑡)=0). Also, for simplicity we take 𝜎_𝑐=1 and 𝑇₀=1. In such a case it is well known [3] that all the probability densities with support in the interval [−1,1] are fixed points. Also, if 𝛼≤0.5, then there no more fixed points, and if 𝛼>0.5, then there exists a unique fixed point not contained inside the interval [−1,1]. It is known that this fixed point is asymptotically stable [6].

In this paper, we prove that if 𝛼<0.5, a subset of probability densities with support in the interval [−1,1] is stable.

Furthermore, based on our numerical simulations we make some more conjectures about the stability of fixed points and also about the asymptotic behavior of solutions as time goes to +∞.

First, it is natural to think that when 𝛼>0.5 the unique fixed point with support not contained inside the interval [−1,1] if globally asymptotically stable, that is, if the initial data is not a fixed point, then the corresponding solution will converge to this point as time goes to +∞ (which would imply that all the other fixed points are unstable). However, it seems that this is not true, as we found numerically a solution converging to a fixed point with support in the interval [−1,1]. Hence, several stable fixed points seem to coexist, although this fact has not been proved yet.

Secondly, the numerical simulations support the conclusion that every solution of the equation converges to a fixed point as time goes to +∞, in other words, that the omega-limit set of every trajectory is a fixed point. Again, there is still no proof of this fact.

Finally, it is natural to think that most fixed points with support in the interval [−1,1] are unstable, since if we take an arbitrary fixed point and a close initial data, usually the numerical solution converges to distant fixed point.

2. Implicit Finite Difference Scheme

We shall consider the simplified model that results from (1.1) after taking 𝑏(𝑡)≡0. Also, we set 𝜎_𝑐=1 for simplicity, although our results can be obtained for an arbitrary 𝜎_𝑐>0. Hence, we consider the equation𝜕_𝑡𝑝=𝐷𝑝(𝑡)𝜕²_𝜎𝜎1𝑝−𝑇₀𝜒_{ℝ⧵[−1,1]}(𝜎)𝐷𝑝(𝑡)𝑝+𝛼𝛿₀(𝜎)𝑝(),𝑝≥0,0,𝜎=𝑝⁰(𝜎), (2.1) with 𝑝=𝑝(𝑡,𝜎),𝑡≥0,𝜎∈ℝ,𝐷(𝑝(𝑡))≡𝐷(𝑝(𝑡,⋅))=𝛼/𝑇₀∫_|𝜎|>1𝑝(𝑡,𝜎)𝑑𝜎. Here, 𝜒_𝐼 is the characteristic function on the interval 𝐼.

Consider an initial condition 𝑝⁰ satisfying 𝑝⁰∈𝐿^∞(ℝ)∩𝐿¹(ℝ),𝑝⁰≥0,_ℝ𝑝⁰𝑑𝜎=1𝐷𝑝⁰>0. (2.2) Then it is proved in [3, Theorem 1.1] that problem (2.1) has a unique solution 𝑝(𝑡,𝜎) satisfying the following properties for all 𝑇>0:𝑝∈𝐿^∞0,𝑇;𝐿²(ℝ)∩𝐿¹(ℝ)∩𝐿²0,𝑇;𝐻¹(ℝ),𝑝∈𝐿^∞((0,𝑇)×ℝ)∩𝐶[0,𝑇],𝐿²(ℝ)∩𝐿¹(ℝ),_ℝ𝑝()(𝑡)𝐷𝑝(𝑡)𝑡,𝜎𝑑𝜎=1,𝑝≥0∀𝑡≥0,∈𝐶([0,𝑇]),min_{0≤𝑡≤𝑇}𝐷𝑝(𝑡)≥𝜈(𝑇)>0,𝜎𝑝∈𝐿^∞0,𝑇;𝐿¹(ℝ), (2.3) where 𝜈(𝑇) exists for any 𝑇>0. The last property implies that the average stress 𝜏(𝑡)=∫_ℝ𝜎𝑝(𝑡,𝜎)𝑑𝜎 belongs to 𝐿^∞(0,𝑇). From 𝛿₀∈𝐻⁻¹(ℝ) and the above properties, it follows that 𝜕_𝑡𝑝∈𝐿²(0,𝑇;𝐻⁻¹(ℝ)).

Our aim is to approximate the solution to the nonlinear partial differential equation given in (2.1) by solving an implicit discrete problem based on finite differences. We begin by discretizing the unbounded spacial domain ℝ via a uniform grid with space step ℎ=Δ𝜎. Let (𝜎_𝑖)_𝑖∈ℤ⊂ℝ be such that 𝜎₀=0 and |𝜎_𝑖+1−𝜎_𝑖|=ℎ. We fix ℎ so as 𝜎_−2𝑛1=−1 and 𝜎_2𝑛1=1. In our numerical simulations, we will always take 𝑝⁰ in such a way that its support is contained in [−𝑆,𝑆] with 𝑆=ℎ𝑁 sufficiently large, and we truncate this unbounded grid for |𝑖|>𝑁. If 𝑝⁰ has an unbounded support, we can take 𝑆 large enough in order to achieve ∫^𝑆_−𝑆𝑝⁰𝑑𝜎≈1. Similarly, we discretize the time domain with a time step Δ𝑡=𝑠. The resulting numerical grid is illustrated in Figure 1.

Figure 1 

Numerical grid.

Denote by ‖Δ‖_ℝ^2𝑁+1 and (Δ,Δ)_ℝ^2𝑁+1 the usual norm and scalar product in ℝ^2𝑁+1. Let us define the following norms in ℝ^2𝑁+1:𝑝_𝑙^1𝑁∶=∑^𝑁_𝑖=−𝑁(|𝑖|+1)𝑝_𝑖𝑝_𝑙𝑆^1𝑁ℎ∶=3_{|𝑖|𝑖≤𝑁iseven}2𝑝_𝑖|𝑖|+≤𝑁_𝑖isodd4𝑝_𝑖, (2.4) for 𝑝∈ℝ^2𝑁+1. Obviously, all these norms are equivalent in ℝ^2𝑁+1.

We shall implement an implicit scheme in finite differences by using the following approximations. Let 𝑝_ℎ^(𝑛)=(𝑝^𝑛_𝑖)_{𝑖∈(−𝑁,𝑁)} denote the vector grid satisfying 𝑝^𝑛_𝑖≈𝑝(𝑡^𝑛,𝜎_𝑖). Since the partial differential equation (2.1) is of evolution type, we approximate 𝜕_𝑡𝑝 at the point (𝑡^𝑛,𝜎_𝑖) by the forward difference 𝜕_𝑡𝑝𝑡^𝑛,𝜎_𝑖≈𝑝_𝑖^𝑛+1−𝑝^𝑛_𝑖𝑠. (2.5)

In order to approximate 𝜕_𝜎𝜎𝑝 at the point (𝑡^𝑛,𝜎_𝑖), we define the operator 𝐴^𝑁_ℎ∶ℝ^2𝑁+1→ℝ^2𝑁+1 by 𝐴^𝑁_ℎ∶=(1/(2ℎ)²)𝐴^𝑁, where𝐴^𝑁∶=10−10⋯000010−10⋯00−1020−10⋯00⋱⋱⋱⋱⋱⋱⋮⋮⋱⋱⋱⋱⋱⋱00⋯0−1020−100⋯0−1010000⋯0−101_{(2𝑁+1)×(2𝑁+1)}, (2.6) and we take −(𝜕²𝑝/𝜕𝜎²)(𝑡^𝑛,𝜎_𝑖)≈(𝐴^𝑁_ℎ𝑝_ℎ^(𝑛+1))_𝑖. This choice is necessary for the matrix 𝐴^𝑁_ℎ to satisfy the condition𝐴^𝑁_ℎ𝑢,𝑣_ℝ2𝑁+1=𝑢,𝐴^𝑁_ℎ𝑣_ℝ2𝑁+1=0, (2.7) where 𝑣=𝑎_{𝑖|𝑖|≤𝑁}=ℎ32^𝛽,…,2,4,2^()𝑖=0,4,2,…,2^𝛽, (2.8) being 𝛽=1, if 𝑁 is even, and being 𝛽=2, if 𝑁 is odd. This follows from (𝐴^𝑁_ℎ𝑣)_𝑖=0, for any 𝑖. Observe also that (𝑢,𝑣)_ℝ^2𝑁+1=‖𝑢‖_𝑙𝑆^1𝑁, which together with equality (2.7) is a key fact to prove Theorem 2.1 (see below).

We define also the matrix 𝐵^𝑁_ℎ∶=(1/(2ℎ))𝐵^𝑁, where𝐵^𝑁∶=0000⋯00000000⋯0010−1000⋯00⋱⋱⋱⋱⋱⋱⋮⋮⋱⋱⋱⋱⋱⋱⋮0⋯010−10000⋯010−10000⋯010−1_{(2𝑁+1)×(2𝑁+1)}, (2.9) which satisfies 𝐴^𝑁_ℎ=(𝐵^𝑁_ℎ)^𝑇𝐵^𝑁_ℎ, where the superscript “𝑇” stands for “transpose”.

The integral (𝛼/𝑇₀)∫_|𝜎|>1𝑝(𝑡,𝜎)𝑑𝜎 is approximated by 𝐷^𝑁_ℎ𝑝^𝑛_ℎ=𝛼𝑇₀ℎ3_{2𝑛1<|𝑘|𝑘≤𝑁isodd}4𝑝^𝑛_𝑘+_{2𝑛1<|𝑘|𝑘≤𝑁iseven}2𝑝^𝑛_𝑘, (2.10) which is a slight variation of the Simpson formula. The only difference is that the weight for the terms 𝑝_2𝑛1 and 𝑝_−2𝑛1 is equal to zero instead of 1. Observe that we can also set 𝐷^𝑁_ℎ(𝑝_ℎ^(𝑛))=𝛼/𝑇₀(𝑝_ℎ^(𝑛),𝜒^𝑁_𝑛1𝑣)_ℝ2𝑁+1, where𝜒^𝑁_𝑛1=1⋱𝑁−2𝑛₁1𝟎𝟎𝟎0⋱4𝑛₁+10𝟎𝟎𝟎1⋱𝑁−2𝑛₁1_{(2𝑁+1)×(2𝑁+1)}. (2.11)

Finally, the 𝛿-function 𝛿₀ will be approximated by the following piecewise parabolic function: 𝛿^𝑐(𝜎)=30for𝜎≤−,42𝑐9𝑐³𝜎²+3𝑐²𝜎+9𝑐43for−42𝑐≤𝜎≤0,9𝑐³𝜎²−3𝑐²𝜎+9𝑐43for0≤𝜎≤,32𝑐0for𝜎≥.2𝑐 (2.12)

This function is continuous, its support is [−3/(2𝑐),3/(2𝑐)], and it attains the maximum value at 𝜎=0, which is equal to 𝑐>0. Also, the integral ∫_ℝ𝛿^𝑐(𝜎)𝑑𝜎 equals 1 (see Figure 2). Let us take 𝑐 such that 3/(2𝑐ℎ)=𝑛_𝑐 is even and let 𝑛_𝑐<2𝑛₁ so as 3/(2𝑐)<1, that is, the support of 𝛿^𝑐 is strictly included in the interval [−1,1]. Then at 𝜎_𝑖=𝑖ℎ this function looks as 𝛿^𝑐𝜎_𝑖=0for𝑖≤−𝑛_𝑐,3𝜎²_𝑖+6ℎ𝑛_𝑐𝜎_𝑖+3ℎ²𝑛²_𝑐2ℎ³𝑛³_𝑐for−𝑛_𝑐≤𝑖≤0,3𝜎²_𝑖−6ℎ𝑛_𝑐𝜎_𝑖+3ℎ²𝑛²_𝑐2ℎ³𝑛³_𝑐for0≤𝑖≤𝑛_𝑐,0for𝑖≥𝑛_𝑐, (2.13) or 𝛿^𝑐𝜎_𝑖=0for𝑖≤−𝑛_𝑐,3𝑖²+6𝑛_𝑐𝑖+3𝑛²_𝑐2ℎ𝑛³_𝑐for−𝑛_𝑐≤𝑖≤0,3𝑖²−6𝑛_𝑐𝑖+3𝑛²_𝑐2ℎ𝑛³_𝑐for0≤𝑖≤𝑛_𝑐,0for𝑖≥𝑛_𝑐. (2.14) Denote 𝛿^𝑐_𝑖=𝛿^𝑐(𝜎_𝑖).

Figure 2 

Approximation of Dirac delta function.

Substituting all the above discretization formulas into the partial differential equation (2.1), we get the following algebraic system for each integer 𝑛≥0:𝑝_𝑖^𝑛+1−𝑝^𝑛_𝑖𝑠=−𝐷^𝑁_ℎ𝑝_ℎ^(𝑛)𝐴^𝑁_ℎ𝑝_ℎ^(𝑛+1)_𝑖−1𝑇₀𝜒_{ℤ⧵[−2𝑛1,2𝑛1]}(𝑖)𝑝^𝑛_𝑖+𝐷^𝑁_ℎ𝑝_ℎ^(𝑛)𝛼𝛿^𝑐_𝑖,|𝑖|≤𝑁,𝑛≥0. (2.15)

Bearing in mind that 𝐴^𝑁_ℎ=(1/(2ℎ)²)𝐴^𝑁 and denoting by 𝜅=𝑠/(2ℎ)²=𝑠/(4ℎ²), the difference scheme system (2.15) can then be equivalently written in the form 𝑝_𝑖^𝑛+1+𝜅𝐷^𝑁_ℎ𝑝_ℎ^(𝑛)⋅𝐴^𝑁𝑝_ℎ^(𝑛+1)_𝑖=𝑝^𝑛_𝑖−𝑠𝑇₀𝜒_{ℤ⧵[−2𝑛1,2𝑛1]}(𝑖)𝑝^𝑛_𝑖+𝑠𝛼𝐷^𝑁_ℎ𝑝_ℎ^(𝑛)𝛿^𝑐_𝑖, (2.16) which can be expressed in matrix form as𝑀_𝑛𝑝_ℎ^(𝑛+1)=𝐶𝑝_ℎ^(𝑛)+𝑠𝛼𝐷_𝑛𝛿^𝑐_ℎ, (2.17) where 𝑀_𝑛=𝐼_2𝑁+1+𝜅𝐷_𝑛𝐴^𝑁,𝐶=𝐼_2𝑁+1−𝑠𝑇₀𝜒^𝑁_𝑛1,𝐷_𝑛=𝐷^𝑁_ℎ𝑝_ℎ^(𝑛),𝛿^𝑐_ℎ=𝛿^𝑐_𝑖. (2.18) Observe that the matrix 𝐶 depends on the time step 𝑠, while 𝑀_𝑛 depends on the steps ℎ and 𝑠. It is easy to check that 𝑀_𝑛 is strictly diagonally dominant for all 𝑛∈ℕ. In particular, this fact guarantees that 𝑀_𝑛 is nonsingular and the system (2.17) has a unique solution for every 𝑛>0. Recall also the following result.

Theorem 2.1 (see [5]). Assume that 0<𝑠≤𝑇₀. If 𝑝⁰_ℎ≥0 and ‖𝑝⁰_ℎ‖_{𝑙𝑆1𝑁}=1, then 𝑝_ℎ^(𝑛)≥0 and ‖𝑝_ℎ^(𝑛)‖_{𝑙𝑆1𝑁}=1, for all 𝑛≥1.

3. Numerical Simulation

We present in this section some numerical solutions of the simplified model (2.1) made with the finite difference scheme (2.15). For simplicity, let us consider the case where 𝑇₀=1. We shall consider two cases depending whether 𝛼>0.5 or 𝛼≤0.5 with different initial conditions 𝑝⁰, which may or not have support intersecting the interval [−1,1]. We shall not plot cases when the support of 𝑝⁰ is totally contained in [−1,1], since in such a case 𝑝⁰ would be a fixed point of the discrete system (2.15). Therefore, we will focus on cases such that supp(𝑝⁰)⋂(ℝ⧵[−1,1])≠∅.

Case 1 (𝛼>0.5). In this case all the probability densities whose support is totally contained in [−1,1] are fixed points of the system. It is known [3] that for the continuous system, the unique fixed point whose support is not totally contained in [−1,1] is given by 𝑝_𝑓(𝜎)=√𝐷𝑒2𝛼^{√(1+𝜎)/𝐷}1for𝜎≤−1,√2𝛼𝜎+𝐷+1−12𝛼for−1≤𝜎≤0,√2𝛼𝜎+𝐷+1√2𝛼for0≤𝜎≤1,𝐷𝑒2𝛼^{√(1−𝜎)/𝐷}for𝜎≥1, (3.1) where 𝐷=𝐷(𝑝_𝑓) satisfies the equation 𝐷+√𝐷=𝛼−(1/2) in order for the normalization constraint ∫_ℝ𝑝_𝑓=1 to hold true. Furthermore, this fixed point is stable (see [6]).
In the next experiments, we shall take the constants 𝛼=0.8 and 𝑇₀=1 and we define the spacial and time mesh size Δ𝜎=ℎ=0.001 and Δ𝑡=𝑠=0.1, respectively. Other choices of 𝛼>0.5 lead to similar results. The numerical simulations will be plotted as a sequence of discrete functions converging to a fixed point. The initial condition and the fixed point will be drawn in black, and between both there are a sequence of plots converging to the fixed point which are drawn in a thinner black. We shall consider different cases depending on the initial conditions, and we shall discuss the results.
Subcase 1. Consider the initial condition 𝑝⁰ with support out of [−1,1] given by 𝑝⁰(𝜎)=0for𝜎<−2.5,0.75for−2.5≤𝜎≤−1.5,0for−1.5<𝑥<1.5,0.25for1.5≤𝜎≤2.5,0for𝜎>2.5. (3.2) The numerical results (see Table 1 and Figure 3) show how the sequence converges to the fixed point (3.1).Subcase 2. Consider the initial condition 𝑝⁰ with support in and out of [−1,1] given by 𝑝⁰(𝜎)=0for𝜎<0,0.5for0≤𝜎≤2,0for𝜎>2. (3.3) We observe (see Table 2 and Figure 4) the same behavior as in Subcase 1. Subcase 3. Consider a parabolic-type initial condition with support in and out of [−1,1] given by 𝑝⁰(𝜎)=30for𝜎<−2,𝜎16²for−2≤𝜎≤2,0for𝜎>2. (3.4) We obtain a similar behavior as in the above cases (see Table 3 and Figure 5).Subcase 4. We start from an initial condition verifying 𝐷(𝑝⁰)≈0 given by 𝑝⁰(𝜎)=0for𝜎<−0.3,1.65for−0.3≤𝜎≤0.3,0for0.3<𝜎<1.5,0.02for1.5≤𝜎≤2,0for𝜎>2. (3.5) The sequence approaches to a fixed point whose support is contained in [−1,1] (see Table 4 and Figure 6).
This fact leads to the hypothesis that besides the fixed point (3.1) there are stable fixed points 𝑝 for (2.1) with support contained in [−1,1].Subcase 5. As in Subcase 4, consider an initial condition verifying 𝐷(𝑝⁰)≈0 given by 𝑝⁰(𝜎)=0for𝜎<−1,0.49995for−1≤𝜎≤1,0for1<𝜎<1.5,0.0001for1.5≤𝜎≤2.5,0for𝜎>2.5. (3.6) On the contrary, this sequence approaches to the fixed point (3.1) (see Table 5 and Figure 7). In fact, this is the usual behavior of the system in the majority of cases.

Table 1 

Numerical results: Subcase 1

Table 2 

Numerical results: Subcase 2.

Table 3 

Numerical results: Subcase 3.

Table 4 

Numerical results: Subcase 4.

Table 5 

Numerical results: Subcase 5.

Figure 3 

Graphics sequence: Subcase 1.

Figure 4 

Graphics sequence: Subcase 2.

Figure 5 

Graphics sequence: Subcase 3.

Figure 6 

Graphics sequence: Subcase 4.

Figure 7 

Graphics sequence: Subcase 5.

We note that, usually, when we take an arbitrary fixed point with support in the interval [−1,1] and an initial data close to it, then the numerical solution converges as time goes +∞ to a fixed point situated far away. Hence, we think that most of the fixed points of such type are unstable. However, we prove in Section 2 that some of them are stable if 𝛼<0.5.

Case 2 (𝛼≤0.5). In this case the probability densities with support totally included in [−1,1] are the unique fixed points. Let us take the constants 𝛼=0.3 and 𝑇₀=1. We define the spacial and time mesh size Δ𝜎=ℎ=0.001 and Δ𝑡=𝑠=0.1, respectively. As in Case 1, the numerical simulations will be plotted by a sequence of plots converging to a fixed point. The initial condition and the fix point are both drawn in a thicker black. Let us take the same first three initial conditions as in Case 1.Subcase 1. Consider an initial condition 𝑝⁰ with support out of [−1,1] given by (3.2). In this case, the sequence converges to a fixed point with support in [−1,1] (see Table 6 and Figure 8).Subcase 2. Consider an initial condition 𝑝⁰ with support in and out of [−1,1] given by (3.3). We also obtain a sequence approaching to a fixed point with support in [−1,1] (see Table 7 and Figure 9).Subcase 3. Consider a parabolic initial condition 𝑝⁰ with support in and out of [−1,1] given by (3.4). A similar behavior is obtained (see Table 8 and Figure 10).
Observe that in the three cases we obtain a sequence converging to different fixed points, although all of them have support [−1,1]. This fact leads to the hypothesis that if 𝛼≤0.5 and is taken an initial condition whose support is not contained in [−1,1], then the iterates of the discrete system converge to a fixed point with support in [−1,1].

Table 6 

Numerical results: Subcase 1.

Table 7 

Numerical results: Subcase 2.

Table 8 

Numerical results: Subcase 3.

Figure 8 

Graphics sequence: Subcase 1.

Figure 9 

Graphics sequence: Subcase 2.

Figure 10 

Graphics sequence: Subcase 3.

4. Stability

Again we take 𝑇₀=1. In this section we prove that if 𝛼<0.5, then some fixed points with support totally included in [−1,1] are stable. In fact, due to the numerical simulations, we believe that this result holds true also for 𝛼≥0.5. Some examples are also exhibited.

We consider the function 𝜃(𝜎)=1(1−𝑎)²(𝜎+𝑎)²1for−1≤𝜎≤−𝑎,0for−𝑎≤𝜎≤𝑎,(1−𝑎)²(𝜎−𝑎)²|𝜎||𝜎|for𝑎≤𝜎≤1,1≤𝜃≤2for1≤≤1+𝑏,2for≥1+𝑏, (4.1) where 0<𝑎<1,𝑏>0, 𝜃(𝜎) is smooth on |𝜎|>𝑎 and |𝜃^′′(𝜎)|≤2/(1−𝑎)² for any 𝜎. Sometimes, total or partial derivatives with respect to 𝜎 be denoted by subindices like, for example, in 𝜃_𝜎𝜎 and 𝑝_𝜎 below.

Lemma 4.1. Denote 𝜂=1/(1−𝑎)² and assume that 𝛼<1/2𝜂. Let 𝑝(⋅) be the unique solution corresponding to the initial data 𝑝⁰ satisfying (2.2). Then ^𝑡_0|𝜎|>1𝑝()1𝑡,𝜎𝑑𝜎𝑑𝑡≤1−2𝜂𝛼_ℝ𝑝⁰(𝜎)𝜃(𝜎)𝑑𝜎,∀𝑡≥0. (4.2)

Proof. Let 𝑝(⋅) be the unique solution corresponding to the initial data 𝑝⁰ such that 𝐷(𝑝⁰)>0. Multiplying (2.1) by 𝜃 we have 𝑑𝑑𝑡_ℝ𝑝(𝑡,𝜎)𝜃(𝜎)𝑑𝜎+_ℝ𝟏_{ℝ⧵[−1,1]}(𝜎)𝑝(𝑡,𝜎)𝜃(𝜎)𝑑𝜎=𝐷𝑝(𝑡)_ℝ𝑝(𝑡,𝜎)𝜃_𝜎𝜎(𝜎)𝑝(𝑡)𝑑𝜎≤2𝜂𝐷=2𝜂𝛼_|𝜎|>1𝑝(𝑡,𝜎)𝑑𝜎. (4.3) Hence, 𝑑𝑑𝑡_ℝ𝑝(𝑡,𝜎)𝜃(𝜎)𝑑𝜎+1−2𝜂𝛼_|𝜎|>1𝑝(𝑡,𝜎)𝑑𝜎≤0,1−2𝜂𝛼^𝑡_0|𝜎|>1𝑝(𝑡,𝜎)𝑑𝜎𝑑𝑡≤_ℝ𝑝(0,𝜎)𝜃(𝜎)𝑑𝜎. (4.4) These calculations are formal, but can be justified by using a cutoff function.

Let 𝑓(𝜎)∶[0,1]→ℝ⁺ be such that 𝑓∈𝐶²([0,1]), 𝑓^′′(𝜎)≥0, 𝑓^′(𝜎)≤0, 𝑓(𝑎)=0, and ∫^𝑎₀𝑓(𝜎)𝑑𝜎=1/2. We define 𝑣(𝜎) by 𝑣(𝜎)=𝑓(𝜎)𝑓()|𝜎|for0≤𝜎≤𝑎,−𝜎for−𝑎≤𝜎≤0,0for≥𝑎. (4.5) We note that 𝑣∈𝐻¹(ℝ). As 𝐷(𝑣)=0 and ∫_ℝ𝑓(𝜎)𝑑𝜎=1, we obtain that 𝑣 is a fixed point.

Denote by ‖⋅‖ the norm in 𝐿²(ℝ).

Theorem 4.2. Let 𝛼<(1−𝑎)²/2. If −2𝛼𝑓_𝜎(0)+1²(0)8𝛼−𝑓≤0, (4.6) then the fixed point 𝑣 is stable with respect to initial data 𝑝⁰ satisfying (2.2) in the following sense: for any 𝜀>0 there exists 𝛿>0 such that if ‖𝑝⁰-𝑣‖<𝛿 and ∫_ℝ|𝑝⁰(𝜎)-𝑣(𝜎)|𝑑𝜎<𝛿, then 𝑝(𝑡)-𝑣<𝜀∀𝑡≥0. (4.7)

Proof. Let 𝑝(⋅) be the unique solution corresponding to the initial data 𝑝⁰ such that 𝐷(𝑝⁰)>0. Denote 𝑧(𝑡)=𝑝(𝑡)-𝑣. Multiplying (2.1) by z, we have 12𝑑‖𝑧‖𝑑𝑡²=𝛼_|𝜎|>1𝑝(𝑡,𝜎)𝑑𝜎−_ℝ𝑝_𝜎²𝑑𝜎+_ℝ𝑝_𝜎𝑣_𝜎−𝑑𝜎_|𝜎|>1𝑝(𝑡,𝜎)²𝑑𝜎+_|𝜎|>1𝑝(𝑡,𝜎)𝑑𝜎𝑝(𝑡,0)−𝑓(0). (4.8)
We use that _ℝ𝑝_𝜎𝑣_𝜎𝑑𝜎=𝑓_𝜎(𝑎)𝑝(𝑡,𝑎)−𝑓_𝜎(0)𝑝(𝑡,0)−^𝑎₀𝑝𝑓_𝜎𝜎𝑑𝜎+𝑓_𝜎(𝑎)𝑝(𝑡,−𝑎)−𝑓_𝜎(0)𝑝(𝑡,0)−⁰_−𝑎𝑝𝑓_𝜎𝜎𝑑𝜎≤𝑓_𝜎(𝑎)𝑝(𝑡,𝑎)+𝑝(𝑡,−𝑎)−2𝑓_𝜎(0)𝑝(𝑡,0). (4.9) Denote 𝑏=−𝑓_𝜎(0), 𝑐=𝑓(0). Thus 12𝑑‖𝑧‖𝑑𝑡²+_|𝜎|>1𝑝(𝑡,𝜎)²≤𝑑𝜎_|𝜎|>1𝑝(𝑡,𝜎)²𝑑𝜎(2𝛼𝑏+1)𝑝(𝑡,0)−𝛼_ℝ𝑝_𝜎²𝑑𝜎+𝛼𝑓_𝜎(𝑎)𝑝(𝑡,𝑎)+𝑝(𝑡,−𝑎)−𝑐. (4.10)
Now 𝑝(𝑡,0)=𝑝(𝑡,−𝑎)+⁰_−𝑎𝑝_𝜎()𝑝()()−𝑡,𝜎𝑑𝜎,𝑡,0=𝑝𝑡,𝑎^𝑎₀𝑝_𝜎(𝑡,𝜎)𝑑𝜎, (4.11) so ()𝑝()≤()2𝛼𝑏+1𝑡,02𝛼𝑏+12𝑝()()+()𝑡,−𝑎+𝑝𝑡,𝑎2𝛼𝑏+12^𝑎_−𝑎𝑝_𝜎(≤()𝑡,𝜎)𝑑𝜎2𝛼𝑏+12𝑝()()+()𝑡,−𝑎+𝑝𝑡,𝑎2𝛼𝑏+12^𝑎_−𝑎𝑝_𝜎(𝑡,𝜎)²𝑑𝜎^1/2(2𝑎)^1/2≤(2𝛼𝑏+1)2𝑝(𝑡,−𝑎)+𝑝(𝑡,𝑎)+𝛼^𝑎_−𝑎𝑝_𝜎(𝑡,𝜎)²()𝑑𝜎+2𝛼𝑏+1².8𝛼 (4.12)
Hence 12𝑑‖𝑧‖𝑑𝑡²+_|𝜎|>1𝑝(𝑡,𝜎)²𝑑𝜎≤_|𝜎|>1𝑝(𝑡,𝜎)²()𝑑𝜎2𝛼𝑏+1²8𝛼−𝑐−𝛼_|𝜎|>𝑎𝑝_𝜎²+𝑑𝜎𝛼𝑓_𝜎(𝑎)+(2𝛼𝑏+1)2𝑝(𝑡,𝑎)+𝑝(𝑡,−𝑎). (4.13) Denote 𝑑=𝛼𝑓_𝜎(𝑎)+(2𝛼𝑏+1)/2. Using (4.6), (4.13), and the inequality |𝑢(𝑥)|²≤2‖𝑢‖_𝐿²_(𝐼)‖𝑢_𝜎‖_𝐿2(𝐼), for any 𝑥∈𝐼 (an interval of ℝ) and 𝑢∈𝐻¹(𝐼) (see [7, page 129]), we have 𝑑𝛼𝑝()≤𝑑𝑡,−𝑎𝛼√2^−𝑎_−∞𝑝²𝑑𝜎^−𝑎_−∞𝑝_𝜎²𝑑𝜎^1/4≤^−𝑎_−∞𝑝_𝜎²𝑑𝜎+34𝑑𝛼^4/3−𝑎_−∞𝑝²𝑑𝜎^1/3,𝑑𝛼𝑝()≤𝑑𝑡,𝑎𝛼√2^∞_𝑎𝑝²𝑑𝜎^∞_𝑎𝑝_𝜎²𝑑𝜎^1/4≤^∞_𝑎𝑝_𝜎²𝑑𝜎+34𝑑𝛼^4/3∞_𝑎𝑝²𝑑𝜎^1/3,12𝑑‖𝑧‖𝑑𝑡²+∫_|𝜎|>1𝑝(𝑡,𝜎)²𝑑𝜎≤3𝛼2𝑑𝛼^4/3∫_|𝜎|>1𝑝(𝑡,𝜎)𝑑𝜎∫_|𝜎|>𝑎𝑝(𝑡,𝜎)²𝑑𝜎^1/3. (4.14) It is known [4] that for any 𝑅>0 there exists 𝑁(𝑅) such that 𝑝(𝑡)≤𝑁(𝑅),∀𝑡≥0, (4.15) if ‖𝑝⁰‖≤𝑅. Hence by Lemma 4.1, we have 𝑝(𝑡)−𝑣²=‖𝑧(𝑡)‖²≤‖𝑧(0)‖²𝑑+3𝛼𝛼^4/3(𝑁(𝑅))^2/3𝑡_0|𝜎|>1𝑝(≤𝑝𝑡,𝜎)𝑑𝜎₀−𝑣²+3𝛼^−1/3𝑑^4/3(𝑁(𝑅))^2/31−2𝜂𝛼_ℝ𝑝⁰(𝜎)𝜃(𝜎)𝑑𝜎. (4.16)
We note that ∫_ℝ𝑣(𝜎)𝜃(𝜎)𝑑𝜎=0. Hence, for any 𝛾>0 there exists 𝛿>0 such that ∫_ℝ|𝑝⁰(𝜎)−𝑣(𝜎)|𝑑𝜎<𝛿 implies _ℝ𝑝⁰(𝜎)𝜃(𝜎)𝑑𝜎≤_ℝ𝑝⁰(𝜎)−𝑣(𝜎)𝜃(𝜎)𝑑𝜎≤2_ℝ𝑝⁰(𝜎)−𝑣(𝜎)𝑑𝜎<𝛾. (4.17) From here the stability follows.

Let us consider now particular examples of stable fixed points.

Example 4.3. Take 𝛼=0.1 and 𝑣(𝜎)=𝐶𝐴+𝐵𝜎+1𝜎−𝐷for0≤𝜎≤2,𝐶𝐴−𝐵𝜎+1−𝜎−𝐷for−2|𝜎|≥1≤𝜎≤0,0for2, (4.18) where 𝐶=0.05,𝐷=−0.01, and 𝐵=8𝐶/(𝐷−1)+8𝐶∗log(|(1−2𝐷)/(2𝐷)|)−4, 𝐴=2𝐶/(2𝐷−1)−𝐵/2. Here 𝑎=0.5 and 𝛼<(1−𝑎)²/2=1/8.
This function satisfies ∫₀^1/2𝑣𝑑𝜎=1/2, 𝑣(1/2)=0, and 𝑣_𝜎0⁺=−6.5536,𝑣(0)=8.1537,𝐾=1+2𝛼−𝑣_𝜎0⁺²−8𝛼𝑣(0)=−1.1836, (4.19) so that (4.6) holds. Hence, by Theorem 4.2, 𝑣 is stable in the above sense.
Taking an initial condition close to 𝑣, for instance, 𝑝⁰(𝜎)=𝑣()𝑣()𝜎−𝑎for𝜎≤0,𝜎+𝑎for0≤𝜎<1.5,𝑐for1.5≤𝜎≤2,0for𝜎>2, (4.20) where 𝑎=0.001and𝑐=0.02238, we obtain the numerical results collected in Table 9 (see also Figure 11).
If we take an initial condition closer to 𝑣 than before, for instance, 𝑝⁰(𝜎)=𝑣()𝑣()𝜎−𝑎for𝜎≤0,𝜎+𝑎for0≤𝜎<1.5,𝑐for1.5≤𝜎≤2,0for𝜎>2, (4.21) where 𝑎=0.0001,and𝑐=0.00057414, then we obtain the numerical results collected in Table 10 (see also Figure 12).

Table 9 

Numerical results: Stable fixes point.

Table 10 

Numerical results: Stable fixed point.

Figure 11 

Graphics sequence: Example 4.3.

Figure 12 

Graphics sequence: Example 4.3.

Hence, the solution remains near this fixed point if the initial data is close to it.

Example 4.4. Take 𝛼=0.1 and 𝑣(𝜎)=𝐶𝐴+𝐵𝜎+1𝜎−𝐷for0≤𝜎≤2,𝐶𝐴−𝐵𝜎+1−𝜎−𝐷for−2|𝜎|≥1≤𝜎≤0,0for2, (4.22) where 𝐶=0.1,𝐷=−0.01 and 𝐵=8(𝐶/(𝐷−1))+8𝐶∗log(|((1−2𝐷)/(2𝐷))|)−4, 𝐴=2𝐶/(2𝐷−1)−𝐵/2. Here 𝑎=0.5 and 𝛼<(1−𝑎)²/2=1/8.
This function satisfies ∫₀^1/2𝑣𝑑𝜎=1/2, 𝑣(1/2)=0 and 𝑣_𝜎0⁺=−9.1071,𝑣(0)=14.3075,𝐾=1+2𝛼−𝑣_𝜎0⁺²−8𝛼𝑣(0)=−3.4855, (4.23) so that (4.6) holds.
Hence, by Theorem 4.2, 𝑣 is stable.

Example 4.5. Take 𝛼=0.3 and 𝑣(𝜎)=𝐶𝐴+𝐵𝜎+𝐶𝜎−𝐷for0≤𝜎≤0.2,𝐴−𝐵𝜎+|𝜎|−𝜎−𝐷for−0.2≤𝜎≤0,0for≥0.2, (4.24) where 𝐶=0.11,𝐷=−0.005, and 𝐵=2(𝐶/(𝑎(𝐷−𝑎)))+2(𝐶/𝑎²)∗log(|(𝑎−𝐷)/𝐷|)−1/𝑎², 𝐴=(𝐶/(𝐷−𝑎)))−𝑎𝐵. Here 𝑎=0.2, and 𝛼<(1−𝑎)²/2=0.32.
This function satisfies ∫₀^0.2𝑣𝑑𝜎=1/2, 𝑣(0.2)=0 and 𝑣_𝜎0⁺=−10.0512,𝑣(0)=23.4517,𝐾=1+2𝛼−𝑣_𝜎0⁺²−8𝛼𝑣(0)=−6.8529, (4.25) so that (4.6) holds. Hence, by Theorem 4.2  𝑣 is stable.
We can see that for a given 𝛼<1/2 more than one stable point exists. Finally, we shall give an example of an unstable fixed point for 𝛼=0.1. Let us consider the fixed point given by 𝑢(𝜎)=0for𝜎<−1,0.5for−1≤𝜎≤1,0for𝜎>1, (4.26) and the initial condition 𝑝⁰(𝜎)=𝑐0for𝜎<−1,₁𝑐for−1≤𝜎≤1,0for1<𝜎<1.5,₂for1.5≤𝜎≤2.5,0for𝜎>2.5, (4.27) where the constants 𝑐₁ and 𝑐₂ are chosen in such a way that 𝑝⁰ is close enough. For 𝑐₁=0.495 and 𝑐₂=0.01 we obtain the numerical results given in Table 11 (see Figure 13). For 𝑐₁=0.49995 and 𝑐₂=0.0001 the numerical results are given in Table 12 (see Figure 14). These numerical results show that the fixed point 𝑝_𝑓 is unstable.

Table 11 

Numerical results: Unstable fixed point.

Table 12 

Numerical results: Unstable fixed point.

Figure 13 

Graphics sequence: Unstable fixed point.

Figure 14 

Graphics sequence: Unstable fixed point.

5. Conclusions

For problem (2.1) with 𝑇₀=1, we have proved in Theorem 4.2 that if 𝛼<0.5, then a class of fixed points with support inside the interval [−1,1] is stable. It would be interesting to prove such result for 𝛼≥0.5, as the numerical simulations suggest that this is true. It follows from such hypothesis that the unique fixed point with support not contained in the interval [−1,1] is not globally asymptotically stable, that is, not all solutions which are not fixed points converge to it as time goes to +∞, although it seems that most of solutions behave in this way.

On the other hand, we have obtained in this paper numerical simulations of the solutions leading to the following conjectures.(1)Every solution of the equation converges to a fixed point as time goes to +∞, that is, the omega-limit set of every trajectory is a fixed point.(2)Most fixed points with support in the interval [−1,1] are unstable.

These are interesting open problems to consider for the near future.

Acknowledgments

This work was partially supported by Spanish Ministerio de Ciencia e Innovación, Projects MTM2008-00088 and MTM2009-11820, the Consejería de Innovación, Ciencia y Empresa (Junta de Andalucía), Proyecto de Excelencia P07-FQM-02468, and the Consejería de Cultura y Educación (Comunidad Autónoma de Murcia), Grant 08667/PI/08.