Table of Contents Author Guidelines Submit a Manuscript
Discrete Dynamics in Nature and Society
Volume 2015 (2015), Article ID 542507, 8 pages
Research Article

Simulation of a Vibrant Membrane Using a 2-Dimensional Cellular Automaton

1Modeling and Simulation Laboratory, Centro de Investigación en Computación, Instituto Politécnico Nacional, Avenida Juan de Dios Bátiz, Esquina Miguel Othón de Mendizábal, Colonia Nueva Industrial Vallejo, 07738 Gustavo A. Madero, DF, Mexico
2Postgraduate and Research Section, Escuela Superior de Cómputo, Instituto Politécnico Nacional, Avenida Juan de Dios Bátiz, Esquina Miguel Othón de Mendizábal, Colonia Lindavista, 07738 Gustavo A. Madero, DF, Mexico
3Instituto Mexicano del Petróleo, Eje Central Lázaro Cárdenas Norte 152, 07730 Gustavo A. Madero, DF, Mexico

Received 21 April 2015; Revised 4 June 2015; Accepted 8 June 2015

Academic Editor: Tetsuji Tokihiro

Copyright © 2015 I. Huerta-Trujillo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


This paper proposes a 2-dimensional cellular automaton (CA) model and how to derive the model evolution rule to simulate a two-dimensional vibrant membrane. The resulting model is compared with the analytical solution of a two-dimensional hyperbolic partial differential equation (PDE), linear and homogeneous. This models a vibrant membrane with specific conditions, initial and boundary. The frequency spectrum is analysed as well as the error between the data produced by the CA model. Then it is compared to the data provided by the solution evaluation to the differential equation. This shows how the CA obtains a behavior similar to the PDE. Moreover, it is possible to simulate nonclassical initial conditions for which there is no exact solution using PDE. Very interesting information could be obtained from the CA model such as the fundamental frequency.

1. Introduction

The infinitesimal calculus and its descendants [1] have been one of the dominant branches in mathematics since it was developed by Newton and Leibniz. Differential equations are the main core of calculus and have been the cornerstone for understanding sciences, particularly physics. It is often necessary to consider more than one variable since it is required to change this in function to another one or the time. Therefore, to model physical systems, the differential equations have had great success due to more than three centuries of experience with methods to give the symbolic solutions. Even so, few partial differential equations have an exact solution [2].

Moreover, the current necessity to experiment with physical systems in order to recognize their behavior make necessary to develop models that simulate the systems and become less complex allowing manipulation and approximation as close as possible to reality. In this order of ideas, the discrete techniques have had more success when they have been implemented for simulation purposes. An example of this is the cellular automata (CA). Wolfram [3, 4] defined the CA as a mathematical idealization of physical systems whose time and space are discrete and in which the physical quantities can be grouped in a finite set of values. The CA are appropriate in physical systems with a highly nonlinear regime such as chemical or biological systems, where there are discrete thresholds [5].

For example, the two-dimensional wave equation is an important one since it represents the hyperbolic partial differential equations. Although it has many analytical solutions, if the initial or boundary conditions are changed, it cannot be solved. There are attempts to simulate the wave equation behavior using a CA as presented by Chopard and Droz [6] and Kawamura et al. [7]; they have only been performed in a one-dimensional automaton. This due to the complexity of applying a suitable rule for the CA evolution.

The evolution rule for CA is a fundamental problem that has no analytical solution [2]. In this paper a methodology for clarifying the CA evolution rule that simulates the behavior of a vibrating membrane is compared with the analytical solution of hyperbolic PDE. This simulates a proposed frequency spectrum observing similarity solutions. This shows that the discrete modeling can be an alternative for systems in which the boundary conditions or other circumstances make the PDE intractable. Section 1 offers a brief introduction; Section 2 presents the two-dimensional continuous wave equation model as well as the solution to the equation given the initial and boundary conditions. Section 3 shows the methodology used to make the problem discrete and obtain the CA evolution rule and the discrete model definition. Section 4 analyses and discusses the results of the continuous model against the proposed CA discrete model. Finally, Section 5 gives the conclusions.

