- Email: [email protected]
Some applications of Filippov’s dynamical systems
Journal of Computational and Applied Mathematics 254 (2013) 132–143
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Some applications of Filippov’s dynamical systems Martin Biák, Tomáš Hanus, Drahoslava Janovská ∗ Department of Mathematics, Institute of Chemical Technology, Technická 5, 16628 Prague 6, Czech Republic
article
info
Article history: Received 28 September 2012 Received in revised form 8 March 2013 Dedicated to Prof. Klaus Böhmer on his 75th birthday.
abstract The Filippov systems theory is applied to selected problems from biology and chemical engineering. In particular, we explore (a) new formulation of Bazykin’s ecological model, (b) a predator–prey model with human intervention, (c) an ideal closed gas–liquid system including DAE formulation of the system with a chemical reaction. We illustrate the theory by simulations of the behavior of the particular systems. © 2013 Elsevier B.V. All rights reserved.
MSC: 37N30 34C23 65P30 Keywords: Filippov systems theory Sliding bifurcations DAE formulation Ecological models Closed gas–liquid system System with a chemical reaction
1. Introduction Discontinuous dynamical systems arise in a large number of applications, [1,2]. We concentrate on applications in biology and chemistry. Predator–prey models govern many phenomena in population dynamics, immunology, medicine etc. We assume that there are only two competing species: one species (predator) feeds on another species (prey), which in turn feeds on other things. The model is a variation of the Lotka–Volterra system, a particular case of the Bazykin model, [3,4]. As a Filippov system, we formulate a population model with a human intervention, in particular we consider an ecosystem that consists of spiders that prey upon insects. If the population of the insects exceeds the given limit although the spiders are active, an insecticide is spread on the vineyard, [5]. The model is analyzed according to the theory described in [2,6]. In the originally introduced model of the ideal closed gas–liquid system, [7], the solution depends on ten parameters, [8]. The dimensionless formulation [9] significantly reduces the number of parameters to only four. Another obvious advantage of the dimensionless formulation is its independence on the model scale. We study and simulate the solution of the dimensionless problem subjected to a given parameter set, [10–13]. Similarly as the system of ODEs, also the system of DAEs can be reformulated as a Filippov system. As an example, we studied the soft-drink system, i.e. the system with a chemical reaction, [14,15]. We illustrate the theory by simulations of the behavior of the dynamical systems. All simulations are performed in Matlab, in particular we use a modified version of the program developed by Petri T. Piiroinen and Yuri A. Kuznetsov, see [13].
∗
Corresponding author. Tel.: +420 731577040; fax: +420 233 323 917. E-mail addresses: [email protected] (M. Biák), [email protected] (T. Hanus), [email protected] (D. Janovská).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.03.034
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
133
2. Preliminaries Let the state space S ⊂ Rn be divided into a system of m subsets Si , i = 1, . . . , m, such that Si ∩ Sj = ∅ for i, j = 1, . . . , m, i ̸= j, m
Si = S ,
Σij = Σji = Si ∩ Sj ,
i ̸= j,
Σ=
i =1
m m
Σij ,
i=1 j=i+1
where Σij is the boundary between sets Si and Sj , Σ is the boundary of the whole state space S. Each set Si is equipped with a continuous and smooth vector field f(i) : Rn → Rn . Let Z be a system of ordinary differential equations dx
Z :
dt
= f(i) (x),
x ∈ Si , i = 1, . . . , m.
(1)
System (1) is called continuous, if lim f(i) (u) = lim f(j) (v),
u→x
v→ x
∀x ∈ Σij , ∀i ∀j, i ̸= j,
and it is called a system with discontinuous right hand side, a Filippov system, if there exists x ∈ Σij for some i ̸= j such that lim f(i) (u) ̸= lim f(j) (v).
u→x
v→ x
Let us restrict ourselves to the state space S = R2 which is divided only into two regions S1 and S2 . Let the boundary Σ := Σ12 consist of parts, which are smooth curves topologically equivalent to a circle or a line. 2.1. Generic Filippov systems Let F be a planar dynamical system f(1) (x), f(2) (x),
F : x˙ =
x ∈ S1 , x ∈ S2 .
(2)
where S1 = {x ∈ S ⊂ R2 : H (x) < 0}, S2 = {x ∈ S ⊂ R2 : H (x) > 0}, and the function H : S → R has a non-vanishing gradient Hx on Σ = {x ∈ S , H (x) = 0}. Let us define the function σ : Σ −→ R as
σ (x) = ⟨Hx (x), f(1) (x)⟩ · ⟨Hx (x), f(2) (x)⟩, where ⟨·, ·⟩ means the dot product in R2 . We will define the vector field g on the boundary Σ , see [2]. 1. At these points x ∈ Σ in which σ (x) > 0 we define the vector field g as an arithmetic mean of the vectors f(1) and f(2) g := gA (x) =
1 (1) (f (x) + f(2) (x)),
2
x ∈ Σ.
In this case we say that x belongs to the crossing set Σc ⊂ Σ , where Σc = {x ∈ Σ : σ (x) > 0}.
The subset Σc of Σ contains these points x ∈ Σ , in which both vectors f(1) and f(2) point into the same region.
134
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
Fig. 1. Stable sliding trajectory on the left. On the right: Time reversed ⇒ the trajectory is unstable, different from the original one.
2. In all other cases of configuration of the vectors f(1) and f(2) we apply the so-called Filippov’s method and define the vector field g on the boundary Σ as a convex combination of the vector fields f(1) and f(2) , g := gF (x) = λ(x)f(1) (x) + (1 − λ(x))f(2) (x),
λ : Σ −→ [0, 1],
λ(x) =
x ∈ Σ , where
(2)
Hx (x) · f Hx (x) ·
(f(2) (x)
(x) , − f(1) (x))
(3)
Hx (x) is the normal vector at the point x ∈ Σ . In this case we say that x belongs to the sliding set Σs . The set Σs is a complement of Σc in Σ , i.e. Σs = {x ∈ Σ : σ (x) ≤ 0}.
On the sliding boundary Σs special points (sliding points) can be detected. Let us classify some of them. 1. Singular sliding point is such point x ∈ Σs that
⟨Hx (x), f(1) (x)⟩ = 0 and also ⟨Hx (x), f(2) (x)⟩ = 0. At these points, both vectors f(1) (x) and f(2) (x) are tangent to Σs . 2. The point x ∈ Σs is a generic pseudoequilibrium if f(1) (x) ̸= 0,
g(x) = 0,
f(2) (x) ̸= 0.
At these points, the vectors f(1) (x) and f(2) (x) are anti-collinear. 3. In a boundary equilibrium x ∈ Σs , one of the vectors f(i) (x) vanishes, f(1) (x) = 0 or f(2) (x) = 0. 4. The point x ∈ Σs is a tangent point if both f(1) (x) ̸= 0, f(2) (x) ̸= 0 and
⟨Hx (x), f(1) (x)⟩ = 0 or ⟨Hx (x), f(2) (x)⟩ = 0. In this case, both vectors f(1) (x), f(2) (x) are nonzero, but one of them is tangent to Σ . The tangent point terminates Σs in Σ . Let us note that in general, in Filippov systems the theorem on existence and uniqueness does not apply (see Fig. 1). Two Filippov systems are called topologically equivalent if there is a homeomorphism h : R2 −→ R2 – that maps the discontinuity boundary of one system onto the discontinuity boundary of the other one – orbits of one system are mapped onto the corresponding orbits of the other one with time orientation preserved – crossing and sliding segments of any orbit are mapped onto the corresponding segments of its image. Remark 1. The another possibility how to define the vector field g on the boundary Σ is to apply the so-called Utkin’s equivalent control method, see [2]. In the Utkin’s method the vector field g on the boundary Σ is defined as g := gU (x) =
1 (2) 1 (1) f (x) + f(2) (x) + f (x) − f(1) (x) β(x), 2
2
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
135
Table 1 Defining systems on the boundary of the discontinuity H (x, α) = 0. Discontinuity boundary Crossing boundary Sliding boundary Tangent points of the vector field fi Pseudoequilibrium Boundary equilibrium
H (x, α) = 0 ⟨Hx , f(1) ⟩⟨Hx , f(2) ⟩ > 0 ⟨Hx , f(1) ⟩⟨Hx , f(2) ⟩ ≤ 0 ⟨Hx (x, α), f(i) (x, α)⟩ = 0 λ1 f(1) (x, α) + λ2 f(2) (x, α) = 0, λ1 + λ2 − 1 = 0 f(i) (x, α) = 0
where the control function β : Σ −→ [−1, 1] has the form
β(x) = −
⟨Hx (x), f1 (x) + f2 (x)⟩ . ⟨Hx (x), f2 (x) − f1 (x)⟩
(4)
If we put β(x) := 1 − 2λ(x) in (4), we obtain the Eq. (3), i.e. the Utkin’s control method and the Filippov’s method are algebraically equivalent.
2.2. Parameter-dependent Filippov systems In piecewise-smooth systems the interaction of invariant sets with the ‘‘switching’’ boundary often produces bifurcations that are forbidden in smooth systems. These are known as discontinuity-induced bifurcations. We will study those bifurcations that are likely to occur as a single parameter is varied. Let us have a piecewise-smooth dynamical system dependent on a parameter α ,
x˙ =
f(1) (x, α), f(2) (x, α),
x ∈ S1 (α), x ∈ S2 (α),
(5)
x ∈ R2 , α ∈ R, and f(i) : R2 × R −→ R2 , i = 1, 2, are smooth functions of (x, α), S1 (α) = x ∈ R2 : H (x, α) < 0 ,
S2 (α) = x ∈ R2 : H (x, α) > 0 ,
Σ (α) = {x ∈ S : H (x, α) = 0},
H (x, α) is a smooth function such that Hx (x, α) ̸= 0 for all (x, α) ∈ Σ (α), where Σ (α) is the boundary of discontinuity (‘‘switching’’ boundary). We say that the system (5) exhibits a bifurcation at the parameter value α = α0 if any arbitrarily small parameter perturbation results in a topologically nonequivalent system. Bifurcations that occur within the single region S1 or S2 for all parameter values of interest can be studied using smooth bifurcation theory. For trajectories crossing the boundary one can adopt the same theory, [2,6]. We will be interested in the discontinuity-induced bifurcations, [2,1]. They must occur if we continuously vary parameters. A piecewise-smooth system is piecewise-structurally stable if there exists an ε > 0 such that all C 1 perturbations of maximum size ε of the vector field f lead to topological equivalent phase portraits. We say that a discontinuity-induced bifurcation occurs at a parameter value at which a piecewise-smooth system is not piecewise-structurally stable, i.e., there exists an arbitrarily small perturbation that leads to a system that is not piecewise-topological equivalent. The simplest kind of the discontinuity-induced bifurcations occurs when an equilibrium point lies precisely on a discontinuity boundary. In Filippov systems there might occur so called pseudoequilibria, which are equilibria of the sliding flow but are not equilibria of any of the vector fields of the original system. For examples of the defining equations, see Table 1. Let us mention only three types of local sliding bifurcations that occur on the sliding set: – Collision of an equilibrium, for example a node, with the boundary Σ . – Collision of tangent points (double tangency). – Collision of pseudoequilibria, for example pseudo-saddle–node bifurcation.
136
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
Fig. 2. Trajectories of the system in dependence on the parameter α .
Example 1. Let us simulate the behavior of the trajectories of the following parameter-dependent piecewise-smooth dynamical system. Under a small perturbation a node collides with the boundary Σ . f(1) (x, α) =
−3x1 − x2 , −x1 − 3x2
f(2) (x, α) =
0 , 1
H (x, α) = x2 + α.
In Fig. 2, for α < 0 there is a node Xα in the region S1 . The sliding region Σs goes from the left and ends at the tangent point Tα ∈ Σ . The red dashed part of the Σ represents the crossing boundary Σc . The trajectories that start on the left of Tα in S1 are heading to Σs , slide on it up to the tangent point and then they are attracted by the node. Trajectories that start in S2 tend towards the boundary Σ . To the left of the tangent point they slide to Tα and end in the node. If they go to the right of the tangent point, they cross the crossing boundary and again end in the node. If we increase α the node moves to the boundary and for α = 0 collides with it, X0 ≡ Tα , X0 becomes a pseudonode Pα . Then if α further increases the pseudonode Pα stays on the sliding boundary Σs and the tangent point Tα secedes from it and moves to the right again marking the end of the sliding set. 3. Ecological applications 3.1. Bazykin’s ecological model Bazykin’s ecological model is the dynamical system defined by the following differential equations: dx1
ax1 x2
= x1 − x21 −
dt dx2
= −dx2 +
dt
b + x1 ax1 x2
b + x1
,
,
x1 ∈ ⟨0, ∞), x2 ∈ ⟨0, ∞), a > 0, b > 0, d > 0.
We reformulate this system as a Filippov system. We divide the state space S = ⟨0, ∞) × ⟨0, ∞) ⊂ R2 into two open regions S1 , S2 separated by a boundary Σ defined by a smooth parameter dependent function H. ax1 x2
dx dt
x1 − x21 −
= f(1) (x) =
ax1 x2 b + x1
dx dt
b + x1 ,
− dx2
x1 − x21 −
x ∈ S1 ,
ax1 x2
b + x1 , = f(2) (x) = ax1 x2 − dx2 − Ex2 b + x1
H (x, α) = x2 − α,
(6)
x ∈ S2 ,
(7)
x ∈ Σ,
(8)
where x1 > 0 is the mass of the prey, x2 > 0 is the mass of the predator, α, a, b, d, E are parameters, and S1 (α) = x ∈ R2 : H (x, α) < 0 ,
S2 (α) = x ∈ R2 : H (x, α) > 0 ,
Σ (α) = {x ∈ S : H (x, α) = 0}.
Let us note that f(1) is a continuous vector field on the open set S1 , f(2) on S2 . On the boundary Σ , we define
x˙ =
gA (x), gF (x),
x ∈ Σc = {x ∈ Σ : σ (x) > 0}, x ∈ Σs = {x ∈ Σ : σ (x) ≤ 0},
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
137
where σ (x) = ⟨Hx (x), f(1) (x)⟩ · ⟨Hx (x), f(2) (x)⟩, ⟨·, ·⟩ is the dot product in R2 and Hx (x) is the gradient of H at the point x ∈ Σ. We performed a one parameter study for a continuation parameter α ∈ R; the values of the other parameters were fixed. We apply the theory of Filippov systems and compute generic pseudoequilibria, i.e. such points on Σs in which the vectors f(1) (x) and f(2) (x) are anticolinear. For λ ∈ (0; 1), f(1) (x, α) ̸= 0, f(2) (x, α) ̸= 0, we have to solve the following system of algebraic equations: H (x, α) = 0,
λf(1) (x, α) + (1 − λ)f(2) (x, α) = 0.
Tangent points, i.e. ending points of the sliding set, are such points in which both vectors f(i) (x) are nonzero, but one of them is tangent to Σ . The defining systems for the tangent points are H (x, α) = 0
⟨Hx (x, α), f(1) (x, α)⟩ = 0 for f(1) tangent to Σ , ⟨Hx (x, α), f(2) (x, α)⟩ = 0 for f(2) tangent to Σ .
and
For more details, see [2,4,16]. 3.2. Population model: spiders on the vineyard Population models with a human intervention can be formulated as Filippov systems and then analyzed according to the theory described in [2]. The population of the spiders that predate on the insects in the vineyards is an example of the situation, when the standard predator–prey model is modified by human intervention if one of the populations exceeds the given limit. The ecosystem we consider consists of the spiders populations s(t ) and the insects on which they prey. Let f (t ) be the insects population living in the woods and the green patches bordering the cultivated land, and let v(t ) denote the insects having the vineyards as habitat. In [5], the following model was proposed f˙ = rf
f
1−
s˙ = s −a +
− csf ,
W kbv
H +v
v˙ = v e −
bs
H +v
+ kcf
,
(9)
,
where k is the conversion factor of prey into new spiders, k < 1, r is the specific growth rate of the prey in the woods, e is the specific growth rate of the prey in the vineyards, W and H are carrying capacities of the woods and the vineyard, respectively. All parameters are positive. Let us note that the first and the third equations describe the evolution of the wood insects and the growth of the insects population in open fields, respectively. The dynamics of spiders is given in the second equation. Now, let us model the human intervention, for example some spraying effect on the ecosystem of the vineyard. We define the function
ϕ(f , s, v) = vm − v,
(10)
where vm is a prescribed limit of the population of the insects on the vineyards. If ϕ(f , s, v) > 0 then there is no spraying needed and the system (9) describes the situation. If the population of the insects exceeds the limit, i.e. ϕ(f , s, v) = vm − v < 0, the spraying starts. Then we have to modify the model accordingly:
f
− csf − h(1 − q)f , s˙ = s −a + + kcf − hKqs, H +v bs v˙ = v e − − hqv, H +v
f˙ = rf
1−
W kbv
(11)
where q is the portion of insecticide sprayed directly on the vineyards, (1 − q) denotes the portion which may accidentally be dispersed in the woods, h stands for effectiveness against the parasites and K models the possibly smaller effect of the spraying on the spiders, 0 < K < 1. The Eqs. (9), (11) together with (10) define the Filippov system F . The boundary Σ of the system F is given by
Σ = {(f , s, v) : ϕ(f , s, v) = 0}.
(12)
The results of the dynamical simulations are summarized in Fig. 3: The trajectory reaches the boundary Σ , slides on it and then again leaves it. The situation can be interpreted as follows. The amount of the insects in the vineyards grows. When the boundary is reached, i.e. vm = v , the insecticide is applied. It takes some time to see the effect of spraying (this is the
138
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
Fig. 3. Spiders on the vineyard. Parameters: s(0) = 0, f (0) = 2, v(0) = 1.5, vm = 3.
time of sliding). But then the amount of insects again decreases. The whole process can be repeated if we set a lower bound to the amount of insects in the vineyards. 4. Closed gas–liquid system A closed gas–liquid system, more precisely an ideal gas–liquid system, see [7], is a simplified model of the slurry reactor. In the slurry reactor, the feed contains a liquid phase (a catalyst in diesel oil) and a gas phase (hydrogen and ethylene), see [17]. The output from the reactor consists of solid HDPE dispersed in diesel oil and unreacted gases. The process runs continuously, because the feed is sent in and also the product is withdrawn continuously. Two different mathematical models are needed to describe this system: one for the liquid output, the other one for the gas output. Only one of these two models takes place at a given time, because the dispersed polymer and gases cannot come out from the reactor simultaneously. Under suitable assumptions, the closed gas–liquid system may serve as a good approximation of the industrial HDPE reactor. Like the slurry reactor, a closed gas–liquid system consists of a closed tank with only one outlet tube. Let us start with a simple system shown in Fig. 4. The feed, that consists of a liquid and a gas, is sent to the closed tank, which has the only one outlet tube. If the level of the liquid is above the outlet tube opening, see Fig. 4(a), the liquid comes out from the tank. Otherwise, the gas comes out, see Fig. 4(b). Let FL and FG be molar inflow rates, G and L be molar outflow rates of the liquid and the gas, respectively. The following assumptions help to simplify the model. 1. 2. 3. 4. 5. 6. 7. 8.
The feed consists of an ideal gas G and an incompressible liquid L. The gas and the liquid do not react. The liquid has negligible vapor pressure at the operating conditions. The gas does not dissolve in the liquid at the operating conditions. In the liquid, there are no gas bubbles that might be leaving the tank. The valve dynamics is ignored. The flow rate through the valve is proportional to the difference of the tank pressure and the outlet pressure. The temperature, the molar inflow rates, the outlet pressure and the valve opening are all constant.
Classically, we need two sets of differential equations to describe the behavior of this system: one for the liquid model, the other one for the gas model. The following model equations of the ideal gas–liquid system were proposed by Moudgalya and Ryali, [7]. Liquid model: ML /ρL > Vd dMG = FG , dt dML = FL − L, dt RT ML + = V, MG P ρL L = kL x(P − Pout ).
Gas model: ML /ρL < Vd dMG = FG − G, dt dML = FL , dt RT ML MG + = V, P ρL G = kG x(P − Pout ).
(13)
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
a
b
139
c
Fig. 4. An ideal gas–liquid system.
Here, ML and MG are the molar hold-ups of the liquid and the gas, ρL is the molar density, T is the absolute temperature, P and Pout are pressures in the tank and in the outlet, respectively. R is the gas constant, V is the volume of the tank and Vd denotes the volume of the liquid below the outlet tube opening in the tank, see Fig. 4(c). The valve opening is denoted by x, 0 < x ≤ 1. The valve coefficients are kL for the liquid and kG for the gas flow. With respect to time, FG and FL have to be positive constant. The valve opening x is a constant, too. Let us point out that the system keeps switching between the liquid and gas models and that the discontinuity is a natural part of the system. Applying the theory of Filippov systems [2,18], both the outgoing liquid phase and the outgoing gas phase can be defined using only one set of ordinary differential equations with a discontinuous right-hand side, [12]. From the third equation in (13) we derive the expression for the pressure P in the tank, P =
MG RT V − ML /ρL
.
(14)
Substituting this P into the fourth equations in (13) we obtain
L = kL x
MG RT V − ML /ρL
− Pout ,
(15)
− Pout .
(16)
and
G = kG x
MG RT V − ML /ρL
The first and second equations in (13) together with (15) and (16) define the equivalent Filippov system F . After a few simplifying assumptions, the behavior of this system depends on ten parameters. We have chosen two of them namely, FG and FL , a molar inflow of both gas and liquid. With respect to time, FG and FL are positive constants. We study the dependence of the solution (ML MG ) (the molar hold-ups of the liquid and gas) on these two parameters. Now we are ready to formulate the system (13) as the following Filippov system F .
F :
d
MG ML
dt
f(1) =
FG − kG x FL
f(1) (MG , ML , FG , FL ), f(2) (MG , ML , FG , FL ),
=
− Pout , V − ML /ρL MG RT
ϕ(MG , ML , FG , FL ) < 0, ϕ(MG , ML , FG , FL ) > 0, FG , MG RT f(2) = FL − kL x − Pout V − ML /ρL
(17)
where
ϕ(MG , ML ) := ϕ(MG , ML , FG , FL ) = ML − ρL Vd , Σ = {(MG , ML ) ∈ R2 : ϕ(MG , ML ) = 0}, 2 S1 = {(MG , ML ) ∈ R : ϕ(MG , ML ) < 0}, S2 = {(MG , ML ) ∈ R2 : ϕ(MG , ML ) > 0}. Here, V is the volume of the tank and Vd , 0 < Vd < V , denotes the volume below the outlet tube opening in the tank. Two nullclines of the system (17) together with the discontinuity boundary Σ split the state space S into four regions with the different qualitative behavior of the system, [7].
140
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
Fig. 5. Dependence of the solution (MG , ML ) on the molar inflow of the gas for FG = 0.1, FL = 2.5. On the right: Zoom of the sliding set.
The following theorem says that the sliding set Σs is an attractor for the whole state space S. Theorem 1. Let the Filippov system be described by Eq. (17). Then for any x ∈ S \ Σs , x = (MG , ML ), every orbit γx ends on the sliding set Σs . For proof, see [7]. We simulated the behavior of the system (17) for different values of the parameters FG , FL , see Fig. 5, and constructed a solution diagram. It shows that the coordinates MG and ML of the pseudoequilibria and tangent points on the sliding and crossing sets depend linearly on the molar inflow of the gas FG and liquid FL . 4.1. Dimensionless closed gas–liquid system It turned out that a transformation of the model equations into a dimensionless form reduces the number of parameters substantially—to only four. It can be proven that the behavior of the dimensionless system is qualitatively the same as the behavior of the system (17), but the reduction of the number of parameters substantially simplifies the bifurcation analysis. The Filippov system (17) depends on 10 parameters, in particular the parameter set is p = (FG , FL , ρL , V , Vd , T , Pout , x, kL , kG ),
p ∈ R10 .
If we define dimensionless state variables, then choose a suitable scaling and simplify the resulting equations we obtain four ˜ K˜ , β˜ L , γ˜ ), q˜ ∈ R4 , see [9] for the detailed transformation. dimensionless parameters q˜ = (δ, The resulting dimensionless Filippov system F˜ then depends only on q˜ . It reads (tildes are omitted): f(1) (MG , ML , q), f(2) (MG , ML , q),
ϕ(MG , ML , q) < 0, ϕ(MG , ML , q) > 0, MG 2 + K β 1 − K δβ L L , f(1) (MG , ML , q) = 1 − ML F (q) :
d
dt
MG ML
=
1
(2)
f
(MG , ML , q) =
1 1 − δβL2
MG 1 − ML
+ βL
,
ϕ(MG , ML , q) = ML − γ . It can be shown that the behavior of the dimensionless system is qualitatively the same as the behavior of the original one, but the reduction of the number of parameters only to four simplifies substantially the bifurcation analysis, [9]. 4.2. Soft drink process The process of manufacturing soft-drink depicted in Fig. 6 is based on the reaction between CO2 and water: CO2 + H2 O → H2 CO3 .
(18)
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
141
To simplify the model, we will suppose that – – – – – –
The system contains only components CO2 , H2 O and H2 CO3 (denoted by indices 1, 2 and 3, respectively). Intermediate ionization reactions and dissociation of H2 CO3 are ignored. In the liquid, there are no gas bubbles. The valve dynamics is ignored. The flow rate through the valve is proportional to the difference of the tank pressure P and the outlet pressure Pout . The temperature T , the molar inflow rates F1 and F2 , the outlet pressure, valve coefficients kG and kL and the valve opening x are all constant.
Let M1 = M1 (t ),
M2 = M2 (t ),
M3 = M3 (t ),
be the total molar hold-ups of CO2 , H2 O and H2 CO3 , respectively. For a fixed t, let us define a scalar function ϕ = ϕ(M1 , M2 , M3 ),
ϕ(M1 , M2 , M3 ) =
M2
ρL
+
M3
ρa
− Vd ,
(19)
where ρL , ρa are molar densities of water and acid, respectively. The volume of the whole tank is equal to V and the part of the volume that is below the opening of the dip tube is denoted as Vd , 0 < Vd < V . Similarly as in [8,9], in the tank two systems take place: the liquid model (the liquid leaves the tank) if ϕ(M1 , M2 , M3 ) > 0 or the gas model (the gas leaves the tank) for ϕ(M1 , M2 , M3 ) < 0. The acid phase consists of H2 CO3 , H2 O and dissolved CO2 while the gas phase contains only CO2 . As a consequence, the acid model is described by 3 ODEs and 6 algebraic equations, the gas model needs also 3 ODEs but only 4 algebraic equations. Let us give the list of these equations. Differential equations: Liquid model : ϕ(M1 , M2 , M3 ) > 0
Gas model : ϕ(M1 , M2 , M3 ) < 0
dM1
dM1
dt dM2 dt dM3 dt
= F1 − L1 − rV ,
dt dM2
= F2 − L2 − rV ,
dt dM3
= −L3 + rV ,
dt
= F1 − G − rV , = F2 − rV , = rV ,
M1 , M2 and M3 are the total molar hold-ups of CO2 , H2 O, H2 CO3 , respectively. F1 and F2 are the molar inflow rates. The molar flow rates of the components through the valve are denoted L1 , L2 and L3 in the liquid model and G in the gas model. The rate r of the reaction (18) is given by r = κc
M1 M2 V2
,
where κc is the rate constant.
(20)
Algebraic equations: Liquid model : ϕ(M1 , M2 , M3 ) > 0
Gas model : ϕ(M1 , M2 , M3 ) < 0
0 = M1 − (Mℓ + Mg ),
0 = M1 − (Mℓ + Mg ),
0=P−
σ Mℓ Mℓ + M2 + M3
0=V− 0= 0=
M1 RT P
+
M2 Mℓ + M2 + M3 M3 Mℓ + M2 + M3
M2
ρL − −
, +
0=P− M3
ρa
,
L2 L1 + L2 + L3 L3 L1 + L2 + L3
σ Mℓ Mℓ + M2 + M3
0=V−
,
M1 RT P
+
M2
ρL
, +
M3
ρa
,
0 = G − kG x(P − Pout ),
,
0 = L1 + L2 + L3 − kL x(P − Pout ), P and T means pressure and temperature in the tank, the hold-ups of CO2 in liquid and gas are denoted Mℓ and Mg , R is a gas constant and σ is Henrys constant for CO2 . The straightforward computation shows that both the system of DAEs for the gas mode and the system of DAEs for the liquid model has index 1, [15].
142
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
a
b
Fig. 6. The soft-drink process.
Similarly as in the system of ODEs for the ideal gas–liquid model, the system of DAEs describing the soft-drink model can be reformulated by making use of Filippov’s method, [14]. As a result we obtain the formulation of the soft-drink process as Filippov system. The corresponding differential equations are: M1 t M2 t M3 t
= λ(G − L1 ) + F1 − G − rV , = −λL2 + F2 − rV , = −λL3 + rV ,
where λ =
F2 ρa + rV (ρL − ρa ) L2 ρa + L3 ρL
.
Similarly, the corresponding algebraic equations are: 0 = M1 − (Mℓ + Mg ),
σ Mℓ
0=P−
Mℓ + M2 + M3
0=V−
0= 0=
M1 RT P
+
M2 Mℓ + M2 + M3 M3 Mℓ + M2 + M3
M2
ρL − −
, +
M3
ρa
,
L2 L1 + L2 + L3 L3 L1 + L2 + L3
, ,
0 = L1 + L2 + L3 − kL x(P − Pout ), 0 = G − kG x(P − Pout ). Let us note that missing equations in both liquid and gas models that should describe the same phenomenon are replaced by equations 0 = 0. For a particular example, the results are depicted in Fig. 7. 5. Conclusions We considered three different parameter-dependent Filippov’s dynamical systems. The objective was to perform a numerical analysis by dynamical simulations. We used the software [13]. Concerning the ecological applications, we have shown that Filippov’s formulation is very natural (e.g. comparing [5] with our approach). As far as the study of the closed gas–liquid system is concerned, it is just the first step towards to a model of the real HDPE reactor. Acknowledgments The research was supported by MSM 6046137306, the project financed by The Ministry of Education, Youth and Sports, Czech Republic. The authors were also supported from specific university research grant (MSMT no. 21/2010).
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
143
Fig. 7. Soft-drink process: (a)–(c) The integral curves of the state variables M1 , M2 and M3 . (d) The trajectory of the soft-drink system starting at the point (0.72, 95, 0).
References [1] Yu.A. Kuznetsov, Elements of Applied Bifurcation Theory, second ed., Springer, New York, 1998. [2] M. di Bernardo, C.J. Budd, A.R. Champneys, P. Kowalczyk, Piecewise-Smooth Dynamical Systems, Theory and Applications, Springer, London, 2008. [3] A.D. Bazykin, A.I. Khibnik, B. Krauskopf, Nonlinear Dynamics of Interacting Populations, in: World Scientific Series on Nonlinear Science, Series A, vol. 11, 1998. [4] W. Govaerts, Numerical Methods for Bifurcations of Dynamical Equilibria, SIAM, Philadelphia, 2000. [5] E. Venturino, M. Isaia, F. Bona, S. Chatterjee, G. Badino, Biological controls of intensive agroecosystems: wanderer spiders in the Langa Astigiana, Ecological Complexity 5 (2008) 157–164. [6] D.J.W. Simpson, Bifurcations in Piecewise-Smooth Continuous Systems, World Scientific, University of British Columbia, Canada, 2010. [7] K.M. Moudgalya, V. Ryali, Chemical Engineering Science 56 (2001) 3595–3609. [8] M. Biák, D. Janovská, Filippov dynamical systems, in: Seminar on Numerical Analysis & Winter School, Proceedings of the Conference, SNA’09, Institute of Geonics AS CR Ostrava, 2009, pp. 1–4. [9] M. Biák, D. Janovská, A parametric study of the dimensionless closed gas–liquid system, in: Nové Hrady, Hana Bílková (Eds.), Seminar on Numerical Analysis & Winter School, Proceedings of the Conference, SNA’10, Institute of Computer Science AS CR Praha, 2010, pp. 26–28. [10] L.L. Böhm, P. Goebel, P.R. Schöneborn, Detailed reaction engineering as a basis of modern slurry technology for PE-HD-production, Die Angewandte Makromolekulare Chemie 174 (1990) 189–203. [11] F. Dercole, Yu.A. Kuznetsov, ACM Transactions on Mathematical Software 31 (2005) 95–119. [12] A.F. Filippov, Differential Equations with Discontinuous Righthand Sides, Kluwer Academic Publishers, Dordrecht, 1988. [13] P.T. Piiroinen, Yu.A. Kuznetsov, ACM Transactions on Mathematical Software 34 (2008) 1–24. [14] J. Agrawal, K.M. Moudgalya, A.K. Pani, Sliding motion of discontinuous dynamical systems described by semi-implicit index one differential algebraic equations, Chemical Engineering Science 61 (2006) 4722–4731. [15] P. Kunkel, V. Mehrmann, Differential–Algebraic Equations: Analysis and Numerical Solution, European Mathematical Society, 2006. [16] T. Hanus, D. Janovská, Discontinuous, piecewise-smooth dynamical systems, in: CHISA 2006, Praha, Czech Republic, CD-ROM of Full Texts CHISA 2006, Title Index: D, Org.No.: F4.8. [17] K.M. Moudgalya, S. Jaguste, A class of discontinuous dynamical systems II., an industrial slurry high density polyethylene reactor, Chemical Engineering Science 56 (2001) 3611–3621. [18] M. di Bernardo, P. Kowalczyk, A. Nordmark, Bifurcations of dynamical systems with sliding: derivation of normal-form mappings, Physica D (2002) 175–205.
Contents lists available at SciVerse ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Some applications of Filippov’s dynamical systems Martin Biák, Tomáš Hanus, Drahoslava Janovská ∗ Department of Mathematics, Institute of Chemical Technology, Technická 5, 16628 Prague 6, Czech Republic
article
info
Article history: Received 28 September 2012 Received in revised form 8 March 2013 Dedicated to Prof. Klaus Böhmer on his 75th birthday.
abstract The Filippov systems theory is applied to selected problems from biology and chemical engineering. In particular, we explore (a) new formulation of Bazykin’s ecological model, (b) a predator–prey model with human intervention, (c) an ideal closed gas–liquid system including DAE formulation of the system with a chemical reaction. We illustrate the theory by simulations of the behavior of the particular systems. © 2013 Elsevier B.V. All rights reserved.
MSC: 37N30 34C23 65P30 Keywords: Filippov systems theory Sliding bifurcations DAE formulation Ecological models Closed gas–liquid system System with a chemical reaction
1. Introduction Discontinuous dynamical systems arise in a large number of applications, [1,2]. We concentrate on applications in biology and chemistry. Predator–prey models govern many phenomena in population dynamics, immunology, medicine etc. We assume that there are only two competing species: one species (predator) feeds on another species (prey), which in turn feeds on other things. The model is a variation of the Lotka–Volterra system, a particular case of the Bazykin model, [3,4]. As a Filippov system, we formulate a population model with a human intervention, in particular we consider an ecosystem that consists of spiders that prey upon insects. If the population of the insects exceeds the given limit although the spiders are active, an insecticide is spread on the vineyard, [5]. The model is analyzed according to the theory described in [2,6]. In the originally introduced model of the ideal closed gas–liquid system, [7], the solution depends on ten parameters, [8]. The dimensionless formulation [9] significantly reduces the number of parameters to only four. Another obvious advantage of the dimensionless formulation is its independence on the model scale. We study and simulate the solution of the dimensionless problem subjected to a given parameter set, [10–13]. Similarly as the system of ODEs, also the system of DAEs can be reformulated as a Filippov system. As an example, we studied the soft-drink system, i.e. the system with a chemical reaction, [14,15]. We illustrate the theory by simulations of the behavior of the dynamical systems. All simulations are performed in Matlab, in particular we use a modified version of the program developed by Petri T. Piiroinen and Yuri A. Kuznetsov, see [13].
∗
Corresponding author. Tel.: +420 731577040; fax: +420 233 323 917. E-mail addresses: [email protected] (M. Biák), [email protected] (T. Hanus), [email protected] (D. Janovská).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.03.034
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
133
2. Preliminaries Let the state space S ⊂ Rn be divided into a system of m subsets Si , i = 1, . . . , m, such that Si ∩ Sj = ∅ for i, j = 1, . . . , m, i ̸= j, m
Si = S ,
Σij = Σji = Si ∩ Sj ,
i ̸= j,
Σ=
i =1
m m
Σij ,
i=1 j=i+1
where Σij is the boundary between sets Si and Sj , Σ is the boundary of the whole state space S. Each set Si is equipped with a continuous and smooth vector field f(i) : Rn → Rn . Let Z be a system of ordinary differential equations dx
Z :
dt
= f(i) (x),
x ∈ Si , i = 1, . . . , m.
(1)
System (1) is called continuous, if lim f(i) (u) = lim f(j) (v),
u→x
v→ x
∀x ∈ Σij , ∀i ∀j, i ̸= j,
and it is called a system with discontinuous right hand side, a Filippov system, if there exists x ∈ Σij for some i ̸= j such that lim f(i) (u) ̸= lim f(j) (v).
u→x
v→ x
Let us restrict ourselves to the state space S = R2 which is divided only into two regions S1 and S2 . Let the boundary Σ := Σ12 consist of parts, which are smooth curves topologically equivalent to a circle or a line. 2.1. Generic Filippov systems Let F be a planar dynamical system f(1) (x), f(2) (x),
F : x˙ =
x ∈ S1 , x ∈ S2 .
(2)
where S1 = {x ∈ S ⊂ R2 : H (x) < 0}, S2 = {x ∈ S ⊂ R2 : H (x) > 0}, and the function H : S → R has a non-vanishing gradient Hx on Σ = {x ∈ S , H (x) = 0}. Let us define the function σ : Σ −→ R as
σ (x) = ⟨Hx (x), f(1) (x)⟩ · ⟨Hx (x), f(2) (x)⟩, where ⟨·, ·⟩ means the dot product in R2 . We will define the vector field g on the boundary Σ , see [2]. 1. At these points x ∈ Σ in which σ (x) > 0 we define the vector field g as an arithmetic mean of the vectors f(1) and f(2) g := gA (x) =
1 (1) (f (x) + f(2) (x)),
2
x ∈ Σ.
In this case we say that x belongs to the crossing set Σc ⊂ Σ , where Σc = {x ∈ Σ : σ (x) > 0}.
The subset Σc of Σ contains these points x ∈ Σ , in which both vectors f(1) and f(2) point into the same region.
134
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
Fig. 1. Stable sliding trajectory on the left. On the right: Time reversed ⇒ the trajectory is unstable, different from the original one.
2. In all other cases of configuration of the vectors f(1) and f(2) we apply the so-called Filippov’s method and define the vector field g on the boundary Σ as a convex combination of the vector fields f(1) and f(2) , g := gF (x) = λ(x)f(1) (x) + (1 − λ(x))f(2) (x),
λ : Σ −→ [0, 1],
λ(x) =
x ∈ Σ , where
(2)
Hx (x) · f Hx (x) ·
(f(2) (x)
(x) , − f(1) (x))
(3)
Hx (x) is the normal vector at the point x ∈ Σ . In this case we say that x belongs to the sliding set Σs . The set Σs is a complement of Σc in Σ , i.e. Σs = {x ∈ Σ : σ (x) ≤ 0}.
On the sliding boundary Σs special points (sliding points) can be detected. Let us classify some of them. 1. Singular sliding point is such point x ∈ Σs that
⟨Hx (x), f(1) (x)⟩ = 0 and also ⟨Hx (x), f(2) (x)⟩ = 0. At these points, both vectors f(1) (x) and f(2) (x) are tangent to Σs . 2. The point x ∈ Σs is a generic pseudoequilibrium if f(1) (x) ̸= 0,
g(x) = 0,
f(2) (x) ̸= 0.
At these points, the vectors f(1) (x) and f(2) (x) are anti-collinear. 3. In a boundary equilibrium x ∈ Σs , one of the vectors f(i) (x) vanishes, f(1) (x) = 0 or f(2) (x) = 0. 4. The point x ∈ Σs is a tangent point if both f(1) (x) ̸= 0, f(2) (x) ̸= 0 and
⟨Hx (x), f(1) (x)⟩ = 0 or ⟨Hx (x), f(2) (x)⟩ = 0. In this case, both vectors f(1) (x), f(2) (x) are nonzero, but one of them is tangent to Σ . The tangent point terminates Σs in Σ . Let us note that in general, in Filippov systems the theorem on existence and uniqueness does not apply (see Fig. 1). Two Filippov systems are called topologically equivalent if there is a homeomorphism h : R2 −→ R2 – that maps the discontinuity boundary of one system onto the discontinuity boundary of the other one – orbits of one system are mapped onto the corresponding orbits of the other one with time orientation preserved – crossing and sliding segments of any orbit are mapped onto the corresponding segments of its image. Remark 1. The another possibility how to define the vector field g on the boundary Σ is to apply the so-called Utkin’s equivalent control method, see [2]. In the Utkin’s method the vector field g on the boundary Σ is defined as g := gU (x) =
1 (2) 1 (1) f (x) + f(2) (x) + f (x) − f(1) (x) β(x), 2
2
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
135
Table 1 Defining systems on the boundary of the discontinuity H (x, α) = 0. Discontinuity boundary Crossing boundary Sliding boundary Tangent points of the vector field fi Pseudoequilibrium Boundary equilibrium
H (x, α) = 0 ⟨Hx , f(1) ⟩⟨Hx , f(2) ⟩ > 0 ⟨Hx , f(1) ⟩⟨Hx , f(2) ⟩ ≤ 0 ⟨Hx (x, α), f(i) (x, α)⟩ = 0 λ1 f(1) (x, α) + λ2 f(2) (x, α) = 0, λ1 + λ2 − 1 = 0 f(i) (x, α) = 0
where the control function β : Σ −→ [−1, 1] has the form
β(x) = −
⟨Hx (x), f1 (x) + f2 (x)⟩ . ⟨Hx (x), f2 (x) − f1 (x)⟩
(4)
If we put β(x) := 1 − 2λ(x) in (4), we obtain the Eq. (3), i.e. the Utkin’s control method and the Filippov’s method are algebraically equivalent.
2.2. Parameter-dependent Filippov systems In piecewise-smooth systems the interaction of invariant sets with the ‘‘switching’’ boundary often produces bifurcations that are forbidden in smooth systems. These are known as discontinuity-induced bifurcations. We will study those bifurcations that are likely to occur as a single parameter is varied. Let us have a piecewise-smooth dynamical system dependent on a parameter α ,
x˙ =
f(1) (x, α), f(2) (x, α),
x ∈ S1 (α), x ∈ S2 (α),
(5)
x ∈ R2 , α ∈ R, and f(i) : R2 × R −→ R2 , i = 1, 2, are smooth functions of (x, α), S1 (α) = x ∈ R2 : H (x, α) < 0 ,
S2 (α) = x ∈ R2 : H (x, α) > 0 ,
Σ (α) = {x ∈ S : H (x, α) = 0},
H (x, α) is a smooth function such that Hx (x, α) ̸= 0 for all (x, α) ∈ Σ (α), where Σ (α) is the boundary of discontinuity (‘‘switching’’ boundary). We say that the system (5) exhibits a bifurcation at the parameter value α = α0 if any arbitrarily small parameter perturbation results in a topologically nonequivalent system. Bifurcations that occur within the single region S1 or S2 for all parameter values of interest can be studied using smooth bifurcation theory. For trajectories crossing the boundary one can adopt the same theory, [2,6]. We will be interested in the discontinuity-induced bifurcations, [2,1]. They must occur if we continuously vary parameters. A piecewise-smooth system is piecewise-structurally stable if there exists an ε > 0 such that all C 1 perturbations of maximum size ε of the vector field f lead to topological equivalent phase portraits. We say that a discontinuity-induced bifurcation occurs at a parameter value at which a piecewise-smooth system is not piecewise-structurally stable, i.e., there exists an arbitrarily small perturbation that leads to a system that is not piecewise-topological equivalent. The simplest kind of the discontinuity-induced bifurcations occurs when an equilibrium point lies precisely on a discontinuity boundary. In Filippov systems there might occur so called pseudoequilibria, which are equilibria of the sliding flow but are not equilibria of any of the vector fields of the original system. For examples of the defining equations, see Table 1. Let us mention only three types of local sliding bifurcations that occur on the sliding set: – Collision of an equilibrium, for example a node, with the boundary Σ . – Collision of tangent points (double tangency). – Collision of pseudoequilibria, for example pseudo-saddle–node bifurcation.
136
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
Fig. 2. Trajectories of the system in dependence on the parameter α .
Example 1. Let us simulate the behavior of the trajectories of the following parameter-dependent piecewise-smooth dynamical system. Under a small perturbation a node collides with the boundary Σ . f(1) (x, α) =
−3x1 − x2 , −x1 − 3x2
f(2) (x, α) =
0 , 1
H (x, α) = x2 + α.
In Fig. 2, for α < 0 there is a node Xα in the region S1 . The sliding region Σs goes from the left and ends at the tangent point Tα ∈ Σ . The red dashed part of the Σ represents the crossing boundary Σc . The trajectories that start on the left of Tα in S1 are heading to Σs , slide on it up to the tangent point and then they are attracted by the node. Trajectories that start in S2 tend towards the boundary Σ . To the left of the tangent point they slide to Tα and end in the node. If they go to the right of the tangent point, they cross the crossing boundary and again end in the node. If we increase α the node moves to the boundary and for α = 0 collides with it, X0 ≡ Tα , X0 becomes a pseudonode Pα . Then if α further increases the pseudonode Pα stays on the sliding boundary Σs and the tangent point Tα secedes from it and moves to the right again marking the end of the sliding set. 3. Ecological applications 3.1. Bazykin’s ecological model Bazykin’s ecological model is the dynamical system defined by the following differential equations: dx1
ax1 x2
= x1 − x21 −
dt dx2
= −dx2 +
dt
b + x1 ax1 x2
b + x1
,
,
x1 ∈ ⟨0, ∞), x2 ∈ ⟨0, ∞), a > 0, b > 0, d > 0.
We reformulate this system as a Filippov system. We divide the state space S = ⟨0, ∞) × ⟨0, ∞) ⊂ R2 into two open regions S1 , S2 separated by a boundary Σ defined by a smooth parameter dependent function H. ax1 x2
dx dt
x1 − x21 −
= f(1) (x) =
ax1 x2 b + x1
dx dt
b + x1 ,
− dx2
x1 − x21 −
x ∈ S1 ,
ax1 x2
b + x1 , = f(2) (x) = ax1 x2 − dx2 − Ex2 b + x1
H (x, α) = x2 − α,
(6)
x ∈ S2 ,
(7)
x ∈ Σ,
(8)
where x1 > 0 is the mass of the prey, x2 > 0 is the mass of the predator, α, a, b, d, E are parameters, and S1 (α) = x ∈ R2 : H (x, α) < 0 ,
S2 (α) = x ∈ R2 : H (x, α) > 0 ,
Σ (α) = {x ∈ S : H (x, α) = 0}.
Let us note that f(1) is a continuous vector field on the open set S1 , f(2) on S2 . On the boundary Σ , we define
x˙ =
gA (x), gF (x),
x ∈ Σc = {x ∈ Σ : σ (x) > 0}, x ∈ Σs = {x ∈ Σ : σ (x) ≤ 0},
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
137
where σ (x) = ⟨Hx (x), f(1) (x)⟩ · ⟨Hx (x), f(2) (x)⟩, ⟨·, ·⟩ is the dot product in R2 and Hx (x) is the gradient of H at the point x ∈ Σ. We performed a one parameter study for a continuation parameter α ∈ R; the values of the other parameters were fixed. We apply the theory of Filippov systems and compute generic pseudoequilibria, i.e. such points on Σs in which the vectors f(1) (x) and f(2) (x) are anticolinear. For λ ∈ (0; 1), f(1) (x, α) ̸= 0, f(2) (x, α) ̸= 0, we have to solve the following system of algebraic equations: H (x, α) = 0,
λf(1) (x, α) + (1 − λ)f(2) (x, α) = 0.
Tangent points, i.e. ending points of the sliding set, are such points in which both vectors f(i) (x) are nonzero, but one of them is tangent to Σ . The defining systems for the tangent points are H (x, α) = 0
⟨Hx (x, α), f(1) (x, α)⟩ = 0 for f(1) tangent to Σ , ⟨Hx (x, α), f(2) (x, α)⟩ = 0 for f(2) tangent to Σ .
and
For more details, see [2,4,16]. 3.2. Population model: spiders on the vineyard Population models with a human intervention can be formulated as Filippov systems and then analyzed according to the theory described in [2]. The population of the spiders that predate on the insects in the vineyards is an example of the situation, when the standard predator–prey model is modified by human intervention if one of the populations exceeds the given limit. The ecosystem we consider consists of the spiders populations s(t ) and the insects on which they prey. Let f (t ) be the insects population living in the woods and the green patches bordering the cultivated land, and let v(t ) denote the insects having the vineyards as habitat. In [5], the following model was proposed f˙ = rf
f
1−
s˙ = s −a +
− csf ,
W kbv
H +v
v˙ = v e −
bs
H +v
+ kcf
,
(9)
,
where k is the conversion factor of prey into new spiders, k < 1, r is the specific growth rate of the prey in the woods, e is the specific growth rate of the prey in the vineyards, W and H are carrying capacities of the woods and the vineyard, respectively. All parameters are positive. Let us note that the first and the third equations describe the evolution of the wood insects and the growth of the insects population in open fields, respectively. The dynamics of spiders is given in the second equation. Now, let us model the human intervention, for example some spraying effect on the ecosystem of the vineyard. We define the function
ϕ(f , s, v) = vm − v,
(10)
where vm is a prescribed limit of the population of the insects on the vineyards. If ϕ(f , s, v) > 0 then there is no spraying needed and the system (9) describes the situation. If the population of the insects exceeds the limit, i.e. ϕ(f , s, v) = vm − v < 0, the spraying starts. Then we have to modify the model accordingly:
f
− csf − h(1 − q)f , s˙ = s −a + + kcf − hKqs, H +v bs v˙ = v e − − hqv, H +v
f˙ = rf
1−
W kbv
(11)
where q is the portion of insecticide sprayed directly on the vineyards, (1 − q) denotes the portion which may accidentally be dispersed in the woods, h stands for effectiveness against the parasites and K models the possibly smaller effect of the spraying on the spiders, 0 < K < 1. The Eqs. (9), (11) together with (10) define the Filippov system F . The boundary Σ of the system F is given by
Σ = {(f , s, v) : ϕ(f , s, v) = 0}.
(12)
The results of the dynamical simulations are summarized in Fig. 3: The trajectory reaches the boundary Σ , slides on it and then again leaves it. The situation can be interpreted as follows. The amount of the insects in the vineyards grows. When the boundary is reached, i.e. vm = v , the insecticide is applied. It takes some time to see the effect of spraying (this is the
138
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
Fig. 3. Spiders on the vineyard. Parameters: s(0) = 0, f (0) = 2, v(0) = 1.5, vm = 3.
time of sliding). But then the amount of insects again decreases. The whole process can be repeated if we set a lower bound to the amount of insects in the vineyards. 4. Closed gas–liquid system A closed gas–liquid system, more precisely an ideal gas–liquid system, see [7], is a simplified model of the slurry reactor. In the slurry reactor, the feed contains a liquid phase (a catalyst in diesel oil) and a gas phase (hydrogen and ethylene), see [17]. The output from the reactor consists of solid HDPE dispersed in diesel oil and unreacted gases. The process runs continuously, because the feed is sent in and also the product is withdrawn continuously. Two different mathematical models are needed to describe this system: one for the liquid output, the other one for the gas output. Only one of these two models takes place at a given time, because the dispersed polymer and gases cannot come out from the reactor simultaneously. Under suitable assumptions, the closed gas–liquid system may serve as a good approximation of the industrial HDPE reactor. Like the slurry reactor, a closed gas–liquid system consists of a closed tank with only one outlet tube. Let us start with a simple system shown in Fig. 4. The feed, that consists of a liquid and a gas, is sent to the closed tank, which has the only one outlet tube. If the level of the liquid is above the outlet tube opening, see Fig. 4(a), the liquid comes out from the tank. Otherwise, the gas comes out, see Fig. 4(b). Let FL and FG be molar inflow rates, G and L be molar outflow rates of the liquid and the gas, respectively. The following assumptions help to simplify the model. 1. 2. 3. 4. 5. 6. 7. 8.
The feed consists of an ideal gas G and an incompressible liquid L. The gas and the liquid do not react. The liquid has negligible vapor pressure at the operating conditions. The gas does not dissolve in the liquid at the operating conditions. In the liquid, there are no gas bubbles that might be leaving the tank. The valve dynamics is ignored. The flow rate through the valve is proportional to the difference of the tank pressure and the outlet pressure. The temperature, the molar inflow rates, the outlet pressure and the valve opening are all constant.
Classically, we need two sets of differential equations to describe the behavior of this system: one for the liquid model, the other one for the gas model. The following model equations of the ideal gas–liquid system were proposed by Moudgalya and Ryali, [7]. Liquid model: ML /ρL > Vd dMG = FG , dt dML = FL − L, dt RT ML + = V, MG P ρL L = kL x(P − Pout ).
Gas model: ML /ρL < Vd dMG = FG − G, dt dML = FL , dt RT ML MG + = V, P ρL G = kG x(P − Pout ).
(13)
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
a
b
139
c
Fig. 4. An ideal gas–liquid system.
Here, ML and MG are the molar hold-ups of the liquid and the gas, ρL is the molar density, T is the absolute temperature, P and Pout are pressures in the tank and in the outlet, respectively. R is the gas constant, V is the volume of the tank and Vd denotes the volume of the liquid below the outlet tube opening in the tank, see Fig. 4(c). The valve opening is denoted by x, 0 < x ≤ 1. The valve coefficients are kL for the liquid and kG for the gas flow. With respect to time, FG and FL have to be positive constant. The valve opening x is a constant, too. Let us point out that the system keeps switching between the liquid and gas models and that the discontinuity is a natural part of the system. Applying the theory of Filippov systems [2,18], both the outgoing liquid phase and the outgoing gas phase can be defined using only one set of ordinary differential equations with a discontinuous right-hand side, [12]. From the third equation in (13) we derive the expression for the pressure P in the tank, P =
MG RT V − ML /ρL
.
(14)
Substituting this P into the fourth equations in (13) we obtain
L = kL x
MG RT V − ML /ρL
− Pout ,
(15)
− Pout .
(16)
and
G = kG x
MG RT V − ML /ρL
The first and second equations in (13) together with (15) and (16) define the equivalent Filippov system F . After a few simplifying assumptions, the behavior of this system depends on ten parameters. We have chosen two of them namely, FG and FL , a molar inflow of both gas and liquid. With respect to time, FG and FL are positive constants. We study the dependence of the solution (ML MG ) (the molar hold-ups of the liquid and gas) on these two parameters. Now we are ready to formulate the system (13) as the following Filippov system F .
F :
d
MG ML
dt
f(1) =
FG − kG x FL
f(1) (MG , ML , FG , FL ), f(2) (MG , ML , FG , FL ),
=
− Pout , V − ML /ρL MG RT
ϕ(MG , ML , FG , FL ) < 0, ϕ(MG , ML , FG , FL ) > 0, FG , MG RT f(2) = FL − kL x − Pout V − ML /ρL
(17)
where
ϕ(MG , ML ) := ϕ(MG , ML , FG , FL ) = ML − ρL Vd , Σ = {(MG , ML ) ∈ R2 : ϕ(MG , ML ) = 0}, 2 S1 = {(MG , ML ) ∈ R : ϕ(MG , ML ) < 0}, S2 = {(MG , ML ) ∈ R2 : ϕ(MG , ML ) > 0}. Here, V is the volume of the tank and Vd , 0 < Vd < V , denotes the volume below the outlet tube opening in the tank. Two nullclines of the system (17) together with the discontinuity boundary Σ split the state space S into four regions with the different qualitative behavior of the system, [7].
140
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
Fig. 5. Dependence of the solution (MG , ML ) on the molar inflow of the gas for FG = 0.1, FL = 2.5. On the right: Zoom of the sliding set.
The following theorem says that the sliding set Σs is an attractor for the whole state space S. Theorem 1. Let the Filippov system be described by Eq. (17). Then for any x ∈ S \ Σs , x = (MG , ML ), every orbit γx ends on the sliding set Σs . For proof, see [7]. We simulated the behavior of the system (17) for different values of the parameters FG , FL , see Fig. 5, and constructed a solution diagram. It shows that the coordinates MG and ML of the pseudoequilibria and tangent points on the sliding and crossing sets depend linearly on the molar inflow of the gas FG and liquid FL . 4.1. Dimensionless closed gas–liquid system It turned out that a transformation of the model equations into a dimensionless form reduces the number of parameters substantially—to only four. It can be proven that the behavior of the dimensionless system is qualitatively the same as the behavior of the system (17), but the reduction of the number of parameters substantially simplifies the bifurcation analysis. The Filippov system (17) depends on 10 parameters, in particular the parameter set is p = (FG , FL , ρL , V , Vd , T , Pout , x, kL , kG ),
p ∈ R10 .
If we define dimensionless state variables, then choose a suitable scaling and simplify the resulting equations we obtain four ˜ K˜ , β˜ L , γ˜ ), q˜ ∈ R4 , see [9] for the detailed transformation. dimensionless parameters q˜ = (δ, The resulting dimensionless Filippov system F˜ then depends only on q˜ . It reads (tildes are omitted): f(1) (MG , ML , q), f(2) (MG , ML , q),
ϕ(MG , ML , q) < 0, ϕ(MG , ML , q) > 0, MG 2 + K β 1 − K δβ L L , f(1) (MG , ML , q) = 1 − ML F (q) :
d
dt
MG ML
=
1
(2)
f
(MG , ML , q) =
1 1 − δβL2
MG 1 − ML
+ βL
,
ϕ(MG , ML , q) = ML − γ . It can be shown that the behavior of the dimensionless system is qualitatively the same as the behavior of the original one, but the reduction of the number of parameters only to four simplifies substantially the bifurcation analysis, [9]. 4.2. Soft drink process The process of manufacturing soft-drink depicted in Fig. 6 is based on the reaction between CO2 and water: CO2 + H2 O → H2 CO3 .
(18)
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
141
To simplify the model, we will suppose that – – – – – –
The system contains only components CO2 , H2 O and H2 CO3 (denoted by indices 1, 2 and 3, respectively). Intermediate ionization reactions and dissociation of H2 CO3 are ignored. In the liquid, there are no gas bubbles. The valve dynamics is ignored. The flow rate through the valve is proportional to the difference of the tank pressure P and the outlet pressure Pout . The temperature T , the molar inflow rates F1 and F2 , the outlet pressure, valve coefficients kG and kL and the valve opening x are all constant.
Let M1 = M1 (t ),
M2 = M2 (t ),
M3 = M3 (t ),
be the total molar hold-ups of CO2 , H2 O and H2 CO3 , respectively. For a fixed t, let us define a scalar function ϕ = ϕ(M1 , M2 , M3 ),
ϕ(M1 , M2 , M3 ) =
M2
ρL
+
M3
ρa
− Vd ,
(19)
where ρL , ρa are molar densities of water and acid, respectively. The volume of the whole tank is equal to V and the part of the volume that is below the opening of the dip tube is denoted as Vd , 0 < Vd < V . Similarly as in [8,9], in the tank two systems take place: the liquid model (the liquid leaves the tank) if ϕ(M1 , M2 , M3 ) > 0 or the gas model (the gas leaves the tank) for ϕ(M1 , M2 , M3 ) < 0. The acid phase consists of H2 CO3 , H2 O and dissolved CO2 while the gas phase contains only CO2 . As a consequence, the acid model is described by 3 ODEs and 6 algebraic equations, the gas model needs also 3 ODEs but only 4 algebraic equations. Let us give the list of these equations. Differential equations: Liquid model : ϕ(M1 , M2 , M3 ) > 0
Gas model : ϕ(M1 , M2 , M3 ) < 0
dM1
dM1
dt dM2 dt dM3 dt
= F1 − L1 − rV ,
dt dM2
= F2 − L2 − rV ,
dt dM3
= −L3 + rV ,
dt
= F1 − G − rV , = F2 − rV , = rV ,
M1 , M2 and M3 are the total molar hold-ups of CO2 , H2 O, H2 CO3 , respectively. F1 and F2 are the molar inflow rates. The molar flow rates of the components through the valve are denoted L1 , L2 and L3 in the liquid model and G in the gas model. The rate r of the reaction (18) is given by r = κc
M1 M2 V2
,
where κc is the rate constant.
(20)
Algebraic equations: Liquid model : ϕ(M1 , M2 , M3 ) > 0
Gas model : ϕ(M1 , M2 , M3 ) < 0
0 = M1 − (Mℓ + Mg ),
0 = M1 − (Mℓ + Mg ),
0=P−
σ Mℓ Mℓ + M2 + M3
0=V− 0= 0=
M1 RT P
+
M2 Mℓ + M2 + M3 M3 Mℓ + M2 + M3
M2
ρL − −
, +
0=P− M3
ρa
,
L2 L1 + L2 + L3 L3 L1 + L2 + L3
σ Mℓ Mℓ + M2 + M3
0=V−
,
M1 RT P
+
M2
ρL
, +
M3
ρa
,
0 = G − kG x(P − Pout ),
,
0 = L1 + L2 + L3 − kL x(P − Pout ), P and T means pressure and temperature in the tank, the hold-ups of CO2 in liquid and gas are denoted Mℓ and Mg , R is a gas constant and σ is Henrys constant for CO2 . The straightforward computation shows that both the system of DAEs for the gas mode and the system of DAEs for the liquid model has index 1, [15].
142
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
a
b
Fig. 6. The soft-drink process.
Similarly as in the system of ODEs for the ideal gas–liquid model, the system of DAEs describing the soft-drink model can be reformulated by making use of Filippov’s method, [14]. As a result we obtain the formulation of the soft-drink process as Filippov system. The corresponding differential equations are: M1 t M2 t M3 t
= λ(G − L1 ) + F1 − G − rV , = −λL2 + F2 − rV , = −λL3 + rV ,
where λ =
F2 ρa + rV (ρL − ρa ) L2 ρa + L3 ρL
.
Similarly, the corresponding algebraic equations are: 0 = M1 − (Mℓ + Mg ),
σ Mℓ
0=P−
Mℓ + M2 + M3
0=V−
0= 0=
M1 RT P
+
M2 Mℓ + M2 + M3 M3 Mℓ + M2 + M3
M2
ρL − −
, +
M3
ρa
,
L2 L1 + L2 + L3 L3 L1 + L2 + L3
, ,
0 = L1 + L2 + L3 − kL x(P − Pout ), 0 = G − kG x(P − Pout ). Let us note that missing equations in both liquid and gas models that should describe the same phenomenon are replaced by equations 0 = 0. For a particular example, the results are depicted in Fig. 7. 5. Conclusions We considered three different parameter-dependent Filippov’s dynamical systems. The objective was to perform a numerical analysis by dynamical simulations. We used the software [13]. Concerning the ecological applications, we have shown that Filippov’s formulation is very natural (e.g. comparing [5] with our approach). As far as the study of the closed gas–liquid system is concerned, it is just the first step towards to a model of the real HDPE reactor. Acknowledgments The research was supported by MSM 6046137306, the project financed by The Ministry of Education, Youth and Sports, Czech Republic. The authors were also supported from specific university research grant (MSMT no. 21/2010).
M. Biák et al. / Journal of Computational and Applied Mathematics 254 (2013) 132–143
143
Fig. 7. Soft-drink process: (a)–(c) The integral curves of the state variables M1 , M2 and M3 . (d) The trajectory of the soft-drink system starting at the point (0.72, 95, 0).
References [1] Yu.A. Kuznetsov, Elements of Applied Bifurcation Theory, second ed., Springer, New York, 1998. [2] M. di Bernardo, C.J. Budd, A.R. Champneys, P. Kowalczyk, Piecewise-Smooth Dynamical Systems, Theory and Applications, Springer, London, 2008. [3] A.D. Bazykin, A.I. Khibnik, B. Krauskopf, Nonlinear Dynamics of Interacting Populations, in: World Scientific Series on Nonlinear Science, Series A, vol. 11, 1998. [4] W. Govaerts, Numerical Methods for Bifurcations of Dynamical Equilibria, SIAM, Philadelphia, 2000. [5] E. Venturino, M. Isaia, F. Bona, S. Chatterjee, G. Badino, Biological controls of intensive agroecosystems: wanderer spiders in the Langa Astigiana, Ecological Complexity 5 (2008) 157–164. [6] D.J.W. Simpson, Bifurcations in Piecewise-Smooth Continuous Systems, World Scientific, University of British Columbia, Canada, 2010. [7] K.M. Moudgalya, V. Ryali, Chemical Engineering Science 56 (2001) 3595–3609. [8] M. Biák, D. Janovská, Filippov dynamical systems, in: Seminar on Numerical Analysis & Winter School, Proceedings of the Conference, SNA’09, Institute of Geonics AS CR Ostrava, 2009, pp. 1–4. [9] M. Biák, D. Janovská, A parametric study of the dimensionless closed gas–liquid system, in: Nové Hrady, Hana Bílková (Eds.), Seminar on Numerical Analysis & Winter School, Proceedings of the Conference, SNA’10, Institute of Computer Science AS CR Praha, 2010, pp. 26–28. [10] L.L. Böhm, P. Goebel, P.R. Schöneborn, Detailed reaction engineering as a basis of modern slurry technology for PE-HD-production, Die Angewandte Makromolekulare Chemie 174 (1990) 189–203. [11] F. Dercole, Yu.A. Kuznetsov, ACM Transactions on Mathematical Software 31 (2005) 95–119. [12] A.F. Filippov, Differential Equations with Discontinuous Righthand Sides, Kluwer Academic Publishers, Dordrecht, 1988. [13] P.T. Piiroinen, Yu.A. Kuznetsov, ACM Transactions on Mathematical Software 34 (2008) 1–24. [14] J. Agrawal, K.M. Moudgalya, A.K. Pani, Sliding motion of discontinuous dynamical systems described by semi-implicit index one differential algebraic equations, Chemical Engineering Science 61 (2006) 4722–4731. [15] P. Kunkel, V. Mehrmann, Differential–Algebraic Equations: Analysis and Numerical Solution, European Mathematical Society, 2006. [16] T. Hanus, D. Janovská, Discontinuous, piecewise-smooth dynamical systems, in: CHISA 2006, Praha, Czech Republic, CD-ROM of Full Texts CHISA 2006, Title Index: D, Org.No.: F4.8. [17] K.M. Moudgalya, S. Jaguste, A class of discontinuous dynamical systems II., an industrial slurry high density polyethylene reactor, Chemical Engineering Science 56 (2001) 3611–3621. [18] M. di Bernardo, P. Kowalczyk, A. Nordmark, Bifurcations of dynamical systems with sliding: derivation of normal-form mappings, Physica D (2002) 175–205.