Chaotic Attractors of Control Systems

Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Chaotic Attractors of Control Systems David N. Cheban1 and Cristiana Mammana2 1 2 State University of Moldova, str. A. Mattevich 60 , MD-2009, Chişinău, Republic of Moldova (E-mail: [email protected]) Institute of Economics and Finances, University of Macerata, str. Crescimbeni 14, I–62100 Macerata, Italy (E-mail: [email protected]) Abstract. The paper is dedicated to the study of the problem of existence of chaotic attractors of discrete control systems and to description of its structure. We consider a family of continuous mappings of metric space W into itself, and (W, fi )i∈I is the family of discrete dynamical systems. On the metric space W we consider a discrete inclusion ut+1 ∈ F (ut ) associated by M := {fi : i ∈ I}, where F (u) = {f (u) : f ∈ M} for all u ∈ W. If the family M consists of a finite number of maps, then the corresponding attractor is chaotic. We study this problem in the framework of non-autonomous dynamical systems (cocyles). Keywords: Chaotic attractor; set-valued dynamical system; control system, collage, cocycle. 1 Introduction The aim of this paper is the study of the problem of existence of compact global chaotic attractors of discrete control systems (see, for example, Cheban [4,5] Bobylev, Emel’yanov and Korovin [2] and the references therein). Let W be a metric space, M := {fi : i ∈ I} be a family of continuous mappings of W into itself and (W, fi )i∈I be the family of discrete dynamical systems, where (W, f ) is a discrete dynamical system generated by positive powers of continuous map f : W → W . On the space W we consider a discrete inclusion ut+1 ∈ F (ut ) associated by M := {fi : i ∈ I} (denotation DI(M)), where F (u) = {f (u) : f ∈ M} for all u ∈ W. A solution of the discrete inclusion DI(M) is called (see, for example, [2,7]) a sequence {{xj } | j ≥ 0} ⊂ W such that xj = fij xj−1 (1) for some fij ∈ M (trajectory of DI(M)). We can consider that it is a discrete control problem, where at each moment of the time j we can apply a control from the set M, and DI(M) is the set of possible trajectories of the system. 131 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece The problem of existence of compact global attractors for a discrete inclusion arise in a number of different areas of mathematics (see, for example, [5,6] and the references therein). A sequence ω : Z+ 7→ {1, 2, . . . , m} is called p ∈ N := {1, 2, . . .} periodic, if ω(n + p) = ω(n) for all n ∈ Z+ . A point x0 ∈ W is said to p-periodic for DI(M) if there exists an pperiodic sequence ω : Z+ 7→ {1, 2, . . . , m} such that solution {x(k)}k∈Z+ of equation (1) (ω(i) = ij for all i ∈ Z and ij ∈ {1, 2, . . . , m}) with initial data x(0) = x0 is mp-periodic, i.e., x(k + p) = x(k) for all k ∈ Z+ . It is well known the following result. Theorem 1. [1, Ch.II,IV] Let M = {f1 , f2 , . . . , fm } be a finite family of continuous mappings from W into itself. If there exists a number q ∈ (0, 1) such that ρ(fi (x1 ), fi (x2 )) ≤ qρ(x1 , x2 ) for all x1 , x2 ∈ W and i ∈ {1, 2, . . . , m}, then the following statement hold: 1. DI(M) admits a compact global attractor I, i.e., (a) I is a nonempty, compact and invariant set (F (I) = I, S where F (x) := {f1 (x), f2 (x), . . . , fm (x)) for all x ∈ W and F (I) := {F (x) : x ∈ I}); (b) lim β(F n (x), I) = 0 for all x ∈ W uniformly with respect to x n→∞ on every compact subset M from W , where β(A, B) := sup ρ(a, B) a∈B (A, B ⊆ W ). 2. I coincides with the closure of the all periodic points of DI(M). In the book [5] (Chapter VI) it was generalized this theorem for the finite family M = {f1 , f2 , . . . , fm } when it is contracting in the generalized sense, i.e., there are two positive numbers N and q ∈ (0, 1) such that ρ(fin fin−1 . . . fi1 (x1 ), fin fin−1 . . . fi1 (x2 )) ≤ N q n ρ(x1 , x2 ) (2) for all x1 , x2 ∈ W and n ∈ N , where i1 , i2 , . . . , in ∈ {1, 2, . . . , m}. 2 Some Notions and Facts from Dynamical and Control Systems In this paper we will use some notions and facts from the theory of dynamical and control systems [5]. In this Section we collect some of them. 2.1 Set-valued dynamical systems and their compact global attractors Let (X, ρ) be a complete metric space, S be a group of real (R) or integer (Z) numbers, T (S+ ⊆ T ) be a semi-group of additive group S. By C(X) we denote the family of all non-empty compact subsets of X. For every point x ∈ X and number t ∈ T we put S in correspondence a closed compact subset π(t, x) ∈ C(X). So, if π(P, A) = {π(t, x) : t ∈ P, x ∈ A}(P ⊆ T ), then 132 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 1. π(0, x) = x for all x ∈ X ; 2. π(t2 , π(t1 , x)) = π(t1 + t2 , x) for all x ∈ X; β(π(t, x), π(t0 , x0 )) = 0 for all x0 ∈ X and t0 ∈ T , where 3. lim x→x0 ,t→t0 β(A, B) = sup{ρ(a, B) : a ∈ A} is a semi-deviation of the set A ⊆ X from the set B ⊆ X. In this case it is said [11] that there is defined a set-valued semi-group dynamical system. Let T ⊂ T ′ ⊂ S. A continuous mapping γx : T ′ → X is called a motion of the set-valued dynamical system (X, T, π) issuing from the point x ∈ X at the initial moment t = 0 and defined on T ′ , if a. γx (0) = x; b. γx (t2 ) ∈ π(t2 − t1 , γx (t1 )) for all t1 , t2 ∈ T ′ (t2 > t1 ). The set of all motions of (X, T, π), passing through S the point x at the initial moment t = 0 is denoted by Fx (π) and F(π) := {Fx (π) | x ∈ X} (or simply F). The trajectory γ ∈ F(π) defined on S is called a full (entire) trajectory of the dynamical system (X, T, π). Denote by Φ(π) the set ofTall full trajectories of the dynamical system (X, T, π) and Φx (π) := Fx (π) Φ(π). A system (X, T, π) is called [3],[5] compactly dissipative, if there exists a nonempty compact K ⊆ X such that lim β(π t M, K) = 0; t→+∞ for all M ∈ C(X). Let (X, T, π) be compactly dissipative and K be a compact set attracting every compact subset of X. Let us set \ [ J := ω(K) := π τ K. (3) t≥0 τ ≥t It can be shown [3],[5] that the set J defined by equality (3) does not depend on the choice of the attractor K, but it is characterized only by the properties of the dynamical system (X, T, π) itself. The set J is called Levinson center of the compact dissipative dynamical system (X, T, π). 2.2 Discrete inclusions, ensemble of dynamical systems (collages) and cocycles Let W be a topological space. Denote by C(W, W ) the space of all continuous operators f : W → W equipped with the compact-open topology. Consider a set of operators M ⊆ C(W, W ) and, respectively, an ensemble (collage) of 133 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece discrete dynamical systems (W, f )f ∈M ((W, f ) is a discrete dynamical system generated by positive powers of map f ). A discrete inclusion DI(M) is called (see, for example, [2,7]) a set of all sequences {xj } ⊂ W (j ∈ Z+ ) such that xj = fij xj−1 (4) for some fij ∈ M (trajectory of DI(M)). A bilateral sequence {xj } ⊂ W (j ∈ Z) is called a full trajectory of DI(M) (entire trajectory or trajectory on Z), if xn+j = fij xn+j−1 for all n ∈ Z. Let us consider the set-valued function F : W → C(W ) defined by the equality F(x) := {f (x) |f ∈ M}. Then the discrete inclusion DI(M) is equivalent to the difference inclusion xj ∈ F(xj−1 ). (5) inclusion (5) (or Denote by Fx0 the set of all trajectories of discrete S DI(M)) issuing from the point x0 ∈ W and F := {Fx0 | x0 ∈ W }. Below we will give a new approach concerning the study of discrete inclusions DI(M) (or difference inclusion (5)). Denote by C(Z+ , W ) the space of all continuous mappings f : Z+ → W equipped with the compact-open topology. Denote by (C(Z+ , X), Z+ , σ) a dynamical system of translations (shifts dynamical system or dynamical system of Bebutov [9,10]) on C(Z+ , W ), i.e. σ(k, f ) := fk and fk is a k ∈ Z+ shift of f (i.e., fk (n) := f (n + k) for all n ∈ Z+ ). We may now rewrite equation (4) in the following way: xj+1 = ω(j)xj , (ω ∈ Ω := C(Z+ , M)) (6) where ω ∈ Ω is the operator-function defined by the equality ω(j) := fij+1 for all j ∈ Z+ . We denote by ϕ(n, x0 , ω) the solution of equation (6) issuing from the point x0 ∈ E at the initial moment n = 0. Note that Fx0 = {ϕ(·, x0 , ω) | ω ∈ Ω} and F = {ϕ(·, x0 , ω) | x0 ∈ W, ω ∈ Ω}, i.e., DI(M) (or inclusion (5)) is equivalent to the family of non-autonomous equations (6) (ω ∈ Ω). From the general properties of difference equations it follows that the mapping ϕ : Z+ × W × Ω → W satisfies the following conditions: 1. ϕ(0, x0 , ω) = x0 for all (x0 , ω) ∈ W × Ω; 2. ϕ(n + τ, x0 , ω) = ϕ(n, ϕ(τ, x0 , ω), σ(τ, ω)) for all n, τ ∈ Z+ and (x0 , ω) ∈ W × Ω; 3. the mapping ϕ is continuous; 4. for any n, τ ∈ Z+ and ω1 , ω2 ∈ Ω there exists ω3 ∈ Ω such that U (n, ω2 )U (τ, ω1 ) = U (n + τ, ω3 ), (7) Qn where ω ∈ Ω, U (n, ω) := ϕ(n, ·, ω) = k=0 ω(k), ω(k) := fik (k = 0, 1, . . . , n) and fi0 := IdW . 134 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Let T1 ⊆ T2 be two sub-semigroups of group S, X, Y be two metric (or topological) spaces and (X, T1 , π) (respectively (Y, T2 , σ)) be a semigroup dynamical system on X (respectively on Y ). A triplet h(X, T1 , π), (Y, T2 , σ), hi is called a non-autonomous dynamical system, where h : X → Y is a homomorphism from (X, T1 , π) onto (Y, T2 , σ), i.e., h(π(t, x)) = σ(t, h(x)) for all x ∈ X and t ∈ T1 . Let W, Ω be two topological spaces and (Ω, T2 , σ) be a semi-group dynamical system on Ω. Recall [9] that a triplet hW, ϕ, (Ω, T2 , σ)i (or briefly ϕ) is called a cocycle over (Ω, T2 , σ) with the fiber W , if ϕ is a mapping from T1 × W × Ω to W satisfying the following conditions: 1. ϕ(0, x, ω) = x for all (x, ω) ∈ W × Ω; 2. ϕ(t+τ, x, ω) = ϕ(t, ϕ(τ, x, ω), σ(τ, ω)) for all t, τ ∈ T1 and (x, ω) ∈ W ×Ω; 3. the mapping ϕ is continuous. Let X := W ×Ω, and define the mapping π : X ×T1 → X by the equality: π((u, ω), t) := (ϕ(t, u, ω), σ(t, ω)) (i.e., π = (ϕ, σ)). Then it is easy to check that (X, T1 , π) is a dynamical system on X, which is called a skew-product dynamical system [9]; but h = pr2 : X → Ω is a homomorphism of (X, T1 , π) onto (Ω, T2 , σ) and hence h(X, T1 , π), (Ω, T2 , σ), hi is a non-autonomous dynamical system. From the presented above it follows that every DI(M) (respectively, inclusion (5)) in a natural way generates a cocycle hW, ϕ, (Ω, Z+ , σ)i, where Ω = C(Z+ , M), (Ω, Z+ , σ) is a dynamical system of shifts on Ω and ϕ(n, x, ω) is the solution of equation (6) issuing from the point x ∈ W at the initial moment n = 0. Thus, we can study inclusion (5) (respectively, DI(M)) in the framework of the theory of cocycles with discrete time. Below we need the following result. Theorem 2. [6] Let M be a compact subset of C(W, W ) and hW, φ, (Ω, Z+ , σ)i be a cocycle generated by DI(M). Then 1. Ω = P er(σ), where P er(σ) is the set of all periodic points of (Ω, Z+ , σ) (i.e. ω ∈ P er(σ), if there exists τ ∈ N such that σ(τ, ω) = ω); 2. the set Ω is compact; 3. Ω is invariant, i.e., σ t Ω = Ω for all t ∈ Z+ ; 4. ϕ satisfies the condition (7). 3 Chaotic attractors of discrete control systems In Section 3 we give the conditions of existence of chaotic attractor for discrete control systems. Denote by A the set of all mapping ψ : Z+ × R+ 7→ R+ possessing the following properties: 135 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece (G1) ψ is continuous; (G2) there exists a positive number t0 such that: (a) ψ(t0 , r) < r for all r > 0; (b) the mapping ψ(t0 , ·) : R+ 7→ R+ is monotone increasing. (G3) ψ(t + τ, r) ≤ ψ(t, ψ(τ, r)) for all t, τ ∈ Z+ and r ∈ R+ . Remark 1. 1. Note that the functions ψ(t, r) = N q t r (N > 0 and q ∈ (0, 1)) r and ψ(t, r) = 1+rt belong to A, where (t, r) ∈ Z+ × R+ . 2. Let f : R+ be a continuous function satisfying the conditions: 1. f (r) < r for all r > 0; 2. f is monotone increasing. Then the mapping ψ : Z+ × R+ 7→ R+ defined by equality ψ(t, r) = x(t) (8) for all (t, r) ∈ Z+ × R+ , where x(t) is a unique solution of difference equation xt+1 = f (xt ) with initial data x0 = r, belongs to A. Let ψ ∈ A. contracting, if A set M of operators from C(W, W ) is said to be ψ- ρ(fit fit−1 . . . fi1 (x1 ), fit fit−1 . . . fi1 (x2 )) ≤ ψ(t, ρ(x1 , x2 )) (9) for all x1 , x2 ∈ W and t ∈ N , where fi1 , fi2 , . . . , fit ∈ C(W ) and i1 , i2 , . . . , it ∈ N. The set S ⊂ W is 1. nowhere dense, provided the interior of the closure of S is empty set, int(cl(S)) = ∅; 2. totally disconnected, provided the connected components are single points; 3. perfect, provided it is closed and every point p ∈ S is the limit of points qn ∈ S with qn 6= p. The set S ⊂ W is called a Cantor set, provided it is totally disconnected, perfect and compact. The subset M of (X, T, π) is called (see, for example, [8]) chaotic, if the following conditions hold: 1. the set M is transitive, i.e. there exists a point x0 ∈ X such that M = H(x0 ) := {π(t, x0 ) : t ∈ T }; 2. M = P er(π), where P er(π) is the set of all periodic points of (X, T, π). Recall that a point x ∈ X of the dynamical system S T, π) is called T (X, Poisson’s stable, if x belongs to its ω-limit set ωx := t≥0 τ ≥t π(τ, x). Theorem 3. Suppose that the following conditions are fulfilled: 136 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece a. M is a finite subset of C(W, W ), i.e., M := {f1 , f2 , . . . , fm } and m ≥ 2; b. M is ψ-contracting for some ψ ∈ A. Then the following statement hold: 1. the cocycle hW, ϕ, (Ω, Z+ , σ)i (Ω := C(Z+ , M)) generated by DI(M) is compactly dissipative; 2. the skew-product dynamical system (X, Z+ , π) generated by DI(M) is compactly dissipative; 3. I = P er(ϕ), where P er(ϕ) := {u ∈ W : ∃τ ∈ N and ω ∈ Ω such that σ(τ, ω) = ω and ϕ(τ, u, ω) = u}; 4. if every map f ∈ M is invertible, then 1. Levinson’s center J of the skew-product dynamical system (X, Z+ , π) is a chaotic Cantor set; 2. there exists a residual subset J0 ⊆ J (large in the sense of Baire category), consisting from Poisson’s stable points, such that the positive semi-trajectory of every point x0 ∈ J0 is dense on J; 3. I = pr1 (J) (pr1 : X → Ω, where I is the Levinson’s center of cocycle ϕ and X := W × Ω), i.e., I is a continuous image of the Cantor set J. Proof. Let Y = Ω := C(Z+ , Q) and (Y, Z+ , σ) be a semi-group dynamical system of shifts on Y . Then Y is compact. By Theorem 2, P er(σ) = Ω and Ω is compact and invariant. Let hW, ϕ, (Ω, Z+ , σ)i be a cocycle generated by DI(M) (i.e., ϕ(n, u, ω) := Qn U (n, ω)u, where U (n, ω) = k=0 ω(k) (ω ∈ Ω)), (X, Z+ , π) be a skewproduct system associated by the cocycle ϕ (i.e., X := W ×Ω and π := (ϕ, σ)) and h(X, Z+ , π), (Y, Z+ , σ), hi (h := pr2 : X → Y ) be a non-autonomous dynamical system generated by cocycle ϕ. Under the conditions of Theorem 2 we have ρ(ϕ(n, u1 , ω), ϕ(n, u2 , ω)) ≤ ψ(n, ρ(u1 , u2 )) (10) for all n ∈ Z+ , u1 , u2 ∈ W and ω ∈ Ω. Now to finish the proof of the theorem it is sufficient to apply Theorem 6.1.3 [5, Ch.VI] and Theorem 2 to the non-autonomous dynamical system h(X, Z+ , π), (Ω, Z+ , σ), hi and denote by Iω := pr1 (JT ω ), where J is Levinson center of the dynamical system (X, Z+ , π), Jω := J h−1 (ω) and h := pr2 . Let M ⊂ C(W ), hW, ϕ, (Ω, Z+ , π)i (respectively (X, Z+ , π) ) be a cocycle (a skew-product dynamical system) generated by DI(M) and let I (J) be Levinson center of the cocycle ϕ (respectively, skew-product dynamical system (X, Z+ , π)). The set I is said to be a chaotic attractor of DI(M), if 1. the set J is chaotic, i.e. J is transitive and J = P er(σ), where J is the Levinson center of the skew-product dynamical system (X, Z+ , π) generated by DI(M); 137 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 2. I = pr1 (J). Remark 2. 1. Theorem 3 it was proved in [6] for the special case, when ψ(t, r) = N q t ((t, r) ∈ Z+ × R+ , N > 0 and q ∈ (0, 1)). 2. The problem of the existence of compact global attractors for DI(M) with finite M (collage or iterated function system (IFS)) was studied before in works [1,2] (see also the bibliography therein). In [1,2] the statement close to Theorem 3 was proved. Namely: 1. in [1] it was announced the first and proved the second statement of Theorem 3, if ψ(t, r) = q t r (t ∈ Z+ and q ∈ (0, 1)); 2. in [2] they considered the case when W is a compact metric space and every map f ∈ M = {f1 , f2 , . . . , fm } (i = 1, . . . , m) is contracting (not obligatory invertible). For this type of DI(M) it was proved the existence of a compact global attractor A such that for all u ∈ A and almost all ω ∈ Ω (with respect toQcertain measure on Ω) the trajectory ϕ(n, u, ω) = n U (n, ω)u (U (n, ω) := k=0 fik , (ik ∈ {1, . . . , m}) and fi0 := IdW ) was dense in A. The problem of description of the structure of the attractor I of DI(M) in general case (when the maps f ∈ M are not invertible) is more complicated. We plan to study this problem in one of our next publication. References 1.M. F. Barnsley. Fractals everywhere. N. Y.: Academic Press, 1988. 2.N. A. Bobylev, S. V. Emel’yanov and S. K. Korovin. Attractors of Discrete Controlled Systems in Metric Spaces. Computational Mathematics and Modeling, v.11, No.4, pp.321-326, 2000 [Translated from Prikladnaya Mathematika i Informatika, No.3, 1999, pp.5-10.] 3.D. N. Cheban. Global Attractors of Nonautonomous Dissipstive Dynamical Systems. Interdisciplinary Mathematical Sciences, vol.1, River Edge, NJ: World Scientific, 2004, xxiii+502 pp. 4.D. N. Cheban. Compact Global Attractors of Control Systems. Journal of Dynamical and Control Systems, vol.16, no.1, pp. 23–44, 2010. 5.D. N. Cheban. Global Attractors of Set-Valued Dynamical and Control Systems. Nova Science Publishers Inc, New York, 2010, xvii+269 pp. 6.D. N. Cheban D. N. and C. Mammana. Global Compact Attractors of Discrete Inclusions. Nonlinear Analyses, v.65 , No.8, pp.1669-1687, 2006. 7.L. Gurvits. Stability of Discrete Linear Inclusion. Linear Algebra Appl., 231, pp.47-85, 1995. 8.C. Robinson C. Dynamical Systems: Stabilty, Symbolic Dynamics and Chaos (Studies in Advanced Mathematics). Boca Raton Florida: CRC Press, 1995. 9.G. R. Sell. Topological Dynamics and Ordinary Differential Equations. Van Nostrand-Reinhold, London, 1971. 10.B. A. Shcherbakov. Topological Dynamics and Poisson’s Stability of Solutions of Differential Equations. Kishinev, Shtiintsa, 1972. (in Russian) 11.K. S. Sibirskii and A. S. Shube. Semidynamical Systems. Kishinev, Shtiintsa,1987. (in Russian) 138 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Implementation of CNN–Based Viterbi Algorithm Using FPGA Thang Van Chu2 , Viet–Thanh Pham1 , Mattia Frasca1 , Thang Manh Hoang2 , and Luigi Fortuna1 1 2 Dipartimento di Ingegneria Elettrica Elettronica e Informatica, Università degli Studi di Catania, 95125 Catania, Italy (E-mail: [email protected]) School of Electronics and Telecommunications, Hanoi University of Science and Technology, Hanoi, Vietnam (E-mail: [email protected]) Abstract. Viterbi decoder is commonly utilized in digital communication because of its effectiveness in error correction. In this paper, we implement Viterbi algorithm based on Cellular Neural Network (CNN) in which cells play roles as nodes in the trellis diagram. The simulation and realization of CNN–based Viterbi decoder show the feasibility and the accuracy of our approach. Moreover, our results are compared with ones of current FPGA–based Viterbi decoders. Keywords: Viterbi decoder, Cellular Neural Network, Field Programmable Gate Array. 1 Introduction Digital communication systems have been widely used because of their better features comparing to analog ones Viterbi and Omura[1]. Generally, quality of digital system is estimated by the bit error rate (BER). Due to the channel impairment as noise and fading, channel coding is required to improve performance of digital signal. There are two common types of channel coding: block coding and convolutional coding. Different from block coding, convolutional coding can be utilized for both block of data and data stream. Hence, convolutional coding has been applied to amount of practical applications, especially in wireless communication. An effective decoding algorithm for convolutional coding was proposed by Viterbi Viterbi[2]. Viterbi algorithm reduced the computational load by taking advantage of the special structure of the trellis, was therefore realized by an application–specific integrated circuit Dawid et al.[3]. Applying the Cellular Neural Network (CNN) to the Viterbi decoder is another approach for achieving full parallelism Kim et al.[4]. Kim et al. proposed a decoding structure with a circularly connected CNN cell array for a fully parallel Viterbi decoder. This structure, corresponding to CNN analog processing array, had specific features as lower power consumption, no path memory requirement and no trace–back requirement. 139 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece With the developing of electronic design automatic tools, Viterbi decoder is possible to be synthesized on an Field Programmable Gate Array (FPGA) chip Yeh et al.[5], Hema et al.[6]. Viterbi decoders have designed by different hardware description languages such as VHDL Zengliang et al.[7] or Verilog Kavinilavu et al.[8]. In addition, the principal parameters of decoders can vary for diverse applications, for example general wireless communication Ali et al.[9], Wong et al.[10]. However, these FPGA–based Viterbi decoders are always required memories to store the metric and previous state information. Therefore, parallel computing ability of decoder is reduced significantly. In this paper, different from conventional Viterbi decoders, we introduce the CNN–based Viterbi decoding which is expected to implement effectively on FPGA platform. The paper is structured as follows. Viterbi decoding algorithm is reviewed in the next Section while our CNN–based architecture for Viterbi decoder is introduced in Section 3. Section 4 presents obtained results. Finally some conclusions are given in the last Section. 2 Viterbi decoding algorithm In order to understand the Viterbi algorithm, we consider firstly the operation of a convolutional encoder. A convolutional code is described by two parameters R and K, the code rate and the constraint length, respectively. Fig. 1 illustrates a convolutional encoder with R = 1/2 and K = 3. The encoder consists of 2–stages shift register and three modulo–2 adders. Two bit stream outputs are created when one bit stream input is applied. Another way to represent a convolutional encoder is with a state diagram as shown in Fig. 2 where the states of register are noted as 00, 01, 10, and 11. There are two transitions from each state, corresponding to zero input bit (dashed path) and one input bit (solid path). Notice that the output word is written next to each path between two states. Viterbi algorithm bases on trellis diagram for decoding as shown in Fig. 3. The algorithm calculates Hamming distances of paths entering each state node and chooses the path with smallest distance. On the other word, Viterbi algorithm performs minimum distance metric Kim et al.[4] as follows Di,j = min {Dk,l + dij,kl , (k, l) ∈ S} , (1) where Di,j , Dk,l and dij,kl are the accumulated minimum distance metric at node (i, j), the minimum total distance metric from the start node to the node (k, l) and the distance metric for a data bit assigned on the path between the two nodes (i, j) and (k, l), respectively. S is the set of the nodes in the neighbour of the node (i, j). A hardware realization of Viterbi decoder [11] as presented in Fig. 4 is often includes: branch metric unit (BMU), add compare and select unit (ACS), trace back unit (TBU), bit error rate monitor (BER), memory (MEM) 140 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Output 1 Input Reg. 2 Reg. 1 Output 2 Fig. 1. Convolutional encoder with R = 1/2 and K = 3. 01 output 11 00 state 10 00 00 11 10 01 01 input: 1 10 11 input: 0 Fig. 2. Encoder state diagram with R = 1/2 and K = 3. and memory management unit (MMU). MEM stores the accumulated metric and the previous state information while MMU generates addresses and read/write signals for MEM during different phases of operation. It is significant that the store requirement of Viterbi decoder increases exponentially with parameter K Sklar[12] as given by u = h2K−1 , (2) where u is the amount of path storage and h is the length of the information bit path history per stage. A novel structure of Viterbi encoder based on CNN is expected to solve this problem. 3 CNN–based architecture for Viterbi decoding algorithm CNN has introduced by Chua and Yang[13], [14] in which cells are connected directly to nearest neighbours. In addition, because each cell is consisted of a capacitor, resistors and voltage–controlled current, CNN has been implemented successfully through VLSI process Roska and Chua[15]. Applications of CNN have been reported in a considerable amount of fields such as image 141 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Time Input 0 1 2 3 4 5 6 1 0 1 1 0 0 0 7 00 01 10 11 Encode: 11 10 00 01 01 11 00 Decode: 11 10 00 01 01 11 00 Fig. 3. Decoder trellis diagram with R = 1/2 and K = 3. MMU din BMU MEM ACS TBU dout BER BER Fig. 4. Block diagram of a classical Viterbi decoder. processing, solving partial differential equation Chua and Yang[13], modelling complex system Arena et al.[16] or realizing FFT/IFFT algorithm Le et al.[17]. CNN is suitable for real–time parallel signal processing, thus CNN is possible in convolutional decoding. By employing CNN cells at nodes on the trellis diagram, a circularly connected two–dimensional analog CNN enabled very high performance Viterbi decoder Kim et al.[4]. In this part, we present a novel CNN structure for Viterbi decoding algorithm as in Fig. 5. Each cell communicates with two cells at the left–hand side and two cells at the right– hand side. This structure is similar to the trellis diagram where there are two paths entering each state. Cell (i, j) (see Fig. 6) calculates the accumulate minimum distance metric P mi,j of distance metrics received from left–hand side cells P mk,l and the distance metric at the cell Dmi,j P mi,j = min{P mk,l + Dmi,j }. (3) Thus, from Eq. (1) and Eq. (3) it can be proved that the new CNN structure can implement Viterbi algorithm. 142 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig. 5. Modified CNN model which is similar to a trellis diagram. Time t-1 00 minPmk ,l  Dmi, j  t Pmk,l  Dmi, j 00 01 01 10 10 11 11 Fig. 6. Detail connections of cells in the modified CNN. 4 Implementation results Firstly basic cell is designed in Matlab/Simulink as illustrated in the Fig. 7. Here accumulation distance metric is calculated and is transmitted to the next right–hand side cells. Moreover the cell also corrects output bits. After constructing basic cells, they are connected together to create a complete Viterbi decoder. The simulation results show the correction of Viterbi decoder. Using the hardware description language VHDL, CNN–based Viterbi decoder is implemented on EP2C35F672C6 FPGA chip. Detail of basic cell 143 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Table 1. Utilized FPGA resources for realization of CNN–based Viterbi decoder FPGA Components Utilized Total logic elements 2% Total registers 23 Total memory bits 0% Table 2. Comparison with different architectures This work Zengliang et al. [7] Kavinilavu et al. [8] Wong et al .[10] R 1/2 1/2 1/2 1/2 K Memory 3 no 5 yes 7 yes 3 yes Fig. 7. Detail of a basic cell constructed in Matlab/Simulink. is displayed in Fig. 8 by Register–Transfer Level (RTL) Viewer tool. Table 1 presents the resource utilization of our design. Note that no memory requirement is obtained in this case. From the simulation of the decoder by Quartus II software as in Fig. 9, it is obvious that it takes two clock durations to complete the calculation. The comparison of our work and other implementations is summarized in Table 2. Obviously, our architecture has the advanced features comparing to conventional others because of non–memory necessity. 5 Conclusion In this paper, we propose a new architecture of Viterbi decoder based on CNN. Unlike current Viterbi decoders, our architecture does not require memories, the vital factor reduces the calculating rate. As the result our CNN–based Viterbi decoder could be implemented easily on programmable 144 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig. 8. RTL representation of a basic cell synthesised on FPGA. Fig. 9. Simulation of the CNN–based Viterbi decoder which is implemented by VHDL. devices. The new Viterbi decoder has been shown to be able to implement in FPGA chip of Altera. Although the parameters of the decoder are chosen relatively small R = 1/2 and K = 3, they could be modified to support diverse applications in communication. Acknowledgement This work is supported by the Vietnam’s National Foundation for Science and Technology Development (NAFOSTED) under Grant No.102.992010.17. 145 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece References 1.A. J. Viterbi and J. K. Omura. Principles of Digital Communication and Coding, 1979. McGraw–Hill. 2.A. J. Viterbi. Error bounds for convolutional codes and an asymptotically optimum decoding algorithm. IEEE Trans. Inf. Theory, IT–13:260–269, 1967. 3.H. Dawid, G. Fettweis, and H. Meyr. A CMOS IC for Gb/s Viterbi decoding: System design and VLSI implementation. IEEE Trans. Very Large Scale Integration, 4:17–31, 1996. 4.H. Kim, H. Son, T. Roska, and L. O. Chua. Highperformance Viterbi decoder with circularly connected 2D CNN unilateral cell array. IEEE Trans. Circuits Systems, 52:2208–2218, 2005. 5.D. Yeh, P. Chow, and G. Feygin. A multiprocessor Viterbi decoder using Xilinx FPGAs. Proc. Canadian Workshop on Field-Programmable Devices, 1–8, 1994. 6.S. Hema, V. S. Babu, and P. Ramesh. FPGA implementation of Viterbi decoder. Proc. 6th WSEAS Int. Conf. on Electronics Hardware Wireless and Optical Communications, 162–167, 2007. 7.Z. Zengliang and X. Cheng. Design method of Viterbi decoding of convolutional code based on VHDL. Proc. 3rd IEEE International Communication Software and Network, 178–181, 2011. 8.V.Kavinilavu, S. Salivahanan, V. S. K. Bhaaskaran, S. Sakthikumaran, B. Brindha, and C. Vinoth. Implementation of convolutional encoder and Viterbi decoder using Verilog HDL. Proc. 3rd Int. Conf. on Electronics Computer Technology, 297–300, 2011. 9.M. L. Ali, M. B. Karim, and S. M. T. Ahmad. Simulation and design of parameterized convolutional encoder and Viterbi decoder for wireless communication. Proc. 2nd Inte. Conf. on Industrial Engineering and Operations Management, 1206–1211, 2011. 10.Y. S. Wong, W. J. Ong, J. H. Chong, C. K. Ng, and N. K. Noordin. Implementation of convolutional encoder and Viterbi decoder using VHDL. Proc. IEEE Student Conference on Research and Development, 22–25, 2009. 11.Viterbi Decoder Users Guide, Lattice Semiconductor, Inc., Hillsboro, OR, 2005, release 4.2. 12.B. Sklar. Digital Communication: Fundamentals and Applications, 2001. Prentice Hall. 13.L. O. Chua and L. Yang. Cellular neural network: Theory. IEEE Trans. Circuits Systems, 35:1257–1272, 1988. 14.L. O. Chua and L. Yang. Cellular neural network: Applications. IEEE Trans. Circuits Systems, 35:1273–1290, 1988. 15.T. Roska and L. O. Chua. The CNN universal machine: An analogic array computer. IEEE Trans. Circuits Systems II, 40:163–173, 1993. 16.P. Arena, R. Caponetto, L. Fortuna, and G. Manganaro. Cellular neural networks to explore complexity. Soft Computing, 01:120–136, 1997. 17.H. H. Le, N. T. Dzung, and T. M. Hoang. Realization of CNN–based FFT/IFFT algorithms on FPGA. Proc. IEICE International Symposium on Nonlinear Theory and its Applications, 277–280, 2010. 146 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece The Effects of Diffusion of Individuals between Two Single-Species Populations Dumitru D. Deleanu Maritime University of Constanta, Romania E-mail: [email protected] Abstract: In the paper we focused on a general model for the growth of a single-species population with non-overlapping generations. The data we have used correspond to Nicholson’s blow-flies population and lie in the chaotic regime. The population was divided in two groups. If these groups evolve in distinct locations, their behavior is chaotic and, after a few generations, the initial small difference in number of individuals becomes big enough and behaves randomly. The question I want to answer in the paper is: What happens with the two populations if the individuals can migrate in both directions within the time intervals between their reproduction and death? The effect of coupling the two groups consisted in a rich dynamic behavior depending on the coupling strength. It was found that there is a consistent region where the coupling brings out the full synchronization of the two chaotic systems, two transition regions where an intermittent behavior was observed and two peripheral regions where control of chaos is shown to coexist with quasi-periodic and chaotic regimes. Keywords: Single-species populations, Synchronization, Intermittent chaos, Control of chaos. 1. Introduction According to May [1], models for population growth in a limited environment are based on two fundamental premises: a) the populations have the potential to increase exponentially; b) there is a density-dependent feedback that progressively reduces the actual rate of increase. By using a variety of data from field and laboratory populations, some researchers have proposed continuous or discrete models of population growth. The most known of these models is the logistic equation (Verhulst, 1838). Other simple models were introduced by May (1974), Li & Yorke (1975), May & Oster (1976), and Hassel et al (1976). Their models, which refer to single-species population with discrete, nonoverlapping generation, predict that most of the populations show monotonic damping back to an equilibrium following a disturbance, with some exceptions of oscillatory damping or some sort of low-order limit cycles. They concluded that high-order limit cycles and chaos appear to be relatively rare phenomena in naturally occurring single-species populations. Guckenheimer et al (1987) have found that more realistic models of population growth, such as these that include overlapping generations, are more likely to exhibit complex behaviors. If data from laboratory population are used, even for these simple models, it was found that some populations will not exhibit stable equilibrium points but stable cycles or chaotic behavior [2]. That is because the laboratory situation (homogeneous environment, constant food supply, no competitors, no predators) make possible 147 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece an exaggerated non-linear behavior. In this paper we focused on a general model for growth of a single-species population with non-overlapping generations, namely N t +1 = f ( N t ) = λ N t (1 + a N t ) −b (1) where N t and N t +1 are the populations in successive generations, λ is the finite rate of increase and a, b are constant defining the density-dependent feedback term. The values for parameters correspond to Nicholson’s blowflies and lie in the chaotic regime [3]. The population of blowflies was divided in two groups. If these groups evolve in distinct locations their behavior is chaotic and, after a few generations, the initial small difference in number of individuals becomes big enough and behaves randomly. The question I want to answer in the paper is: What happens with the two populations if the individuals can migrate in both directions within the time intervals between their reproduction and death? Divergence of two nearby orbits 1500 1000 DNt 500 0 -500 -1000 -1500 0 5 10 15 20 25 t 30 35 40 45 50 Fig. 1. Divergence of the two isolated populations versus time 2. The Model of Two-coupled Single-species Populations To answer the above question let us hereafter turn our attention towards the following system of two-coupled single-species populations: N t +1 = f (N t ) + c [ f (M t ) − f (N t )] , M t +1 = f (M t ) + c [ f (N t ) − f (M t )] (2) where the coupling parameter c can be thought as the fraction of the two populations which migrate to the neighboring location. Throughout the paper I used the fixed parameter values λ = 60, a = 0.003 and b = 6. The total population N t = 3950 was divided in two unequal groups, N t = 1950 and 148 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece M t = 2000. If no change between the groups was permitted, the initial small difference in number of individuals, ∆N t = 50 , increased quickly and behaved chaotically (see Figure 1). The effect of coupling consisted in a rich dynamic behavior having the main features as follows. 2.1. Complete synchronization If two or more chaotic systems are couple, it is possible that the attractive effect of a suitable coupling to counterbalance the trend of the trajectories to separate due to chaotic dynamics. Synchronization of chaotic systems can be explained by the suppression of expanding dynamics in the state space transversal to the synchronization manifold (here M t = N t ). It is natural then to ask for which values of coupling strength c the two systems will oscillate in a coherent and synchronized way. Laureano et al [4] have demonstrated that, for this kind of coupling, the range of synchronization (in the linear approximation) is given by 0.5 (1 − exp(−λ u ) ) < c < 0.5 (1 _ exp(−λ u ) ) (3) where λ u is the Lyapunov exponent for the uncoupled map f. For our data it was found that λ u ≅ 0.35 , so c ∈ (0.15 ; 0.85) = S .As an example, let consider c = 0.16 . The synchronization takes place after 200 generations (see figure 2). Synchronization in time 10 8 6 4 DNt 2 0 -2 -4 -6 -8 -10 0 50 100 150 200 t 250 300 350 400 Fig. 2. Evolution to synchronous state for c = 0.16 Each of the systems shows chaos and their states are identical at each moment in time (full synchronization). To verify that the synchronous state is chaotic, a Lyapunov exponent versus coupling strength diagram was considered (see Figure 3. 149 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Dependence λ= λ(c) 0.8 0.6 0.4 0.2 0 λ -0.2 -0.4 -0.6 -0.8 -1 -1.2 0 0.1 0.2 0.3 0.4 0.5 c 0.6 0.7 0.8 0.9 1 Fig. 3. Lyapunov exponent versus coupling strength If c is chosen deep inside the interval S, the synchronous state is reached after only few steps (see Figure 4). Otherwise, if c is chosen near the borders of S the synchronization is hard to obtain, a lot of steps being necessary (e.g. 2000 steps for c = 0.15 . Synchronization in time 10 8 6 4 DNt 2 0 -2 -4 -6 -8 -10 0 5 10 15 20 25 t 30 35 40 45 Fig. 4. Evolution to synchronous state for c = 0.25 150 50 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 2.2. Intermittent chaos If the coupling strength c falls short of the critical value c crit = 0.15 the synchronized state M t = N t becomes unstable and an intermittent dynamics is observed. Figure 5 shows the time evolution of the transverse coordinate DN t = N t − M t for c = 0.1467 . The time periods of synchronicity are interrupted by aperiodic chaotic bursts. Synchronization in time 800 600 400 DNt 200 0 -200 -400 -600 0 1000 2000 3000 4000 5000 t 6000 7000 8000 9000 10000 Fig. 5. Time periods of synchronicity interrupted by aperiodic chaotic bursts ( c = 0.1467 ) Synchronization in time 1000 800 600 400 DNt 200 0 -200 -400 -600 -800 -1000 0 0.2 0.4 0.6 0.8 1 t 1.2 1.4 1.6 1.8 2 Fig. 6. A completely erratic state for c = 0.14 151 4 x 10 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece The basic intermittency mechanism comes from the competition between the trajectory instability of chaotic elements and the synchronization tendency due to the diffusion-type coupling [8]. For c = 0.14 the chaotic bursts were already merged so the synchronization started to dissolve into a completely erratic state (see Figure 6). 2.3. Stabilization to an ordered state Outside the interval of synchronization the dynamics is quite complicated. For very small values of c (weak coupling) the system behaves chaotically, the N t values being distributed over an entire interval. By increasing c the chaotic distribution of N t comes undone in strips, thinner and thinner (see Figure 7). Bifurcation diagram 1400 1200 1000 Nt 800 600 400 200 0 -200 6.8 7 7.2 7.4 c 7.6 7.8 8 -3 x 10 Fig. 7. A part of the bifurcation diagram N t = N t (c) At c ≅ 0.007 the system entered a periodic regime, and was subjected to a sequence of changes from a 2 n - period cycle to a 2 n −1 - period cycle. A 8period cycle was obtained for c ∈ (0.0072 ; 0.0080) (see Figure 8). Then, a quasi-periodic regime with two strips appeared (Figure 9) which, in its turn, was changed by a 2-period cycle for c ∈ (0.013 ; 0.11) . This periodic regime is interrupted by windows corresponding to a 4-period cycle or even to thin windows of chaotic regime. Beginning with c ≅ 0.1 the number of steps required for stabilization to the 2period cycle became bigger and bigger so, finally, the chaotic regime was reached. An analogous discussion can be done for c ∈ (0.85 ; 1) . 152 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Time series 1800 Nt 1600 Mt 1400 1200 Nt, Mt 1000 800 600 400 200 0 -200 1980 1982 1984 1986 1988 1990 t Fig. 8. Time evolution of N t , M t 1992 1994 1996 1998 2000 for c = 0.0079 (8-period cycle) Time series 1250 1200 1150 1100 Nt 1050 1000 950 900 850 800 1500 1550 1600 1650 1700 1750 t 1800 1850 1900 1950 2000 Fig. 9. Time evolution of N t for c = 0.99 (quasi-period regime) 153 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 3. Conclusions The dynamics for many biological populations, which breed seasonally and have non-overlapping generations, are described by a density-dependent relation of the form N t +1 = f (N t ) = λ N t (1 + a N t )−b . If data from laboratory tests are used it was found that populations can exhibit even a chaotic behavior. Two almost identical populations, living in distinct locations, evolved so that the initial small difference in number of individuals became big enough and behaved randomly. If the individuals representing the two populations could migrate in both directions within the time intervals between their reproduction and death then a rich dynamic behavior depending on the coupling strength was observed. It was found that there is a consistent region where the coupling brings out the full synchronization of the two chaotic systems, two transition regions where an intermittent behavior was observed and two peripheral regions where control of chaos is shown to coexist with quasi-periodic and chaotic regimes. References 1. M. Doebeli. Intermittent chaos in populations dynamics, J. Theor. Biol., vol. 166, 325330, 1994. 2. M. P. Hassell, J.H. Lawton and R.M.May. Patterns of Dynamical Behavior in singlespecies populations. The Journal of Animal Ecology, vol. 45, no. 2, 471–486, 1976. 3. M. P. Hassell. Density-dependence in single-species populations. The Journal of Animal Ecology, vol. 49, no.1, 283-295, 1975. 4. R. Laureano, D.A. Menales and M.A.M. Ferreira. Efficient synchronization of onedimensional chaotic quadratic maps by different coupling terms. Journal of Mathematics and Technology, 5–12, 2010. 5. R. M. May. Bifurcation and dynamic complexity in simple ecological models. The American naturalist, vol. 110, no. 974, 573–599, 1976. 6. L. D. Mueller and H. J. Ayala. Dynamics of single-species population growth: stability or chaos? Ecology, vol. 62, no.5, 1148–1154, 1981. 7. L. M. Pecora and T.L. Carrol. Synchronization in chaotic systems. Phs. Rev. Lett., vol. 64, no. 8, 821-828, 1990. 8. H.G. Schuster and W. Just. Deterministic chaos: An introduction. Wiley-VCH Verlag. 2005. 154 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Thermal Emission, Density and Light Scattering in water Dimitrios Dellaportas Mechanical Engineering Technology, Pedagogics, Master of Physics, Principal of 1st EPAL School of Pireas e-mail:[email protected] Abstract: At temperatures above absolute zero, particles can emit as well as absorb and scatter electromagnetic radiation. Emission does not strictly fall within the bounds. It is more akin to such phenomena as luminescence than to elastic scattering. Consider an enclosure of dimensions large compared with any wavelengths under consideration, which is opaque but otherwise arbitrary in shape and composition. If the enclosure is maintained at a constant absolute temperature T, the equilibrium radiation field will be isotropic, homogenous, and unpolarized. At any point the amount of radiant energy per unit frequency interval, confined to a unit solid angle about any direction, which crosses a unit area normal to this direction in unit time is given by the Plank function Pe= ћω³/4π³c²· 1/ exp(ћω/kBT)-1. Optical properties of liquids are similar in many ways to those of solids. Electrically, there are metallic liquids such as mercyry and molten iron, but the majority of common liquids are nonmetallic. As an illustration of a liquid we have chosen H2O, a ubiquitious substance on our planet; water dominates not atmospheric processes but the chemistry of life. The optical properties of water had been studied for centuries; the modern results are scattered throughout the literature of many scientific fields. Dielectric functions calculated from refractive indices tabulated by Hale and Querry. Optical constants at wavelengths shorter than 0.2 µm were taken from the work of Kerr et al, where electronic absorption by H2O begins. Dipole relaxation in H2O maintains a high level of absorption from the microwave well in to the infrared; this also is the cause of the major differences between the optical properties of liquid and solid water. Introduction Let us contemplate the sea water near the costs. This water includes some particles of sodium chloride, the quantity of them is increasing or decreasing according to meteorological conditions and to electromagnetical moon’s functions. So the thermal emission and light scattering is up to particles assemple, phenomenon that affects the life of the aquatic organisms. Scattering or extinction measurements allow to determine size and distribution of sizes, shape and orientation, and chemical composition of the particles. If we consider that the salt (sodium chloride) consists of small scaterring and absorbing particles, we write the complex refractive index of the salt as m=ñ-iñ'. The light scattered from a fluid containing optically isotropic molecules. Many new features arise in the spectrum of light scattered from fluids containing optically anisotropic molecules, of course. In the dense water of these arises as 155 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece to how the collective motions of the fluid can be described by the rotation and how the thermal emission is changed. Let us note that the depolarization spectrum Ivн(ω) consists of Lorentzian bads, all centered at zero frequency. If we have symmetric top rotors, the spectrum consists of a single band with a width [q2D+6Θ] which depends only on the translational self- diffusion D and on the rotational diffusion coefficient Θ. Then, the spectrum appears to be split. So, it is presented a spectrum theory from a molecular point of view. 0PTICAL PROPERTIES OF WATER: Optical properties of water are similar in many ways to this of solid. Electrically, there are metallic liquids such as mercury and molten iron, but the majority of common liquids are non metallic. As liquid we have chosen the water because it is much important for the chemistry of life. During the centuries, the optical properties of water have been studied by the scientists. Irvine and Pollack (1968), Hale and Querry (1973) have calculated n and k for small particle calculations .Dielectric functions calculated by Hale and Querry. Instant to other scientists who use electron volts, we use a wavelength scale. As was done for ClNa, e' is plotted linearly and e'' is plotted logarithmically. We also include on the e'' plot some results for solid H2O from Irvine and Pollack. The comparison of ice and water gives the field of the similarities and differencies between the solid and the liquid phase. The insulating liquid H2O is not unlike the insulating solid ClNa. The electrons and vibrational excitation regions for both materials are well separated by a highly transparent region. Electronic absorption by H2O begins at about 0.2 µm. The major difference between clear water and salt water occurs longward of the vibrational absorption region. Dipole relaxation of H2O maintains a high level of absorption from the microwave well well into the infrared. This is the cause that makes the differencies between the optical properties of liquid and solid water. 156 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece THERMAL EMISSION In our case, let us have a spherical particle that is placed in the enclosure, then in equilibrium the distribution of radiation is uncharged. A spherical surface of radius R is surrounding the particle that has radius a, each element has entropydS of which is the source of a nearly plane wave that illuminates the 157 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece particle with irradiance Pe.dS/R2. The amount of energy absorbed per unit time by the particle is; ∫₀∞ ∫s PeCabs .d/S. dω=4π∫₀∞ PeCabc.dω. CALCULATION OF THE ENTROPY PRODUCTION The sunrays are scattering on the water particles and this scattering is related to the water density, because of the salt (ClNa). By Gibbs equation we have: n TdS=dE + PdV-∑µidni where S,E,P,V,µi,di, are the total entropy, internal ene i=1 rgy, pressure, volume, chemical potential of species I and the number of moles of salt respectively. The integrated form of this equation: n TS=E+PV-∑µini. i=1 Dividing this later equation by V, we get an expression for the local densities s,e,ci, corresponding to S,E,ni. So we have: Ts=e+P-∑µici. i 158 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Then, substituting sV to S, eV to E, ciV to ni we have: Td(Vs)=d(Ve)+PdV-∑µid(Vci) i So; V(TdS-de+∑µidci)=dV(-Ts+e+p-∑µici). The right part of this equation is zero and we obtain: n Tds=de-∑ -µidci i=1 This equation is a view of entropy density, internal energy density, and the molecules concentration. Alternatively, this equation can be expressed: TdS=dq-∑µidci i=1 dq is the increment of pure heat per unit volume. CONCLUSION These equations make evidence how much serious and important for the seawater life is the influence of the light scattering, the density of the water, and the thermal emission. The density of the water can be changes because of the pouring rain, the snow’s dissolving, the atmosphere thermal changes. If the seawater has high density of sodium chloride, then, it has the possibility to gather much more energy by the sunshine. It creates the presuppositions for growth of phytoplangton, zooplangton and fish. References: 1. Asano S., 1979. Light Scattering properties of spheroidal particles 2. Crosswhite, H. M. and H. W. Moos 1967. Optical Properties of Ions in Crystals, Interscience, New York. 3. Faxvog, F. R., and D. M. Roessler, 1981. Optical absorption in thin slabs and spherical particles. 4. Bruce J. Berne and Robert Pecora 2000. Dynamic Light Scattering, Dover Publications New York. 5. Gerald C. Pomraning 2005. The Equations of Radiation Hydrodynamics, Dover Publications New York. 6. H. C. van de Hulst 1981. Light Scattering by Small Particles, Dover Publicatiions, New York 159 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 160 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Bifurcation Analysis of Nociceptive Neurons* Olga E. Dick and Boris V. Krylov Pavlov Institute of Physiology of Russian Academy of Science, St. Petersburg, Russia E-mail: [email protected] Abstract: As known, modification of specified slow sodium channels (Na1,8) in the membrane of nociceptive neurons is the basis of the pain perception by the human brain. The work is devoted to determination of parameters of the channels most sensitive to perceiving the painful signals. Using the bifurcation analysis of the model system describing the impulse activity of the membrane of mammalian nociceptive neuron we partition the parameter planes into the regions corresponding to stable and unstable periodic solutions. The left boundary of the region corresponds to subcritical Hopf bifurcation and emergence of the rough excitation in the form of large amplitude oscillations. The right boundary relates to supercritical Hopf bifurcation and appearance of the smooth excitation in the form of small large amplitude. Integrating inside the region of stable solutions we obtain the relationship between the parameter and frequency values. Bifurcation parameters such as the effective charge transfer of the activation gating system of the sodium channels and the maximal conductance of the channels play the main role in increasing the frequency and, hence, in transformation of the unpainful stimulus into the painful one. The results explain ionic mechanisms of action of analgesic drugs having high selectivity to NaV1,8 channels independently of the primary target of action. Keywords Hopf bifurcation, Membrane model, Sodium channels, Nociceptive neuron. 1. Introduction It is known that in response to injury of nervous system nociceptive neurons can become hyper-excitable and generate spontaneous impulse activity of unusual frequency [1]. Perception of painful feeling is connected with activation of peripheral nociceptors recording painful signals and transmitting them by afferent nerve fibers to nociceptive neurons soma of which are in spinal ganglia. Low frequency of nerve impulses carries information about adequate tactile action and rise of the frequency for amplification of signal testifies about possible injury [2]. Slow sodium NaV1,8 channels are considered significant in generation of painful feeling since the enhancement of synthesis and functional activity of these channels is related to hyper-excitability of nociceptive neurons and high frequency neurophatic pain [3, 4]. The failure in the synthesis of the channels causes the reduce of neurophatic pain [5]. Modulation of activity of the channels by mediators of inflammation can lead to pathological state such as hyperalgesia (an increase of painful sensitivity) [6]. Hyperalgesia is removed by agents descending impulse activity of NaV1,8 channels [7].That is why these * Paper included in Chaotic Systems: Theory and Applications, 161 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece agents are believed as the analgesic highly selective drugs [2].The aim of the work is to answer the question: what parameters of the slow sodium NaV1,8 channels do maximal influence on pain signaling transduction?To answer the question it is necessary 1) to study relations between these parameters, an applied external stimulus and a type of stable solution of the model system describing the impulse activity of the nociceptive neuron;2) to clarify what parameters do determine the possibility of the nociceptive neuron to generate spontaneously a signal of a painful range frequency? 2. The model We have used the space-clamped Hodgkin-Huxley type model: cm dE = I − g Naf m 3 h( E − E Na ) − g K n 4 ( E − E K ) − dt g L ( E − E L ) − g Nas ms3 hs ( E − E Na ), dm = α m ( E )(1 − m) − β m ( E )m, dt dh = α h ( E )(1 − h) − β h ( E )h, dt dn = α n ( E )(1 − n) − β n ( E )n, dt dms = α m s ( E )(1 − ms ) − β ms ( E )ms , dt dhs = α h s ( E )(1 − hs ) − β hs ( E )hs , dt (1) where E is the membrane potential, the variables m, h, n, ms, hs represent the probabilities of activation and inactivation of fast sodium, potassium and slow sodium channels, respectively. gNaf =40 nS, gK = 20 nS, gL = 5 nS, ENa =55 mV, The constants cm=20 pF, EK= - 85 mV, EL= -70 mV are the membrane capacitance, the maximal conductance of the fast sodium, potassium and leakage ions channels and the reversal potentials for Na+, K+ and leakage ions. The voltage-dependent expressions α m (E) = 0.115(1 + e( E + 70) / 10 ) , 1 + e ( E + 40) / 42) α h ( E ) = 0.012(1 + e − ( E + 43) / 10 ), α n (E) = 0.006( E + 45) , 1 − e − ( E + 45) / 12) β m ( E ) = 0.015(1 + e ( E + 25) / 8 ), 1.32 , 1 + 0.2e ( E + 10) / 7) β n ( E ) = 0.13e − ( E + 45) / 30 ), βh (E) = 162 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece α m s ( E ) = e k1( E + G ) + d1 , β ms ( E ) = e k 2 ( E + G ) + d 2 , α h s ( E ) = 0.0015e − ( E + 4) / 30 ), β hs ( E ) = 0.01 − 1 + 0.2e ( E + 10) / 7) describe rates of transfer of the activation and inactivation gating structures of ionic channels between the closed and open states. According to the Boltzmann’s principle for the channel with the two-state openclosed structure the ratio of the number of open channels (N0) to the number of closed channels (NC) is determined by Z e ( E − E ) / kT No ms , = = e eff N c 1 − ms where Zeff is the effective charge of the activation gating structure (in electron units) coupled with conformational change of the gating structure during the ion transfer through the membrane, k is the Boltzmann’s constant, T is the absolute temperature, e is the electron charge, E is the membrane potential such that N0=NC. Then at E = E for the activation gating structure of the slow sodium channels one can write α m = β m , whence it follows that the effective charge value of S S the activation gating structure can be gained as Z eff = kT (k1 + k2 ) . e 3. Partition of the model parameter space into regions of qualitatively different solutions To obtain relationship between the type of stable solution of the system, its parameters and an applied external stimulus it is sufficient to find points belonging to the boundary partitioning the parameter space of the model system into the regions of the qualitatively different types of stable solutions (steady states and stable periodic oscillations). For constructing the boundary the method of bifurcation analysis is applied. On the I axis there are at least 3 bifurcation points (I0<I1<I2) [8]. For I<I0 and I>I1 there is a one-to-one correspondence between the type of steady state (unstable or stable) and the presence or absence of a stable periodic solution. For I≤I0 and I≥I2 the steady state is stable and a limit cycle does not exist. While the bifurcation parameter I increases in interval (I0<I≤I1) the steady state is stable and a stable and unstable periodic solutions coexist appearing via fold limit cycle bifurcation. The unstable periodic solution shrinks down to the rest state and makes it lose stability via subcritical Andronov-Hopf bifurcation. 163 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Therefore, for I>I1 the stable periodic oscillations of large amplitude exist both with decreasing and increasing I value. But for I0<I≤I1 the stable limit cycle of large amplitude is exhibited only with decreasing I Since for the Hodgkin-Huxley type system I0≈I1 for all the physiologically possible parameter values [9], the value of I1 can be used as an approximate value of I0. That is why the task of finding the boundary of qualitatively different types of stable solutions can be reduced to the more simple numerical task of constructing the boundary of various steady states (stable and unstable). We write system (1) in the form dx = F ( x, p, I ), dt (2) where x=(E, m, h, n, mS, hS) is a vector of the phase coordinates, p=(gNaS, k1, k2) is a vector of parameters which can be considered as bifurcation ones. The method for determining the boundary points of the region of stable periodic solutions is reduced to the sequence of operations: 1) finding the equilibrium state of system (2) as a unique solution x0 (p, I) of the equation F(x, p, I)=0, 2) calculating the eigenvalues {λi ( p, I )}16 of the Jacobian matrix     ∂F i J ( p, I ) =  , i, j = 1,...,6  ,   ∂x j x = x0( p , I )   3) finding the parameter values satisfying the Hopf bifurcation, namely, arising of a pair of purely imaginary eigenvalues λ1 = iω , λ2 = −iω , λ3 < 0, λ4 < 0, λ5 < 0, λ6 < 0 . To determine values I1 (p) and I2 (p) at which maximal real part (λm ( p, I )) becomes equal to zero the following algorithm is used. 1) The interval [I0, IK] of possible values is discretized with k consecutive subintervals of length ∆. In the case of existence even though one solution I1 (p) of equation λm ( p, I ) = 0 , 164 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece the i - subinterval involving the solution, is determined by the consecutive search beginning from the left side of the interval. The value of this solution is determined by linear interpolation. 2) The value of I2 (p) is calculated in the interval [i*∆, IK] by the method of bisection followed by linear interpolation. The numerical solution of system (1) inside the obtained region of stable periodic solutions is found by a fourth-order Runge-Kutta method with a modified variable step size and Gear algorithm. The frequency of the periodic solution is calculated by the time values corresponding to local maxima. 4. Results and discussion To elucidate the role of slow sodium channels in generation of the painful stimulus the maximal conductance of the slow sodium channels (gNaS), the effective charge transfer of the activation gating system of the channels (Zeff) and the shift (G) of the activation curve along the membrane potential axis have been used as variable parameters. The family of the plane sections of the boundary partitioning the parameter space (gNaS, Zeff, I) into the regions of stable and unstable steady states are given in Fig.1 a, b. Inside the each found region the steady state is unstable and there is a stable limit cycle corresponding to stable periodic solution. Stain-stepping effect of the left boundary of the region is related to features of arising limit cycles on the left and right sides of the boundary. The left boundary of the region corresponds to subcritical Hopf bifurcation and emergence of the rough excitation in the form of large amplitude oscillations. The right boundary relates to supercritical Hopf bifurcation and appearance of the smooth excitation in the form of small large amplitude. As is seen, if gNaS =0, periodic oscillations are absent at any stimulus value. The minimal value of gNaS such that the oscillations emerge is equal to 14,9 nS at the stimulus – 142,5 pA and the value grows when Zeff increases. When the effective charge is less than 5e the periodic oscillations arise only by hyperpolarizing stimulus (I<0 pA). With increasing Zeff the steady periodic solutions region extends significantly and shifts in direction of depolarizing stimulus (I>0 pA). Integrating inside the constructed regions we obtain the relationship between the model parameter and frequency values. The examples of steady periodic solutions are represented in Fig.2 and Fig.3. The periodic oscillations emerging on the left boundary of the region have large amplitude and small frequency. When moving inside the region from left to right an amplification of the external stimulus tends to change in amplitude and frequency of the nociceptive neuron. In other words, for the constant maximal conductance of the slow sodium channels and effective charge transfer of the activation gating system an enhancement of the external stimulus leads to the increase of the frequency of periodic oscillations and then their disruption. 165 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 400 a 350 gnas (nS) 300 Zeff=6.5 Zeff=6 250 Zeff=5.5 200 Zeff=5 150 Zeff=4.5 100 50 Zeff=3.5 0 -300 Zeff=4 -200 -100 0 current I (pA) 100 200 600 Zeff=10 gnas (nS) 400 300 b Zeff=9 500 Zeff=8 Zeff=7 200 100 0 -100 0 100 200 current I(pA) 300 400 Fig. 1. The examples of the plane sections of the boundary partitioning the parameter space (gNaS, Zeff, I) into the regions of stable and unstable steady states. Each section is constructed with 800 points on the (gNaS, I) plane corresponding to 700 net values of the parameter gNaS .Values G=10 mV, Z eff 11 = {3.5e , 4e , 4.5e , 5e , 5.5e , 6e , 6.5e , 7e , 8e , 9e , 10e } . 1 { } The periodic oscillations emerging on the left boundary of the region have large amplitude and small frequency. When moving inside the region from left to right an amplification of the external stimulus tends to change in amplitude and frequency of the nociceptive neuron. In other words, for the constant maximal conductance of the slow sodium channels and effective charge transfer of the activation gating system an enhancement of the external stimulus leads to the increase of the frequency of periodic oscillations and then their disruption. 166 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Zeff=6.5, I=0 pA, f=7 Hz 50 E(mV) E(mV) 50 0 -50 0.5 1 0 50 I=75 pA, f=44 Hz E(mV) E(mV) 0 -50 0 50 I=50 pA, f=38 Hz 0 -50 0.5 time(sec) 1 I=100 pA, f=46 Hz 0 -50 0 0.5 time(sec) 1 0 0.5 time(sec) Fig. 2. The examples of steady solutions for Z eff = 6.5e , 1 g NaS = 100nS and various values of stimulus. 50 E(mV) E(mV) 50 I=0 pA, gnas=100 nS, f=7 Hz 0 -50 0.5 time(sec) 1 0 gnas=200 nS 50 E(mV) E(mV) 0 -50 0 50 gnas=150 nS, f=11 Hz 0 -50 0 0.5 time(sec) 1 gnas=250 nS 0 -50 0.5 time(sec) 1 0 0.5 time(sec) Fig.3. The examples of steady solutions Z eff = 6.5e , 1 I = 0 and different values of g NaS . The periodic oscillations emerging on the left boundary of the region have large amplitude and small frequency. When moving inside the region from left to right an amplification of the external stimulus tends to change in amplitude and frequency of the nociceptive neuron. In other words, for the constant maximal conductance of the slow sodium channels and effective charge transfer of the activation gating system an enhancement of the external stimulus leads to the increase of the frequency of periodic oscillations and then their disruption. 167 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece When moving inside the region from bottom to top a growth of the maximal conductance of the slow sodium channels for the constant stimulus value results in increase of the frequency, period doubling and also failure of periodic oscillations (Fig.3). Thus, both factors, namely, decrease of effective charge transfer of the activation gating system of the slow sodium channels for the constant maximal conductance of the channels and decrease of the maximal conductance of the slow sodium channels for the constant stimulus decline the frequency of impulse activity. Since an increase in the frequency of impulse activity of nociceptive neurons is related to the emergence of neuropathic pain, our findings indicate the direction of looking for chemical agents possessing analgesic properties. 3. Conclusion The form of the constructed regions demonstrates that ability of each parameter to be bifurcation one significantly depends on the other parameter values. Thus, the conclusions about bifurcation properties of the system parameters are determined by the investigated point in the parameter space. The character of changes in the system solutions and in the frequency of periodic solutions can be used in searching of chemical agents aimed for selective removal of neuropathic pain. References 1. J.K. Diss, S. P. Fraser, M. B. Diamoz Voltage-gated Na channels: 2. 3. 4. 5. 6. 7. 8. multiplicity of expression, plasticity, functional implications and pathophysiological aspects. Eur. Biophys. J. 33: 180-193, 2004. E.A. Karymova, I. E. Katina, V. B. Plakhova, et. al. Possible coding mechanism of nociceptive signals: role of slow sodium channels. Sensory systems, 22: 257-270, 2008. A.E. Goldin Evolution of voltage-gated Na channels. J. Exper. Biol. 205: 575-584, 2001. J. Lai, F. Porreca, J.C. Hunter, et al. Voltage-gated sodium channels and hyperalgesia. Ann. Rev. Pharmacol. Toxicol. 44: 371-397, 2004. S.G. Waxman The molecular pathophysiology of pain: abnormal expression of sodium channel genes and its contrubutions to hyperexcitability of primary sensory neurons. Pain Aug. Suppl. 6: 133-140, 1999. S.G. Waxman, T.R. Cummins, S. D. Dib-Hajj, et.al. Voltage-gated sodium channels and the molecular pathogenesis of pain. J. Rehabil. Res. Dev. 37: 517-528, 2000. N. Ogata, Y. Ohishi The molecular diversity of structure and function of the voltage-gated Na channels. Jpn. J. Pharmacol. 88: 365-377, 2002. B.Hassard Bifurcation of periodic solutions of the Hodgkin-Huxley model for the squid giant axon. J. Theor. Biol. 71: 401-420, 1978. 168 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 9. Y.A. Bedrov, G. N. Akoev, O.E. Dick. Partition of the Hodgkin-Huxley type model parameter space into the regions of qualitatively different solutions. Biol. Cybern. 66: 413-418, 1992. Acknowledgements We thank E.A. Karymova for her help with experimental recordings of ionic channels. 169 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 170 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Modeling of a Scenario of Transition to Chaos in Plasma through Sub-harmonic Bifurcation Dan G. Dimitriu1,a, Silviu Gurlui1,b, Maricel Agop2,c 1 Faculty of Physics, “Alexandru Ioan Cuza” University, Iasi, Romania Department of Physics, “Gh. Asachi” Technical University, Iasi, Romania a E-mail: [email protected] b E-mail: [email protected] c E-mail: [email protected] 2 Abstract: We report on theoretical modeling of a scenario of transition to chaos in plasma through sub-harmonic bifurcations, based on scale relativity theory. Experimental results show that a discharge plasma which self-structures as concentrically multiple double layers is able to evolve towards chaos through a cascade of spatio-temporal subharmonic bifurcations when the applied constrain is gradually increased. By considering that the plasma particles (electrons, ions, neutrals) move on continuous but nondifferentiable curves (because of the collisions), i.e. on fractal curves, we develop a model based on scale relativity theory capable to explain the self-structuring of the electrically conductive plasma as multiple double layers, as well as the transition to chaos through sub-harmonic bifurcations of these structures when the applied constrain increases. Thus, by considering that the complexity of the physical processes in plasma is replaced by fractality, it is no longer necessary to use the whole classical “arsenal” of quantities from the standard physics (differentiable physics). The plasma system will behave as a special interaction-less “fluid” by means of geodesics in a fractal space. In such context, by using the fractional revival formalism, a Reynold’s fractional criterion of evolution to chaos through cascade of spatio-temporal sub-harmonic bifurcations was stated. A good agreement between the experimental results and those provided by the theoretical model was obtained. Keywords: Chaotic modeling, Fractal, Scale relativity theory, Pattern formation, Subharmonic bifurcations. 1. Introduction Concentrically multiple double layers are complex nonlinear potential structures in plasmas, consisting of two or more double layers attached to the anode of a dc glow discharge [1] or to a positively biased electrode immersed into plasma [2]. They appear as several bright and concentric shells attached to the electrode. The number of layers depends on the background gas, on its pressure, on the electrode voltage and on the discharge current. The axial profile of the plasma potential has a stair step shape, with potential jumps close to the ionization potential of the working gas [2]. At high values of the voltage applied to the electrode, the multiple double layer structure passes into a dynamic state which consists of periodic disruptions and re-aggregations of the constituent double layers [2]. In many systems where deterministic chaos arises, spatial and temporal structures were also experimentally observed. For time scales large with respect 171 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece to the inverse of maximum Lyapunov exponent, deterministic trajectories can be replaced by families of potential trajectories and the concept of definite positions by that of probability density. This allows the description of chaos in a stochastic way by a diffusion process as described in [3]. By considering that the particles movement takes place on continuous but non-differentiable curves, i.e. on fractal curves, the scale relativity theory approaches the chaotic effects in the same way as in [3], except diffusion, which becomes a spatio-temporal scale dependent process [4]. The theory of scale relativity (SR) applies the principle of relativity to scale transformations. In such conjecture, the principle of SR requires that the fundamental laws of nature apply whatever the state of scale of the coordinate system. In order to include the non-differentiable quantum motion, the quantum space-time must be considered relative and fractal, as described in [5]. In this theoretical framework, it is not necessary to endow a point particle with mass, energy, momentum or velocity. The particle may be reduced to and identified with its own trajectory [4]. Here we report on experimental results revealing how the electrically conductive plasma can self-structure in the form of concentrically multiple double layers. We show that these new structures evolve towards chaos through a cascade of spatio-temporal sub-harmonic bifurcations when the applied constraint in form of the voltage applied on a supplementary electrode increases. A theoretical model based on the scale relativity theory is developed, and its results are found to be in good agreement with the experimental data. 2. Experimental Results The experiments were performed in a plasma diode, schematically presented in Fig. 1. Plasma is created by an electrical discharge between the hot filament (marked by F in Fig. 1) as cathode and the grounded metallic tube as anode. The plasma was pulled away from equilibrium by gradually increasing the voltage applied to a tantalum disk electrode (marked by E in Fig. 1) with 1 cm diameter, under the following experimental conditions: argon pressure p = 10-2 mbar, plasma density npl = 109 cm-3 and electron temperature kTe ≅ 2 eV. When the voltage on the electrode reaches VE ≅ 55 V, a double layer structure appears in front of the electrode (see photo in Fig. 2(a)). The structure is in dynamic state, phenomenon extensively described in [6]. The FFT of the oscillations of the current collected by the electrode is shown in Fig. 3(a). By a further increase of the voltage on the electrode, new double layers develop in front of the electrode, giving rise to a dynamic multiple double layer structure (see photos in Fig. 2(b)2(e)). Simultaneously with every new double layer formation, a new subharmonic appears in the FFT spectrum of the current oscillations (see Fig. 3(b)3(f). Thus, we recorded, in fact, spatio-temporal bifurcations in the plasma system (sudden changes in the spatial symmetry and in the temporal dynamics of the plasma system). At high values of the applied potential, the plasma system passes into a chaotic state, characterized by uncorrelated and intermittent oscillations (see FFT’s of these oscillations in Fig. 3(g)). 172 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece U1 U2 R2 F E A R1 PS Fig. 1. Schematic of the experimental setup (F – filament, A – anode, E – supplementary electrode, U1 – power supply for heating the filament, U2 – power supply for discharge, PS – power supply for supplementary electrode bias, R1, R2 – load resistors). (a) (b) (c) (d) (e) Fig. 2. Photos of the multiple double layers in different stages (at different increasing value of the potential applied on the electrode) of its formation. 3. Theoretical Model Assuming that the discharge plasma particles (electrons, ions, neutrals) motions take place on continuous but non-differentiable curves, i.e. on fractal curves of fractal dimension DF, the dynamics of such a system can be described by the complex operator (see [7]): 173 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece (a) (b) (d) (e) (c) (f) (g) Fig. 3. FFT of the oscillations of the current collected by the supplementary electrode, at different increasing values of the voltage applied on it. dˆ ∂ ˆ (2 D )−1 = + V ⋅ ∇ − iD (dt ) F ∆ dt ∂t (1) Vˆ = V − iU (2) where The real part V of the complex field Vˆ represents the standard classical velocity, which is differentiable and independent of resolution dt , while the imaginary part U is a new quantity arising from fractality, which is nondifferentiable and resolution-dependent. The quantity D is Nottale’s coefficient and corresponds to the fractal – non-fractal transition. The geodesics equation (a generalization of the first Newton’s principle) can be written in the form: dˆVˆ ∂Vˆ + Vˆ ⋅ ∇ Vˆ − η ∆Vˆ = dt ∂t ( ) 174 (3) Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece which is a Navier-Stokes type equation with the imaginary “viscosity coefficient”: η = iD (dt )(2 D F )−1 (4) If the motions of the fractal fluid are irrotational, i.e. ∇ × Vˆ = 0 , we can choose V̂ on the form: (2 D )−1 Vˆ = −2iD (dt ) F ∇ lnψ (5) with ψ being the scalar potential of the complex velocity. Then, by substituting (5) in (3) and using the method described in [9] it results: dVˆ (2 D )−1  ∂ lnψ (2 D )−1 ∇ψ = −2iD (dt ) F ∇ − 2iD (dt ) F dt ψ  ∂t   = 0  (6) This equation can be integrated in a universal way and yields: Lˆ ψ = 0 (4 D )− 2 (2 D )−1 ∂ Lˆ = 4D 2 (dt ) F ∆ + 2iD (dt ) F ∂t (7a,b) up to an arbitrary phase factor which may be set to zero by a suitable choice of the phase of ψ . Equation (7a) is of Schrödinger type. If the fractal fluid is placed into the external field W, then (7a,b) become: Lˆ 'ψ = 0 W (4 D )− 2 (2 D )−1 ∂ Lˆ ' = 4D 2 (dt ) F ∆ + 2iD (dt ) F − ∂t m0 (7c,d) where m0 is the rest mass of the fractal fluid particle. In our experiment the particle interacts with an external potential which results in a discontinuous change in momentum. This interaction can be theoretically modeled by means of a one-dimensional infinite square-well potential that confines a particle to a box of given width. Localized time-dependent solutions for this problem can be constructed in a very straightforward way from solutions of the free-particle problem [9]. The one-dimensional infinite square-well potential confines a particle to a box of width a and is described by 175 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 0 , W (x ) =  ∞ , x ≤a 2 (8) x >a 2 A particle of mass m0 is placed in this potential. The energy eigenvalues and eigenstates are found by solving the time-independent Schrödinger-type equation (as (7c,d)). The discrete energy eigenvalues are 2  nπ  (2 D )−1 En = 2m0 D 2   , D = D (dt ) F  a  (9a,b) for positive integer n . The scalar potential eigenvalues corresponding to the energy eigenvalues (9a,b) are  2  nπ x  sin  , x ≤ a 2 , n even   a   a  2  nπ x  φn ( x ) =  cos  , x ≤ a 2 , n odd a  a    0, x >a 2  (10) Let’s consider the expression kn2 = n 2 k12 induced by the generalized coherence (in the present context, the physical mean of the generalized coherence refers to the generation of multiple double layers) (see for details [10]). Then, the relation En = E1n 2 , E1 = 2m0 D 2 k12 , k1 = π (11a-c) a expanded around n either in the form En = E1n 2 + 2 E1n (n − n ) + E1 (n − n ) 2 (12) or in the form  (n − n ) (n − n )2  2 + En = En + 4πm0 D   , En = E1n Tβ   Tα induces the characteristic times: 176 (13a,b) Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Tα = 2π m0 D 4π m0 D , Tβ = nE1 E1 (14a,b) Because Tβ is independent of n , expression (14b) defines a universal time scale. By means of relations (14b) and (11b), a characteristic frequency can be associated: k  E  1 f0 = = 1  1  Tβ 2π  2m0  1 2 (15) Let us evaluate expression (15) with respect to the experimental results. Thus, by identifying a with a characteristic length of the double layer, namely the width of the double layer, a = ( 2ε 0V en0 ) , expression (15) takes the form 12 1  en E f 0 =  0 1 2  m0ε 0V    1 2 (16) where E1 is the ion energy, m0 is the ion mass, n0 is the ion density in the double layer and V is the voltage on the electrode. For our experimental conditions f0 ≅ 150 kHz. This value is in a good agreement with those experimentally obtained by us (see Fig. 3). Now, the fractional revival formalism may be applied. A fractional revival of a physical function occurs when a physical function evolves in time to a state that can be described as a collection of spatially distributed physical sub-functions that each closely reproduces the initial physical function shape – see for details [11]. In this context, we write the particle’s initial ( t = 0 ) scalar potential velocity function in the square well as: ψ ( x, t = 0 ) = ψ i ( x ) (17) We expand this scalar potential velocity function by using the relation (10) ∞ ψ i ( x ) = ∑ cnφn ( x ) (18) n =1 with +∞ cn = ∫ φn ( x )ψ i ( x )dx −∞ 177 (19) Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece By using the time scale Tβ , the time evolution of the system is found from the Schrödinger’s type equation (the charge transport takes place on fractal curves – see the operator L̂' from (7c,d)) (4 4D 2 (dt ) DF )− 2 W ∂ 2ψ (2 D )−1 ∂ψ − ψ ≡0 + 2iD (dt ) F 2 ∂t m0 ∂x (20) to be ψ (x, t ) = ∑ exp[−i 2π (t Tβ )n 2 ]cnψ n (x ) (21) n A function F ( n ) whose domain is restricted to the integers ( n ∈ Z ) can be written as a finite sum of exponentials if and only if it is r periodic, that is, if there is an integer r such that F ( n ) = F ( n + r ) for all n . Such a finite sum is called the finite Fourier series. In our case, we identify F ( n ) = exp  −i 2π ( t Tβ ) n 2  . The necessary and sufficient condition for this exponential to be a periodic function of the quantum number n is that the time ratio t Tβ must be rational, and we write T ( p, q ) = p Tβ q (22) for relatively prime integers p and q (that is, p q forms a simplified fraction). In terms of frequency, expression (48) takes the form: f ( p, q ) = 1 q = f0 , T ( p, q ) p f0 = 1 Tβ (23a,b) The condition p > q specifies the presence of sub-harmonics in the Fourier spectra. Now, through the fractal expressions E1 = m0V 2 p = 4πm0 D f ( p, q ) 2 q we can introduce the Reynolds’s “fractional” criterion 178 (24) Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece VL p Re( p, q ) = c~ c = ν q (25) where we used the substitutions: (2 D )−1 Vc = V , Lc = Vf −1 ( p, q ) , ν~ = 8πD (dt ) F (26a-c) From (25) and [12] it results a critical value for the Reynolds number, Re c . Up to this value the fractal fluid flow becomes turbulent. Because from (23a) with p > q results sub-harmonics for Rec , according to [13-15] a criterion of evolution to chaos through a cascade of spatio-temporal sub-harmonic bifurcations is stated. We call this criterion the “fractional” criterion of transition towards chaos. 3. Conclusions A scenario of transition to chaos through cascade of sub-harmonic bifurcations was experimentally evidenced in connection to the generation and dynamics of concentric multiple double layers in plasma. By considering that the particles movements in the gas discharge plasma take place on fractal curves, a mathematical model according with the scale relativity theory was developed in order to describe its dynamics. In this model, the selfstructuring of a discharge plasma as concentric multiple double layers was analyzed, a good similarity between the experimental and theoretical dependences being obtained. By using the fractional revival formalism, a Reynold’s fractional criterion of evolution to chaos through cascade of spatiotemporal sub-harmonic bifurcations was stated. Acknowledgment This work was supported by the Romanian National Authority for Scientific Research, CNCS-UEFISCDI, project number PN-II-ID-PCE-2011-3-0650. References 1. L. Conde and L. Leon. Multiple double layers in a glow discharge. Phys. Plasmas 1:2441-2447, 1994. 2. D. G. Dimitriu, M. Aflori, L. M. Ivan, C. Ionita and R. W. Schrittwieser. Common physical mechanism for concentric and non-concentric multiple double layers in plasma. Plasma Phys. Control. Fusion 49:237–248, 2007. 3. A. J. Lichtenberg and M. A. Lieberman. Regular and Stochastic Motion, SpringerVerlag, New York, 1983. 4. L. Nottale. Fractal Space-Time and Microphysics: Towards a Theory of Scale Relativity, World Scientific, Singapore, 1993. 5. L. Nottale. Fractals and the quantum theory of spacetime. Int. J. Mod. Phys. A 4:5047– 5117, 1989. 179 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 6. R. L. Stenzel, C. Ionita and R. Schrittwieser. Dynamics of fireball. Plasma Source Sci. Technol. 17:035006, 2008. 7. M. Agop, O. Niculescu, A. Timofte, L. Bibire, A. S. Ghenadi, A. Nicuta, C. Nejneru and G. V. Munceleanu. Non-differentiable mechanical model and its implications. Int. J. Theor. Phys. 49:1489-1506, 2010. 8. O. Niculescu, D. G. Dimitriu, V. P. Paun, P. D. Matasaru, D. Scurtu and M. Agop. Experimental and theoretical investigations of a plasma fireball dynamics. Phys. Plasmas 17:042305, 2010. 9. M. Agop, G. V. Munceleanu, O. Niculescu and T. Dandu-Bibire. Static and free timedependent fractal systems through an extended hydrodynamic model of the scale relativity theory. Phys. Scripta 82:015010, 2010. 10. D. Stoler. Equivalence classes of minimum-uncertainty packets. II. Phys. Rev. D 4:1925-1926, 1971. 11. D. L. Aronstein and C. R. Stroud. Fractional wave-function revivals in the infinite square well. Phys. Rev. A 55: 4526-4537, 1997. 12. L. D. Landau and E. M. Lifshitz. Fluid Mechanics 2nd ed., Butterworth-Heinemann, Oxford, 1987. 13. F. T. Arecchi, R. Meucci, G. Puccioni and J. Tredicce. Experimental evidence of subharmonic bifurcations, multistability, and turbulence in a Q-switched gas laser. Phys. Rev. Lett. 49:1217-1220, 1982. 14. M. Dubois, M. A. Rubio and P. Berge. Experimental evidence of intermittencies associated with a subharmonic bifurcation. Phys. Rev. Lett. 51:1446-1449, 1983. 15. A. Atipo, G. Bonhomme and T. Pierre. Ionization waves: from stability to chaos and turbulence. Eur. Phys. J. D 19:79-87, 2002. 180 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Height Field Representation and Compression Using Fractal Interpolation Surfaces Vassileios Drakopoulos1 and Polychronis Manousopoulos2 1 2 National and Kapodistrian University of Athens, Panepistimioupolis, 15784 Athens, Greece (E-mail: [email protected]) National and Kapodistrian University of Athens, Panepistimioupolis, 15784 Athens, Greece (E-mail: [email protected]) Abstract. Height fields provide efficient means for representing surface elevation data which can be used for rendering 3D terrains or landscapes. In this paper, a novel method for representing height fields using fractal interpolation techniques is presented. The proposed methodology allows describing natural surfaces with an intrinsic fractal structure in a more convenient manner. Specifically, fractal interpolation surfaces constructed on rectangular domains have been used. Results indicate the feasibility and advantages of the proposed method in terms of quality of representation as well as compression ratios. Keywords: Fractals, IFS, Interpolation, Surfaces. 1 Introduction Height fields provide an efficient tool for representing surface elevation data and are often used, among other applications, in 3D computer graphics for rendering 3D terrains or landscapes. A height field is essentially a 2D array of height values and is usually stored as a raster image; the pixel intensity corresponds to the height at the location defined by the pixel coordinates. An example of landscape rendering based on a height field is shown in Figure 1. Fig. 1. Rendering of a landscape defined by a height field. 181 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fractal interpolation as defined by Barnsley[1] and other researchers is based on the theory of iterated function systems. It provides an alternative to traditional interpolation techniques, aiming mainly at data which present detail at different scales or possess some degree of self-similarity. These properties define an irregular, non-smooth, structure which is inconvenient to be described by using elementary functions such as polynomials. A fractal interpolation function is a continuous function whose graph is the attractor of an appropriately chosen iterated function system. In case this graph, usually of non-integral dimension, belongs to the three-dimensional space and has Hausdorff - Besicovitch dimension between 2 and 3, the resulting attractor is called fractal interpolation surface. In this paper, a novel methodology for representing height fields using fractal interpolation techniques is introduced. Our motivation stems from the fact that natural surfaces, such as earth terrains, present an intrinsic fractal structure, i.e. detail at multiple scales and some degree of self-similarity. The most important and non-trivial part for constructing fractal interpolation surfaces on rectangular domains involves ensuring their continuity. We also present the application of the proposed methodology to terrain data, indicating its feasibility and advantages in terms of quality of representation as well as compression ratios. The paper is structured as follows. In Section 2 we briefly review the theory of iterated function systems. The necessary concepts of fractal interpolation surfaces, focusing on the proposed construction on rectangular domains, are introduced in Section 3. In Section 4, we present the proposed methodology for height field representation and compression using the surfaces of the previous section. Section 5 contains the result of the application of our method to terrain data, in terms of quality of representation as well as compression ratios. Finally, Section 6 summarizes our conclusions and indicates areas of further work. 2 Iterated function systems Let X, Y ⊂ Rn . A function f : X → Y is called a Hölder function of exponent a if |f (x) − f (y)| ≤ c |x − y|a for x, y ∈ X, a ≥ 0 and for some constant c. Note that, if a > 1, the functions are constants. Obviously, c ≥ 0. The function f is called a Lipschitz function if a may be taken to be equal to 1. A Lipschitz function is a contraction with contractivity factor c, if c < 1. An iterated function system, or IFS for short, is a collection of a complete metric space (X, ρ) together with a finite set of continuous mappings wn : X → X, n = 1, 2, . . . , N , where ρ is a distance between elements of X. It is often convenient to write an IFS formally as {X; w1 , w2 , . . . , wN } or, somewhat more briefly, as {X; w1−N }. 182 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece The associated map of subsets W : H(X) → H(X) is given by W (E) = N  wn (E) for all E ∈ H(X), n=1 where H(X) is the metric space of all non-empty, compact subsets of X with respect to some metric, e.g. the Hausdorff metric. The map W is called the collage map to alert us to the fact that W (E) is formed as a union or ‘collage’ of sets. Sometimes H(X) is referred to as the “space of fractals in X” (but note that not all members of H(X) are fractals). If wn are contractions with corresponding contractivity factors sn for n = 1, 2, . . . , N , the IFS is termed hyperbolic and the map W itself is then a contraction with contractivity factor s = max{s1 , s2 , . . . , sN } (Barnsley[1], Theorem 7.1, p. 81). In what follows we abbreviate by f k the k-fold composition f ◦ f ◦ · · · ◦ f . The graph of f is the set of points G(f ) = {(x, f (x)) : x ∈ X}. The attractor of a hyperbolic IFS is the unique set A∞ for which limk→∞ W k (E0 ) = A∞ for every starting set E0 . The term attractor is chosen to suggest the movement of E0 towards A∞ under successive applications of W . A∞ is also the unique set in H(X) which is not changed by W , so W (A∞ ) = A∞ , and from this important perspective it is often called the invariant set of the IFS. 3 Rectangular subdomain fractal interpolation surfaces Fractal interpolation surfaces constructed as attractors of iterated function systems were first proposed by Peter R. Massopust[4], where he considered the case of a triangular domain with coplanar boundary data. A slightly more general construction of such fractal surfaces was later presented by Jeffrey S. Geronimo and Douglas Hardin[3], where the domain used was a polygonal region with arbitrary interpolation points but same contractivity factors. Here, we focus on fractal interpolation surfaces constructed on rectangular domains with arbitrary boundary data and same contractivity factors. Let D be a closed nondegenerate rectangular region in R2 and let S = {q0 , q1 , . . . , qm−1 } be m, with m > 4, distinct points in D such that {q0 , q1 , q2 , q3 } are the vertices of D. Given real numbers z0 , z1 , . . . , zm−1 we wish to construct a function f such that f (qj ) = zj , j = 0, 1, . . . , m − 1 and whose graph is self-similar. Notice that the constructed FIS in the present work is not self-affine since it is resulting from bivariate functions. Let us denote by C(D) the linear space of all real-valued continuous functions defined on D, i.e. C(D) = {f : D → R | f continuous}. The basic idea is to decompose D into N subrectangles R1 , R2 , . . . , RN with vertices chosen from S and define affine maps Li : D → Ri and contractions Fi : D × R → R, i = 1, 2, . . . , N such that Φ defined by −1 (Φf )(x, y) = Fi (L−1 i (x, y), f (Li (x, y))) 183 (1) Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece maps an appropriate subset of C(D) onto itself. If Li is invertible, G(f ) is mapped onto G (Φ(f )|Ri ) by (Li (x, y), Fi (x, y, f (x)). Henceforth we assume that {Ri }N i=1 consists of nondegenerate rectangles whose interiors are nonN intersecting, L−1 i (Ri ) = D and that the set of vertices of {Ri }i=1 equals S. Let k: {1, 2, . . . , N } × {0, 1, 2, 3} → {0, 1, . . . , m − 1} be such that {qk(i,j) }3j=0 gives the vertices of {Ri }N i=1 . Let i ∈ {1, 2, . . . , N }. Since D and Ri are nondegenerate, there is a unique invertible affine map Li : R2 → R2 satisfying Li (qj ) = qk(i,j) , j = 0, 1, 2, 3. (2) Let si be given such that |si | < 1 and Fi : R3 → R be defined by Fi (x, y, z) = ai x + bi y + gi xy + si z + ci , (3) where ai , bi and ci are uniquely determined by Fi (qj , zj ) = zk(i,j) , j = 0, 1, 2, 3. (4) With these definitions for Li and Fi we have Φ(f )|Ri ∈ C(Ri ) and (Φf )(qk(i,j) ) = zk(i,j) , j = 0, 1, 2, 3, whenever f ∈ C(D) and f (qj ) = zj , j = 0, 1, 2, 3. If Ri and Ri′ are adjacent rectangles with common edge qj qj ′ , it remains to be determined whether Φf is well-defined along qj qj ′ , i.e., whether Φf satisfies the “join-up” condition −1 −1 −1 Fi (L−1 i (x, y), f (Li (x, y))) = Fi′ (Li′ (x, y), f (Li′ (x, y))), for all (x, y) ∈ qj qj ′ . We consider the case where the graph associated with the tesselation {Ri }N i=1 has chromatic number 4. The chromatic number of a graph is the least number of symbols required to label the vertices of the graph so that any two adjacent vertices (i.e., joined by an edge) have distinct m−1 labels. Since each edge is part of some Ri this implies the vertices {qj }j=0 can be labelled with l = l(j) ∈ {0, 1, 2, 3} such that the vertices of each Ri have distinct labels. For i = 1, 2, . . . , N and j = 0, 1, 2, 3 let k(i, j) be determined by the condition k(i, l(j ′ )) = j ′ for all vertices qj ′ of Ri . Then, Eqs. (2) and (4) become Li (ql(j) ) = qj , Fi (ql(j) , zl(j) ) = zj (5) for each of the vertices qj of Ri . Let CB (D) denote the collection of continuous functions f such that f (qj ) = zj , qj ∈ ∂D. Theorem 31 Suppose the graph associated with {Ri }N i=1 has chromatic number 4. Let Li and Fi , i = 1, 2, . . . , N , be determined by (3) and (5) with si = s(|s| < 1). Let Φ be defined by (1). Then Φ: CB (D) → CB (D) is welldefined and contractive in the sup-norm with contractivity s. Furthermore (Φf )(qj ) = zj , j = 0, 1, . . . , m − 1 and f ∈ CB (D). 184 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Proof. See Drakopoulos and Manousopoulos[2] ⊓ ⊔ Then the corresponding IFS is of the form {R3 ; w1−N }, where wi (x, y, z) = (Li (x, y), Fi (x, y, z)). An illustration of this is shown in Fig. 2, where the left part indicates the vertices and connecting edges of D and the middle and right parts of the figure indicate where these vertices are mapped by the domain contractions. For larger data sets, this pattern is repeated as necessary. 3 2 1 0 0 1 0 3 2 3 0 0 1 3 2 3 2 0 1 0 1 3 2 3 2 0 1 0 1 Fig. 2. Rectangular domain contractions to satisfy join-up conditions. 4 Height fields A height field or relief map (see e.g. Theoharis et al.[5], p. 505) is defined as a 2D array of height values: H = {(xi , yj , zij : i = 0, 1, . . . , M and j = 0, 1, . . . , N }, where the x, y coordinates define a rectangular grid on the plane and the z coordinate defines the height. The underlying grid is usually regular, i.e. xi = x0 + i∆x , yj = y0 + j∆y , for every i = 0, 1, . . . , M and j = 0, 1, . . . , N , where ∆x = (xM − x0 )/M and ∆y = (yN − y0 )/N . From the above definition, it is clear that a height field can be directly represented by a fractal interpolation surface of Section 3. The only issue to be determined is whether all of the height field data will be used in the construction of the surface or only a subset of them. This can be achieved by regularly sampling the height field along the x, y dimensions. The sampling frequency defines a trade-off between quality of representation and compression ratio. 185 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece The proposed representation is expected to be especially fruitful for height fields defining natural surfaces, such as terrains. These often possess an intrinsic fractal structure which is conveniently described by fractal interpolation models. The idea of representing natural surfaces using fractal interpolation has also been suggested in Xie et al.[6], where the generation of rock fracture surfaces using fractal interpolation was examined. 5 Results A height field of resolution 257 × 257 is presented in Figures 3(a) and 3(b). Specifically, the former figure depicts the height field as a 2D raster image where brighter areas indicate greater height; the latter figure contains its 3D depiction. This height field, which was created using TerragenTM Classic, contains a total of 257 × 257 = 66049 points. y x (a) (b) Fig. 3. The original height field depicted (a) as a 2D image and (b) as a 3D surface. Figures 4(a) and 4(b) show the 2D and 3D representation of this height field, respectively, using the proposed method. Specifically, it has been represented by a fractal interpolation surface constructed on a subset of the original data with s = 0.02; every 8th point along each dimension of the height field has been chosen as interpolation point. This results in a rectangular grid of resolution 33 × 33, containing 1089 points in total, i.e. about 1.65% of the original points. Despite the significant sparsity of the interpolation points, the quality of the reconstructed height field is satisfactory. Another example is given in Figures 5(a) and 5(b), where the same height field has been represented by a fractal interpolation surface using even fewer interpolation points. Specifically, every 16th point along each dimension has been chosen as interpolation point. This results in a rectangular grid of resolution 17 × 17, containing 289 points in total, i.e. about 0.44% of the original points. Also in this case, the results are satisfactory despite the even smaller number of interpolation points. These results indicate that fractal interpolation surfaces are indeed capable of describing natural surfaces, such 186 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece y x (a) (b) Fig. 4. The reconstructed height field, using every 8th data point as interpolation point, depicted (a) as a 2D image and (b) as a 3D surface. as terrains, with considerable quality even when high compression ratios are involved. y x (a) (b) Fig. 5. The reconstructed height field, using every 16th data point as interpolation point, depicted (a) as a 2D image and (b) as a 3D surface. Figures 6–7 depict an artistic rendering of the original height field as well as its two aforementioned reconstructions; these figures were created using TerragenTM Classic. As shown in the figures, the reconstructed height fields produce equivalent results compared to the original one, even though the significant sparsity of the interpolation points. 6 Conclusions and future work We have presented a novel method for the representation and compression of height fields using fractal interpolation techniques. Specifically, we have represented a height field as a fractal interpolation surface constructed on the rectangular domain defined by a subset of the original data. The results indicate that the proposed methodology is feasible, while achieving satisfactory results in terms of quality of representation as well as compression ratios. 187 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig. 6. Artistic rendering of the original height field. (a) (b) Fig. 7. Artistic rendering of the reconstructed height field, using as interpolation point (a) every 8th data point and (b) every 16th data point. Further work will focus on the calculation of the optimal values of the vertical scaling factors in order to achieve increased localized accuracy, as well as on the exploration of alternative fractal interpolation surface models, affine or bivariate, including recurrent fractal interpolation surfaces. References 1.M. F. Barnsley. Fractals everywhere, 2nd ed., San Diego, 1993. Academic Press Professional. 2.V. Drakopoulos and P. Manousopoulos. Bivariate fractal interpolation surfaces: Theory and applications. International Journal of Bifurcation and Chaos, to appear, 2012. 3.J. S. Geronimo and D. Hardin. Fractal interpolation surfaces and a related 2-D multiresolution analysis. Journal of Mathematical Analysis and Applications, 176:561-586, 1993. 4.P. R. Massopust. Fractal surfaces. Journal of Mathematical Analysis and Applications, 151:275-290, 1990. 5.T. Theoharis, G. Papaioannou, N. Platis, and N. M. Patrikalakis. Graphics & Visualization: Principles and Algorithms, Massachusetts, 2008. AK Peters. 6.H. Xie, H. Sun, Y. Ju and Z. Feng. Study on generation of rock fracture surfaces by using fractal interpolation. International Journal of Solids and Structures, 38:5765-5787, 2001. 188 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece The Baltic Sea coupled ice-ocean model Lidia Dzierzbicka, Jaromir Jakacki, Maciej Janecki, Artur Nowicki Institute of Oceanology Polish Academy of Sciences, Sopot, Poland E-mail: [email protected] Abstract: Community Earth System Model (CESM) is a coupled climate model that consist of five separate components with additional central coupler which controls time, exchanging data, domains, grids and other model data. Here CESM was adopted for the Baltic Sea and called 3D-CEMBS. This is not fully coupled configuration. In our case we have been taken into account ocean (POP model, version 2.1) and ice (CICE model, varsion 4.0) models which are forced by atmospheric data model (datm7). Other models are turned off in our configuration. However the main task of the datm7 is to interpolate data on model domain, we have interpolated data outside of the models. As an external forces daily-averaged forty-years reanalysis data derived from European Centre for Medium-range Weather Forecasts (ECMWF, ERA-40 reanalisys) has beed used. Currently in pre-operational state 3D-CEMBS is using 48-hour atmospheric forcing data from ICM-UM model (University of Warsaw). 3D-CEMBS model has also ecosystem part (currently work in progress). The study was financially supported by the Polish State Committee of Scientific Research (grants: No N N305 111636, N N306 353239). Keywords: Baltic Sea, 3D model, hydrodynamic variables 1. Introduction In 2011, the operational ecohydrodynamic model (3D CEMBSv2 - new version) was launched at the Institute of Oceanology PAS in the parallel version on the 2 km grid with rivers and the open boundary for the hydrodynamic module http://deep.iopan.gda.pl/CEMBaltic/new_lay/index.php). Therefore, the following works were completed: - a new version of the model was launched (3D CEMBSv2), which includes the ocean-ice module POP-CICE (POP, version 2.1, the ice module CICE, version 4.0) and the ecosystem module (Fig.1): • compilation of the 2 km bathymetric grid, • adaptation of the model to meteorological data UM (ICM UW), • compilation of a data set with meteorological forecasts (UM) from 2010 and 2011 in the required format, preparation of the initial fields for the Baltic area in the 2 km grid, - development of the operational system for data retrieval necessary for the model 3D CEMBSv2, - development of the operational ecohydrodynamic model (3D CEMBSv2) in the parallel version on the 2 km grid with rivers and the open boundary: • compilation of climatic river discharge , • compilation of climatic conditions along the open boundary. This paper presents an integrated, operational model of the Baltic ecosystem – the hydrodynamic part, the ocean-ice module POPCICE (version 2.1 and 4.0) 189 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece and the initial verification, the comparison of the obtained results with data from other models and measurements in situ. Fig. 1. Schematic presentation of the Baltic model 2. 3D-CEMBSv2 MODEL CCSM4.0/CESM1.0 (the Community Climate System Model/Community Earth System Model) consists of five separate elements with an additional coupler (CPL7), which monitors the time, exciting forces, domains, grids and exchange of information between models. The central part of the model is based on MCT (The Model Coupling Toolkit) – which is a parallel tool providing a number of services, such as a register of particular components of the model, descriptors of the domain distribution into processes, communication, redistribution of data and other very helpful tools. For our purposes, CESM was adapted for the Baltic Sea and was called 3D-CEMBS. However, it is not an entirely active configuration. The ocean model (POP, version 2.1) and the ice model (CICE, version 4.0) work in the active mode, and they are imposed by the model of atmospheric data (datm7). Other models are excluded from this configuration (the stub mode). The main task of datm7 is interpolation of atmospheric data into the domain of the model. External forces are daily-averaged for the period of forty years, which come from the ECMWF re-analysis (ERA-40). At present, in the operational mode, 48-hour atmospheric forecasts are used, which are supplied by the UM model of the Interdisciplinary Centre for Mathematical and Computational Modelling of the Warsaw University. 3D-CEMBS Model is also equipped with the ecosystem module, on which works are being currently carried out in order to incorporate it into the operational mode. POP The ocean model is based on the Parallel Ocean Program (POP, [7]) from the National Laboratory in Los Alamos (LANL), which is derived from the global ocean model ([6]) with additional conditions for free surface. This is a model of ‘z’ type (identical thickness of layers for every cell); the three-dimensional 190 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece equations describing the behaviour of the stratified ocean are solved by parametrization. Numerically the model defines spatial derivatives in the spherical coordinates using the method of finite elements. Physical quantities of the model are embedded in the spherical grid of Arakaw B ([1]). The barotropic equation is solved by ’preconditioned conjugate gradient solver’ (PCG), and advection is determined by a centred differential. Parametrization of horizontal mixing is accomplished by a biharmonic operator, and vertical turbulence is determined by KPP parametrization. The equation of state, introduced by McDougall et al. ([5]), is also used. POP is a three-dimensional hydrodynamic model derived from the ocean model created in the late 1960s by Kirk Bryan ([3]) and Michael Cox from the NOAA Geophysical Fluid Dynamics Laboratory in Princeton. The model was later modified and adapted by Semtner ([6]) for vector processors. The whole class of models from which POP is derived are of Brayan-Cox-Semtner type (B-C-S). The code of the model is adapted inter alia for supercomputers, but is well adaptable (porting) also for machines of different architecture, such as cluster types. The code of the POP model is characterised by good numerical performance and is well scalable on a large number of O(1000) processors. A special feature of the POP code (written in Fortranie 90) is orientation to parallel machines and the use of MPI (Message Passing Interface) or SHMEM (Shared Memory) libraries to communicate between processors. And therefore, the most technologically advanced machines can be used to perform calculations and to solve major computational problems. CICE CICE (Community Ice CodE) is based on elastic-visco-plastic (EVP) rheology ([4]). It is designed to work in accordance with the POP ocean model using the parallel computing machines. It consists of several interactive elements: the thermodynamic model, which computes local growth rates of snow and ice owing to vertical conduction of energy and momentum fluxes. It also defines velocity of each ice cell based on wind and ocean velocity. It has a few vertical categories, so that the stress distribution is much closer to the real one. CICE was applied worldwide in different configurations. It was used to study the variability and the impact of ice on the processes occurring in the atmosphere and the ocean in time scales ranging from decades to hundreds of years ([2]). It is used in regional models and global applications. Incorporation of the CICE model (version 4.0) into the regional POP model (version 2.1), completed within the scope of the research grant (No NN305111636 − the Polish State Committee for Scientific Research), is the first application of such an advanced ice model for the Baltic Sea. Ecosystem The ecosystem model consists of 11 main components: zooplankton, small phytoplankton, diatoms, cyanobacteria, two types of detritus deposits (fraction of dissolved DOM and molecular POM detritus), dissolved oxygen, and components of nutrients, such as: nitrates, ammonia, phosphates and silicates. The class of small phytoplankton should reflect nano- and picophytoplankton, and can be limited by nitrates and phosphates. The class of larger phytoplankton 191 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece is represented by diatoms and can be limited by the aforementioned factors, as well as by silicates. The growth rate of cyanobacteria may depend on phosphates and the available light. Many biotic and detrital compartments contain multiple elemental pools as we track carbon, nitrogen, phosphorus and silicon through the ecosystem. 3D-CEMBS model is currently configured at two horizontal resolutions of 9 km and approximately 2km (1/12° and 1/48° respectively). The model bathymetry is represented by 21 vertical levels and the thickness of the first four surface layers was five metres. The bottom topography is based on ETOPO1 1 arcminute global relief model (http://www.ngdc.noaa.gov/mgg/global/global.html, National Geophysical Data Center). The bathymetric data were interpolated to the model grid using the kriging method. The initial state of the ocean model was prepared using temperature and salinity climatic data. The ocean surface level (5m deep) is restored based on the monthly timescale to the monthly average T and S climatology, as a correction term to the explicitly calculated fluxes and overlying atmosphere or sea ice. The restoring time was set to 30 days at the surface and 10 days at the domain boundary. 3D-CEMBS model domain is based on stereographic coordinates, but the equator of these coordinates is in the centre of the Baltic Sea (so we actually use rotated stereographic coordinates) and we can assume that shape of the cells is square and they are identical. The current calculations are performed on supercomputers of cluster type, Galera, which is located at TASK (the Academic Computer Centre in Gdańsk). The time needed for computing 1 model year, for the ecohydrodynamic model, for resolution of 9km, is 30 h on 16 processors, and for resolution 2km – 120 h on 256 processors. 3. Results In the first half of 2011, the new CESM model (Community Earth System Model - UCAR/NCAR USA) was adapted for the Baltic domain. Models have been adapted and work properly (see the website). The results of simulations for 48-hour forecast are presented for the model with resolution of 9 km and 2 km. Select 2km forecast for the area or the point (Fig. 2). 192 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig. 2. Select 2km forecast for the area. Simulations for 48-hour forecast for 2 km: Forecast for the area: • select the time of forecast (48-hour forecast) (Fig. 3a), • select a required variable (temperature, salinity, currents, sea surface height, ice cover area) (Fig. 3b), • select the depth, one of the ten layers for which you would like to see the model results (Fig. 3c), • optionally you may provide boundaries to narrow the area of your query by selecting two random points (corners of a rectangle) on the map (Fig. 3d), • to change coordinates x and y of selected points, press: Reset Coordinates (Fig. 3e), • after all parameters are selected, press: Submit (Fig. 3f). 193 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig. 3. Forecast for the area. Choose parameters of forecast for 05.03.2012. The example is presented below for the following situations: 1. Forecast start: 2012.03.05, 18:00 UTC (Fig. 4) Select forecast parameters: hour: +6 (forecast is for 2012.03.06, 00:00 UTC), variable: currents (cm s-1); depth: 20-26 m (5th layer), coordinate X: 187 to 425, coordinate Y: 32 to 168. After selection of parameters, press: Submit and the screen will show a drawing with results of the model. 194 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig. 4. Forecast for the area. Forecast for the point : • select the depth, • select a point on the map to get 48 h time series, • to change coordinates x and y of selected points, press: Reset Coordinates, • after selection of all parameters, press: Submit The example is presented for the following situation: forecast start: 2012.03.05, 12:00 UTC (Fig. 5). 4. Conclusion The work presents the idea of a three-dimensional Coupled Ecosystem Model Baltic Sea version 2 (3D CEMBSv2) – a hydrodynamic part, which determines the main physical parameters of the environment: horizontal components of the velocity u,v and the vertical component of velocity w, pressure p, density ρ , temperature T and salinity S of water. Variables presented on the website for 48hour forecast are as follows: temperature, salinity, currents, sea surface height and ice area cover. The 3D CEMBSv2 model (at present – the hydrodynamic module) is a suitable tool for studying the annual, seasonal, monthly and daily variability of environmental parameters in the southern Baltic Sea. It can therefore be applied in the forecasting of ecological changes in the Baltic. 195 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig. 5. Forecast for the point.. References 1. A. Arakawa and V.R. Lamb. Computational design of the basic dynamic processes of the UCLA general circulation model. Methods Comput Phys 17: 173–265, 1977. 2. C.M. Bitz and W.H.Lipscomb. An energy-conserving thermodynamic sea ice model for climate study. Journal of Geophysical Research, Ocean 104: 15669–15677, 1999. 3. K.A. Bryan. Numerical Method for the Study of the Circulation of the World Ocean. J Comput Phys 4(3): 47–376, 1969. 4. E.C. Hunke and J.K. Dukowicz. An elastic-viscous-plastic model for sea ice dynamics. J Phys Oceanogr 27:1849–1867, 1997. 5. T.J. McDougall, D.R. Jackett, D.G. Wright and R. Feistel. Accurate and computationally efficient algorithms for potential temperature and density of seawater. J Atm Ocean Technol 20: 730–741, 2003. 6.. A.J. Semtner. A general circulation model for the World Ocean. UCLA Dept. of Meteorology Tech.Rep No.8: 99 pp., 1974. 7. R. Smith and P. Gent. Reference manual for the Parallel Ocean Program (POP), Los Alamos National Lab., New Mexico, 75 pp., 2004. 196 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Extraction of the Dominant Features of Complex Dynamics in Experimental Air-Water Two-Phase Flows Alberto Fichera, Arturo Pagano University of Catania, Department of Industrial and Mechanical Engineering, Viale Andrea Doria, 6, I-95125 Catania, Italy E-mail: [email protected] Abstract: This paper aims at describing a novel approach to the analysis of the experimental void fraction time series detected from air-water upward two phase flows in a vertical channel. This system can express a great variety of different flow patterns, whose characterisation and classification strongly depends on the approach used for feature extraction. Phase space analysis in a traditional delayed embedding has allowed for the observation of the complex dynamics of the system. Nonetheless, the attractors obtained in a delayed embedding, though characterised by a regular complex structure, appear partly folded and are affected by noisy hydrodynamic high order dynamics. The present paper proposes the application of the SVD approach with the aim of assessing a more appropriate embedding into the phase space spanned by the principal vectors. In this way the dominant features of the system dynamics can be separated from noise-like dynamics, i.e. unfolded and noise-free attractors can be obtained. Keywords: Feature extraction, Two-phase flows, Experimental nonlinear dynamics, SVD analysis. 1. Introduction Several basic industrial processes, ranging from power generation, chemical and processing plants to oil pipelines, present heat and mass transfer applications of two phase flows. When two phase flows occur, very different flow patterns can be observed as well as transitions from a flow pattern to another. Indeed, the dynamical behaviours associated to the various types of flow pattern established in the system represent critical factors for the performances of such industrial systems. This explains the great efforts that have been and are still devoted to flow patterns identification, which represents a fundamental basis for appropriate characterisation of two phase flow systems. The dynamics of two phase flows are typically of highly complex pulsating nature, under the effect of several nonlinearities deriving from the strong coupling of different mechanisms and of the dependence on various factors. Among the others, the most important factors are the differential action of gravity on the two phases and the effect of shear and surface tension forces at their interface. As a consequence, several different flow patterns can be identified, each of which can be characterised in terms of the dynamical behaviour of the void fraction time series. 197 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Among the other, two phase flows of air-water mixtures are often theoretical and experimental analysed with the aim of achieving a reference perspective on the general dynamical behaviours, often valid also for more complex flows, such as those arising in presence of phase changes. In particular, the present study aims at analysing the behaviour of ascending air-water two phase flows in vertical pipes. For this kind of flows heat transfer phenomena connected to phase change are not involved, so that the flow pattern established in the system mainly depends on the mass flow rates of the two phases. By varying the mass flow rate of the two phases, in fact, bubbly, slug, churn and annular flows can be identified as the main flow patterns typical of several classifications [1-3]. The bubbly flow exists for low values of the gas mass flow rate and consists in the motion of dispersed and small gas bubbles in the liquid phase. Coalescence phenomena are at the basis of the transition from bubbly to slug flow, which can be observed by increasing the gas mass flow rate. Slug flow is characterised by gas bubbles, namely Taylor bubbles, enveloped by a liquid film separating them from the pipe walls, alternated to liquid slugs. In the class of slug flow, it is possible to distinguish between: cap flow, with short air bubbles (with the head approximately connected to the tail) separated by long liquid slugs; plug flow, with gas bubbles and liquid slugs of comparable length; proper slug flow, characterised by elongated gas bubbles separated from relatively short liquid slugs, often aerated for the presence of small dispersed air bubbles. For growing gas mass flow rate, bubble coalescence and increasing aeration of the liquid slug leads to a highly unstable flow pattern addressed as churn flow, characterised by waves propagating through the liquid film enveloping the bubbles and occasionally falling within the tube, so to form a short, unstable and highly aerated liquid slug. Finally, the annular flow consists of a thin annular liquid film at the tube wall on which small ripples, interspersed occasionally with large disturbance waves, flow in a regular manner up the tube. It is usual practice to perform flow pattern identification on the basis of the differences of the dynamical behaviour of the time series of the local void fraction. Therefore, the reliability of the identification approach is highly dependent on the accuracy of the technique adopted to measure the void fraction. Several techniques have been proposed [4-9] and impedance measurements seem to be recognized as the most reliable [6]. At the same time, the performances of flow patter identification approaches depend also on the techniques adopted for time series analysis and feature extraction. Statistical [1, 2, 6] or spectral [9-12] techniques indeed represent the typical approach for flow patterns identification on the basis of the analysis of the experimental void fraction time series. Nonlinear techniques have been also adopted, among the others see [10, 13-16], but a main drawback has been represented by the relatively poor spatial and temporal resolution of the experimental time series. In order to address this problem, the experimental time series considered in the present study have been detected by means of a resistive probe characterised by high temporal and spatial resolution, which has been appositely set-up as described in [17]. The preliminary analysis in a delayed embedding of the void 198 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece fraction time series detected by means of this sensor has shown the existence of strange attractors of interesting morphology for the various flow patterns [18]. Nonetheless, attractors obtained in this way are somewhat noisy as a consequence of the superposition of high order dynamics to the dominant dynamics characterizing the flow pattern. Among the others, the most important high order “noisy” dynamics are those of hydrodynamic nature associated to small diameter bubbles dispersed in the liquid slugs and to disturbances on the liquid film enveloping the Taylor bubbles. Therefore, the present study aims at extracting the dominant features of the flow dynamics under various flow pattern conditions so to separate the dominant features of the system dynamics from noise-like dynamics. The proposed approach is analogous to that proposed in [19] and is based on the calculation of the singular vectors of a n-dimensional delayed embedding, through the application of the technique known as Singular Value Decomposition (SVD) [20], and in the analysis of the restricted portion of the dynamics that is obtained by projecting the attractor onto the phase space spanned by the singular vectors corresponding to the three highest singular values. Reported results, show that the attractors described in the new embedding present a well defined and regular structure, indicating the existence of a low order source of the system dynamics, which will be analysed in future studies. Fig.2 Experimental apparatus 2. Experimental Apparatus The experimental apparatus reported in Fig. 1 has been built and tested in order to study the dynamics of two-phase flow in vertical pipes. The test section is a vertical pipe of diameter 0.026 m diameter and length 3 m. The apparatus is equipped by an electromagnetic flowmeter and three air flow metres, 199 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece respectively used for the measure of the water velocity and mass flow rate and for the regulation of the air flow rate in the range between 10 and 210 l/min. The air is supplied to the mixing section by a pressurised tank fed by a compressor, whereas the water flow rate can be varied in the range 0-150 l/min by means of a series of valves and bypasses placed at the pump outlet. A resistive probe for the measure of the void fraction is placed at a distance of over 100 times the diameter of the pipe from the mixing section, i.e. over the required entry region for two phase flows, in order to ensure a well established flow regime. In particular, the void fraction probe has been designed and realised for the experimental campaign and operates in the resistive range (carrier frequency of 20 kHz). The sampling frequency was set at 1 kHz with a cut-off frequency of 200 Hz. A detailed description of the experimental probe and on the wide set of experimental tests performed is reported in [17]. 3. Dynamical Feature Extraction The results of preliminary linear analyses of the experimental time series have been shown to be unable to deal with the intrinsic complexity of two phase flows dynamics. Hence, in [18] a morphological analysis of the threedimensional attractors has been proposed in a classical Takens’ delayed embedding of the experimental void fraction time series [21]. In particular, it has been observed that the attractors obtained for some of the flow patterns are characterised by a regular fractal structure, which is indeed one of the most important evidences of deterministic chaotic behaviour. In the present study, the aim is to improve the dynamical representation by adopting a new embedding, derived through the application of Singular Value Decomposition technique, SVD [20], to the classical delayed embedding based on Takens’ theorem, similarly to the approach proposed in [19]. The new representation is characterised by a drastic reduction of noisy dynamics and, above all, a sensitive improvement of the attractor unfolding, so that the dominant morphological characteristic can be fully exploited. As a first step, the phase space reconstruction consists in the creation of a n × w matrix, S, where n is the length of a window moving through the data and w are the independent variables defining the phase space, i.e. delayed version of the experimental void fraction time series s(t)=(s0, s1, s2, …, si, …), with each column delayed τ time steps from the previous. The condition w>2d+1 for an appropriate embedding is implicitly respected if w is set much greater than the unknown fractal dimension d on the basis of a preliminary estimation. The second step consists in the application of the SVD approach to matrix S. This is done through the calculation of a new diagonal matrix, equivalent to the original one, i.e. with identical singular values but in decreasing order. In particular, S is factorized into its singular values according to equation: Λ = MT S C (1) In (1) Λ is the diagonal matrix containing the singular values of S in decreasing order and M and C are the matrices of the singular vectors associated with Λ. 200 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Details on the factorization can be found in [20]; what is interesting for the scopes of the present study is that the high level singular values in Λ are associated to the dominant singular vectors, i.e. those representing the dominant features of the system dynamics, whereas the low level ones correspond to local behaviors or noise-like components. Therefore, the system can be virtually partitioned into two subsystem: the first deriving from noise free data (i.e. the main features and the relevant details) and the second from noisy dynamical behaviours, which can be considered superimposed and then eliminated. 4. Results and Discussion The described approach has been used in the present study in order to obtain a denoised and unfolded representation of the experimental dynamics. The SVD technique has been applied to the delayed embedding S of the experimental void fraction time series, created considering τ=1 and w=40 in order to ensure that w is sufficiently greater than m, i.e. greater than the (unknown) system dimension. The length n of the observation window has been set at 10000 data samples in order to be wide enough to obtain a well defined attractor in phase space, i.e. an attractor whose morphology does not change if further data samples are added. The claimed advantages of the proposed methodology can be observed in the results reported in the Fig. 3 to 8, which report the attractors of the same operating condition in two different embeddings. In particular, the phase space adopted for the plots on the left hand side of each figure is the basic three dimensional Takens’ delayed embedding, whereas the projections on the pseudo-phase space spanned by the three dominant principal vectors of the improved embedding obtained through application of SVD are those reported on the right hand side of each figure. By comparing the two methods of representation it is possible to observe that the attractors in the delayed phase space are in all cases sensibly affected by a higher noise level and are not sufficiently unfolded with respect to the corresponding attractors in the principal component embeddings, the last being characterised by a very low level of noise and a satisfactory unfolding. It is worth to remind that, even if the two attractors of each flow pattern appear different, they are, nonetheless, expressions of the same dynamical behaviour. In fact, they are morphologically equivalent and, therefore, characterised by the same invariants of the dynamics, such as fractal dimension and Lyapunov exponents [22-25]. The successful unfolding contributes to the achievement of a clear and well defined morphology of the attractors. This is a main advantage for the distinction of different flow patterns through a comparison of the representation of their dynamics in the phase space spanned by the principal components. Moreover, in some cases the proposed embedding amplifies important characteristics of the system dynamics. For example, the right hand cap flow attractor in Fig. 4 shows a clear distribution of the trajectory in alternated bands, which is a hint of the fractal (i.e. chaotic) nature of the system dynamics. Finally, the representation in the principal component phase space is very effective in underlining the differences between the various flow patterns. 201 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig.3 Attractors in the delayed and principal component embeddings for the bubbly flow; air flow rate 2 lit/min - water flow rate 32.4 lit/min. Fig.4 Attractors in the delayed and principal component embeddings for the cap flow; air flow rate 5 lit/min - water flow rate 20.28 lit/min. Fig.5 Attractors in the delayed and principal component embeddings for the plug flow; air flow rate 10 lit/min - water flow rate 9.06 lit/min. Each type of flow pattern is, in fact, characterised by a specific morphology, sufficiently different from that of the other flow patterns. In particular: - each flow pattern attractor occupies a different phase space region; - each attractor differently “fills” its own region of phase space; for example, the cap flow attractor (properly 3-D) has a higher filling rate than that of the plug flow (which moves around a sort of 2-D limit cycle); 202 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece - the attracting region is progressively shifted, with a continuous trend from bubbly to annular flow, with respect to the first principal component. These differences are very important as the morphological considerations drawn insofar are related to the fractal nature and to the stretch and folding behaviour of the attractors [24], which can be considered as the topological expressions of the mentioned invariants of the dynamics, whose calculation is behind the scope of the present study and will be the object of future studies. 5. Conclusions This study proposes a phase space approach for the description of typical complex dynamics of two-phase flow. At first the singular vectors of the classical delayed embedding are calculated and the attractors of the system dynamics are projected on the state space spanned by these eigenvectors. In this way the dominant feature of the dynamics, corresponding to a subset of the highest singular values, are separated from noisy dynamics in the time series, corresponding to the remaining lower singular values. The morphology of the attractors in the obtained unfolded and noise-free representation are analysed. Reported results demonstrate that the proposed approach represent a powerful tool for the identification of two-phase flow patterns. Fig.6 Attractors in the delayed and principal component embeddings for the slug flow; air flow rate 40 lit/min - water flow rate 16.80 lit/min. Fig.7 Attractors in the delayed and principal component embeddings for the churn flow; air flow rate 80 lit/min - water flow rate 9.01 lit/min. 203 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece Fig.8 Attractors in the delayed and principal component embeddings for the annular flow; air flow rate 80 lit/min - water flow rate 5.58 lit/min. References 1. G. Costigan, P. B. Whalley, Slug flow regime identification from dynamic void fraction measurement in vertical air-water flows, International Journal of Multiphase Flow 23, pp 263–282, 1997. 2. Y. Mi, M. Ishii, L. H. Tsoukalas, Vertical two-phase flow identification using advanced instrumentation and neural networks, Nuclear Eng. Design 184, pp 409420, 1998. 3. F. S. Qi, G. H. Yeoh, S.C.P. Cheung, J. Y. Tu, E. Krepper, D. Lucas, Classification of bubbles in vertical gas–liquid flow: Part 1 – An analysis of experimental data, International Journal of Multiphase Flow 39, pp 121–134, 2012. 4. F. Devia, M. Fossa, Design and optimisation of impedance probes for void fraction measurements, Flow Measurement and Instrumentation 14, pp 139-149, 2003. 5. D. Lowe, K. S. Rezkallah, A capacitance sensor for the characterization of microgravity two-phase liquid-gas flows, Measurement Science and Technology 10, pp 965-975, 1999. 6. J. K. Keska, B.E. Williams, Experimental comparison of flow pattern detection techniques for air-water mixture flow, Exp. Thermal Fluid Science 19, pp 1-12, 1999. 7. P. Andreussi, A. Di Donfrancesco, M. Messia, An impedance method for the measurement of liquid hold-up in two phase flow, International Journal of Multiphase Flow 14, pp 777–785, 1988. 8. C. Vial, E. Camarasa, S. Poncin, G. Wild, N. Midoux, J. Bouillard, Study of hydrodynamic behaviour in bubble columns and external loop airlift reactors through analysis of pressure fluctuations, Chemical Eng. Science 55, pp 2957-2973, 2000. 9. C. H. Song, M. K. Chung, H. C. No, Measurements of void fraction by an improved multi-channel conductance void meter, Nuclear Eng. Design 184, pp 269-285, 1998. 10. N. D. Jin, X. B. Nie, Y. Y. Ren, X. B. Liu, Characterization of oil/water two-phase flow patterns based on nonlinear time series analysis, Flow Measurement and Instrumentation 14, pp 169-175, 2003. 11. G. P. Lucas, I. C. Walton, Flow rate measurement by kinematic wave detection in vertically upward, bubbly two-phase flows, Flow Measurement and Instrumentation 8(3-4), pp 133-143, 1997. 12. M. J. Watson, G. F. Hewitt, Pressure effects on the slug to churn transition, International Journal of Multiphase Flow 25, pp 1225–1241, 1999. 13. J. Drahos, J. Zahradnik, M. Puncochar, M. Fialova, K. Chen, F. Bradka, Effect of operating conditions of the characteristics of pressure fluctuations in a bubble column, Chemical Engineering Processing 29, pp 107-105, 1991. 204 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 14. H. Letzel, J. Schouten, R. Krishna, C. M. van den Bleek, Characterization of regimes and regime transitions in bubble columns by chaos analysis of pressure signals, Chemical Engineering Science 52, pp 4447-4459, 1997. 15. J. Drahos, J. Tihon, C. Serio, A. Lubbert, Deterministic chaos analysis of pressure fluctuations in a horizontal pipe at intermittent flow regime. Chemical Engineering Journal 64, pp 149–156, 1996. 16. J. Y. Lee, N. S. Kim, M. Ishii, Flow regime identification using chaotic characteristics of two-phase flow, Nuclear Eng. Design 238, pp 945–957, 2008. 17. L Cantelli, A. Fichera, A. Pagano, A high-resolution resistive probe for nonlinear analysis of two-phase flows, Journal of Thermodynamics 2011, Article ID 491350, p. 10, 2011. 18. L. Cantelli, A. Fichera, I. D. Guglielmino, A. Pagano, Non-Linear Dynamics of AirWater Mixtures in Vertical Pipes: Experimental Trends, International Journal of Bifurcation and Chaos 16-9, 2006. 19. G. Cammarata, A. Fichera, A. Pagano, Nonlinear analysis of a rectangular natural circulation loop, Int. Comm. Heat Mass Transfer 27 (8), pp 1077-1089, 2000. 20. D.S. Broomhead, G. P. King, Extracting Qualitative Dynamics from Experimental Data, Physica D 20, 217-236 1986. 21. F. Takens, Lecture Notes in Mathematics, Dynamical System and Turbulence, D. A. Rand & L.S. Young Editions, Springer, New York, 1981. 22. J. M. T. Thompson, H. B. Stewart, Nonlinear Dynamics and Chaos, J. Wiley & Sons, New York, 1986. 23. S. N. Rasband, Chaotic Dynamics of Nonlinear Systems, J. Wiley & Sons, New York, 1990. 24. S. H. Strogatz, Nonlinear Dynamics and Chaos, with Applications to Physics, Biology, Chemistry, and Engineering, Perseus Books, Cambridge, 1994. 25. H. Kantz, T., Schreiber, Nonlinear time series analysis, Cambridge University Press, 1999. 205 Proceedings, 5th Chaotic Modeling and Simulation International Conference, 12 – 15 June 2012, Athens Greece 206