2. Vibrating Membrane, Classical Model

The motion equation for a membrane is based on the assumption that it is thin and uniform with negligible stiffness that is perfectly elastic and without spring vibrating with small amplitude movements. The equation governing the transverse membrane vibrations is given bywhere is the membrane deflection and , where is the membrane mass density [8] and is membrane tension per unit length. Equation (1) is called wave equation in two dimensions.

Considered as a square membrane (see Figure 1), the boundary conditions are defined as follows:

Figure 1: Rectangular membrane of length.

The initial conditions are defined aswhere and are the membrane displacement and initial speed, respectively.

The initial conditions are defined as

Then, the repose membrane starts with a spatial distribution that can be seen in Figure 2.

Figure 2: The membrane system initial and boundary conditions.

Differential equation solution (1) with the initial conditions and boundary shown in (4) and (5), a length , and a constant is

The dynamic response of the central membrane point according to (6) can be observed in Figures 3(a) and 3(b).

Figure 3: PDE graph simulation (a) presents the cell displacement on the -axis; (b) zoom-in of the PDE simulation graph.

3. Discrete Membrane Model

This paper considers a synchronous CA [9] where the underlying topology is a rectangular two-dimensional finite lattice. For this type of CA there are two classic neighborhoods; the case of the CA is considered a Von Neumann neighborhood type [10] consisting of a central cell and its four geographical neighbors to the north, south, east, and west (). In this sense, the importance of extending the CA to two dimensions lies in that [11] the extension is the emergence of two-dimensional patterns that have no simple analogue in one dimension, and two-dimensional dynamic sometimes allows a direct physical systems comparison.

The methodology used for the discretization of the PDE on a -dimensional CA is presented.

3.1. Formal Definition of Cellular Automata

There are different authors [4, 9, 11, 12] that have defined the CA, but all coincide in four elements. Taken as a base a CA is as follows.

Definition 1. A CA is a 4-tuple whereis a regular lattice for a -dimensional lattice,is a finite set of all possible cell states, ,is a finite set of cells that define the cell neighborhood,, is a transition function applied simultaneously to the cells that conform the lattice [9].

The objective cell state actualization requires knowing the neighborhood cell states [6, 9].

3.2. Analytical Model Discretization

Suppose there is a membrane which can be thought of as a succession of specific mass points connected by springs (see Figure 4). Each point of the membrane is attached to the four orthogonal neighbors, where the membrane mass is spread over the attachment points and not on the springs, and the membrane boundary is subject to a fixed surface.

Figure 4: Representation of the membrane as a succession of masses and springs, each mass joined to its neighbors as north, south, east, and west.

Let us call (or equilibrium length) the spring length that joins two masses in an equilibrium state.

For a vibrating membrane system with initial position conditions described in (4) and starting in repose. Each membrane junction point is subjected to four forces acting toward each orthogonal neighbor, called , , , , and , for the central cell, north, south, east, and west, respectively (see Figure 5).

Figure 5: Representation of the central cell and its four orthogonal neighbors forming a Von Neumann CA neighborhood, as well as the direction in which the neighbor force acts, represented by the direction of the vectors.

Suppose it is necessary to know the force each neighbor applied over . Following the defined process by Huerta-Trujillo et al. [13], the force that a neighbor exercises over the central cell is computed (see Figure 6).

Figure 6: Representation of the vectors associated with the central cell, a neighbor, and the direction of the force exercised by the spring in the neighbor cell direction; , , and represent the central cell vectors, the west neighbor, and the spring vector that join them, respectively.

Now, taking into account Figure 6 gives

Knowing that represents the separation length between masses, it is possible to write this scalar as the sum of the spring length increasing the deformation undergone by the same spring due to the central mass position change. Thus,where is the spring deformation value; therefore it is necessary to compute this with the purpose of knowing the force increase from the equilibrium to the analysis point. Thus,

