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

Authors: Ángel Giménez; Francisco Morillas; José Valero; José María Amigó
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

## 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.


---

## Body

## 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.1)∂tp=-b(t)∂σp+D(p(t))∂σσ2p-1T0χR∖[-σc,σc](σ)p+D(p(t))αδ0(σ),p≥0,p(0,σ)=p0(σ),
with p=p(t,σ),t≥0,σ∈ℝ, and D(p(t))≡D(p(t,·))=α/T0∫|σ|>σcp(t,σ)dσ. Furthermore, χI denotes the characteristic function of the interval I and δ0 is the Dirac delta function supported at the origin.We point out the important fact that the functionp(t,σ) is a probability density function for every fixed t, so it must satisfy that p(t,σ)≥0 and ∫ℝp(t,σ)dσ=1, for any t≥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 b(t)∂σp 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 ∂σσ2p and using a Simpson's quadrature for the integral D(p). 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.,b(t)=0). Also, for simplicity we take σc=1 and T0=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 b(t)≡0. Also, we set σc=1 for simplicity, although our results can be obtained for an arbitrary σc>0. Hence, we consider the equation(2.1)∂tp=D(p(t))∂σσ2p-1T0χR∖[-1,1](σ)p+D(p(t))αδ0(σ),p≥0,p(0,σ)=p0(σ),
with p=p(t,σ),t≥0,σ∈ℝ,D(p(t))≡D(p(t,·))=α/T0∫|σ|>1p(t,σ)dσ. Here, χI is the characteristic function on the interval I.Consider an initial conditionp0 satisfying (2.2)p0∈L∞(R)∩L1(R),p0≥0,∫Rp0dσ=1D(p0)>0.
Then it is proved in [3, Theorem 1.1] that problem (2.1) has a unique solution p(t,σ) satisfying the following properties for all T>0:(2.3)p∈L∞(0,T;L2(R)∩L1(R))∩L2(0,T;H1(R)),p∈L∞((0,T)×R)∩C([0,T],L2(R)∩L1(R)),∫Rp(t,σ)dσ=1,p(t)≥0∀t≥0,D(p(t))∈C([0,T]),min⁡0≤t≤TD(p(t))≥ν(T)>0,σp∈L∞(0,T;L1(R)),
where ν(T) exists for any T>0. The last property implies that the average stress τ(t)=∫ℝσp(t,σ)dσ belongs to L∞(0,T). From δ0∈H-1(ℝ) and the above properties, it follows that ∂tp∈L2(0,T;H-1(ℝ)).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 h=Δσ. Let (σi)i∈ℤ⊂ℝ be such that σ0=0 and |σi+1-σi|=h. We fix h so as σ-2n1=-1 and σ2n1=1. In our numerical simulations, we will always take p0 in such a way that its support is contained in [-S,S] with S=hN sufficiently large, and we truncate this unbounded grid for |i|>N. If p0 has an unbounded support, we can take S large enough in order to achieve ∫-SSp0dσ≈1. Similarly, we discretize the time domain with a time step Δt=s. The resulting numerical grid is illustrated in Figure 1.Figure 1
Numerical grid.Denote by∥Δ∥ℝ2N+1 and (Δ,Δ)ℝ2N+1 the usual norm and scalar product in ℝ2N+1. Let us define the following norms in ℝ2N+1:(2.4)‖p‖l¯1N:=∑i=-NN(|i|+1)|pi|‖p‖lS1N∶=h3[∑|i|≤Niis even2|pi|+∑|i|≤Niis odd4|pi|],
for p∈ℝ2N+1. Obviously, all these norms are equivalent in ℝ2N+1.We shall implement an implicit scheme in finite differences by using the following approximations. Letph(n)=(pin)i∈(-N,N) denote the vector grid satisfying pin≈p(tn,σi). Since the partial differential equation (2.1) is of evolution type, we approximate ∂tp at the point (tn,σi) by the forward difference (2.5)∂tp(tn,σi)≈pin+1-pins.In order to approximate∂σσp at the point (tn,σi), we define the operator AhN:ℝ2N+1→ℝ2N+1 by AhN:=(1/(2h)2)AN, where(2.6)AN:=[10-10⋯000010-10⋯00-1020-10⋯00⋱⋱⋱⋱⋱⋱⋮⋮⋱⋱⋱⋱⋱⋱00⋯0-1020-100⋯0-1010000⋯0-101](2N+1)×(2N+1),
and we take -(∂2p/∂σ2)(tn,σi)≈(AhNph(n+1))i. This choice is necessary for the matrix AhN to satisfy the condition(2.7)(AhNu,v)R2N+1=(u,AhNv)R2N+1=0,
where (2.8)v=(ai)|i|≤N=h3(2β,…,2,4,2(i=0),4,2,…,2β),
being β=1, if N is even, and being β=2, if N is odd. This follows from (AhNv)i=0, for any i. Observe also that (u,v)ℝ2N+1=∥u∥lS1N, which together with equality (2.7) is a key fact to prove Theorem 2.1 (see below).We define also the matrixBhN:=(1/(2h))BN, where(2.9)BN:=[0000⋯00000000⋯0010-1000⋯00⋱⋱⋱⋱⋱⋱⋮⋮⋱⋱⋱⋱⋱⋱⋮0⋯010-10000⋯010-10000⋯010-1](2N+1)×(2N+1),
which satisfies AhN=(BhN)TBhN, where the superscript “T” stands for “transpose”.The integral(α/T0)∫|σ|>1p(t,σ)dσ is approximated by (2.10)DhN(phn)=αT0h3(∑2n1<|k|≤Nkis odd4pkn+∑2n1<|k|≤Nkis even2pkn),
which is a slight variation of the Simpson formula. The only difference is that the weight for the terms p2n1 and p-2n1 is equal to zero instead of 1. Observe that we can also set DhN(ph(n))=α/T0(ph(n),χn1Nv)ℝ2N+1, where(2.11)χn1N=[1⋱N-2n1︸10000⋱4n1+1︸00001⋱N-2n1︸1](2N+1)×(2N+1).Finally, theδ-function δ0 will be approximated by the following piecewise parabolic function: (2.12)δc(σ)={0forσ≤-32c,49(c3σ2+3c2σ+9c4)for-32c≤σ≤0,49(c3σ2-3c2σ+9c4)for0≤σ≤32c,0forσ≥32c.This function is continuous, its support is[-3/(2c),3/(2c)], and it attains the maximum value at σ=0, which is equal to c>0. Also, the integral ∫ℝδc(σ)dσ equals 1 (see Figure 2). Let us take c such that 3/(2ch)=nc is even and let nc<2n1 so as 3/(2c)<1, that is, the support of δc is strictly included in the interval [-1,1]. Then at σi=ih this function looks as (2.13)δc(σi)={0fori≤-nc,3σi2+6hncσi+3h2nc22h3nc3for-nc≤i≤0,3σi2-6hncσi+3h2nc22h3nc3for0≤i≤nc,0fori≥nc,
or
(2.14)δc(σi)={0fori≤-nc,3i2+6nci+3nc22hnc3for-nc≤i≤0,3i2-6nci+3nc22hnc3for0≤i≤nc,0fori≥nc.
Denote δic=δc(σi).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 n≥0:(2.15)pin+1-pins=-DhN(ph(n))(AhNph(n+1))i-1T0χZ∖[-2n1,2n1](i)pin+DhN(ph(n))αδic,|i|≤N,n≥0.Bearing in mind thatAhN=(1/(2h)2)AN and denoting by κ=s/(2h)2=s/(4h2), the difference scheme system (2.15) can then be equivalently written in the form (2.16)pin+1+κDhN(ph(n))⋅(ANph(n+1))i=pin-sT0χZ∖[-2n1,2n1](i)pin+sαDhN(ph(n))δic,
which can be expressed in matrix form as(2.17)Mnph(n+1)=Cph(n)+sαDnδhc,
where (2.18)Mn=I2N+1+κDnAN,C=I2N+1-sT0χn1N,Dn=DhN(ph(n)),δhc=(δic).
Observe that the matrix C depends on the time step s, while Mn depends on the steps h and s. It is easy to check that Mn is strictly diagonally dominant for all n∈ℕ. In particular, this fact guarantees that Mn is nonsingular and the system (2.17) has a unique solution for every n>0. Recall also the following result.Theorem 2.1 (see [5]).
Assume that0<s≤T0. If ph0≥0 and ∥ph0∥lS1N=1, then ph(n)≥0 and ∥ph(n)∥lS1N=1, for all n≥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 T0=1. We shall consider two cases depending whether α>0.5 or α≤0.5 with different initial conditions p0, which may or not have support intersecting the interval [-1,1]. We shall not plot cases when the support of p0 is totally contained in [-1,1], since in such a case p0 would be a fixed point of the discrete system (2.15). Therefore, we will focus on cases such that supp(p0)⋂(ℝ∖[-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
(3.1)pf(σ)={D2αe(1+σ)/Dforσ≤-1,12ασ+D+12αfor-1≤σ≤0,-12ασ+D+12αfor0≤σ≤1,D2αe(1-σ)/Dforσ≥1,
where D=D(pf) satisfies the equation D+D=α-(1/2) in order for the normalization constraint ∫ℝpf=1 to hold true. Furthermore, this fixed point is stable (see [6]).
In the next experiments, we shall take the constantsα=0.8 and T0=1 and we define the spacial and time mesh size Δσ=h=0.001 and Δt=s=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 conditionp0 with support out of [-1,1] given by
(3.2)p0(σ)={0forσ<-2.5,0.75for-2.5≤σ≤-1.5,0for-1.5<x<1.5,0.25for1.5≤σ≤2.5,0forσ>2.5.
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 conditionp0 with support in and out of [-1,1] given by
(3.3)p0(σ)={0forσ<0,0.5for0≤σ≤2,0forσ>2.
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
(3.4)p0(σ)={0forσ<-2,316σ2for-2≤σ≤2,0forσ>2.
We obtain a similar behavior as in the above cases (see Table 3 and Figure 5).Subcase 4.
We start from an initial condition verifyingD(p0)≈0 given by
(3.5)p0(σ)={0forσ<-0.3,1.65for-0.3≤σ≤0.3,0for0.3<σ<1.5,0.02for1.5≤σ≤2,0forσ>2.
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 p for (2.1) with support contained in [-1,1].Subcase 5.
As in Subcase4, consider an initial condition verifying D(p0)≈0 given by
(3.6)p0(σ)={0forσ<-1,0.49995for-1≤σ≤1,0for1<σ<1.5,0.0001for1.5≤σ≤2.5,0forσ>2.5.
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: Subcase1
Time(=tn)Dn∥pf–ph(n)∥lS1N∥ph(n)∥lS1N00.81.0565210.10.7054580.81650510.20.6227590.68483510.30.5547320.59054310.50.4529660.460410.80.3531250.33999211.30.2558660.2282812.30.1634150.12748416.30.07781010.02935511500.0582590.001745921Table 2
Numerical results: Subcase2.
Time(=tn)Dn∥pf–ph(n)∥lS1N∥ph(n)∥lS1N00.3997330.70745810.10.3642710.59509410.30.3094020.46901110.60.2529060.35034411.10.194970.234926120.1400350.13347515.20.08014110.03539121500.0582590.001744381Table 3
Numerical results: Subcase3.
Time(=tn)Dn∥pf–ph(n)∥lS1N∥ph(n)∥lS1N00.69990.91209610.10.6003630.69584910.20.5253820.58137610.40.4215930.43918710.70.3266350.31887811.20.237970.21357312.30.1499040.11472918.20.06809220.01541311500.05825910.001745391Table 4
Numerical results: Subcase4.
Time(=tn)Dn∥ph(500)–ph(n)∥lS1N∥ph(n)∥lS1N00.0080.39705610.10.00720.26250510.30.0058320.17232910.90.003099360.07183971501.18676e-2501Table 5
Numerical results: Subcase5.
Time(=tn)Dn∥pf–ph(n)∥lS1N∥ph(n)∥lS1N08e−050.27068210.60.0276150.21778211.30.05890460.12957913.10.07392960.03684121500.05825890.001736311Figure 3
Graphics sequence: Subcase1.Figure 4
Graphics sequence: Subcase2.Figure 5
Graphics sequence: Subcase3.Figure 6
Graphics sequence: Subcase4.Figure 7
Graphics sequence: Subcase5.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 T0=1. We define the spacial and time mesh size Δσ=h=0.001 and Δt=s=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 conditionp0 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 conditionp0 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 conditionp0 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: Subcase1.
Time(=tn)Dn∥ph(500)–ph(n)∥lS1N∥ph(n)∥lS1N00.31.1409510.10.2687920.94535210.20.2402640.81759110.30.2148260.71676410.50.1727950.56454110.70.1404510.454826110.1048490.33921911.40.07304290.24002612.10.04111770.14477114.30.009298530.05143361509.89956e−0601Table 7
Numerical results: Subcase2.
Time(=tn)Dn∥ph(500)–ph(n)∥lS1N∥ph(n)∥lS1N00.14990.7816710.10.1358350.68245610.30.11260.55450510.60.08641830.424413110.06222930.31118811.60.03961640.20914612.80.01790620.1122251110.0006976780.01474861501.42e−0501Table 8
Numerical results: Subcase3.
Time(=tn)Dn∥ph(500)–ph(n)∥lS1N∥ph(n)∥lS1N00.2624631.0287510.10.2316990.85225710.20.2055210.73888510.40.1639060.57704310.60.1328270.4644710.90.09924770.34840211.30.06957070.24995812.10.03693940.1457414.50.008060680.05133281501.18855e−0501Figure 8
Graphics sequence: Subcase1.Figure 9
Graphics sequence: Subcase2.Figure 10
Graphics sequence: Subcase3.


## 4. Stability

Again we takeT0=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(4.1)θ(σ)={1(1-a)2(σ+a)2for-1≤σ≤-a,0for-a≤σ≤a,1(1-a)2(σ-a)2fora≤σ≤1,1≤θ≤2for1≤|σ|≤1+b,2for|σ|≥1+b,
where 0<a<1,b>0, θ(σ) is smooth on |σ|>a and |θ′′(σ)|≤2/(1-a)2 for any σ. Sometimes, total or partial derivatives with respect to σ be denoted by subindices like, for example, in θσσ and pσ below.Lemma 4.1.
Denoteη=1/(1-a)2 and assume that α<1/2η. Let p(·) be the unique solution corresponding to the initial data p0 satisfying (2.2). Then
(4.2)∫0t∫|σ|>1p(t,σ)dσdt≤11-2ηα∫Rp0(σ)θ(σ)dσ,∀t≥0.Proof.
Letp(·) be the unique solution corresponding to the initial data p0 such that D(p0)>0. Multiplying (2.1) by θ we have
(4.3)ddt∫Rp(t,σ)θ(σ)dσ+∫R1R∖[-1,1](σ)p(t,σ)θ(σ)dσ=D(p(t))∫Rp(t,σ)θσσ(σ)dσ≤2ηD(p(t))=2ηα∫|σ|>1p(t,σ)dσ.
Hence,
(4.4)ddt∫Rp(t,σ)θ(σ)dσ+(1-2ηα)∫|σ|>1p(t,σ)dσ≤0,(1-2ηα)∫0t∫|σ|>1p(t,σ)dσdt≤∫Rp(0,σ)θ(σ)dσ.
These calculations are formal, but can be justified by using a cutoff function.Letf(σ):[0,1]→ℝ+ be such that f∈C2([0,1]), f′′(σ)≥0, f′(σ)≤0, f(a)=0, and ∫0af(σ)dσ=1/2. We define v(σ) by (4.5)v(σ)={f(σ)for0≤σ≤a,f(-σ)for-a≤σ≤0,0for|σ|≥a.
We note that v∈H1(ℝ). As D(v)=0 and ∫ℝf(σ)dσ=1, we obtain that v is a fixed point.Denote by∥·∥ the norm in L2(ℝ).Theorem 4.2.
Letα<(1-a)2/2. If
(4.6)(-2αfσ(0)+1)28α-f(0)≤0,
then the fixed point v is stable with respect to initial data p0 satisfying (2.2) in the following sense: for any ɛ>0 there exists δ>0 such that if ∥p0–v∥<δ and ∫ℝ|p0(σ)–v(σ)|dσ<δ, then
(4.7)‖p(t)–v‖<ɛ∀t≥0.Proof.
Letp(·) be the unique solution corresponding to the initial data p0 such that D(p0)>0. Denote z(t)=p(t)–v. Multiplying (2.1) by z, we have
(4.8)12ddt‖z‖2=α∫|σ|>1p(t,σ)dσ(-∫R|pσ|2dσ+∫Rpσvσdσ)-∫|σ|>1|p(t,σ)|2dσ+∫|σ|>1p(t,σ)dσ(p(t,0)-f(0)).
We use that(4.9)∫Rpσvσdσ=fσ(a)p(t,a)-fσ(0)p(t,0)-∫0apfσσdσ+fσ(a)p(t,-a)-fσ(0)p(t,0)-∫-a0pfσσdσ≤fσ(a)(p(t,a)+p(t,-a))-2fσ(0)p(t,0).
Denote b=-fσ(0), c=f(0). Thus
(4.10)12ddt‖z‖2+∫|σ|>1|p(t,σ)|2dσ≤∫|σ|>1|p(t,σ)|2dσ((2αb+1)p(t,0)-α∫R|pσ|2dσ+αfσ(a)(p(t,a)+p(t,-a))-c).
Now(4.11)p(t,0)=p(t,-a)+∫-a0pσ(t,σ)dσ,p(t,0)=p(t,a)-∫0apσ(t,σ)dσ,
so
(4.12)(2αb+1)p(t,0)≤(2αb+1)2(p(t,-a)+p(t,a))+(2αb+1)2∫-aa|pσ(t,σ)|dσ≤(2αb+1)2(p(t,-a)+p(t,a))+(2αb+1)2(∫-aa|pσ(t,σ)|2dσ)1/2(2a)1/2≤(2αb+1)2(p(t,-a)+p(t,a))+α∫-aa|pσ(t,σ)|2dσ+(2αb+1)28α.
Hence(4.13)12ddt‖z‖2+∫|σ|>1|p(t,σ)|2dσ≤∫|σ|>1|p(t,σ)|2dσ((2αb+1)28α-c-α∫|σ|>a|pσ|2dσ+(αfσ(a)+(2αb+1)2)(p(t,a)+p(t,-a))).
Denote d=αfσ(a)+(2αb+1)/2. Using (4.6), (4.13), and the inequality |u(x)|2≤2∥u∥L2(I)∥uσ∥L2(I), for any x∈I (an interval of ℝ) and u∈H1(I) (see [7, page 129]), we have
(4.14)dαp(t,-a)≤dα2(∫-∞-a|p|2dσ∫-∞-a|pσ|2dσ)1/4≤∫-∞-a|pσ|2dσ+34(dα)4/3(∫-∞-a|p|2dσ)1/3,dαp(t,a)≤dα2(∫a∞|p|2dσ∫a∞|pσ|2dσ)1/4≤∫a∞|pσ|2dσ+34(dα)4/3(∫a∞|p|2dσ)1/3,12ddt‖z‖2+∫|σ|>1|p(t,σ)|2dσ≤3α2(dα)4/3∫|σ|>1|p(t,σ)|dσ(∫|σ|>a|p(t,σ)|2dσ)1/3.
It is known [4] that for any R>0 there exists N(R) such that
(4.15)‖p(t)‖≤N(R),∀t≥0,
if ∥p0∥≤R. Hence by Lemma 4.1, we have
(4.16)‖p(t)-v‖2=‖z(t)‖2≤‖z(0)‖2+3α(dα)4/3(N(R))2/3∫0t∫|σ|>1|p(t,σ)|dσ≤‖p0-v‖2+3α-1/3d4/3(N(R))2/31-2ηα∫Rp0(σ)θ(σ)dσ.
We note that∫ℝv(σ)θ(σ)dσ=0. Hence, for any γ>0 there exists δ>0 such that ∫ℝ|p0(σ)-v(σ)|dσ<δ implies
(4.17)∫Rp0(σ)θ(σ)dσ≤∫R|p0(σ)-v(σ)|θ(σ)dσ≤2∫R|p0(σ)-v(σ)|dσ<γ.
From here the stability follows.Let us consider now particular examples of stable fixed points.Example 4.3.
Takeα=0.1 and
(4.18)v(σ)={A+Bσ+Cσ-Dfor0≤σ≤12,A-Bσ+C-σ-Dfor-12≤σ≤0,0for|σ|≥12,
where C=0.05,D=-0.01, and B=8C/(D-1)+8C*log(|(1-2D)/(2D)|)-4, A=2C/(2D-1)-B/2. Here a=0.5 and α<(1-a)2/2=1/8.
This function satisfies∫01/2vdσ=1/2, v(1/2)=0, and
(4.19)vσ(0+)=-6.5536,v(0)=8.1537,K=(1+2α(-vσ(0+)))2-8αv(0)=-1.1836,
so that (4.6) holds. Hence, by Theorem 4.2,  v is stable in the above sense.
Taking an initial condition close tov, for instance,
(4.20)p0(σ)={v(σ-a)forσ≤0,v(σ+a)for0≤σ<1.5,cfor1.5≤σ≤2,0forσ>2,
where a=0.001andc=0.02238, we obtain the numerical results collected in Table 9 (see also Figure 11).
If we take an initial condition closer tov than before, for instance,
(4.21)p0(σ)={v(σ-a)forσ≤0,v(σ+a)for0≤σ<1.5,cfor1.5≤σ≤2,0forσ>2,
where a=0.0001,andc=0.00057414, then we obtain the numerical results collected in Table 10 (see also Figure 12).Table 9
Numerical results: Stable fixes point.
Time(=tn)Dn∥v–ph(n)∥lS1N∥ph(n)∥lS1N00.0011190.042046610.10.00100710.15299710.40.0007341760.24172714.11.48858e−050.315771501.64377e−260.3166671Table 10
Numerical results: Stable fixed point.
Time(=tn)Dn∥v–ph(n)∥lS1N∥ph(n)∥lS1N02.87e−050.004168631504.21594e−280.09492171Figure 11
Graphics sequence: Example4.3.Figure 12
Graphics sequence: Example4.3.Hence, the solution remains near this fixed point if the initial data is close to it.Example 4.4.
Takeα=0.1 and
(4.22)v(σ)={A+Bσ+Cσ-Dfor0≤σ≤12,A-Bσ+C-σ-Dfor-12≤σ≤0,0for|σ|≥12,
where C=0.1,D=-0.01 and B=8(C/(D-1))+8C*log(|((1-2D)/(2D))|)-4, A=2C/(2D-1)-B/2. Here a=0.5 and α<(1-a)2/2=1/8.
This function satisfies∫01/2vdσ=1/2, v(1/2)=0 and
(4.23)vσ(0+)=-9.1071,v(0)=14.3075,K=(1+2α(-vσ(0+)))2-8αv(0)=-3.4855,
so that (4.6) holds.
Hence, by Theorem4.2,  v is stable.Example 4.5.
Takeα=0.3 and
(4.24)v(σ)={A+Bσ+Cσ-Dfor0≤σ≤0.2,A-Bσ+C-σ-Dfor-0.2≤σ≤0,0for|σ|≥0.2,
where C=0.11,D=-0.005, and B=2(C/(a(D-a)))+2(C/a2)*log(|(a-D)/D|)-1/a2, A=(C/(D-a)))-aB. Here a=0.2, and α<(1-a)2/2=0.32.
This function satisfies∫00.2vdσ=1/2, v(0.2)=0 and
(4.25)vσ(0+)=-10.0512,v(0)=23.4517,K=(1+2α(-vσ(0+)))2-8αv(0)=-6.8529,
so that (4.6) holds. Hence, by Theorem 4.2  v 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
(4.26)u(σ)={0forσ<-1,0.5for-1≤σ≤1,0forσ>1,
and the initial condition
(4.27)p0(σ)={0forσ<-1,c1for-1≤σ≤1,0for1<σ<1.5,c2for1.5≤σ≤2.5,0forσ>2.5,
where the constants c1 and c2 are chosen in such a way that p0 is close enough. For c1=0.495 and c2=0.01 we obtain the numerical results given in Table 11 (see Figure 13). For c1=0.49995 and c2=0.0001 the numerical results are given in Table 12 (see Figure 14). These numerical results show that the fixed point pf is unstable.Table 11
Numerical results: Unstable fixed point.
Time(=tn)Dn∥u–ph(n)∥lS1N∥ph(n)∥lS1N00.0010.012247411.40.001755080.10296815.80.0005450870.1831221503.77079e−060.2301021Table 12
Numerical results: Unstable fixed point.
Time(=tn)Dn∥u–ph(n)∥lS1N∥ph(n)∥lS1N01e−050.00012247412.80.001109660.10146119.90.0002748490.1837571505.14826e-060.216861Figure 13
Graphics sequence: Unstable fixed point.Figure 14
Graphics sequence: Unstable fixed point.


## 5. Conclusions

For problem (2.1) with T0=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.



---

*Source: 415921-2011-05-29.xml*