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

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.


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.

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 by where (, , ) is the membrane deflection and  2 = /, 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: (, , ) = 0 for  = 0,  = ,  = 0,  =  ∀.(2) The initial conditions are defined as where (, ) and (, ) are the membrane displacement and initial speed, respectively.The initial conditions are defined as (, , )          =0 = 0.
Then, the repose membrane starts with a spatial distribution that can be seen in Figure 2.
The dynamic response of the central membrane point  = (/2, /2) according to (6) can be observed in Figures 3(a

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 2-dimensional CA is presented.

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 CA = (, , , Φ) where 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].

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.
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).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).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: Equating component to component of ( 10) and (11) gives 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, yields Accordingly where  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.

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 AC = (, , , Φ), 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.

󳨀 → 𝑎 𝑡
V  is the acceleration that the neighbors of   exercise over the cell in time ;   →  +1   is the final cell position in space, composed of three variables   ,   , and   , and according to its rectangular coordinates, each one is represented by 32 bits; →  +1   is the final speed in time  + 1, composed of three variables V  , V  , and V  , 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 2 193 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  + 1, taking the speed at time .This position is updated, to be the new starting position for  + 2 and so on.Similarly, for (b) the final speed for time  + 1 is updated, with the initial speed for time  + 2.
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].

𝑘 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 2 × 2 nodes, thereupon, following the reduction springs connected in series and in parallel obtains Following the process for a CA consisting of 3 × 3 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.

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,  = 51, based on the necessary number of cells to simulate a linear system [13].For CA the basis is a membrane 10 × 10 cm, stretched by 20% 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 (/2, /2) to match the point (/2, /2) 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.
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 1 and took 1 × 10 4 presented samples.The displacement that was obtained from the CA proposed the same cell (in this case the cell located at (/2, /2)) 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 as where Vac  is the th value estimated by the CA and Vr  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 MSE = 4.4189 × 10 −11 , which ensures that the CA simulates the PDE response.
In contrast to other models [7,15], the proposed CA model uses 50% 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] 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 0.995%, which is smaller than that found in other models [7,14,18].
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 0 <   ≪ , 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 (/2, /2) 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 85 Hz such as that obtained in the previous cases.

Conclusions
The 2-dimensional CA model to a vibrant membrane system shown in this paper was derived from a 1-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 ( 2 ) 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 ( 2 /) 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.

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

Figure 5 :Figure 6 :
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.

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

Figure 9 :
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.

Figure 10 :
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.
So, the rectangular components Δ  , Δ  , and Δ  are obtained corresponding to the displacement increase of the axis coordinates , , and  for the mass   for