Given that both sections of (9) are vectors, their components are represented as follows:and the second vector components

Equating component to component of (10) and (11) gives

So, the rectangular components , , and are obtained corresponding to the displacement increase of the axis coordinates , , and for the mass for . The component values for the vectors , , and are calculated in the same manner.

Following the analysis, for Hooke’s law, given for the particle, there are three forces exercised for in the direction of due to the vector components , , and and similar to each of the particles , , and in the directions , , and , respectively. Thus, the force applied by the neighbors in the component is

Considering that the springs that join the masses to the membrane are equal gives so this yields

In the same manner

These are the forces acting over , which define that the acceleration at the time is oscillating. Since the force is known and the acceleration stays constant in each time step, using the straight-line motion equation with constant acceleration, in terms of the initial speed of for its rectangular components, yieldsAccordinglywhere can be computed using second Newton’s law.

The previous equations define the speed of after time . This permits calculating the new particle position for the same instant time using the following equations:

Speed equations (16), (17a), and (17b) and position (18a), (18b), and (18c) are those used in the evolution function for the proposed CA.

3.3. Vibrant Membrane Model Using a 2-Dimensional CA

Using the developed analysis defines the CA model for a vibrant membrane system, fixed at the borders, of an area as a 4-tuple , where each cell is defined by its mass, initial position, and speed. When the membrane is found in rest state, this gives the following.: this is a regular 2-dimensional lattice.::  : .,where is the acceleration that the neighbors of exercise over the cell in time ; is the final cell position in space, composed of three variables , , and , and according to its rectangular coordinates, each one is represented by bits; is the final speed in time , composed of three variables , , and , and according to its rectangular coordinates, each one is represented by 32 bits. This implementation allows having real states without contradicting the definition of cellular automaton [4] (moreover, the total states that each cell can take are restricted to states).

The evolution function is composed of two rules, both applied simultaneously to all the cells that conform the lattice. This is different from the model proposed by Glabisz [14, 15] where the rules are applied in an asymmetric form.

The rule (a) defines the cell position at time , taking the speed at time . This position is updated, to be the new starting position for and so on. Similarly, for (b) the final speed for time is updated, with the initial speed for time .

An important point in the model definition is the spring restitution coefficient. Unlike other models [13] where the constant restitution is relatively simple to calculate, this model faces a problem: The arrangement of masses and springs has no serial pattern but shows an array of grids. To obtain this coefficient, it follows the process defined by Huerta-Trujillo et al. [16].

3.4. Restitution Coefficient Calculation

To calculate the restitution coefficient value of the springs that bind the corresponding mass, the process begins by assuming that the membrane has a spring constant with a value equal to .

Starting with a membrane represented by the proposed CA consisting of nodes, thereupon, following the reduction springs connected in series and in parallel obtains

Following the process for a CA consisting of nodes gives carrying on the process, the relation for springs has the form where is the spring constant, is the membrane elasticity constant, and is the number of nodes for a CA of nodes.

4. Simulation and Results

CA was implemented applying the object-oriented paradigm using the C++ language in which each cell is implemented as a Bean object which has as attributes the cell spatial location, the value of its mass, its instantaneous speed, or if it is a fixed cell. CA space represented by the square lattice is regularly implemented as an object composed of cell arrangements. The CA as such is composed of a grid object and one that implements the evolution rule and the definition of the CA neighborhood.

The CA dynamic response is given by defining a membrane represented by a lattice of nodes, , based on the necessary number of cells to simulate a linear system [13]. For CA the basis is a membrane cm, stretched by of its length in the direction of its axes and for a length equal to that used in the PDE taken as reference node located at to match the point of the membrane. The mass of cells is proportional to the density used to solve the PDE as shown in Section 2. The CA response is shown in Figures 7(a) and 7(b). Qualitatively, the CA response is relevant to the response of the PDE, taking the same initial and boundary conditions for both models.

Figure 7: Response graph (a) of the simulation made by the CA, presenting the cell on the -axis, (b) zoom-in of the simulation graphic.

