Provided for non-commercial research and educational use. Not for reproduction, distribution or commercial use. The attached copy is furnished for non-commercial research and education use only. Authors are permitted to post this version of the article to their personal websites or institutional repositories and to share with other researchers in the form of electronic reprints. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to third party websites are prohibited. For further information on Pliska Studia Mathematica visit the website of the journal http://www.math.bas.bg/~pliska/ or contact: Editorial Office Pliska Studia Mathematica Institute of Mathematics and Informatics Bulgarian Academy of Sciences Telephone: (+359-2)9792818, FAX:(+359-2)971-36-49 e-mail: pliska@math.bas.bg Pliska Stud. Math. 30 (2019), 171–184 STUDIA MATHEMATICA CNN MODELING OF A CLASS OF INTEGRO-DIFFERENTIAL EQUATIONS Angela Slavova, Zoya Zafirova, Pietro Zecca In this paper we study a class of integro-differential equations. First we present Cellular Nonlinear Networks (CNN) modeling of integro-differential equations. We map the equation under consideration into CNN architecture and study its dynamics. Applying harmonic balance method we obtain periodic solutions of the model. We provide numerical simulations which illustrate the obtained theoretical results. 1. Introduction Recently the most important microprocessor manufactures realized that one of the main challenges for the near future is to build efficient processors and infrastructure for the real time handling of images videos or for general time signals coming from space distributed sources. Because both of these tasks are strictly related to spatio-temporal computing, a great effort is then performed to devise supercomputers able to perform spatio-temporal calculations of integrodifferential equations in real time. From this perspective the possibility to exploit the capabilities of analog computation on signal flows instead of traditional digital computation on bits arises. Integro-differential type of equations are widely used to describe phenomena in different fields, as biology, neuroscience, population dynamics, propellant 2010 Mathematics Subject Classification: 92B20, 35B10, 35Q68, 35Q92. Key words: integro-differential equation, Cellular Nonlinear Networks (CNN), dynamic behaviour, harmonic balance method. 172 A. Slavova, Z. Zafirova, P. Zecca rocket motors and in networked control systems. In this paper we consider a class of integro-differential equations (IDE) described by: (1) d2 u = −b1 dt2 Z t −a1 (t−s) e 0 u(s)ds − b2 Z t e−a2 (t−s) u(s)ds 0 where a1 , a2 , b1 , b2 are constants. In order to explain a possible scheme of appearing of integro-differential equations in mathematical models let us consider a simplified model for motion of a single mass point: y ′′ = f (t) + ∆f (t), t ∈ [0, +∞), where vector y(t) describes the motion, and ∆f (t) is an unknown error in defining the right-hand side f (t) which could be estimated ||∆f (t)|| ≤ δ. We can consider m Z t X the error function to be in the form ∆f (t) = − Kr (t, s)y(s)ds, where the r=1 0 −αr (t−s) kernels Kr (t, s) = e βr . This form of the kernels plays an important role in the applications of integro-differential equations. Without loss of generality we can reduce (1) into a system of first order differential equations: x′1 = −x3 − x4 x′2 = x1 (2) x′3 = b1 x2 − a1 x3 x′4 = b2 x2 − a2 x4 , ′ Z where x1 = u (t), x2 = u(t), x3 = b1 t −a1 (t−s) e Z u(s)ds, x4 = b2 0 t e−a2 (t−s) u(s)ds. 0 System (2) is equivalent to (1) and we shall use it later for our numerical simulations. In Section 2 we shall present Cellular Nonlinear Networks (CNN) and how we can model partial differential equations by this approach. We consider more general form of integro-differential equation (1) and obtain its CNN model. In Section 3 we apply harmonic balance technique for studying the dynamics of our CNN model. Computer simulations are provided as well. 2. Cellular Nonlinear Networks modeling Meanwhile, the electronic industry has developed to a stage where billion transistors, nanosecond time scale, nanometer scale size and millions of pixels, etc. are CNN modeling 173 becoming possible, in laboratories and manufacturing alike. A few characteristic directions are as follows. – The essentially 2D structure of silicon and optical technology, as well as the latest sub-100 nm feature size with its power dissipation and wire length constraints push the architectures to an essentially parallel processor, 2 D cellular, sparsely connected nature. – The “sensory revolution”, a third wave in electronics industry, and its impact on data provides for thousands and millions of analog sensory signals organized in arrays to be processed and embedded into sensory computing devices. Typical data are becoming image flows of various kinds: multi-spectral visual, auditory scene, tactile and somatosensory, etc. Multi-modal sensor fusion is becoming a natural task. – Bionic array interfaces lead us to be interconnected to living organisms with complexities far exceeding our man-made word. The wide use of cochlear prostheses, the first pioneering examples in muscle prosthesis and even an experiment in constructing a retinal prosthesis, are now show the tip of the iceberg. – So far we were able to build systems with high complexity in either space or time, now the spatial-temporal arena is opening. – Software systems around high performance processors and grids of them have become feasible. Moreover, mobile communications with standardized interfaces and huge centers developed successfully. However, thousands of dynamically reprogramable tiny sensory processors, either in a fixed array or in a form of ad-hoc mobile platforms, are emerging around us. Their spatial-temporal control and processing software need a fresh look. Cellular Nonlinear Networks (CNN) [1, 2] are simply an analogue dynamic processor array, made of cells, which contain linear capacitors, linear resistors, linear and nonlinear controlled sources. Many complex computational problems can be formulated as well-defined tasks where the signal values are placed on a regular geometric 2-D or 3-D grid, and the direct interactions between signal values are limited within a finite local neighborhood [3]. CNN is an analog dynamic processor array which reflects just this property: the processing elements interact directly within a finite local neighborhood (see Fig. 1). The concept of CNN is based on some aspects of neurobiology and adapted to integrated circuits. For example, in the brain, the active medium is provided by a sheet-like array of massively interconnected excitable neurons whose energy comes from the burning of glucose with oxygen. In CNN the active medium is provided by the local interconnections of active cells, whose building blocks include active nonlinear devices (e.g., CMOS transistors) powered by dc batteries. 174 A. Slavova, Z. Zafirova, P. Zecca - 6 I 6 I 6 ? R ? R ? (1; 1) (1; 2) (1; 3) (2; 1) (2; 2) (2; 3) (3; 1) (3; 2) (3; 3) 6 I 6 I 6 ? R ? R ? - - Figure 1: 2-D CNN grid 2.1. CNN modelling We shall present CNN modeling on the following example of integro-differential equation. It is a more general form of the Hodgkin–Huxley model for the propagation of the voltage pulse through a nerve axon which is referred to as the FitzHugh–Nagumo equation [9, 10, 11]: (3) ut − uxx = u(u − Θ)(1 − u) − b Z t u(s, x)ds, 0 0 < x, t < 1, 0 < Θ < 1/2, b ≥ 0. The proposed equation (3) is nonlinear parabolic integro-differential equation, in which ut is the first partial derivative of u(t, x) with respect to t, uxx is the second derivative of u with respect to x, u is a membrane potential in a nerve axon, the steady state u = 0 represents the resting state of the nerve. Now if we map u(x, t) into a CNN layer such that the state voltage of a CNN cell vxkl (t) at a grid point (k, l) is associated with u(kh, t), h = ∆x and using the one-dimensional discretized Laplacian template A1 [9], it is easy to design the CNN model of the proposed FitzHugh-Nagumo equation (3): (1) CNN cell dynamics: (4) duj − Ijs = uj (uj − Θ)(1 − uj ) − b dt Z 0 t uj (s)ds. 175 CNN modeling (2) CNN synaptic law: Ijs = (5) 1 (uj−1 − 2uj + uj+1 ). h2 Let us assume for simplicity that the grid size of our CNN model is h = 1 and let us denote the nonlinearity n(uj ) = uj (uj − Θ)(1 − uj ). Substituting (5) into (4) we obtain: (6) duj − (uj−1 − 2uj + uj+1 ) = n(uj ) − dt Z t uj (s)ds, 1 ≤ j ≤ N. −b 0 Equation (6) is actually integro-differential equation which is identified as the state equation of an autonomous CNN made of N × N cells [1,2]. In this section we shall consider more general form of (1) in which u is a function of (t, x): (7) ∂2u ∂2u = − b1 ∂t2 ∂x2 Z 0 t −a1 (t−s) e u(s, x)ds − b2 Z t e−a2 (t−s) u(s, x)ds. 0 Now we shall map the integro-differential equation (7) into CNN architecture made of 2-d grid of L × L cells: (8) duj = vj (t), j = 1, . . . , n = L.L dt dvj = (uj−1 − 2uj + uj+1 )− dt Z t Z t −a1 (t−s) −b1 e uj (s)ds − b2 e−a2 (t−s) uj (s)ds. 0 0 For analytical investigations, it is often necessary to assume an autonomous CNN of infinite size, i.e. L → ∞. In this case, the boundary conditions are replaced by the prescribed behavior of the solution at infinity [9]. We simulate our CNN model (8) [3] and obtain the following figure: 176 A. Slavova, Z. Zafirova, P. Zecca 6 x 10 15 10 5 0 −5 2 1.5 2 1.5 1 1 0.5 0.5 0 v 0 u Figure 2: Simulation of the CNN model 3. Dynamics of the CNN model (8) 3.1. Application of harmonic balance technique for studying dynamics of CNN We shall consider the following CNN system with N cells: (9) ẋi (t) = −xi (t) − syi−1 (t) + pyi (t) + syi+1 (t) yi (t) = f (xi (t)), where the nonlinearity f (.) is the piecewise-linear function f (x) = (|x + 1| − |x − p−1 1|)/2 and s > so that the network operates in a global propagation mode. 2 p−1 In fact it has been shown in [4,5,6] that if the condition s > is not satisfied, 2 every trajectory converges to an equilibrium point, which excludes the existence of any periodic solution. Moreover, let us consider a Continuous-Time DiscreteSpace Fourier Transform (CTDSFT) [8] G(ω, Ω) = F{gn (t)} from continuous time t and discrete space n to continuous temporal frequency ω and continuous spatial frequency Ω, defined by Z ∞ ∞ X −jΩn e gn (t)e−jωt dt. (10) G(ω, Ω) = n=−∞ −∞ By taking the CTDSFT of both members of the first equation in (9), which is linear, we obtain that the transfer function of the linear part of the system can CNN modeling 177 be expressed by (11) H(ω0 , Ω0 ) = X(ω0 , Ω0 ) re−jΩ0 + p + sejΩ0 = , Y (ω0 , Ω0 ) 1 + jω0 where X(ω0 , Ω0 ) = F{xi (t)} and Y (ω0 , Ω0 ) = F{yi (t)}. Since we have a finite array of N cells with periodic boundary conditions, the spatial frequency Ω0 assumes only a finite set of values. The fundamental assumptions made in applying the harmonic balance technique are the following: 1) One supposes that (9) has a periodic solution in the form (12) xi (t) = Xm sin(Ω0 i + ω0 t) which amounts to specify an ansatz for ξ(.) in the general form of the system solution, which is ξ(ψ) = Xm sin ψ. The amplitude, the temporal frequency ω0 , and the spatial frequency Ω0 are therefore the unknowns to be determined. 2) The periodic output yi (t) = f (xi (t)) of the nonlinear block corresponding to the sinusoidal input is approximated only by the fundamental component of its Fourier series expansion yi (t) ≈ Ym sin(Ω0 i + ω0 t) where (13) Ym 1 = π Z π f (Xm sin ψ) sin ψ dψ. −π By considering the Lur’e representation [8] of (9), one can easily determine from (11) that the unknown parameters in (12) have to satisfy the following equations: (14) (15) ℜ{H(ω0 , Ω0 )} = Xm p + 2sω0 sin Ω0 = , Ym 1 + ω02 ℑ{H(ω0 , Ω0 )} = 2s sin ω0 − pω0 = 0. 1 + ω02 The second constraint (15) links the spatial and temporal frequencies by (16) ω0 = 2s sin Ω0 p 178 A. Slavova, Z. Zafirova, P. Zecca whereas using (13) and (14), one obtains that the amplitude Xm of the oscillations is given by the root of the following implicit equation: s " # 2p 1 1 (17) Xm = Xm arcsin + 1− 2 . π Xm Xm The conclusion is summarized in the following proposition [6, 8]. Proposition 1. A CNN (9) with periodic boundary conditions, made of N p−1 cells with s > , possesses at least (N − 1)/2 different nontrivial periodic 2 2πk , 1 ≤ k ≤ (N − 1)/2. solutions, whose spatial frequencies are Ω0 = N In order to validate the accuracy of the achieved result we present numerical simulations for N = 20, p = 2, and s = 3. According to the previous proposition, the network possesses at least nine periodic solutions, moreover the corresponding values of the minimal period T0 and Ω0 , that can be computed are reported in Table 1. Table 1: Comparison between the predicted and the simulated values of T0 and Ω0 Simulation T0 Ω0 7.354 0.314 3.883 0.628 2.837 0.943 2.418 1.257 2.300 1.571 2.418 1.884 2.837 2.119 3.883 2.514 7.353 2.827 Prediction (Prop.1) T0 Ω0 6.778 π/10 3.563 π/5 2.589 3π/10 2.202 2π/5 2.094 π/2 2.202 3π/5 2.589 7π/10 3.563 4π/5 6.778 9π/10 3.2. Periodic solutions of our CNN model In this section we shall apply an approximative method, based on a special Fourier transform in order to obtain periodic solutions of our CNN model (8). The idea of CNN modeling 179 using Fourier expansion for finding the solution of a partial differential equation is well known in physics. It is used to predict what spatial frequences or modes will dominate in nonlinear PDEs. In CNN literature this approach has been developed for analyzing the dynamics of CNNs with symmetric templates [1, 2]. It is known that in the frequency domain, the transformation of the spatial index i into a continuous spatial frequency Ω defines the discrete space Fourier transform (DSFT) [7] by: (18) X xΩ (t) = xi (t)e−jΩi , i √ j = −1, if the array of cells in CNN has an infinite length. Whereas its transformation into a discrete spatial frequency Ω = 2πK/N , 0 ≤ K ≤ N − 1, (N -number of cells), defines a discrete Fourier transform (DFT), if the array of cells is circular and made of a finite number N of cells. We can define the dual continuous time Fourier transform (CTFT) [7] from continuous time t to a continuous temporal frequency ω by Z ∞ xi (t)e−jωt dt. (19) Xi (ω) = −∞ Finally, we can combine both transforms to define a continuous time, discrete space Fourier transform (CTDSFT) from discrete space i and continuous time t to continuous spatial frequency Ω and continuous temporal frequency ω by XZ ∞ (20) XΩ (ω) = xi (t)e−j(Ωi+ωt) dt. −∞ i If the array of cells is circular, this transform can still be applied at the finite set of spatial frequencies: (21) Ω= 2πK , N 0 ≤ K ≤ N − 1. It then becomes a continuous time discrete Fourier transform (CTDFT) from finite discrete space i and continuous time t to discrete spatial frequency Ω and continuous frequency ω: (22) XΩ (ω) = XK (ω) = N Z X i=1 ∞ xi (t)e−j((2πKi/N )+ωt) dt. −∞ 180 A. Slavova, Z. Zafirova, P. Zecca In our case we shall apply double Fourier transform F (s, z) of functions continuous in time and discrete in space [9]: (23) F (s, z) = k=∞ X z −k k=−∞ Z ∞ fk (t) exp(−st)dt. −∞ This is Continuous-Time Discrete-Space Fourier Transform (CTDSFT) from continuous time t and discrete space k to continuous temporal frequency ω, and continuous spatial frequency Ω, such that z = exp(iΩ), s = iω, i is the imaginary identity. The method is based on the fact that all cells in CNN are identical, and therefore by introducing a suitable double transform (23), the network can be reduced to a Lur’e system [8] to which the harmonic balance technique is applied for discovering the existence and characteristics of periodic solutions. We develop the following algorithm based on harmonic balance technique described above: 1. Apply transform (23) to the CNN model (8), so we obtain: (24) sU = V sV = (z −1 − 2z + z)U − b1 s−1 G1 (U ) − b2 s−1 G2 (U ) = (z −1 − 2z + z) − F (U (s, z)), where G1 (U ) and G2 (U ) are the Fourier transforms of e−a1 (t−s) uj (s) and e−a2 (t−s) uj (s). 2. Calculate the transfer function H(s, z): (25) H(s, z) = s(z −1 1 . − 2z + z − s2 ) According to the harmonic balance technique [7,8], the transfer function can be expressed in terms of temporal frequency ω and spatial frequency Ω: (26) HΩ (ω) = 1 . iω(2 − 2 cos Ω + ω 2 ) 3. We are looking for possible periodic solutions of our CNN model (8) in the form: (27) uj (t) = ξ(jΩ + ωt), 1 ≤ j ≤ n = L.L for some function ξ : R → R and for some spatial frequency 0 ≤ Ω ≤ 2π and 2π temporal frequency ω = , where T > 0 is the minimal period of (27). T CNN modeling 181 5. Now according to (27) we shall suppose that uj has the form: (28) uj (t) = Um sin(ωt + jΩ), 1 ≤ j ≤ n, which amounts to specify an ansatz for (28) ξ(ψ) = Um sin ψ. 6. We shall approximate the output of our CNN model (8) by the fundamental component of its Fourier expansion: (29) y(uj (t)) = yj (t) = Ym sin(ωt + jΩ). 7. The ratio of the CTDSFT of these periodic solutions is: (30) HΩ (ω) = Um . Ym 8. According to the harmonic balance technique [8] the following constraints hold: (31) Um , Ym I(HΩ (ω)) = 0. R(HΩ (ω)) = 10. Thus (28), (29) and (31) give us necessary set of equations for finding the unknowns Um , Ω and ω. As we mentioned before we are looking for a periodic wave solution of (8), therefore Um will determine approximate amplitude of the 2π wave, an T = will determine the wave speed. Now according to the harmonic ω balance technique, if for a given value of Ω we can find the unknowns (Um , ω), then we can predict the existence of a periodic solution of our CNN model (8) 2π . The following with an amplitude Um and period of approximately T = ω theorem hold: Theorem 1. CNN model (8) of the integro-differential equation (7) with circular array of n = L × L cells has periodic solutions uj (t) with a finite set of 2π 2πk , 0 ≤ k ≤ n − 1 and a period T = . spatial frequencies Ω = n ω P r o o f. We shall consider circular array of cells for our CNN model (8) [9], for which the possible values of Ω can be easily obtained. As uj (t) is assumed to be periodic with minimal period T , one has (32) ξ(jΩ + ωt) = ξ(jΩ + ωt + kωT ), 182 A. Slavova, Z. Zafirova, P. Zecca for any k ∈ N. On the other hand, the periodic boundary conditions impose that (33) ξ(ωt) = ξ(Ωn + ωt). Combining (32) with j = 0 and (33), we get (34) Ω= k 2πk ωT = , 0 ≤ k ≤ n − 1, N n where the range of k is determined by the condition 0 ≤ Ω ≤ 2π. From (29), (31) and (34) we obtain a system of algebraic equations for the unknowns Um , Ω and ω. We solve this system and obtain the unknowns (Um , ω) for a given value of Ω. Therefore according to harmonic balance technique [8] we have proved the existence of periodic solution uj (t) of our CNN model (8) with 2π an amplitude Um and minimal period T = . Thus the theorem is proved.  ω Based on the above algorithm extensive computations were made [3] and we obtain the following periodic solutions of CNN model (8): Figure 3: Periodic solution of CNN model (8) CNN modeling 183 REFERENCES [1] L. O. Chua, L. Yang. Cellular neural networks: Theory. IEEE Trans. Circuits Syst. 35, 10 (1988), 1257–1272. [2] L. O. Chua, L. Yang. CNN: Applications. IEEE Trans. Circuits Syst. 35, 10 (1988), 1273–1299. [3] CNN Software Library, Ver. 1.1, Analogical and Neural Computing Laboratory. Budapest, 2000. [4] R. Genesio, A. Tesi. A harmonic balance approach for chaos prediction: the Chua’s circuit. Int.J.of Bifurcation and Chaos 2, 1 (1992), 61–79. [5] R. Genesio, A. Tesi. Harmonic balance methods for the analysis of chaotic dynamics in nonlinear systems. Automatica 28, 3 (1992), 531–548. [6] R. Genesio, A. Tesi, F. Villoresi. A frequency approach for analyzing and controlling chaos in nonlinear circuits. IEEE Trans. Circuits Syst. – I: Fund. Theory Appl. 40, 11, (1993), 819–827. [7] J. Macki, P. Nistri, P.Zecca. A theoretical justification of the method of harmonic balance for systems with discontinuities. Geoffrey J. Butler Memorial Conference in Differential Equations and Mathematical Biology (Edmonton, AB, 1988). Rocky Mountain J. of Math. 20, 4, (1990), 10791099. [8] A. I. Mees. Dynamics of feedback systems. Chichester, John Wiley & Sons, Ltd., 1981. [9] A. Slavova. Cellular Neural Networks:Dynamics and Modeling. Kluwer Academic Publishers, 2003. [10] A. Slavova, P.Zecca. CNN model for studying FitzHugh-Nagumo equation. C. R. Acad. Bulg. Sci. 53, 6 (2000), 31–34. [11] A. Slavova, P. Zecca. CNN model for studying dynamics and travelling wave solutions of FitzHugh-Nagumo equation. J. Comput. Appl. Math. 151, 1 (2003), 13–24. 184 A. Slavova, Z. Zafirova, P. Zecca Angela Slavova Institute of Mathematics and Informatics Bulgarian Academy of Sciences 1113 Sofia, Bulgaria e-mail: slavova@math.bas.bg Zoya Zafirova Faculty of Applied Mathematics and Informatics Technical University of Sofia 1000 Sofia, Bulgaria e-mail: zafirova@tu-sofia.bg Pietro Zecca Department of Mathematics and Informatics University of Florence 50139 Florence, Italy e-mail: pzecca@unifi.it