On The Theory of Partial Difference Equations: From Numerical Methods to Language of Complexity

On the Theory of Partial Difference Equations: From Numerical Methods to Language of Complexity Bik Kuang Min August 2025 Abstract This work develops a theoretical framework for Partial Difference Equations (P∆E) as a natural mathematical language for modeling discretetime, discrete-space systems. Motivated by the limitations of continuous partial differential equations (PDE) in representing inherently discrete phenomena, we begin by defining P∆E in terms of discrete function spaces and shift operators, contrasting them with ordinary difference equations (O∆E) and PDE, and clarifying the scope of our study. We then examine linear P∆E, outlining their main types, providing formal definitions, and presenting selected analytic solutions in simple cases. Building on this, we introduce the discrete functional analytic setting: discrete function spaces, Hilbert space structure, and discrete operators, including difference and shift operators, and study their algebraic and adjoint properties. The discrete Green’s function is also defined within this framework. As a demonstration of the framework’s unifying power, we reformulate a wide range of well-known discrete models, including elementary cellular automata, coupled map lattices, Conway’s Game of Life, the Abelian sandpile model, the Olami–Feder–Christensen earthquake model, forest fire models, the Ising model, the Kuramoto Firefly Model, the Greenberg–Hastings Model, and the Langton’s ant as explicit P∆E. For each case, we focus on obtaining a concise and mathematically elegant formulation rather than detailed dynamical analysis. Finally, we compare the “discrete universe” of P∆E with the continuous universe of PDE, highlighting their structural parallels and their respective connections to discrete mathematics and continuous analysis. This reveals P∆E and PDE as mathematical “twins”, analogous in form yet rooted in fundamentally different underlying mathematics. The generality of the P∆E formalism suggests broad applicability, from modeling biological and ecological processes to analyzing complex networks, emergent computation, and other spatiotemporally extended systems. 1 Contents Preface 5 1 Introduction to Partial Difference Equations 7 1.1 Motivation and Background . . . . . . . . . . . . . . . . . . . . . 7 1.2 DiscreteFunction........................... 7 1.3 ShiftOperator ............................ 8 1.4 Partial Shift Operator . . . . . . . . . . . . . . . . . . . . . . . . 8 1.5 Ordinary Difference Equations . . . . . . . . . . . . . . . . . . . 8 1.6 Partial Difference Equation . . . . . . . . . . . . . . . . . . . . . 9 1.7 The Order of a Partial Difference Equation . . . . . . . . . . . . 9 1.8 Comparison with Ordinary Difference Equations . . . . . . . . . 10 1.9 Comparison with Partial Differential Equations . . . . . . . . . . 11 1.10 Advantages of Rewriting Complex Systems as Partial Difference Equations............................... 11 1.11 Scope and Limitations . . . . . . . . . . . . . . . . . . . . . . . . 12 1.12 Classification and Examples . . . . . . . . . . . . . . . . . . . . . 13 1.13Notations ............................... 14 2 Linear Partial Difference Equations 16 2.1 Motivation .............................. 16 2.2 Definition ............................... 16 2.3 Discrete 1D Transport Equation . . . . . . . . . . . . . . . . . . 17 2.4 Discrete 1D Diffusion Equation . . . . . . . . . . . . . . . . . . . 20 2.5 Discrete 1D Wave Equation . . . . . . . . . . . . . . . . . . . . . 22 3 Discrete Function Spaces and Operators 24 3.1 Motivation .............................. 24 3.2 Discrete Function Spaces . . . . . . . . . . . . . . . . . . . . . . . 24 3.3 Operators............................... 25 3.4 The Kronecker Delta Function . . . . . . . . . . . . . . . . . . . 32 3.5 The Discrete Green’s Function . . . . . . . . . . . . . . . . . . . 33 3.6 Boundary Value Problems . . . . . . . . . . . . . . . . . . . . . . 35 4 Elementary Cellular Automata 37 4.1 Introduction.............................. 37 4.2 Notation Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.3 TheRule90.............................. 39 4.4 Rule30 ................................ 43 4.5 TheRule153 ............................. 44 4.6 Rule45 ................................ 45 4.7 Rule105................................ 46 4.8 Rule169................................ 47 4.9 Rule210................................ 48 4.10Rule137................................ 49 2 4.11Rule110................................ 51 5 Other 1D Nonlinear Partial Difference Equations 52 5.1 Introduction.............................. 52 5.2 Sierpinski Fan Equation . . . . . . . . . . . . . . . . . . . . . . . 52 5.3 Mod 5 Sierpi´nski Equation . . . . . . . . . . . . . . . . . . . . . . 53 5.4 Mod 4 Sierpinski Fractal Equation . . . . . . . . . . . . . . . . . 54 5.5 Sierpinski Carpet Equation . . . . . . . . . . . . . . . . . . . . . 55 6 Coupled Map Lattice 57 6.1 1D Coupled Map Lattice . . . . . . . . . . . . . . . . . . . . . . . 57 6.2 2D Coupled Map Lattice . . . . . . . . . . . . . . . . . . . . . . . 57 6.3 Explanation.............................. 58 7 The Conway’s Game of Life 59 7.1 Conway’sEquation.......................... 59 7.2 Interpretation............................. 59 8 Abelian Sandpile Model 60 8.1 Abelian Sandpile Equation . . . . . . . . . . . . . . . . . . . . . 60 8.2 Remarks................................ 60 8.3 The Sandpile Fractal . . . . . . . . . . . . . . . . . . . . . . . . . 61 8.4 OctaSandpileModel......................... 62 8.5 Explanation for the Self-Organized Criticality . . . . . . . . . . . 64 8.6 Explanation for the Power Law Distribution . . . . . . . . . . . . 65 9 OFC Earthquake Model 67 9.1 OFC Earthquake Equation . . . . . . . . . . . . . . . . . . . . . 67 9.2 Observations and Phenomena . . . . . . . . . . . . . . . . . . . . 67 10 Forest Fire Model 68 10.1 Forest Fire Equation . . . . . . . . . . . . . . . . . . . . . . . . . 68 10.2Definitions............................... 68 10.3 Transition Logic (via mod 3) . . . . . . . . . . . . . . . . . . . . 68 10.4Interpretation............................. 69 10.5 Observations and Phenomena . . . . . . . . . . . . . . . . . . . . 69 11 Kuramoto Firefly Model 70 11.1 Kuramoto Equation . . . . . . . . . . . . . . . . . . . . . . . . . 70 11.2 Principle of Locality . . . . . . . . . . . . . . . . . . . . . . . . . 70 11.3 Observations and Phenomena . . . . . . . . . . . . . . . . . . . . 70 12 Ising Model 72 12.1 Governing Equation . . . . . . . . . . . . . . . . . . . . . . . . . 72 12.2 Physical Interpretation . . . . . . . . . . . . . . . . . . . . . . . . 72 3 13 Greenberg–Hastings Model 73 13.1 Greenberg-Hastings Equation . . . . . . . . . . . . . . . . . . . . 73 13.2Explanation.............................. 73 14 The Langton’s Ant 74 14.1 The Langton’s Ant Equations . . . . . . . . . . . . . . . . . . . . 74 14.2 Observations and Phenomena . . . . . . . . . . . . . . . . . . . . 75 15 The Discrete Parallel Universe of PDE 76 15.1Motivation .............................. 76 15.2 Calculus and Discrete Calculus . . . . . . . . . . . . . . . . . . . 76 15.3ODEandO∆E............................ 77 15.4 Initial and Boundary Conditions . . . . . . . . . . . . . . . . . . 77 15.5 Laplace Transform and Z Transform . . . . . . . . . . . . . . . . 78 15.6 Fourier Transform and Discrete Fourier Transform . . . . . . . . 79 15.7 Discrete and Continuous Green’s Function . . . . . . . . . . . . . 80 15.8 Gaussian Distribution and Combinatorics . . . . . . . . . . . . . 80 15.9 Discrete and Continuous Dynamical Systems . . . . . . . . . . . 81 15.10Manifolds and Networks . . . . . . . . . . . . . . . . . . . . . . . 82 15.11Boolean Algebra, Theoretical Computer Science and Number Theory................................. 83 15.12Summary ............................... 84 16 Summary 85 16.1 From Numerical Methods to Language of Complexity . . . . . . 85 4 Preface Ever since entering university, I have been trying to model biological evolution. At first, I attempted to use graph theory, but I failed... graph theory could not describe ring species. Then I tried partial differential equations (PDE), but again I failed, as PDE could not describe surface splitting; I wanted to model cellular division, but it was simply beyond reach. One night, a sudden image flashed in my mind: countless organic molecules colliding near a hydrothermal vent on the seafloor... Ah... this is a dynamical system! I suddenly realized: the initial concentration and distribution of molecules are initial conditions, the size of the environment is the boundary condition, and the temperature, pH, salinity, etc., are non-autonomous terms. This is a dynamical system! Then I understood that biology is inherently a discrete system. Molecules, cells, individuals—they are all discrete units (discrete space). Reproduction happens across generations (discrete time). So we should not use differential equations. We should use difference equations. But I couldn’t figure out what partial difference equations were supposed to be. What do they look like? Then it suddenly hit me: the Game of Life, the Sandpile Model... I realized these are, in fact, partial difference equations! They all share one essential feature: discrete space and discrete time. Thus, I inadvertently unified all complex system models. Then I asked myself: partial differential equations describe the evolution of continuous fields. Therefore, partial difference equations describe the evolution of discrete fields! This led me to systematically replace cellular automata, classical evolutionary theory, and many biological models with a new framework I now call Discrete Field Theory. I began incorporating rich terminology from theoretical physics and classical field theory. I realized this theory is far too vast for one person to complete alone. Thus, I sincerely invite all scientists to participate in the development and refinement of this framework. Your contributions will be deeply appreciated. This work is the result of countless nights of failure, reflection, and sudden clarity. I hope this story helps the reader understand not only the theory... but the journey that led to its birth. 5 Prologue: An Invitation, Not a Declaration This paper represents only a fragment of a broader and more ambitious framework that I call Discrete Field Theory. I do not claim to have discovered a final theory for complex systems, rather, I am aware that this framework is vast, unfinished, and in need of many more concrete models and rigorous theorems. As such, I humbly and sincerely invite mathematicians, scientists, and engineers around the world to join in the effort of developing this theory further. If there are any errors in this work, I welcome constructive criticism and corrections. If you find value in this paper and wish to contribute or discuss further, please feel free to contact me. This is my e-mail, b[email protected] 6 1 Introduction to Partial Difference Equations 1.1 Motivation and Background Many systems in modern biology exhibit inherently discrete organization: molecular interaction networks, cellular assemblies, and ecological populations are composed of countable entities; cells occupy spatially distinct locations; reproduction and generational turnover occur in discrete steps. Even gene expression and mutation often manifest as abrupt, stochastic events rather than smooth, continuous processes. Nevertheless, the predominant mathematical frameworks for modeling such systems, particularly partial differential equations (PDE) are formulated in terms of continuous variables and differentiable fields [12,13]. While PDE-based models have been highly successful in physics and engineering, their reliance on smoothness and continuity can make them less well-suited for describing discrete interactions, sudden transitions, or spatially distributed logical rules that frequently arise in biological and other complex systems [2,9]. This motivates the question of whether a discrete analogue of PDE, which we term a Partial Difference Equations (P∆E), that is, a discrete-space, discretetime equation governing the evolution of a field defined on a lattice, could serve as a more natural modeling framework for complex, emergent systems. The aim of this work is to develop a functional-analytic foundation for P∆E and to demonstrate their potential for representing systems in which local, discrete interactions give rise to large-scale collective phenomena. We propose that such a formalism offers a mathematically coherent, computationally grounded, and potentially more faithful alternative to continuous PDE in domains where discreteness is an essential feature. 1.2 Discrete Function In this section, we introduce the concept of a discrete function, which serves as the fundamental object in our framework for discrete dynamical systems. Definition 1.1 (Discrete Function).Adiscrete function is a mapping f: Ω ⊆Zn→C, where Znis the discrete n-dimensional integer lattice, and Cis the codomain of the function (real, integer, or complex values depending on context). In this paper, we primarily focus on the case where f:Zn→R, i.e., realvalued discrete functions. Example 1.1 (Discrete Single-variable Function).Let f:Z→Rbe defined as f(x) = sin(x). This is a real-valued function defined on the integer lattice Z. 7 Example 1.2 (Discrete Multivariable Function).Let u:Z2→Zbe defined as u(t, x) = ⌊x+ sin(t)⌋. Here, umaps discrete spacetime coordinates (t, x)to an integer value via the floor function applied to a real expression. 1.3 Shift Operator We define the shift operator Ekacting on a discrete function x(t) as Ekx:= x(t+k) for any integer k∈Z. Examples. Ex =x(t+ 1) E2x=x(t+ 2) E−1x=x(t−1) 1.4 Partial Shift Operator Definition. Given a discrete scalar field u(x1, x2, . . . , xn), the partial shift operator Eki xiacts on uby shifting the i-th coordinate: Eki xiu:= u(x1, . . . , xi+ki, . . . , xn) Examples. Etu=u(t+ 1, x, y, z) E2 xu=u(t, x + 2, y, z) In this paper, we will write all the Partial Difference Equations in terms of Shift Operators. 1.5 Ordinary Difference Equations Definition 1.2 (Ordinary Difference Equation of Order k).An ordinary difference equation of order kis an equation involving a single-variable function x:Z→C, expressed in terms of shift operators: F(Enx, En−1x, . . . , Ex, x, E−1x, . . . , E−mx, t) = 0, where Eix:= x(t+i),n, m ∈Z+,k=m+nis the total order, and Fis a (possibly nonlinear) function. 8 Definition 1.3 (Definition: Linear Ordinary Difference Equation of Order k). Alinear ordinary difference equation of order kis an equation of the form: an(t)Enx+an−1(t)En−1x+···+a0(t)x+a−1(t)E−1x+···+a−m(t)E−mx=g(t), where: •x=x(t)is the unknown function defined on Z, •Ejx:= x(t+j)is the shift operator, •aj(t)are known coefficient functions, •g(t)is the known input or forcing term, •m, n ∈Nand k=m+nis the total order. 1.6 Partial Difference Equation Definition 1.4 (Partial Difference Equation).Let u:Zn→R(or C) be a scalar function defined on the discrete lattice. A partial difference equation (P∆E) with constant order is an equation of the form: FEk1 x1Ek2 x2···Ekn xnu(k1,k2,...,kn)∈S, x1, x2, ..., xn= 0, where: •u=u(x1, x2, ..., xn) •Fis a given function, which may be linear or nonlinear; •Exidenotes the shift operator in the xidirection, defined by Eki xiu=u(x1, x2, . . . , xi+ki, . . . , xn); •S⊂Znis a finite index set determining the set of applied shifts. 1.7 The Order of a Partial Difference Equation We propose two types of order to classify the structure of partial difference equations with constant order: structural order and absolute order. Definition 1.5 (Structural Order of Partial Difference Equation).Let u= u(x1, x2, . . . , xn)be a scalar function defined on a discrete lattice. A partial difference equation is of the form FEk1 x1Ek2 x2···Ekn xnu(k1,k2,...,kn)∈S, x1, x2, ..., xn= 0, where S⊂Znis a finite index set of shifts. 9 2 Linear Partial Difference Equations 2.1 Motivation Although most well-known cellular automata are nonlinear partial difference equations, we believe it is inappropriate to skip the study of linear partial difference equations and proceed directly to the nonlinear case. However, due to limitations in our current capabilities, this section will be limited to introducing the formal definition of linear partial difference equations, along with three classical examples: the discrete one-dimensional transport equation, the discrete one-dimensional diffusion equation, and the discrete one-dimensional wave equation. We will present a few particular solutions for each equation, rather than attempting to develop a general solution method. 2.2 Definition Definition 2.1 (Linear Partial Difference Equation).Let u=u(x1, x2, . . . , xn) be a scalar function defined on the discrete lattice Zn. A linear partial difference equation is an equation of the form X (k1,...,kn)∈S ak1k2...kn(x1, x2, . . . , xn)Ek1 x1Ek2 x2···Ekn xnu(x1, x2, . . . , xn) = f(x1, x2, . . . , xn), where: •S⊂Znis a finite index set, •Exiis the shift operator in the xidirection: Exiu(...,xi, . . . ) = u(...,xi+ 1, . . . ), •ak1k2...knand fare known functions on Zn. Definition 2.2 (Linear Partial Difference Equation of Absolute Order k).Let u=u(x1, x2, . . . , xn)be a scalar function on a discrete lattice Zn. A linear partial difference equation of absolute order kis an equation of the form: X |k1|+···+|kn|≤k ak1k2...kn(x1, x2, . . . , xn)Ek1 x1Ek2 x2···Ekn xnu=f(x1, x2, . . . , xn), •(k1, k2, ..., kn)∈Zn. •Exiis the shift operator in the xidirection: Exiu(...,xi, . . . ) = u(...,xi+ 1, . . . ), •ak1k2...knand fare known functions on Zn. 16 2.3 Discrete 1D Transport Equation We consider the Discrete Transport Equation: Etu=Ek xu where u:Z2→Ris defined on discrete spacetime. Initial condition: u(0, x) = f(x) We impose no boundary condition. Interpretation: Each value is transported rigidly along a discrete path. The effective velocity is −k, since data moves from x+kto x. Characteristic curve: x(t) = x0+kt Solution We consider the partial difference equation Etu=Ek xu, k ∈Z, with no boundary conditions. Let us define a new variable: ξ=x+kt, and define a new function v(t, ξ) := u(t, x). Then we compute the time-shift of v: Etv=v(t+1, ξ) = v(t+1, x+k(t+1)) = v(t+1, x+kt+k) = v(t+1, ξ+k) = EtEk ξv=EtEk xu Etv=EtEk xu v=Ek xu On the other hand, from the original equation: Etu=u(t+ 1, x) = Ek xu(t, x) = u(t, x +k). Substitute u(t, x) = v(t, ξ), we get: Etu=u(t+ 1, x) = v(t+ 1, ξ) = Etv. 17 And Ek xu=v Thus, Etv=v, or equivalently, Etv=v. This is just an Ordinary Difference Equation. Now we define a new function of ξonly: v(t, ξ) = f(ξ), which solves the equation since Etv=v. Therefore, the general solution is v(t, ξ) = f(ξ) = f(x+kt), so that u(t, x) = f(x+kt), where f(x) = u(0, x) is the initial condition. Conclusion: The solution of the partial difference equation Etu=Ek xu is given by u(t, x) = f(x+kt), where fis the initial profile of uat time t= 0. 18 A General Form of Discrete 1D Linear Transport Equation Etu=p(t, x)Ek(t,x) xu+q(t, x)u+g(t, x) The operator Ek(t,x) xis linear because it does not depend on the function u. Although this may seem counterintuitive at first, we can rigorously verify its linearity: Ek(t,x) x(au +bv) = a u(t, x +k(t, x)) + b v(t, x +k(t, x)) = a Ek(t,x) xu+b Ek(t,x) xv Thus, Ek(t,x) xis a linear operator. Furthermore, we observe that the order of a difference operator does not necessarily need to be a constant; it can also be a function. To avoid confusion, we refer to the classical case where the order is constant as difference equations with constant order. When the order is a function, we refer to them as difference equations with variable order. This is a type of Functional Partial Difference Equations. In this paper, we will only discuss about Partial Difference Equations with constant order. Unless otherwise specified, the term ”difference equation” typically assumes constant order by default. 19 2.4 Discrete 1D Diffusion Equation In this section, we propose a discrete analogue of the classical one-dimensional diffusion equation. The equation is given by ∆tu=α∇2u, where u=u(t, x), u :Z2→R, ∆tu:= Etu−u, ∇2u:= Exu+E−1 xu−2u, and α∈Ris a constant representing the diffusion coefficient. Physical Interpretation. This equation models the diffusion of some quantity (e.g., heat, particles, or chemical concentration) across a one-dimensional discrete lattice over discrete time steps. The Laplacian term ∇2ucaptures the spatial flow of the quantity due to differences between neighbouring sites, while ∆tudescribes the temporal evolution of the system. The Discrete 1D Diffusion Equation models symmetric spreading on a 1D lattice. Special Case We consider the case where the diffusion coefficient α=1 2. Then the equation becomes the discrete heat equation: Etu=1 2Exu+E−1 xu. We first consider the initial condition: u(0, x) = δ(x), with no boundary condition. Then the solution is given by the discrete heat kernel: u(t, x) = δ(mod2(x+t)) Ct, x+t 21 2t , where the binomial coefficient C(x, y) is defined as: C(x, y) =    x! (x−y)!y!,if 0 ≤y≤x, 0,otherwise, with C:Z2→Z. We omit the proof since this solution is well known as the discrete heat kernel, or Green’s function for the discrete heat equation. 20 If the initial condition is instead: u(0, x) = f(x), with no boundary condition, then the solution is given by: u(t, x) = X s∈Z f(s)δ(mod2(x−s+t)) Ct, x−s+t 21 2t . Here, mod2(x) := x−2x 2is the modulo function. 21 2.5 Discrete 1D Wave Equation We define the discrete 1D wave equation as δt∆tu=c2∇2u where u=u(t, x), with u:Z2→R, and the operators are defined via shift operators: δt∆tu=Etu+E−1 tu−2u, δtu=u−E−1 tu(backward difference), ∇2u=Exu+E−1 xu−2u(discrete Laplacian). Physical Interpretation. This equation models a system of coupled oscillators arranged on a one-dimensional lattice. Each oscillator interacts with its nearest neighbours, and the term ∇2urepresents the net restoring force. The parameter cis a constant that governs the wave propagation speed in the discrete medium. Special Case We consider the case where c= 1. Then the equation becomes: Etu=E−1 xu+Exu−E−1 tu with initial conditions: u(0, x) = δ(x), u(1, x) = δ(x) and no boundary condition (i.e., defined on the full lattice Z). Exact Solution: u(t, x) = δ(t)δ(x) + θ(t−1)(−1)|x|+t−1·com(x, t −1) where com(x, t) is a compact support indicator function defined by: com(x, a) := (1,if |x| ≤ a 0,otherwise This can be generalized to be centered at an arbitrary point x0: com(x−x0, a) := (1,if |x−x0| ≤ a 0,otherwise The alternating sign (−1)|x|+t−1introduces a checkerboard-like structure in the wave propagation, while the compact support ensures finite spread with speed 1. The following picture is the spatiotemporal plot of the Discrete 1D Wave Equation. 22 Figure 1: Spatiotemporal plot of the 1D Wave Equation 23 3 Discrete Function Spaces and Operators 3.1 Motivation This section lays the theoretical foundation for the analysis of partial difference equations (P∆E). We systematically introduce the structure and properties of discrete function spaces, difference and shift operators, discrete analogues of Hilbert spaces, the Kronecker delta function, discrete Green’s functions, and boundary value problems. The notation and operator definitions adopted here largely follow established conventions in the literature (see, e.g., [6]). The formulation and conceptual framework are inspired by the classical theory of partial differential equations (see, e.g., [15]), with appropriate adaptations to the discrete setting. This section provides a rigorous functional-analytic formulation of shift operators. Readers primarily interested in applications may proceed directly to Section 4. 3.2 Discrete Function Spaces Definition 3.1 (Discrete Function Space).Let Ω⊆Znbe a discrete domain. We define the discrete function space over Ωas F(Ω) := {f: Ω →C}, i.e., the set of all functions from Ωto the complex numbers. Definition 3.2 (Discrete LpSpace).Let Ω⊆Zn. We define the discrete Lp space for 1≤p < ∞as Lp(Ω) := (f: Ω →CX x∈Ω|f(x)|p<∞), with the corresponding Lpnorm defined by ∥f∥p:= X x∈Ω|f(x)|p!1/p . For p=∞, we define L∞(Ω) := f: Ω →C sup x∈Ω|f(x)|<∞, with the corresponding L∞norm given by ∥f∥∞:= sup x∈Ω|f(x)|. 24 Definition 3.3 (Discrete Hilbert Space).Let Ω⊆Zn. The discrete Hilbert space is the space L2(Ω) defined as L2(Ω) := (f: Ω →CX x∈Ω|f(x)|2<∞), equipped with the inner product ⟨f, g⟩:= X x∈Ω f(x)g(x),for all f, g ∈ L2(Ω). This inner product satisfies the following properties for all f, g, h ∈ L2(Ω) and all scalars α∈C: 1. Conjugate Symmetry: ⟨f, g⟩=⟨g, f⟩. 2. Linearity in the First Argument: ⟨αf +h, g⟩=α⟨f, g⟩+⟨h, g⟩. 3. Positive Definiteness: ⟨f, f⟩ ≥ 0,and ⟨f, f⟩= 0 ⇐⇒ f= 0. The norm induced by this inner product is ∥f∥2:= p⟨f, f⟩= X x∈Ω|f(x)|2!1/2 . Moreover, L2(Ω) is complete with respect to this norm, and hence forms a Hilbert space. 3.3 Operators Definition 3.4 (Discrete Operator).Let V, W be discrete function spaces, i.e., subspaces of F(Ω), where Ω⊆Zn. Adiscrete operator is a map T:V→W, which assigns to each function f∈Va function Tf ∈W. Definition 3.5 (Shift Operator).Let x:Z→C. We define the shift operator Ekacting on a discrete function x(t)by Ekx:= x(t+k) for any integer k∈Z. 25 3.4 The Kronecker Delta Function Definition 3.24 (Kronecker Delta Function).The Kronecker delta function is defined as δ(x−a) = (1,if x=a 0,otherwise for any x, a ∈Z. Definition 3.25 (Multivariable Kronecker Delta Function).Let δ:Zn→Rbe defined, for x= (x1, . . . , xn)and s= (s1, . . . , sn)∈Zn, by δ(x1−s1, . . . , xn−sn) := n Y i=1 δ(xi−si), where δ(x)is the one-dimensional Kronecker delta function Proposition 3.26 (Delta Representation of Discrete Functions).Let Ω = Zn and consider the Hilbert space L2(Ω) with inner product ⟨f, g⟩:= X x∈Ω f(x)g(x). For each s∈Zn, define the multivariable Kronecker delta δ(x−s) := n Y i=1 δ(xi−si). Then every f∈ L2(Zn)can be represented as f(x) = X s∈Zn f(s)δ(x−s), where the series converges in L2. Moreover, the family {δ(x−s) : s∈Zn} forms an orthonormal basis of L2(Zn). Proof. For m,n∈Zn, ⟨δ(x−m), δ(x−n)⟩=X x∈Zn δ(x−m)δ(x−n). The product is nonzero only when x=m=n, hence ⟨δ(x−m), δ(x−n)⟩=(1,m=n, 0,m=n. Thus they are mutually orthogonal and each satisfies ∥δ(x−s)∥2 2=⟨δ(x−s), δ(x−s)⟩= 1, 32 so the set is orthonormal. Finally, for any f∈ L2(Zn), f(x) = X s∈Zn⟨f, δ(x−s)⟩δ(x−s), and ⟨f, δ(x−s)⟩=f(s), which completes the proof. Example 3.4 (One-dimensional case).When n= 1, the multivariable Kronecker delta reduces to the standard one-dimensional Kronecker delta δ(x−s) := (1, x =s, 0, x =s, x, s ∈Z. In L2(Z)with inner product ⟨f, g⟩:= X x∈Z f(x)g(x), the family {δ(x−s) : s∈Z}satisfies ⟨δ(x−m), δ(x−n)⟩=(1, m =n, 0, m =n, so it forms an orthonormal basis of L2(Z). Every f∈ L2(Z)admits the representation f(x) = X s∈Z f(s)δ(x−s), where the series converges in L2. 3.5 The Discrete Green’s Function Definition 3.27 (Discrete Green’s Function).Let Ω⊆Zn, and let x= (x1, x2, . . . , xn)∈Ω. Consider a linear partial difference equation L(u) = f(t, x), where t∈Zis the discrete time variable, Lis a linear shift operator acting on the spatial variables, and the equation is subject to zero initial condition and certain boundary conditions on Ω. The discrete Green’s function G(t, x)is defined as the unique solution of L(G) = δ(t)δ(x), subject to the causality condition G(t, x) = 0, t < 0, where δ(t)and δ(x)are Kronecker delta functions. 33 Proposition 3.28 (Solution as Convolution of the Green’s Function).Let Ω⊆ Znand consider the linear partial difference equation L(u) = f(t, x),x∈Ω, t ≥0, where Lis a linear shift operator in the spatial variables. Let G(t, x)be the discrete Green’s function, i.e., the solution of L(G) = 0 with G(0,x) = δ(x). Then the unique solution u(t, x)subject to zero initial condition is given by u(t, x) = t X τ=0 X s∈Ω G(t−τ, x−s)f(τ, s). Proof of “Solution as Convolution of the Green’s Function”. Fix Ω ⊆Zn. Let Lbe a linear, causal evolution operator that is shift–invariant in space and time, and let G(t, x) denote the discrete Green’s function, i.e. the response to a unit space–time impulse at the origin: input δ(t)δ(x)7→ output G(t, x), t ≥0. Step 1 (delta expansion of the forcing). For any forcing f:Z≥0×Ω→Cwe have the (space–time) delta representation f(t, x) = ∞ X τ=0 X s∈Ω f(τ, s)δ(t−τ)δ(x−s), with convergence in L2(Ω) for each fixed t(or in a suitable function space depending on the problem). Step 2 (responses to elementary impulses). Let eτ,s(t, x) := δ(t−τ)δ(x−s). By timeand space-shift invariance of L, the response to eτ,sis the shifted Green’s function uτ,s(t, x) = G(t−τ, x−s) (understood as 0 for t<τ). Step 3 (linear superposition). By linearity of L, the solution to L(u) = fwith zero initial condition is the superposition of the elementary responses weighted by the coefficients f(τ, s) from Step 1: u(t, x) = t X τ=0 X s∈Ω f(τ, s)G(t−τ, x−s). (The upper limit treflects causality; terms with τ > t vanish.) Thus u=G∗fis the discrete space–time convolution of the Green’s function with the forcing, which proves the stated formula. 34 3.6 Boundary Value Problems We restrict to the one-dimensional case with spatial index x∈ {L, L+1, . . . , M} and unit grid spacing. In the discrete setting, the following are direct analogues of the classical boundary conditions for continuous PDE, obtained by replacing derivatives with finite difference operators. Dirichlet boundary condition. Given constants a, b ∈C, the Dirichlet condition fixes the value of uat both ends: u(t, L) = a, u(t, M) = b. Physically, this corresponds to prescribing a fixed state (e.g. temperature, displacement) at the boundary points. Neumann boundary condition. In the continuous PDE setting, a Neumann boundary condition prescribes the normal derivative of the solution at the boundary: ∂xux=L=a, ∂xux=M=b, where a, b ∈Care prescribed constants. Physically, the Neumann condition represents a prescribed flux across the boundary. For example, in the heat equation, ∂xucorresponds to the heat flux according to Fourier’s law. In the context of partial difference equations, we replace the spatial derivative with the discrete central difference operator ∆c x(unit grid spacing) and introduce ghost points: ∆c xux=L=u(t, L + 1) −u(t, L −1) 2=a, ∆c xux=M=u(t, M + 1) −u(t, M −1) 2=b, where u(t, L−1) and u(t, M +1) lie outside the computational domain. Solving these equations for the ghost points yields: u(t, L −1) = u(t, L + 1) −2a, u(t, M + 1) = u(t, M −1) + 2b. This “ghost point” technique allows the central difference stencil to be applied at the boundary while preserving second-order accuracy. Robin boundary condition. In the continuous PDE setting, a Robin boundary condition is a linear combination of the normal derivative and the function value: ∂u ∂x +p(x)u=q(x), x ∈ {L, M}, where p(x) and q(x) are prescribed functions on the boundary (possibly constants). Physically, Robin conditions often model convective exchange or springtype constraints. 35 In the discrete setting, we replace ∂xwith the central difference operator ∆c x: ∆c xu+p(x)u=q(x), x ∈ {L, M}. Using ghost points, for example at x=L: u(t, L + 1) −u(t, L −1) 2+p(L)u(t, L) = q(L), so the ghost point u(t, L −1) can be eliminated as u(t, L −1) = u(t, L + 1) −2q(L)−p(L)u(t, L). A similar formula holds at x=M. This preserves the central difference stencil at the boundary while enforcing the Robin condition. 36 4 Elementary Cellular Automata 4.1 Introduction The study of one-dimensional cellular automata (1D CA) was significantly advanced by Stephen Wolfram [22], who systematically explored the behaviour of all 256 elementary rules. These models, despite their simplicity, can generate a wide range of complex spatiotemporal patterns, including periodic structures, nested fractals, localized structures, and even chaotic behaviour. All 1D cellular automata can be expressed as autonomous partial difference equations, meaning that the update rule does not include any external forcing term. Moreover, most of these equations are inherently nonlinear, due to the logical (Boolean) nature of the local interactions. In this section, we formulate several 1D cellular automata as partial difference equations. Some of these formulations are based on known Boolean algebra transformations, converting logical rules into Boolean polynomials, and then into partial difference equations. Others are derived heuristically by observing the pattern of evolution. This reformulation has several important advantages: •It provides a unified mathematical framework to study discrete systems. •It allows the application of analytical tools such as operator theory, stability analysis, and even Green’s functions. •It enables the exploration of different types of initial and boundary conditions in a more formalized way. •It makes it possible to introduce non-autonomous terms and study how external forcing influences the evolution. This approach transforms cellular automata from mere computational toys into objects of rigorous mathematical investigation within the broader context of discrete dynamical systems. 4.2 Notation Evolution In the work by Itoh and Chua [8], one-dimensional cellular automata were expressed in the form of difference equations using the notation xn i. For example, Rule 90 can be written as: xn+1 i=xn i−1+xn i+1 (mod 2). However, they did not explicitly introduce the concept of partial difference equations, and continued to refer to such formulations simply as “difference equations.” In contrast, in the present work we explicitly adopt the notion of partial difference equations. My earlier notation was, for example, in the case of Rule 90: u(t+ 1, x) = u(t, x −1) + u(t, x + 1) (mod 2). 37 Although clear for low-dimensional cases, this style becomes cumbersome when multiple spatial variables are involved. For example, with u(t, x, y, z), the equations become long, cluttered, and visually unappealing. Moreover, in more complex models such as the Game of Life or the forest fire model, I initially used custom functions, resulting in update rules that were messy and syntactically similar to programming code rather than mathematical expressions. To address this, I spent several months carefully defining a set of reusable building blocks, including the Kronecker delta function, Heaviside step function, modulo function, shift operators, difference operators, and the discrete Laplacian. With these tools, the notation becomes significantly more compact and expressive. For example, Rule 90 can now be written in operator form as: Etu= mod2(Exu+E−1 xu), u :Z2→Z2. The advantages are even more apparent in higher dimensions. Consider the two-dimensional discrete diffusion equation in traditional notation: u(t+1, x, y) = u(t, x, y)+u(t, x−1, y)+u(t, x+1, y)+u(t, x, y−1)+u(t, x, y+1)−4u(t, x, y), which is long and cumbersome. In operator notation, this simply becomes: Etu=u+∇2u, or equivalently, ∆tu=∇2u. This operator-based notation is not only more concise, but also opens the door to applying discrete functional analysis, operator theory, and other mathematical tools to the study of partial difference equations. 38 4.3 The Rule 90 The Rule 90 cellular automaton can be written as a nonlinear partial difference equation of the form: Etu= mod2E−1 xu+Exu where mod2(x) is defined as: mod2(x) := (1,if xis odd 0,if xis even This function is clearly non-linear, making the equation itself non-linear. We refer to this equation as the Sierpi´nski Equation, due to its deep connection with the Sierpi´nski triangle. Define a discrete delta function as: δ(x−a) := (1,if x=a 0,otherwise Given the initial condition u(0, x) = δ(x) and no boundary condition, the exact solution to the equation is: u(t, x) = mod2(C(2t, x +t)) where C(x, y) is the binomial coefficient, defined as: C(x, y) := (x! (x−y)! y!,if 0 ≤y≤x 0,otherwise For a general initial condition u(0, x) = f(x), the solution becomes: u(t, x) = mod2 X s∈Z f(s)·C(2t, x +t−s)! This system exhibits chaotic behaviour in a discrete binary field. Remarkably, despite its chaotic dynamics, we have found a closed-form analytical solution. This provides a foundation to define and analyze concepts like chaotic functions and quasichaotic functions within the framework of discrete functional analysis. 39 Example 1 This is the spatiotemporal plot of Rule 90 with the initial condition u(0, x) = δ(x) and no boundary conditions. Figure 2: Spatiotemporal plot of Rule 90 with u(0, x) = δ(x) and no boundary conditions. It clearly forms a Sierpi´nski triangle. Note that each cell’s state depends only on its left and right neighbours. Thus, this fractal is not globally planned, but emerges from local interactions. 40 Example 2 This is the spatiotemporal plot of Rule 90 with the initial condition u(0, x) = δ(x) and boundary conditions u(t, −100) = u(t, 100) = 0. Figure 3: Spatiotemporal plot of Rule 90 with u(0, x) = δ(x) and boundary conditions u(t, −100) = u(t, 100) = 0. Initially, we observe a perfect Sierpi´nski triangle. However, once the pattern reaches the boundary, it breaks down into spatiotemporal chaos and becomes unpredictable. 41 4.9 Rule 210 The Boolean update rule for Rule 210 is given by: Etu= mod2E−1 xu·u·Exu+E−1 xu·u·δ(Exu) + E−1 xu·δ(u)·δ(Exu) + δ(E−1 xu)·δ(u)·Exu Example Initial condition: u(0, x) is random. Boundary conditions: u(t, −200) = u(t, 200) = 0. Figure 8: Spatiotemporal plot of Rule 210 with random initial condition and fixed boundaries. As shown in the figure, the system exhibits a highly unusual structure: one region resembles a Sierpi´nski triangle, while the other shows diagonally arranged motifs. This dual pattern reflects complex local interactions and potentially bifurcating spatial dynamics. 48 4.10 Rule 137 We define the evolution equation as follows: Etu= mod2E−1 xu−πu(π+u+Exu)+ 1−u 1 + u+3E−1 xu−4u each state value satisfies u(t, x)∈ {0,1}. This rule is an example of a nonlinear partial difference equation, operating over a binary state space. Interestingly, the system exhibits a mixture of order and apparent randomness. Some regions evolve into highly regular patterns, while others appear chaotic and unpredictable. This suggests that the system may operate near the edge of chaos, a regime known for its potential computational richness and complexity. Although we do not currently possess a closed-form analytic solution to this equation, the intricate structure of its evolving patterns can still be studied using tools such as combinatorics, discrete functional analysis, and spatiotemporal statistics. Example Initial condition: u(0, x) is random. Boundary conditions: u(t, −200) = u(t, 200) = 0. 49 Figure 9: Spatiotemporal plot of Rule 137 with random initial condition and fixed boundaries. The pattern generated by Rule 137 closely resembles that of Rule 110. It exhibits complex behaviour with both periodic islands and chaotic seas. There are large regions of crystalline structure, along with distinct triangular motifs, indicating rich dynamics and spatial self-organization. 50 4.11 Rule 110 Using Boolean algebra, we derive the evolution equation for Rule 110 as follows: Etu= mod2u+Exu+u·Exu+u·Exu·E−1 xu, where Etdenotes the time shift operator, Exis the spatial right-shift operator, and mod2denotes addition modulo 2. Example We consider the following setup: •Initial condition: u(0, x) is chosen randomly with values in {0,1}. •Boundary condition: u(t, −200) = u(t, 200) = 0 for all t. The figure below illustrates the space-time evolution of the system under Rule 110: 51 5 Other 1D Nonlinear Partial Difference Equations 5.1 Introduction Beyond classical two-state cellular automata such as Rule 90 or Rule 137, we have discovered that a wide variety of other one dimensional nonlinear partial difference equations (P∆E) can be constructed. These equations may involve arithmetic operations, logical operations, or floor and absolute value functions. They do not necessarily follow the strict rule-based framework of cellular automata, yet they remain fully deterministic and discrete in both space and time. Such equations often exhibit rich and unexpected dynamical behaviour. Some generate spatially localized chaos, while others display periodic or quasiperiodic structures. This diversity suggests that one-dimensional nonlinear P∆E form a broader and more general class of systems than traditionally studied automata. In the following subsections, we will explore several such examples and discuss their qualitative features. 5.2 Sierpinski Fan Equation Consider the following nonlinear partial difference equation: Etu= mod3E−1 xu+u+Exu with initial condition: u(0, x) = δ(x),and no boundary condition. u=u(t, x)∈ {0,1,2}Here, the delta function δ(x) is defined as: δ(x) = (1, x = 0 0, x = 0 This system produces a fascinating fractal structure resembling a ”fan” or ”tree”, hence we refer to it as the Sierpinski Fan Equation. Although we currently do not have an analytic solution to this equation, the emergent pattern clearly exhibits self-similarity and modular periodicity, which are typical features of modular arithmetic-driven cellular automata. This suggests a rich underlying structure that deserves further investigation. 52 Figure 10: Spatiotemporal evolution of the Sierpinski Fan Equation 5.3 Mod 5 Sierpi´nski Equation We define a new nonlinear partial difference equation as follows: Etu= mod5E−2 xu+E−1 xu+u+E1 xu+E2 xu with the initial condition: u(0, x) = δ(x) where: •δ(x) is the discrete delta function defined by δ(x) = 1 if x= 0, and 0 otherwise; •u=u(t, x)∈ {0,1,2,3,4}for all t, x; •The system is updated on a discrete 1D lattice with no boundary condition. This equation can be seen as a generalization of the Rule 90 or Sierpi´nski triangle in modulo 2. The modulo 5 version creates a far richer fractal structure. From the spatiotemporal plot, we observe that instead of the classic single upside-down triangle in each upright triangle (as in the binary Sierpi´nski triangle), this pattern exhibits ten upside-down triangles within each upright triangle. This suggests a more intricate self-similar structure with higher-order modular symmetry. We call this the Mod 5 Sierpi´nski Equation. 53 Figure 11: Spatiotemporal evolution of the Mod 5 Sierpi´nski Equation 5.4 Mod 4 Sierpinski Fractal Equation We consider the following nonlinear partial difference equation: Etu= mod4E−1 xu+u+Exu+δ(x) with initial condition: u(0, x) = δ(x) and no boundary condition. u=u(t, x)∈ {0,1,2,3} Here, δ(x) is the discrete delta function defined by: δ(x) = (1,if x= 0 0,otherwise The values of u(t, x) are restricted to {0,1,2,3}, and we visualize them using the following colour scheme: Value Colour 0 White 1 Green 2 Purple 3 Yellow The spatiotemporal plot of this equation is shown below: 54 We observe that this system generates a strikingly complicated fractal structure, reminiscent of a Sierpinski triangle, but with richer periodic colour layers. Despite its simplicity, the system exhibits highly nontrivial dynamics that may reflect deep combinatorial or algebraic patterns. 5.5 Sierpinski Carpet Equation We consider the following discrete evolution equation on the integer lattice: Etu= mod3Exu+E−1 xu+E−1 tu, where the state variable is u:Z2→R, and the modulo function is defined by modn(x) = x−nx/n. The initial condition is given by a discrete delta function: u(0, x) = u(1, x) = δ(x), with no imposed boundary conditions. 55 Numerical experiments indicate that the space-time evolution of this system generates a striking self-similar pattern. In particular, the resulting structure coincides with the classical Sierpi´nski carpet fractal. Interestingly, if the modulo operation is removed (i.e. considered over the integers), the recurrence resembles a generalized Pascal-type relation, but with an additional memory term u(t− 1, x). This observation suggests that the solution might be interpreted as a novel form of combinatorial numbers, whose exact combinatorial interpretation remains an open question. This equation was first observed by Wu Han (private communication, 2025), who noticed that it generates the Sierpi´nski carpet. 56 6 Coupled Map Lattice 6.1 1D Coupled Map Lattice The one-dimensional Coupled Map Lattice (CML), originally proposed by Kaneko [24], is a prototypical model for spatiotemporal chaos. It is defined by the following recurrence relation: ut+1 s= (1 −ε)f(ut s) + ε 2f(ut s+1) + f(ut s−1), where ut s∈R,f:R→Ris a nonlinear local map (e.g., logistic map), and ε∈[0,1] is the coupling strength. We reformulate this equation using shift operators to obtain a compact operator-based representation: Etu= (1 −ε)f(u) + ε 2f(Exu) + f(E−1 xu), where u=u(t, x), u :Z2→R. This equation is usually nonlinear autonomous Partial Difference Equation that captures spatiotemporal dynamics and is frequently used to model chaotic behaviour on discrete lattices. 6.2 2D Coupled Map Lattice The two-dimensional Coupled Map Lattice (2D CML) extends the 1D case to a square lattice, where each site interacts with its four nearest neighbours. The classical form is given by: ut+1 i,j = (1 −ε)f(ut i,j) + ε 4f(ut i+1,j) + f(ut i−1,j) + f(ut i,j+1) + f(ut i,j−1), where ut i,j ∈R,f:R→R, and ε∈[0,1] as before. Using shift operators, the system can be rewritten in a compact operatorbased form: Etu= (1 −ε)f(u) + ε 4f(Exu) + f(E−1 xu) + f(Eyu) + f(E−1 yu), where u=u(t, x, y), u :Z3→R. This operator representation highlights the spatial symmetry and facilitates generalization to higher dimensions, anisotropic coupling, or graph-based topologies. 57 Figure 12: Octa Sandpile Model with Moore neighborhood (8-neighbour interaction) and threshold zc= 8. The system consists of 107sand grains evolved over many toppling steps. The resulting pattern exhibits eightfold rotational symmetry (C8), fractal self-similarity, and discrete energy layering. 8.5 Explanation for the Self-Organized Criticality Consider the comparison between Rule 90 and the Sandpile model. Rule 90 is defined by the update rule: Etu= mod2E−1 xu+Exu Both Rule 90 and the Sandpile model are based on local interactions and exhibit nonlinear dynamics. However, they display drastically different macroscopic behaviours: Rule 90 leads to global chaos, while the Sandpile model exhibits self-organized criticality. In Rule 90, any small perturbation is exponentially amplified due to the absence of thresholds, every update directly sums the neighbouring values mod64 ulo 2, regardless of the current state. As a result, the system exhibits global sensitivity to initial conditions, a hallmark of chaos. In contrast, the Sandpile model incorporates a threshold-based collapse mechanism. Before any site reaches the threshold (typically 4 in 2D), the system behaves linearly: each grain of sand simply increases the local state by 1. Once the threshold is reached, a toppling event occurs, distributing sand to neighbours and potentially triggering further topplings. This collapse phase introduces nonlinearity and chaotic propagation, but only in localized regions that have reached the critical state. Therefore, the Sandpile model alternates between two distinct regimes: •Linear growth phase: Sites far below the threshold evolve smoothly and predictably. •Chaotic collapse phase: Sites at the threshold trigger avalanches with unpredictable, nonlinear dynamics. Because the initial configuration is typically random, some sites are nearzero (linear regime), while others are near the threshold (nonlinear regime). This coexistence of order and chaos in space and time allows the Sandpile model to maintain a dynamically balanced state without the need for fine-tuning of parameters. This explains why the Sandpile model self-organizes into a critical state: the rules themselves inherently support a mixture of linear and chaotic zones, allowing the system to spontaneously hover at the edge of chaos. Unlike systems that require tuning to a critical point, the Sandpile model exhibits critical behaviour as an emergent property of its local update rules. 8.6 Explanation for the Power Law Distribution The avalanche size distribution in the Sandpile model follows a power law, typically of the form: P(s)∝s−τ where P(s) is the probability of an avalanche of size s, and τis a constant exponent. This behaviour arises from the interplay between local dynamics, chaotic propagation, and external driving. We explain it in three parts: 1. Why are small avalanches more common than large ones? Because the dynamics are based on local interactions, triggering a large avalanche requires a long chain of causally connected sites to be near the threshold, which is statistically less likely. In contrast, small avalanches require only a short chain and are thus much more probable. 2. Why do large avalanches still occur frequently (not exponentially rare)? The collapse dynamics are chaotic: once a critical site topples, small perturbations can lead to unpredictable and explosive cascades, this is the 65 butterfly effect. Therefore, even though rare, large avalanches are possible due to the system’s inherent sensitivity at the critical state. 3. Why is there no characteristic scale in avalanche size? The update rules impose no intrinsic limit on the spatial extent of an avalanche. The avalanche size depends entirely on the current configuration of the system. If a small localized region is near-critical, only a small avalanche occurs. But if a large domain is close to the threshold, a large avalanche may follow. Furthermore, the external driving is random and slow, allowing the system to evolve into configurations where large avalanches are statistically allowed. Thus, the system exhibits scaleinvariant behaviour. 66 9 OFC Earthquake Model 9.1 OFC Earthquake Equation The Olami-Feder-Christensen (OFC) Model [10] is a discrete-time cellular automaton designed to simulate earthquake dynamics using a self-organized criticality (SOC) framework. The model evolves according to the following partial difference equation: Etu=u−θ(u−1) + αX (i,j)∈V θ(Ei xEj yu−1) + G(t, x, y) where: •u=u(t, x, y)∈[0,1] represents the local shear stress (or load) at site (x, y) at time t,u:Z3→R. •Etu=u(t+ 1, x, y) denotes the stress at the next time step. •θ(x) is the Heaviside function, defined as: θ(x) = (1,if x≥0 0,otherwise •αP(i,j)∈Vθ(Ei xEj yu−1) represents the redistributed stress received from neighbouring sites (i, j)∈Vthat have exceeded the failure threshold. α is a constant. •G(t, x, y) is an external driving term, typically very small, which represents slow tectonic loading. In most implementations, G(t, x, y) = εat a randomly selected site and zero elsewhere, modeling slow, local buildup of stress. The stress uis usually constrained within the interval [0,1), with the threshold for failure set at u= 1. When u≥1, the site fails, redistributes stress to its neighbours, and resets or reduces its own stress, initiating potential cascades (earthquake avalanches) through the system. This model is analogous to the sandpile model but includes continuousvalued stress and dissipation, making it more suitable for simulating earthquakelike events in tectonic systems. 9.2 Observations and Phenomena The Olami–Feder–Christensen (OFC) earthquake model exhibits dynamics similar to the sandpile model. It demonstrates a power-law distribution in event sizes, where small avalanches occur frequently, while large avalanches are rare. This behaviour is a hallmark of self-organized criticality (SOC), indicating that the system naturally evolves toward a critical state without fine-tuning of parameters. 67 10 Forest Fire Model 10.1 Forest Fire Equation We define the forest fire dynamics on a discrete lattice as a Partial Difference Equation with mod-3 state cycling [21]: Etu= mod3(2 δ(u−1) θ(S−1) + G(t, x, y)) S=X (i,j)∈M δEi xEj yu−2 10.2 Definitions •u=u(t, x, y)∈ {0,1,2}is the state of the cell at time tand spatial coordinates (x, y). –0: empty ground (bare land) –1: tree –2: burning tree •Etis the forward time shift operator: Etu:= u(t+ 1, x, y) •Ei xEj yu:= u(t, x +i, y +j) are spatial shift operators in the Moore neighbourhood M. •δ(n) is the Kronecker delta function: δ(n) = 1 if n= 0, else 0. •θ(x) is the Heaviside step function: θ(x) = 1 if x≥0, else 0. •G(t, x, y)∈ {0,1}is the non-autonomous external driving force: –G= 0: no external change –G= 1: apply one-step growth/fire transition (see below) 10.3 Transition Logic (via mod 3) The entire system evolves via modular arithmetic: New state = mod3(current state + external or fire-induced increment) •0+1 −−→ 1: tree grows on empty land •1+1 −−→ 2: tree catches fire 68 •2+1 −−→ 0: burning tree becomes ash (empty) This simple cyclic structure captures the dynamics of vegetation growth, fire propagation, and destruction. 10.4 Interpretation A tree at (x, y) will ignite if: - Its current state is 1 (tree), and - There is at least one burning tree in its neighbourhood (S≥1) Otherwise, only external driving force G(t, x, y) determines its evolution. 10.5 Observations and Phenomena We consider the autonomous equation of the Forest Fire Model, and using Moore Neighbourhood. Numerical simulations reveal that the global behaviour of the system is highly sensitive to the initial condition, particularly the forest density. Assume the initial condition where the entire leftmost column of the grid is set on fire, and all other cells are either trees or empty ground according to a prescribed forest density. We observe the following distinct behaviours: •Low density (e.g., 10%): The fire fails to propagate. Due to the sparsity of trees, the flame quickly extinguishes after burning only a few connected trees. •High density (e.g., 80%): The fire spreads rapidly and extensively. Almost the entire forest is consumed in a short period, with the flame front advancing in a nearly uniform and predictable wave. •Intermediate density (e.g., around 40%): A highly nontrivial behaviour emerges. The fire propagates in an irregular, branching manner, forming a complex fractal-like structure. The fire path becomes winding, unpredictable, and sensitive to small variations in the initial configuration. This suggests the existence of a critical threshold of forest density that separates the extinguishing, fractal, and total-burning regimes. The system exhibits characteristics of critical phenomena and self-organized complexity at intermediate densities. 69 11 Kuramoto Firefly Model 11.1 Kuramoto Equation We define the Kuramoto model as a discrete phase-coupled oscillator network evolving on a 2D grid: Etu= mod2π u+w(x, y) + K 8X (i,j)∈M sin(Ei xEj yu−u)  •u=u(t, x, y)∈[0,2π] is the phase at grid point (x, y) and time t. •w(x, y) is a fixed intrinsic frequency assigned randomly to each point. •Kis the coupling constant. In our simulations, we choose K= 1. •Mis the Moore neighbourhood: M={(±1,0),(0,±1),(±1,±1)} 11.2 Principle of Locality In the original Kuramoto model, oscillators are globally coupled, each unit interacts with every other oscillator in the system. [20]. While this assumption simplifies mathematical analysis, it is biologically implausible. For example, it is unreasonable to assume that a firefly can perceive and compute the phases of all other fireflies in the environment at each time step. To address this limitation, we propose the Kuramoto Firefly Model, which enforces the Principle of Locality. In this model, each oscillator interacts only with its immediate neighbours (e.g., a Moore neighbourhood in two dimensions). This localized interaction is more consistent with real world biological systems and gives rise to rich emergent phenomena such as synchronization, topological defects, and vortex-like phase structures. 11.3 Observations and Phenomena In the simulation, we use periodic boundary condition. From the spatiotemporal plots of the system, we observe the following: •The system exhibits topological defects: points around which the phase (colour) continuously wraps from 0 to 2π. •These defects appear in pairs, with opposite rotational direction: –One clockwise –One counterclockwise •When these defects move and eventually annihilate each other, a large region of synchronization emerges suddenly. 70 This model shows interesting connections to: •Phase synchronization phenomena •Topological vortex dynamics •Pattern formation and self-organization Figure 13: Spatiotemporal evolution of the Kuramoto model. Topological defects (vortices) emerge and move, and their annihilation leads to large-scale synchronization. 71 12 Ising Model 12.1 Governing Equation We propose the following discrete evolution equation for the well-known Ising Model [19]: Etu= (1 −δ(sign(JSVu)−sign(u))) θ(F−ε) sign(JSVu) +δ(sign(JSVu)−sign(u)) θ(F−F0) (−sign(JSVu)) (1) where SVu=X (i,j)∈V Ei xEj yu, •Vis the von Neumann neighbourhood. •u=u(t, x, y)∈ {−1,+1}is the binary state at discrete time tand spatial coordinate (x, y). •Etu:= u(t+ 1, x, y) is the next time step value. •F=F(t, x, y)∼ N(µ, c2) is a random driving field at site (x, y), drawn from a normal distribution centered at µ=kT, where Tis the temperature and kis a constant. •εis a small positive number (baseline activation threshold). •F0is a larger threshold required for flipping a stable site. •Jis the coupling constant determining interaction strength between neighbouring sites. 12.2 Physical Interpretation The function u(t, x, y) represents the spin (or state) of a particle or site at position (x, y) and time t. The variable Etudenotes the updated state at time t+1. The update is governed by local interactions and stochastic environmental drive F. •The first term activates when the site’s current state is different from the local neighbourhood majority (sign(JSVu)= sign(u)), allowing it to flip easily when F > ε. •The second term activates when the site’s current state matches the neighbourhood majority, but may still flip if the driving force exceeds a higher threshold F0, modeling thermal noise or instability. This model expresses the Ising-like behaviour in the language of partial difference equations (P∆E), linking local deterministic update rules with stochastic driving force under thermal control. 72 13 Greenberg–Hastings Model 13.1 Greenberg-Hastings Equation We propose the following formulation of the Greenberg–Hastings model [23] as a partial difference equation (P∆E): Etu=θ(u−1) ·modn(u+ 1) + δ(u)·θ(S−b) where S=X (i,j)∈N δEi xEj yu−1 with the following definitions: •u=u(t, x, y)∈Z, a discrete function defined on space-time: u:Z3→Z. •θ(x) = (1, x ≥0 0, x < 0is the Heaviside step function. •δ(x) = (1, x = 0 0, x = 0 is the Kronecker delta function. •modn(x) = x−n·x nis the modulo-noperation. •n∈Nis the number of discrete states (typically n≥3). •b∈Nis the excitation threshold. •Nis the set of relative spatial coordinates in the neighbourhood (e.g., von Neumann or Moore). 13.2 Explanation The Greenberg–Hastings model is a classical excitable medium model capable of producing complex spatiotemporal patterns such as spiral waves. Traditionally, it is expressed using state variables xn i,j updated via nested if-else rules. In this work, we reformulate the model as a partial difference equation (P∆E) using evolution operators and logical functions. This approach provides a more elegant, symbolic, and generalizable representation that integrates naturally into our proposed framework of discrete dynamical systems. 73 15.7 Discrete and Continuous Green’s Function For simplicity, we only consider the one-dimensional case. Let u=u(t, x), where u:R2→R, and consider a linear partial differential equation (PDE) of the form L(u)=0, where Lis a linear operator. Suppose the initial condition is given by u(0, x) = δ(x), where δ(x) is the Dirac delta function, and assume there are no boundary conditions. If the solution to this PDE is denoted by G(t, x), the Green’s function [15], then the solution corresponding to an arbitrary initial condition u(0, x) = f(x) is given by the convolution: u(t, x) = Z∞ −∞ f(s)G(t, x −s)ds. Now, let u=u(t, x), where u:Z2→R, and consider a linear partial difference equation (P∆E) of the form L(u) = 0, where Lis again a linear operator. Suppose the initial condition is u(0, x) = δ(x), where δ(x) is now the Kronecker delta function, and assume there are no boundary conditions. If the solution is denoted by G(t, x), the discrete Green’s function, then for an arbitrary initial condition u(0, x) = f(x), the solution is given by the discrete convolution: u(t, x) = X s∈Z f(s)G(t, x −s). In the continuous setting, the Green’s function representation of the solution involves an integral over the spatial domain. In the discrete setting, the corresponding representation is given by a summation. These two formulations are structurally analogous, differing only in the underlying measure, Lebesgue measure for the continuous case, and counting measure for the discrete case. 15.8 Gaussian Distribution and Combinatorics In the continuous heat equation, the equation is given by ∂u ∂t =α∂2u ∂x2. Consider the initial condition u(0, x) = δ(x), 80 where δ(x) is the Dirac delta function. With no boundary condition, the solution is u(t, x) = 1 √4παt exp −x2 4αt, which is the Gaussian distribution. In the discrete heat equation, the equation becomes ∆tu=α∇2u, and choosing α=1 2, we write the equation as Etu=1 2Exu+E−1 xu. Consider the initial condition u(0, x) = δ(x), where δ(x) is the Kronecker delta function, and again, no boundary condition is imposed. Then the solution is given by u(t, x) = δ(mod2(x+t)) ·Ct, x+t 2·1 2t , where C(n, k) =    n! k!(n−k)! if 0 ≤k≤n, 0 otherwise. We can see that the Gaussian distribution is closely related to combinatorics. In fact, the Gaussian distribution emerges as the continuous limit of the binomial distribution as the number of trials tends to infinity, with appropriate scaling and centering, a consequence of the Central Limit Theorem. 15.9 Discrete and Continuous Dynamical Systems In the framework of topological dynamics, a dynamical system is defined as a pair (X, φt), where Xis a topological space called the state space, and φt:X→Xis a family of time-evolution operators indexed by t∈T, where T=R(continuous time) or T=Z(discrete time). The map φtsatisfies the semigroup property φt+s=φt◦φsand φ0= id [4]. Continuous Dynamical Systems. Partial Differential Equations (PDE) define continuous dynamical systems, typically of the form: u:Rn→C,e.g., u(t, x, y, z),with t∈R,(x, y, z)∈Ω⊆R3. 81 The state space Xis a continuous function space, such as L2(Ω), Ck(Ω), or Sobolev spaces Hk(Ω). The evolution of the system is governed by differential operators that act on these function spaces. In such systems, chaotic behaviour is rare and usually manifests as smooth strange attractors—continuous curves resembling vortices or spirals. The attractors tend to be geometrically regular and topologically simple. Discrete Dynamical Systems. Partial Difference Equations (P∆E) define discrete dynamical systems, where u:Zn→C,e.g., u(t, x, y, z),with t∈Z,(x, y, z)∈Ω⊆Z3. The state space Xis a discrete function space, i.e., a space of functions defined on a lattice, often taking values in C,Z, or even a finite set (such as Boolean values), for example, discrete Hilbert Space: L2(Ω). Such systems are capable of generating highly complex spatiotemporal patterns and exhibit rich chaotic behaviours. The associated strange attractors can have fractal geometries far more exotic than in continuous systems: they may resemble turbulent flows, ribbons, spider webs, tunnels, wings, bottles, or other intricate structures. Fractals and Complexity. Fractals are rarely seen in PDE due to the smoothness constraints imposed by the function space and differential operators. In contrast, P∆E can easily generate fractal structures. For example: •Rule 90 cellular automaton generates the Sierpi´nski triangle. •The Abelian Sandpile Model exhibits discrete self-organized criticality with fractal avalanche clusters. •Coupled Map Lattices generate chaotic spatial patterns with complex attractors. These observations suggest that discrete dynamical systems possess a fundamentally different character compared to continuous ones: their state spaces are inherently combinatorial, their attractors more irregular, and their emergent behaviours often more computationally universal and structurally complex. 15.10 Manifolds and Networks Partial Differential Equations (PDE) are traditionally defined on smooth manifolds [16], such as curves, surfaces, or general Riemannian manifolds. The underlying space is continuous, and thus the study of PDE often requires tools from differential geometry. PDE can be formulated in various coordinate systems, including rectangular, cylindrical, and spherical coordinates, depending on the geometry of the manifold. In contrast, we observe that Partial Difference Equations (P∆E) need not be confined to rectangular grids. They can also be defined on triangular lattices, 82 hexagonal lattices, and other non-rectangular discretizations [5]. Furthermore, P∆E can be formulated on polyhedral surfaces such as Goldberg polyhedra, zonohedra, and even higher-dimensional polytopes [7]. In this case, the underlying structure is a discrete network, and the evolution takes place over this network. The study of such discrete geometric networks naturally connects to graph theory and discrete differential geometry. Moreover, P∆E exhibit deep relationships with discrete and computational geometry. Although PDE and P∆E belong to continuous and discrete domains respectively, the geometric structures underlying them exhibit a fascinating duality. This suggests a possible formal correspondence between smooth manifolds and discrete networks, and opens a path to unify the continuous and discrete geometric theories under a common framework. 15.11 Boolean Algebra, Theoretical Computer Science and Number Theory In certain classes of Partial Difference Equations (P∆E), such as those arising from elementary cellular automata, the state space is restricted to two values, typically {0,1}. This naturally connects such systems to Boolean algebra, where logical operations govern the update rules. The evolution of these discrete systems is thus governed by logical expressions, and the resulting dynamics can often be described entirely in terms of Boolean functions. Moreover, the well-known Conway’s Game of Life, which can also be viewed as a particular P∆E defined over a two-dimensional grid with binary states, is known to be Turing complete [17]. This highlights a deep and important connection between P∆E and theoretical computer science [18], as such systems are capable of universal computation. Cellular automata and related P∆E systems therefore serve as bridges between dynamical systems, logic, and computation. Additionally, difference equations have strong ties to number theory [18]. For example, the Fibonacci sequence satisfies the recurrence relation x(t+ 1) = x(t) + x(t−1), which is a second-order linear ordinary difference equation. More generally, integer-valued solutions to linear P∆E often exhibit rich combinatorial and arithmetic structures. In many cases, the solutions involve binomial coefficients, discrete convolutions, or integer sequences catalogued in number theory. The fact that P∆E naturally accommodate integer-valued functions u(t, x)∈Z reinforces their relevance to number-theoretic investigations. Hence, P∆E lie at the intersection of Boolean algebra, number theory, and computational theory, offering a unified framework for exploring discrete structures and their evolution. 83 15.12 Summary We observe that although partial differential equations (PDE) and partial difference equations (P∆E) differ fundamentally in nature, one being continuous and the other discrete, they exhibit profound structural similarities, often appearing as duals of one another. Moreover, we have seen that P∆E are deeply interconnected with various domains of discrete mathematics, including combinatorics, number theory, graph theory, and discrete and computational geometry. Partial Difference Equations are not merely exotic multivariable recurrence relations; rather, they serve as a unifying framework that bridges numerous areas in the discrete mathematical landscape. 84 16 Summary 16.1 From Numerical Methods to Language of Complexity We began this work by introducing Partial Difference Equations (P∆E) as a formal mathematical framework for modeling systems that are discrete in both time and space. Starting from their motivation and definitions, we established their relationship to ordinary difference equations (O∆E) and to their continuous counterparts, partial differential equations (PDE), thereby clarifying the scope and intent of this study. We then examined the linear theory, presenting the main types of linear P∆E and illustrating simple analytic solutions without pursuing a full solution theory. This led naturally to the development of discrete functional analysis tools: discrete function spaces, Hilbert space structures, and discrete operators, including shift and difference operators, along with their adjoint properties. Within this setting, we defined the discrete Green’s function and demonstrated how it parallels its continuous analogue. The unifying strength of the P∆E formalism was then illustrated by reformulating a broad spectrum of well-known discrete models, ranging from cellular automata and coupled map lattices to sandpile models, the Ising model, and beyond, into explicit P∆E form. These examples demonstrate that systems traditionally studied in isolation can be expressed within a single coherent mathematical language. Finally, we compared the “discrete universe” of P∆E with the continuous universe of PDE, revealing deep structural parallels while emphasizing their distinct mathematical foundations. In this view, P∆E and PDE emerge as mathematical twins: discrete–continuous analogues that share a common formal architecture yet belong to fundamentally different domains of mathematics. From their origins as tools in numerical methods, difference equations have evolved here into a general language of complexity, capable of describing, unifying, and extending a wide variety of models in physics, biology, and computation. This perspective not only enriches the theoretical foundations of discrete dynamical systems but also opens pathways for future research across the sciences. 85 References [1] William A. Adkins and Mark G. Davidson. Ordinary Differential Equations. Undergraduate Texts in Mathematics. Springer, 2012. [2] Per Bak, Chao Tang, and Kurt Wiesenfeld. Self-organized criticality: An explanation of 1/f noise. Physical Review Letters, 59(4):381–384, 1987. [3] Jean Pierre Boon. Langton’s ant as an elementary turing machine. In Nonequilibrium Thermodynamics and Fluctuation Kinetics, volume 208 of Fundamental Theories of Physics, pages 135–140. Springer, Cham, 2022. [4] Michael Brin and Garrett Stuck. Introduction to Dynamical Systems. Cambridge University Press, 1st edition, 2002. Graduate-level introduction to dynamical systems theory; Cambridge, 1st ed. [5] Satyan L. Devadoss and Joseph O’Rourke. Discrete and Computational Geometry. Princeton University Press, 1st edition, 2011. [6] Saber Elaydi. An Introduction to Difference Equations. Undergraduate Texts in Mathematics. Springer, New York, 3rd edition, 2005. [7] Jacob E. Goodman, Joseph O’Rourke, and Csaba D. T´oth, editors. Handbook of Discrete and Computational Geometry. Discrete Mathematics and Its Applications. Chapman & Hall/CRC, Boca Raton, FL, 3rd edition, 2017. [8] Makoto Itoh and Leon O. Chua. Difference equations for cellular automata. International Journal of Bifurcation and Chaos, 19(3):805–830, 2009. [9] Nathaniel Johnston and Dave Greene. Conway’s Game of Life: Mathematics and Construction. Self-published, 2022. Hardcover, colour printing; US letter size (8.5 ×11 in). [10] Bin-Quan Li and Sheng-Jun Wang. Self-organized criticality in an anisotropic earthquake model. Communications in Theoretical Physics, 69(3):280–284, 2018. [11] Carlo Mariconda and Alberto Tonolo. Discrete Calculus: Methods for Counting, volume 103 of UNITEXT - La Matematica per il 3+2. Springer International Publishing, Switzerland, 2016. [12] J.D. Murray. Mathematical Biology: I. An Introduction, volume 17 of Interdisciplinary Applied Mathematics. Springer-Verlag, New York, Berlin, Heidelberg, third edition, 2002. Printed in the United States of America. [13] J.D. Murray. Mathematical Biology II: Spatial Models and Biomedical Applications, volume 18 of Interdisciplinary Applied Mathematics. SpringerVerlag, New York, Berlin, Heidelberg, third edition, 2003. Printed in the United States of America. 86 [14] R.˜ N. Mutagi. Understanding the discrete fourier transform. PDF lecture notes, January 2004. Uploaded by the author on ResearchGate, Indus University. [15] Peter J. Olver. Introduction to Partial Differential Equations. Undergraduate Texts in Mathematics. Springer, Cham, Heidelberg, New York, Dordrecht, London, 2014. [16] Peter Petersen. Riemannian Geometry, volume 171 of Graduate Texts in Mathematics. Springer, Cham, 3rd edition, 2016. [17] Paul Rendell. Game of life turing machine. In Andrew Adamatzky and Genaro J. Mart´ınez, editors, Turing Machine Universality of the Game of Life, volume 18 of Emergence, Complexity and Computation, pages 45–70. Springer International Publishing, Cham, 2016. [18] Kenneth H. Rosen, editor. Handbook of Discrete and Combinatorial Mathematics. Discrete Mathematics and Its Applications. Chapman & Hall/CRC, 2nd edition, 2017. [19] Satya Pal Singh. The ising model: Brief introduction and its application. In Solid State Physics. IntechOpen, 2020. CC BY 3.0 license. [20] Steven H. Strogatz. From kuramoto to crawford: Exploring the onset of synchronization in populations of coupled oscillators. Physica D: Nonlinear Phenomena, 143(1–4):1–20, 2000. [21] Xuan Sun, Ning Li, Duoqi Chen, Guang Chen, Changjun Sun, Mulin Shi, Xuehong Gao, Kuo Wang, and Ibrahim M. Hezam. A forest fire prediction model based on cellular automata and machine learning. IEEE Access, 12:55389–55403, 2024. CC BY-NC-ND 4.0. [22] Stephen Wolfram. A New Kind of Science. Wolfram Media, Champaign, Illinois, 1st edition, 2002. [23] An-Cai Wu, Xin-Jian Xu, and Ying-Hai Wang. Excitable greenberghastings cellular automaton model on scale-free networks. arXiv preprint cond-mat/0701248, 2007. arXiv:cond-mat/0701248, submitted 11 Jan 2007. [24] Tatsuo Yanagita and Kunihiko Kaneko. Coupled map lattice model for convection. Physics Letters A, 175(6):415–420, 1993. 87