The initial CA conditions are the same as those described for the PDE model ((4) and (5)). The displacement was obtained by the PDE membrane ranging and took presented samples. The displacement that was obtained from the CA proposed the same cell (in this case the cell located at ) at the same time and with the same number of samples. Fifty simulations were performed with the same averaged data for comparison against the obtained PDE. Figures 8(a) and 8(b) show overlapping graphs obtaining oscillations by CA (in solid line) and the PDE (in dotted line) as well as a zoom-in to the same graph. The mean square error is defined aswhere is the th value estimated by the CA and is the reference value taken by the PDE. The error between the measurements generated by the CA and the reference values given by the PDE is , which ensures that the CA simulates the PDE response.

Figure 8: Overlay graphic (a) of the simulation realised by the CA and the PDE shows the displacement of the cell on the -axis in both cases; (b) zoom-in graphic overlay.

In contrast to other models [7, 15], the proposed CA model uses less cells in the lattice to get the PDE response.

In a quantitative form there is a corresponding phase between the two models for the same cell. Analytically, the fundamental frequency is limited by the constants values and . The frequency is defined as [17]

So, a fundamental frequency is expected approximately as

Obtaining frequency spectra for both signals and plotting them observe that there is a congruence between the frequency spectra given by the CA and PDE (see Figure 9(a)). Figure 9(b) shows the fundamental frequency for both models around 85 Hz. The fundamental frequency error between the analytical solution and that simulated by the proposed CA is , which is smaller than that found in other models [7, 14, 18].

Figure 9: Graphic of frequency spectra overlay (a) of the simulation realized by the CA and the PDE, (b) zoom-in of the graphic overlay.

Furthermore, the robustness of the CA model can be tested given a type of irregular initial condition, for example, a random initial position for all cells, with , where is the cell position on the -axis. The density, tension, and boundary conditions are the same as described in Section 2.

In this case the PDE does not have an analytical solution since it is not possible to define a function to describe the initial conditions and to be derivable, but the proposed CA model endures this kind of initial conditions. Simulating the conditions described and plotting the same point as in the previous cases, the CA model provides the central cell response. The results can be seen in Figure 10(a); the frequency spectra of this evolution are shown in Figure 10(b) where the fundamental frequency is approximately such as that obtained in the previous cases.

Figure 10: Graphic of the simulation of lattice’s central cell for the CA model (a); the graphic (b) shows frequency spectra for the central cell.

5. Conclusions

The -dimensional CA model to a vibrant membrane system shown in this paper was derived from a -dimensional CA model that simulates a vibrant string [13]. The CA membrane model is released from the initial conditions so it is not necessary to redefine it; thus it is possible to define the initial conditions and simulate the system behavior immediately. This is different to the 2-dimensional wave PDE that is sensitive to these conditions. Its solution could only be found if the initial conditions are represented by a derivative function in all space; otherwise the PDE would not have an analytical solution. The proposed CA model can be set to almost any initial condition to acquire the response due to that definition.

The model expansion from one to two dimensions entails increasing the model computational complexity, linear in the case of one dimension, to complexity because CA space now has an array of . In this sense, the parallel model is justified if is defined as the number of threads that help the development of CA per unit time. Then, the complexity of the CA will be and if is made large enough such that , then the complexity is reduced to .

On the other hand, in our experience it is possible to extend this model to a three-dimensional CA; the challenge is to obtain the relationship between the elasticity of the material being modeled and the restitution constant for the CA internal springs.

Finally, it is worth mentioning that, in this work, there is not a straight relationship between the CA model and the PDE. From this paper it is clear that the results between the PDE and CA are in excellent agreement. Moreover, the cellular automata could simulate systems which are simulated by PDE under conditions that these equations could not. The latter suggest that perhaps it is possible to find a mathematical transformation from PDE to CA.

Conflict of Interests

All the authors of the paper declare that there is no conflict of interests regarding the publication of this paper.


The authors wish to thank CONACYT, COFAA-IPN (Project no. 20151704), and EDI-IPN for the support given for this work.


  1. J. Paulos, Beyond Numeracy, Vintage Series, Vintage Books, 1992.
  2. T. Toffoli, “Cellular automata as an alternative to (rather than an approximation of) differential equations in modeling physics,” Physica D. Nonlinear Phenomena, vol. 10, no. 1-2, pp. 117–127, 1984. View at Publisher · View at Google Scholar · View at MathSciNet
  3. S. Wolfram, “Statistical mechanics of cellular automata,” Reviews of Modern Physics, vol. 55, no. 3, pp. 601–644, 1983. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  4. S. Wolfram, A New Kind of Science, Wolfram Media, Champaign, Ill, USA, 2002.
  5. S. Wolfram, “Cellular automata as models of complexity,” Nature, vol. 311, no. 5985, pp. 419–424, 1984. View at Publisher · View at Google Scholar · View at Scopus
  6. B. Chopard and M. Droz, Cellular Automata Modeling of Physical Systems, Cambridge University Press, New York, NY, USA, 1998. View at Publisher · View at Google Scholar · View at MathSciNet
  7. S. Kawamura, T. Yoshida, H. Minamoto, and Z. Hossain, “Simulation of the nonlinear vibration of a string using the Cellular Automata based on the reflection rule,” Applied Acoustics, vol. 67, no. 2, pp. 93–105, 2006. View at Publisher · View at Google Scholar · View at Scopus
  8. M. V. Walstijn and E. Mullan, “Time-domain simulation of rectangular membrane vibrations with 1-d digital waveguides,” in Proceedings of the 2011 Forum Acusticum, pp. 449–454, Aalborg, Denmark, 2011.
  9. J. Kari, “Theory of cellular automata: a survey,” Theoretical Computer Science, vol. 334, no. 1–3, pp. 3–33, 2005. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  10. S. H. White, A. M. del Rey, and G. R. Sánchez, “Modeling epidemics using cellular automata,” Applied Mathematics and Computation, vol. 186, no. 1, pp. 193–202, 2007. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  11. A. Ilachinski, Cellular Automata: A Discrete Universe, World Scientific, 2001. View at Publisher · View at Google Scholar · View at MathSciNet
  12. M. Kutrib, R. Vollmar, and T. Worsch, “Introduction to the special issue on cellular automata,” Parallel Computing, vol. 23, no. 11, pp. 1567–1576, 1997. View at Google Scholar
  13. I. Huerta-Trujillo, J. Chimal-Eguía, N. Sánchez-Salas, and J. Martínez-Nuño, “Modelo de autómata celular 1-dimensional para una EDP hiperbólica,” Research in Computing Science, vol. 58, pp. 407–423, 2012. View at Google Scholar
  14. W. Glabisz, “Cellular automata in nonlinear string vibration,” Archives of Civil and Mechanical Engineering, vol. 10, no. 1, pp. 27–41, 2010. View at Publisher · View at Google Scholar · View at Scopus
  15. W. Glabisz, “Cellular automata in nonlinear vibration problems of two-parameter elastic foundation,” Archives of Civil and Mechanical Engineering, vol. 11, no. 2, pp. 285–299, 2011. View at Publisher · View at Google Scholar · View at Scopus
  16. I. Huerta-Trujillo, E. Castillo-Montiel, J. Chimal-Eguía, N. Sánchez-Salas, and J. Martínez-Nuño, “AC 2-dimensional como modelo de una membrana vibrante,” Research in Computing Science, vol. 83, pp. 117–130, 2014. View at Google Scholar
  17. L. Kinsler, Fundamentals of Acoustics, Wiley, 2000.
  18. S. Kawamura, M. Shirashige, and T. Iwatsubo, “Simulation of the nonlinear vibration of a string using the cellular automation method,” Applied Acoustics, vol. 66, no. 1, pp. 77–87, 2005. View at Publisher · View at Google Scholar · View at Scopus