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

On the Theory of Partial Difference Equations: From Numerical Methods to the Language of Complexity Bik Kuang Min Department of Mathematical Sciences, Faculty of Science and Technology, Universiti Kebangsaan Malaysia (UKM), Bangi, Selangor, Malaysia Email: [email protected] ORCID: 0009-0009-8368-3858 12 September 2025 1 Abstract This monograph develops a comprehensive and unified theory of linear and nonlinear partial difference equations (P∆E), extending them far beyond their traditional role in numerical analysis. We establish a rigorous mathematical framework that connects discrete analysis, combinatorics, operator theory, Fourier methods, functional analysis, and the theory of dynamical systems. We begin by introducing the foundations of P∆Es: discrete functions and shifts, the classification of equations (linear, semilinear, quasilinear, fully nonlinear), and the basic regularity properties of discrete solutions. A full operator–theoretic setting is developed through discrete function spaces (Banach, Hilbert, Lp, Schwartz, and finite–support spaces), adjoint theory, and spectral properties of discrete Laplacians—including the Moore Laplacian in higher dimensions. Fundamental tools of functional analysis such as the Hahn–Banach theorem, the Riesz representation theorem, compactness principles, and Banach fixed point theory are established in the discrete setting. Green’s functions are constructed systematically, showing how classical combinatorial structures—binomial, multinomial, and Stirling numbers—arise naturally as fundamental solutions of P∆Es. Discrete Fourier series and transforms are developed in detail, including Parseval’s identity and the use of discrete symbols (Trigonometric polynomials) to analyze well–posedness and classify second–order equations. Explicit classes of equations are solved using discrete analogues of separation of variables, Fourier integrals, and semigroup methods. These include first–order evolution equations, discrete heat and wave equations, and a wide family of higher–dimensional linear models whose solutions are expressed through multinomial evolution kernels. Nonlinear P∆Es are then introduced, with particular attention to mod-nnonlinearities. These equations produce exact fractal solutions such as the Sierpinski triangle, carpet, and pyramid, and lead to the conceptual proposal that fractals can be regarded as solutions of evolution equations. The theory is completed with a systematic treatment of systems of P∆Es, written in evolution form , and with an atlas of nonlinear models including cellular automata, sandpile dynamics, discrete diffusion–reaction systems, and discrete analogues of fluid equations. Overall, this monograph reframes partial difference equations as a fundamental language for discrete dynamics, complex systems, and fractal evolution, advancing them from their numerical origins to a broad mathematical theory capable of describing self–organization and complexity. Keywords: partial difference equations; difference equations; Green’s function; functional analysis; discrete dynamical systems; Fourier transform; combinatorics; fractals 2 Contents Preface 7 1 Introduction to Partial Difference Equations 8 1.1 Motivation and Background . . . . . . . . . . . . . . . . . . . . . 8 1.2 Relation to Existing Literature . . . . . . . . . . . . . . . . . . . 9 1.3 DiscreteFunction........................... 10 1.4 Operators............................... 10 1.5 Ordinary Difference Equations . . . . . . . . . . . . . . . . . . . 11 1.6 Partial Difference Equations . . . . . . . . . . . . . . . . . . . . . 12 1.7 The Order of a Partial Difference Equation . . . . . . . . . . . . 12 1.8 Classification and Examples . . . . . . . . . . . . . . . . . . . . . 13 1.9 Notations ............................... 14 1.10 Notation Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2 Introduction to Linear Partial Difference Equations 17 2.1 Definition ............................... 17 2.2 1D Discrete Transport Equation . . . . . . . . . . . . . . . . . . 18 2.3 3D Discrete Transport Equation . . . . . . . . . . . . . . . . . . 21 3 Discrete Function Spaces and Operators 22 3.1 Introduction.............................. 22 3.2 Discrete Function Spaces . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 TypesofOperators.......................... 27 3.4 OperatorTheory ........................... 32 4 Discrete Functionals and Convergence 41 4.1 Introduction.............................. 41 4.2 Discrete Functionals . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.3 Compactness of the Unit Ball in Discrete Function Spaces . . . . 42 4.4 Hahn–Banach Theorem . . . . . . . . . . . . . . . . . . . . . . . 43 4.5 Riesz Representation Theorem . . . . . . . . . . . . . . . . . . . 44 4.6 Banach Fixed Point Theorem . . . . . . . . . . . . . . . . . . . . 45 4.7 Types of Convergence . . . . . . . . . . . . . . . . . . . . . . . . 45 5 The Discrete Green’s Function 47 5.1 The Kronecker Delta Function . . . . . . . . . . . . . . . . . . . 47 5.2 Discrete Convolution . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.3 BasicProperties ........................... 48 5.4 Fundamental Solution . . . . . . . . . . . . . . . . . . . . . . . . 50 5.5 The Discrete Green’s Function . . . . . . . . . . . . . . . . . . . 50 3 6 Discrete Fourier Transform 52 6.1 Motivation and Background . . . . . . . . . . . . . . . . . . . . . 52 6.2 Discrete Fourier Series . . . . . . . . . . . . . . . . . . . . . . . . 52 6.3 Inner Product and Orthogonality . . . . . . . . . . . . . . . . . . 52 6.4 Discrete Fourier Transform . . . . . . . . . . . . . . . . . . . . . 53 6.5 Parseval’s Identity . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.6 Fourier Multiplier Operators . . . . . . . . . . . . . . . . . . . . 55 6.7 FourierAnsatz ............................ 55 7 Symbolic Analysis and Classification 57 7.1 Introduction.............................. 57 7.2 CauchyProblem ........................... 57 7.3 PrincipalSymbol........................... 59 7.4 Classification of Linear P∆E . . . . . . . . . . . . . . . . . . . . 61 7.5 Trigonometric Polynomial and Quadric Surface . . . . . . . . . . 64 7.6 Hadamard’s Well-Posedness . . . . . . . . . . . . . . . . . . . . . 69 7.7 Stability of Partial Difference Equations . . . . . . . . . . . . . . 69 8 First Order Equations in Time 73 8.1 Introduction.............................. 73 8.2 1D Pascal Evolution Equation . . . . . . . . . . . . . . . . . . . . 73 8.3 1D Discrete Heat Equation . . . . . . . . . . . . . . . . . . . . . 76 8.4 Stirling Second Equation . . . . . . . . . . . . . . . . . . . . . . . 79 9 Second Order Equations in Time 81 9.1 Introduction.............................. 81 9.2 Second–Order Pascal Evolution Equation . . . . . . . . . . . . . 81 9.3 Discrete 1D Wave Equation . . . . . . . . . . . . . . . . . . . . . 84 10 Steady State Problems 87 10.1Introduction.............................. 87 10.2 2D Discrete Laplace Equation . . . . . . . . . . . . . . . . . . . . 88 10.3 2D Discrete Poisson Equation . . . . . . . . . . . . . . . . . . . . 90 10.4 Two-Dimensional Moore Laplace Equation . . . . . . . . . . . . . 92 10.5 List of Discrete Evolution Equations . . . . . . . . . . . . . . . . 95 10.6 List of Steady State Problems . . . . . . . . . . . . . . . . . . . . 97 11 Discrete Evolution Equations 99 11.1Introduction.............................. 99 11.2Definitions............................... 99 11.3 Semigroup Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 103 11.4 Initial Value Problems . . . . . . . . . . . . . . . . . . . . . . . . 104 11.5 Boundary Value Problems . . . . . . . . . . . . . . . . . . . . . . 105 11.6 Initial-Boundary Value Problems . . . . . . . . . . . . . . . . . . 107 11.7 Autonomous and Non–Autonomous Systems . . . . . . . . . . . . 108 11.8 Spatial Homogeneity and Non–Homogeneity . . . . . . . . . . . . 109 4 12 Higher Dimensional Problems 110 12.1Introduction..............................110 12.2 2D Pascal Evolution Equation . . . . . . . . . . . . . . . . . . . . 110 12.3 3D Pascal Evolution Equation . . . . . . . . . . . . . . . . . . . . 111 12.4 n-Dimensional Pascal Evolution Equation . . . . . . . . . . . . . 113 13 Nonlinear Equations with mod n Nonlinearity 115 13.1Introduction..............................115 13.2 Right–Side Sierpinski Triangle Equation . . . . . . . . . . . . . . 115 13.3 Sierpinski Carpet Equation . . . . . . . . . . . . . . . . . . . . . 119 13.4 Sierpinski Pyramid Equation . . . . . . . . . . . . . . . . . . . . 120 13.5 Sierpinski Fan Equation . . . . . . . . . . . . . . . . . . . . . . . 125 13.6 Mod 5 Sierpi´nski Equation . . . . . . . . . . . . . . . . . . . . . . 126 13.7 Mod 4 Sierpinski Fractal Equation . . . . . . . . . . . . . . . . . 127 13.8 Fractals as Solutions to Evolution Equations . . . . . . . . . . . 129 14 Oscillation Theory 130 14.1Motivation ..............................130 14.2 Brillouin Zone and Nyquist Frequency . . . . . . . . . . . . . . . 130 14.3FourierAnalysis ...........................131 14.4WaveletAnalysis...........................131 14.5 Discrete Conservation Law . . . . . . . . . . . . . . . . . . . . . . 131 15 System of Partial Difference Equations 132 15.1Introduction..............................132 15.2 Classification of P∆E Systems . . . . . . . . . . . . . . . . . . . 132 15.3LinearSystems ............................136 15.4 Constant Coefficient Linear P∆E System . . . . . . . . . . . . . 139 15.5 Non-homogeneous Linear System . . . . . . . . . . . . . . . . . . 142 15.6 Semilinear Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 143 16 Atlas of Nonlinear Partial Difference Equations 145 16.1Introduction..............................145 16.2 Coupled Map Lattice . . . . . . . . . . . . . . . . . . . . . . . . . 146 16.3 Elementary Cellular Automata . . . . . . . . . . . . . . . . . . . 148 16.4 The Conway’s Game of Life . . . . . . . . . . . . . . . . . . . . . 160 16.5 Abelian Sandpile Model . . . . . . . . . . . . . . . . . . . . . . . 164 16.6ForestFireModel...........................167 16.7 The Langton’s Ant . . . . . . . . . . . . . . . . . . . . . . . . . . 172 16.8 Kuramoto Firefly Model . . . . . . . . . . . . . . . . . . . . . . . 174 16.9IsingModel..............................176 16.10Discrete Logistic Diffusion Equation . . . . . . . . . . . . . . . . 177 16.11Bacterial Colony Model . . . . . . . . . . . . . . . . . . . . . . . 186 16.12Discrete Navier–Stokes Equations . . . . . . . . . . . . . . . . . . 192 5 17 Summary 196 17.1 From Numerical Methods to the Language of Complexity . . . . 196 17.2 A Glance at the Discrete Field Theory . . . . . . . . . . . . . . . 197 17.3 Limitations and Future Work . . . . . . . . . . . . . . . . . . . . 197 Acknowledgements 198 References 199 6 Preface Since entering university, I have been deeply interested in modeling biological evolution through mathematics. My initial attempts focused on graph theory, but I soon encountered limitations—for instance, classical graph-theoretic approaches could not adequately capture phenomena such as ring species. I then turned to partial differential equations (PDE), hoping to describe continuous processes such as cellular division, but again found that PDEs were insufficient to model discrete structural transformations such as surface splitting. These failures gradually led me to a new perspective. Biological systems are, at their core, discrete in both space and time: molecules, cells, and individuals are discrete entities, and reproduction occurs across generations. This suggested that differential equations might not be the most natural framework for such systems. Instead, one should turn to difference equations. At first, I was unsure how to formulate partial difference equations in a systematic way. The breakthrough came when I recognized that classical models from complex systems science—such as the Game of Life and the Sandpile Model—can, in fact, be naturally expressed as partial difference equations. These models share the same essential structure: dynamics defined on discrete space and discrete time. From this perspective, just as PDEs describe the evolution of continuous fields, partial difference equations can be viewed as describing the evolution of discrete fields. This realization motivated me to systematically reinterpret a wide variety of discrete models—including cellular automata, combinatorial recursions, and evolutionary dynamics—within a single analytic framework that I call Discrete Field Theory. The aim of this framework is not merely to provide another modeling technique, but to unify diverse phenomena under a language inspired by the analysis of PDE, enriched by tools from functional analysis, operator theory, and spectral theory. This work represents only a first step. The theory is far from complete, and much remains to be developed and refined. I hope that by presenting these ideas, I can invite further exploration and collaboration. My goal is to provide not a finished solution, but rather a starting point for a broader research program in the analysis of discrete dynamical systems. 7 1 Introduction to Partial Difference Equations 1.1 Motivation and Background Over the past year, the author has released several preprints [4–6] aiming to describe a wide range of complex systems, including cellular automata, sandpile models, and other discrete dynamical processes, within a single unified framework based on partial difference equations (P∆Es). These works demonstrated that many classical models of complexity, often studied independently in computer science, combinatorics, physics, or nonlinear dynamics, may in fact be expressed as explicit evolution equations built from shift and difference operators. However, despite these applications, a complete and rigorous mathematical foundation for such equations was still missing. The present monograph is therefore devoted to establishing a comprehensive theory of partial difference equations: including their analytic framework, function spaces, operator theory, Fourier analysis, Green’s functions, classification, well-posedness, and methods of solution. Our goal is to create a systematic and self-contained treatment comparable in scope and structure to the classical theory of partial differential equations. The author initially arrived at this subject independently. By analyzing wellknown cellular automata such as Conway’s Game of Life, the Abelian Sandpile Model, and various coupled-map lattice systems, it became evident that each could be written naturally as a partial difference equation in space and time. Searching the literature for “Partial Difference Equations” suggested that the term existed, but almost exclusively in the numerical analysis community, where it referred to finite-difference schemes for approximating PDEs. This was entirely different from the higher-dimensional intrinsic difference equations considered here. For a period of time the author believed that this viewpoint—treating difference equations as genuine multivariable analytic objects, rather than as numerical discretizations, might be entirely new. It was only later, after receiving correspondence from Professor Alexander Lyapin (Siberian Federal University), that the author learned of an existing, though extremely small, research community working on what they call multidimensional difference equations. The approach developed by Lyapin, Leinartas, Apanovich and others is rooted in combinatorics and complex-analytic methods [1,23,25], frequently involving generating functions and functional equations. Although these works share similar motivations, they differ substantially in notation, terminology, and methodology. In this monograph we continue to employ the author’s terminology: •Ordinary Difference Equations (O∆E) refer to one-dimensional discrete equations. •Partial Difference Equations (P∆E) refer to genuinely multivariable discrete equations, defined on Zn. 8 Our notation, operator calculus, and analytical style follow more closely the tradition of partial differential equations: we work extensively with shift operators, adjoint operators, discrete analogues of divergence and Laplacian, spectral theory, eigenfunction expansions, Fourier transforms, Green’s functions, and semigroup methods. Thus, while related to the combinatorial literature, the present monograph develops a fundamentally different and more PDE-oriented viewpoint. The purpose of this volume is to present, for the first time, a unified and systematically developed theory of partial difference equations as an independent mathematical discipline, standing alongside the classical theory of partial differential equations and discrete dynamical systems. 1.2 Relation to Existing Literature Although the title of the present volume coincides with that of Partial Difference Equations by Sui Sun Cheng [11], the mathematical content and methodological framework of the two works are fundamentally distinct. Cheng’s monograph approaches discrete equations primarily from the perspective of classical recurrence relations, combinatorial structures, and finitedifference schemes arising in the numerical analysis of partial differential equations. His treatment is rooted in discrete combinatorics and algorithmic manipulations of recurrence formulas, with an emphasis on explicit computations, special identities, and PDE-inspired discretizations. In contrast, the present monograph develops a self-contained, operatortheoretic and functional-analytic foundation for the theory of partial difference equations (P∆Es). Here, discrete dynamical systems are formulated entirely in terms of shift operators acting on Banach and Hilbert spaces of discrete functions, and the resulting P∆Es are studied as genuine operator equations: F{Ek1 x1···Ekn xnu}(k1,...,kn)∈S, x= 0, u :Zn→C. This functional-analytic setting enables the development of symbolic calculus, spectral theory, Green’s functions, semigroup methods, and Fourier analysis for discrete operators in a manner directly parallel to the modern theory of partial differential equations. Thus, while Cheng’s book belongs to the combinatorial and numerical tradition of discrete equations, the present work aims to establish an operatortheoretic and structural framework for P∆Es, elevating them to the level of a general mathematical theory analogous to that of classical PDEs. The two approaches are therefore complementary but address fundamentally different questions and methodologies within the broader study of discrete dynamical systems. 9 1.10 Notation Evolution In the work by Itoh and Chua [19], 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). 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. 16 2 Introduction to Linear Partial Difference Equations 2.1 Definition In this section we introduce the notion of linear partial difference equations. Definition 2.1 (Linear Partial Difference Equation).Let u: Ω ⊆Zn→R(or C) be a scalar function on the lattice. A linear partial difference equation is an equation of the form X k∈S ak(x)Eku(x) = f(x), where •x= (x1, . . . , xn)∈Zn, •k= (k1, . . . , kn)∈Zn, •Ek=Ek1 x1···Ekn xnare shift operators, •S⊂Znis a finite index set, •ak(x)are given coefficient functions, •f(x)is a prescribed source term. Remark 2.2. •If all coefficients ak(x)are constants, the equation is called a linear constantcoefficient partial difference equation. •If f(x)= 0, the equation is called non-homogeneous. If f(x)=0, it is called homogeneous. 17 2.2 1D Discrete 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. 18 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. Example 2.1 (Discrete Transport Equation).Consider the discrete transport equation Etu=E−1 xu, defined on the discrete spacetime lattice (t, x)∈N0×Z. Initial condition. We prescribe the initial profile u(0, x) = cosπx 60 . Boundary conditions. We impose absorbing (Dirichlet) boundary conditions at the spatial boundaries: u(t, −30) = u(t, 30) = 0, t ≥0. 19 Analytical solution. In the absence of boundary effects, the solution of the discrete transport equation is given by u(t, x) = f(x−t), where f(x) = u(0, x)is the initial profile. Thus, the initial data is transported rigidly to the right with unit speed along discrete characteristic lines. Numerical illustration. Figure 1 shows the spacetime evolution of the solution under the above initial and boundary conditions. Figure 1: Spacetime diagram of the discrete transport equation. The figure clearly exhibits straight characteristic lines, confirming that the discrete transport equation propagates the initial profile rigidly with unit speed, in agreement with the analytical solution u(t, x) = f(x−t) under absorbing boundary conditions. 20 2.3 3D Discrete Transport Equation We consider the three-dimensional discrete transport equation Etu=Ea xEb yEc zu, a, b, c ∈Z, for a function u:N0×Z3→R, u =u(t, x, y, z), together with the initial condition u(0, x, y, z) = f(x, y, z), and no boundary conditions. Solution Introduce the new coordinates ξ=x+at, η =y+bt, ω =z+ct, and define the transformed function v(t, ξ, η, ω) := u(t, x, y, z). We compute the time-shift of v: Etv=v(t+ 1, ξ, η, ω) =vt+ 1, x +at, y +bt, z +ct =vt+ 1, x +at +a, y +bt +b, z +ct +c =v(t+ 1, ξ +a, η +b, ω +c) =EtEa ξEb ηEc ωv =EtEa xEb yEc zu Thus, Etv=EtEa xEb yEc zu, v=Ea xEb yEc zu, and substituting u(t, x, y, z) = v(t, ξ, η, ω) gives Etv=v. Thus vsatisfies the ordinary difference equation Etv=v, whose general solution is v(t, ξ, η, ω) = f(ξ, η, ω). Returning to the original variables, we obtain the solution of the 3D discrete transport equation: u(t, x, y, z) = f(x+at, y +bt, z +ct). 21 3 Discrete Function Spaces and Operators 3.1 Introduction In this section, we introduce and rigorously define several discrete function spaces, such as the discrete Lpspaces, discrete Hilbert spaces, and the discrete Schwartz space. We then proceed to define a variety of operators in the discrete setting, including the shift operator, the difference operator, and the discrete Laplacian. Afterward, we discuss operator theory in this context, covering fundamental notions such as bounded linear operators and adjoints. In this work, we adapt a number of classical theorems from functional analysis (such as the Hahn–Banach theorem, the Riesz representation theorem, and the Banach–Alaoglu theorem) to the discrete framework over Zn. We then illustrate these results through explicit constructions in discrete function spaces and operator theory, thereby establishing an analytic foundation for the study of partial difference equations. Most of these concepts are natural adaptations of existing frameworks from the literature, drawing inspiration from classical references in difference equations, partial differential equations, and functional analysis [8,13,22,29]. 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 (Topological Vector Space).Atopological vector space (TVS) over a field K(where K=Ror C) is a vector space Xtogether with a topology τon Xsuch that: 1. (X, τ)is a topological space. 2. The vector addition map + : X×X→X, (x, y)7→ x+y is continuous with respect to the product topology on X×X. 3. The scalar multiplication map ·:K×X→X, (λ, x)7→ λx is continuous with respect to the product topology on K×X. Definition 3.3 (Normed Vector Space).Let Vbe a vector space over the field Ror C. A norm on Vis a mapping ∥·∥:V→[0,∞) satisfying, for all u, v ∈Vand all scalars α: 22 1. Positive definiteness: ∥u∥= 0 ⇐⇒ u= 0. 2. Homogeneity: ∥αu∥=|α|∥u∥. 3. Triangle inequality: ∥u+v∥ ≤ ∥u∥+∥v∥. A vector space Vequipped with a norm ∥·∥is called a normed vector space, denoted (V, ∥·∥). Definition 3.4 (Banach Space).Let (X, ∥·∥)be a normed vector space over R or C. We say that (X, ∥·∥)is a Banach space if Xis complete with respect to the metric induced by the norm, i.e. d(x, y) = ∥x−y∥, x, y ∈X, meaning that every Cauchy sequence {xn}in Xconverges to some limit x∈X. Definition 3.5 (Measure Space).Ameasure space is a triple (X, A, µ)where 1. Xis a nonempty set (called the underlying set). 2. Ais a σ-algebra of subsets of X; that is: (a) ∅∈ A, (b) If A∈ A, then X\A∈ A, (c) If {Ai}∞ i=1 ⊆ A, then S∞ i=1 Ai∈ A. 3. µ:A → [0,∞]is a measure, i.e. (a) µ(∅)=0, (b) (Countable additivity) For any countable collection of pairwise disjoint sets {Ai}∞ i=1 ⊆ A, µ ∞ [ i=1 Ai!= ∞ X i=1 µ(Ai). Example 3.1 (Counting Measure).Let Xbe any set and let A= 2Xbe the collection of all subsets of X. The counting measure µc:A → [0,∞]is defined by µc(A) = (|A|,if Ais finite, ∞,if Ais infinite. Then (X, A, µc)is a measure space. [21] 23 Definition 3.6 (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)|. [17] Definition 3.7 (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. 24 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. Definition 3.8 (Support of a Function).Let Xbe a set and f:X→Ca function. The support of f, denoted by supp(f), is defined as supp(f) = {x∈X|f(x)= 0}. In words, the support of fis the set of all points where fdoes not vanish. Definition 3.9 (Finite Support Function Space).The Finite Support Function Space on Znis defined as F0(Zn) := {f:Zn→C|supp(f)is finite}. That is, f∈F0(Zn)if and only if there exists a finite set A⊂Znsuch that f(x) = 0 for all x /∈A. Definition 3.10 (Discrete Schwartz Space).Let α= (α1, . . . , αn)∈Nnbe a multi-index and, for x= (x1, . . . , xn)∈Zn, write xα:= n Y i=1 |xi|αi. Define the seminorms pα(f) := sup x∈Znxαf(x)∈[0,∞]. The discrete Schwartz space on Znis S(Zn) := f:Zn→C:pα(f)<∞for all α∈Nn. Equivalently, letting |x|:= px2 1+···+x2 n, S(Zn) = nf:Zn→C: sup x∈Zn (1 + |x|)k|f(x)|<∞for all k∈No. Example 3.2. (1) If fhas finite support, then f∈ S(Zn). (2) The Gaussian restriction f(x) = e−|x|2,x∈Zn, satisfies sup x∈Zn (1 + |x|)ke−|x|2<∞for all k, hence f∈ S(Zn). (3) The polynomially decaying sequence f(x) = (1 + |x|)−mbelongs to S(Zn) iff mcan be taken arbitrarily large; for fixed mit is not in S(Zn). 25 3.4 Operator Theory Definition 3.33 (Linear Operator).Let Ω⊆Zn, and let T:D(T)⊆ L2(Ω) → L2(Ω) be a map with domain D(T). We say that Tis a linear operator if for all u, v ∈ D(T)and all scalars α, β ∈C, the following holds: T(αu +βv) = αT(u) + βT(v). Definition 3.34 (Linear Shift Operator).Let Ω⊆Znand m∈N. A linear shift operator Lof absolute order at most mis any map of the form L(u(x)) = X ∥α∥1≤m aα(x)Eαu(x),x∈Ω, where ∥α∥1:= |α1|+···+|αn|, Eαu(x) := u(x+α), and the coefficient functions aα: Ω →Care given. Here Ek xidenotes the one–step shift in the i-th coordinate, Ek xiu(x1, . . . , xi, . . . , xn) := u(x1, . . . , xi+k,...,xn), so that Eα=Eα1 x1···Eαn xn, where α= (α1, α2, ..., αn)∈Zn. We take D(L)to be the set of all u∈ L2(Ω) for which the right-hand side is well-defined and belongs to L2(Ω). Proposition 3.35. Let Lbe a linear shift operator on a function space F(Ω), where Ω⊆Zn. Then the solution set of the homogeneous partial difference equation L(u)=0 is a function subspace(which is also a vector space) of F(Ω). Proof. If u, v ∈ker Land α, β ∈C, then L(αu +βv) = αL(u) + βL(v) = 0, so αu +βv ∈ker L. Definition 3.36 (Bounded Linear Operator).Let (X, ∥·∥X)and (Y, ∥·∥Y)be normed vector spaces. A linear operator T:X→Yis called a bounded linear operator if there exists a constant C≥0such that ∥T(x)∥Y≤C∥x∥Xfor all x∈X. The smallest such constant Cis called the operator norm of T, denoted by ∥T∥= sup ∥x∥X=1 ∥T(x)∥Y. 32 Example 3.6 (Difference operator is a bounded linear operator).Let X= Lp(Z)for 1≤p < ∞, equipped with the norm ∥u∥Lp=X x∈Z|u(n)|p1/p. Define the difference operator ∆ : Lp(Z)→ Lp(Z)by ∆u(x) = u(x+ 1) −u(x). Proof. First, ∆is linear since ∆(αu(x)+βv(x)) = (αu(x+1)+βv(x+1))−(αu(x)+βv(x)) = α∆u(x)+β∆v(x). Next, by the triangle inequality, |u(x+ 1) −u(x)|p≤2p−1|u(x+ 1)|p+|u(x)|p. Summing over x∈Zgives ∥∆u∥p Lp≤2p∥u∥p Lp. Hence ∥∆u∥Lp≤2∥u∥Lp. Therefore ∆is a bounded linear operator with operator norm ∥∆∥ ≤ 2. Definition 3.37 (Compact Operator).Let Xand Ybe Banach spaces. A bounded linear operator T:X→Yis called a compact operator if for every bounded sequence {xn} ⊂ X, the sequence {Txn}has a convergent subsequence in Y. Equivalently, Tmaps every bounded set in Xinto a relatively compact set in Y. Definition 3.38 (Spectrum of an Operator).Let Xbe a Banach space and let T:X→Xbe a bounded linear operator. The spectrum of Tis defined as σ(T) := {λ∈C|(T−λI)is not invertible as a bounded operator on X}. In other words, λ∈σ(T)if the operator T−λI does not have a bounded inverse on X. Definition 3.39 (Spectral Radius).Let Tbe a bounded linear operator on a Banach space X. The spectral radius of T, denoted by ρ(T), is defined as ρ(T) := sup{|λ|:λ∈σ(T)}, where σ(T)denotes the spectrum of T. 33 Theorem 3.40 (Abelian Group of Partial Shift Operators).Let u:Zn→Cbe a discrete function. For k= (k1, . . . , kn)∈Zn, define the partial shift operator Ek:= Ek1 x1Ek2 x2···Ekn xn,(Eki xiu)(x1, . . . , xn) = u(x1, . . . , xi+ki, . . . , xn). Then the set E={Ek:k∈Zn} forms an abelian group under composition. Moreover, EkEm=Ek+m,(Ek)−1=E−k, and therefore E∼ =(Zn,+). Proof. Closure. For any k,m∈Zn, EkEm=Ek1 x1Em1 x1···Ekn xnEmn xn=Ek1+m1 x1···Ekn+mn xn=Ek+m, hence Eis closed under composition. Identity. The neutral element is the zero shift E0=I. Inverse. For any k∈Zn, EkE−k=E0=I. Associativity. This follows from associativity of operator composition. Commutativity. Since partial shifts in different coordinates commute, Eki xiEmj xj=Emj xjEki xi, i =j, we obtain EkEm=EmEk. Thus the group is abelian. All group axioms are satisfied, and the mapping k7→ Ekis a group isomorphism from (Zn,+) to E. Proposition 3.41 (Lie Group of the Continuous Partial Shift Operators).Let f:Rn→Rbe any function, and for each k∈Rndefine the continuous partial shift operator (Ekf)(x) := f(x+k), x ∈Rn. Then the family of operators G:= {Ek:k∈Rn} forms an n-dimensional abelian Lie group under composition. Moreover, this Lie group is smoothly isomorphic to (Rn,+). 34 Proof. (1) Group structure. For any k, h ∈Rnand any function f, (EkEhf)(x) = Ehf(x+k) = f(x+k+h)=(Ek+hf)(x). Thus EkEh=Ek+h, which shows closure and defines the group law. The identity element is E0, since E0f(x) = f(x). For each k∈Rn, the inverse is (Ek)−1=E−k, since EkE−k=E0. Hence Gis a group. (2) Smooth manifold structure. Define the map Φ : Rn→ G, k 7→ Ek. This map is bijective, with inverse given by Φ−1(Ek) = k. Thus Ginherits a smooth manifold structure from Rnvia Φ. In particular, Gis an n-dimensional smooth manifold. (3) Smoothness of group operations. Under the identification Φ, the group law becomes k·h=k+h, which is smooth on Rn. Similarly, inversion corresponds to k−1=−k, which is also smooth. Hence both multiplication and inversion on Gare smooth with respect to the inherited manifold structure. (4) Conclusion. The set of continuous partial shift operators Gis therefore an n-dimensional abelian Lie group, smoothly isomorphic (indeed isomorphic as Lie groups) to (Rn,+). Definition 3.42 (Adjoint Operator).Let L2(Ω) be a discrete Hilbert space with inner product ⟨f, g⟩:= X x∈Ω f(x)g(x). Let T:D(T)⊆ L2(Ω) → L2(Ω) be a linear operator. We say that a linear operator T∗:D(T∗)⊆ L2(Ω) → L2(Ω) is the adjoint of Tif: ⟨Tf, g⟩=⟨f, T∗g⟩for all f∈ D(T), g ∈ D(T∗). Here, D(T∗)consists of all g∈ L2(Ω) such that the map f7→ ⟨Tf, g⟩ is continuous (i.e., bounded) on D(T). 35 Example 3.7. Let Ω⊆Znbe shift-invariant under α∈Zn(e.g. Ω = Zn, or Ω with periodic boundary conditions), and consider the Hilbert space L2(Ω) with inner product ⟨f, g⟩:= X x∈Ω f(x)g(x). For α∈Zn, define the (multi-dimensional) shift (Eαu)(x) := u(x+α), where Eα=Eα1 x1···Eαn xn, α = (α1, . . . , αn)∈Zn is the vector representing the direction (and magnitude) of the shift. Then for all u, v ∈ L2(Ω), ⟨Eαu, v⟩=X x∈Ω u(x+α)v(x) = X y∈Ω u(y)v(y−α) =⟨u, E−αv⟩. Hence the adjoint of the shift is the opposite shift: (Eα)∗=E−α. Definition 3.43 (Unitary Operator).Let Ω⊆Znand let U:D(U)⊆ L2(Ω) → L2(Ω) be a linear operator with domain D(U). We say that Uis unitary if: U∗U=UU∗=I, that is, U∗=U−1, where U∗is the adjoint of Uand U−1is the inverse operator of U. Equivalently, Uis unitary if it preserves the inner product: ⟨Uf, Ug⟩=⟨f, g⟩for all f, g ∈ L2(Ω). Proposition 3.44. Let Ω⊆Znbe shift-invariant under ±α∈Zn(e.g. Ω = Zn or Ωwith periodic boundary conditions). Consider the Hilbert space L2(Ω) with inner product ⟨f, g⟩:= X x∈Ω f(x)g(x). For α∈Zn, define the shift operator Eαu(x) := u(x+α). Then Eαis a unitary operator on L2(Ω). 36 Proof. From the calculation in the previous example, we have ⟨Eαu, v⟩=⟨u, E−αv⟩, which shows (Eα)∗=E−α. Since EαE−α=E−αEα=I, it follows that (Eα)∗= (Eα)−1. Thus Eαis unitary by definition. Equivalently, for all u∈ L2(Ω), ∥Eαu∥2 2=X x∈Ω|u(x+α)|2=X y∈Ω|u(y)|2=∥u∥2 2, so Eαpreserves the norm and is surjective. Definition 3.45 (Self-adjoint Operator).Let Ω⊆Zn, and let T:D(T)⊆ L2(Ω) → L2(Ω) be a linear operator with domain D(T). We say that Tis self-adjoint if it satisfies: T=T∗and D(T) = D(T∗), where T∗denotes the adjoint operator of T, defined by ⟨Tu, v⟩=⟨u, T ∗v⟩for all u∈ D(T), v ∈ D(T∗). Theorem 3.46 (Self-adjointness of sums and products of commuting operators).Let Ω⊆Zn, and let L1, L2, . . . , Lm:L2(Ω) → L2(Ω) be linear operators satisfying: (1) Self-adjointness: L∗ i=Li, i = 1, . . . , m. (2) Pairwise commutativity: LiLj=LjLi,∀i, j. Then: (a) The sum L= m X i=1 Li is self-adjoint. 37 (b) The product P= m Y i=1 Li is also self-adjoint. Proof. (a) Sum of self-adjoint operators. Using linearity of the adjoint, m X i=1 Li!∗ = m X i=1 L∗ i= m X i=1 Li, so the sum is self-adjoint. (b) Product of commuting self-adjoint operators. Let P=L1L2···Lm. Using the adjoint rule (AB)∗=B∗A∗, P∗= (L1L2···Lm)∗=L∗ mL∗ m−1···L∗ 1. Since each Liis self-adjoint, P∗=LmLm−1···L1. Commutativity gives LmLm−1···L1=L1L2···Lm=P. Thus P∗=P, so Pis self-adjoint. Example 3.8 (Discrete Laplacian is self-adjoint).Let Ω⊆Zbe shift-invariant under ±1(e.g. Ω = Zor periodic boundary conditions on a finite lattice), and consider the Hilbert space L2(Ω) with inner product ⟨f, g⟩:= X x∈Ω f(x)g(x). Define the 1D discrete Laplacian by ∇2u:= Eu +E−1u−2u, i.e. ∇2=E+E−1−2I, where Eku=u(x+k). Since E∗=E−1(and hence (E−1)∗=E), we have (∇2)∗= (E+E−1−2I)∗=E∗+ (E−1)∗−2I=E−1+E−2I=∇2. Equivalently, for all u, v ∈ L2(Ω), ⟨∇2u, v⟩=⟨Eu, v⟩+⟨E−1u, v⟩−2⟨u, v⟩ =⟨u, E−1v⟩+⟨u, Ev⟩−2⟨u, v⟩=⟨u, (E−1+E−2I)v⟩=⟨u, ∇2v⟩. Hence ∇2is self-adjoint on L2(Ω). 38 Example 3.9 (Discrete n-Dimensional Laplacian is self-adjoint).Let Ω⊆Zn be shift-invariant under ±eifor i= 1, . . . , n (e.g. Ω = Znor periodic boundary conditions on a finite lattice), and consider the Hilbert space L2(Ω) with inner product ⟨f, g⟩:= X x∈Ω f(x)g(x). For each coordinate i, let Exidenote the unit shift in the i-th direction: Exiu:= u(x+ei), E−1 xiu:= u(x−ei). Define the discrete Laplacian by ∇2u:= n X i=1 Exiu+E−1 xiu−2u⇐⇒ ∇2= n X i=1 Exi+E−1 xi−2I. Using (Exi)∗=E−1 xi(hence (E−1 xi)∗=Exi), we have (∇2)∗= n X i=1 (Exi)∗+ (E−1 xi)∗−2I= n X i=1 E−1 xi+Exi−2I=∇2. Equivalently, for all u, v ∈ L2(Ω), ⟨∇2u, v⟩= n X i=1 ⟨Exiu, v⟩+⟨E−1 xiu, v⟩−2⟨u, v⟩= n X i=1 ⟨u, E−1 xiv⟩+⟨u, Exiv⟩−2⟨u, v⟩=⟨u, ∇2v⟩. Hence ∇2is self-adjoint on L2(Ω). Example 3.10 (Spectrum of the Discrete Laplacian).Consider the discrete Laplacian defined on a uniform grid (∇2u)(x) = u(x+ 1) + u(x−1) −2u(x). We study the eigenvalue problem −∇2u=λu, u(0) = u(L) = 0. The solutions are given by the discrete sine functions uk(j) = sinkπj L, j = 0,1, . . . , L, with corresponding eigenvalues λk= 4 sin2kπ 2L, k = 1,2, . . . , L −1. Therefore, the spectrum of the discrete Laplacian with Dirichlet boundary conditions is σ(−∇2) = n4 sin2kπ 2L:k= 1,2, . . . , L −1o. 39 Proposition 3.47 (Self-adjointness of the Moore Laplacian).Let Ω⊆Znbe shift-invariant under ±eifor i= 1, . . . , n. The Moore Laplacian ∇2 Mu=X ∅=S⊆{1,...,n} Y i∈S δxi∆xi!u is a self-adjoint operator on L2(Ω). Proof. From the earlier sections, each one-dimensional second-order difference operator δxi∆xisatisfies: (δxi∆xi)∗=δxi∆xi, i = 1, . . . , n, that is, each is self-adjoint on L2(Ω). Moreover, these operators commute pairwise: δxi∆xiδxj∆xj=δxj∆xjδxi∆xi,∀i, j, because all partial shift operators commute. By the theorem on sums and products of commuting self-adjoint operators, any finite product Y i∈S δxi∆xi,∅ =S⊆ {1, . . . , n}, is self-adjoint. Finally, the Moore Laplacian is a finite sum of such products: ∇2 M=X ∅=S⊆{1,...,n}Y i∈S δxi∆xi, hence it is self-adjoint. This completes the proof. 40 4 Discrete Functionals and Convergence 4.1 Introduction In this section, we introduce the notion of discrete functionals, and study the compactness of the unit ball in discrete function spaces. We then recall some fundamental theorems from functional analysis, such as the Hahn–Banach theorem and the Riesz representation theorem, and explain how these ideas extend naturally to the discrete setting. Finally, we present several different notions of convergence for discrete functions, including pointwise convergence, uniform convergence, and weak convergence. [8,22] 4.2 Discrete Functionals Definition 4.1 (Discrete Functional).Let Ω⊆Znand let Vbe a vector space of functions u: Ω →C(or R). A discrete functional on Vis a mapping F:V→C, that assigns to each function u∈Va scalar F(u). Example 4.1. 1. Dirac functional: For a fixed x0∈Ω, define F(u) = u(x0). This functional simply evaluates the function uat the point x0. 2. Summation functional: Define F(u) = X x∈Ω u(x), whenever the sum converges. More generally, with weights w: Ω →C, F(u) = X x∈Ω w(x)u(x). Definition 4.2 (Bounded Linear Functional).Let (X, ∥·∥)be a normed vector space over Ror C. A mapping F:X→R(or C) is called a linear functional if F(αx +βy) = αF(x) + βF(y),∀x, y ∈X, ∀α, β ∈R(or C). The functional Fis said to be bounded (or continuous) if there exists a constant C≥0such that |F(x)| ≤ C∥x∥,∀x∈X. The smallest such constant Cis called the operator norm of F, denoted ∥F∥:= sup x∈X, x=0 |F(x)| ∥x∥. 41 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 5.1 (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. 5.2 Discrete Convolution Let f, g :Z→Cbe discrete functions. The discrete convolution of fand g is defined by (f∗g)(x) = X s∈Z f(x−s)g(s). 5.3 Basic Properties 1. Commutativity f∗g=g∗f. 2. Associativity (f∗g)∗h=f∗(g∗h). 3. Distributivity f∗(g+h) = f∗g+f∗h. 48 4. Shift identity Let Ebe the shift operator E f =f(x+ 1). Then E(f∗g)=(Ef)∗g=f∗(Eg). 5. Difference identity Let the forward difference be ∆f=f(x+ 1) −f(x). Then ∆(f∗g) = (∆f)∗g=f∗(∆g). 6. Shift–invariant operators commute with convolution If Lis a linear operator satisfying L(Ekf) = Ek(Lf) for all k∈Z, then L(f∗g)=(Lf)∗g=f∗(Lg). 7. Fourier transform turns convolution into multiplication The discrete Fourier transform (DFT) is b f(ω) = X x∈Z f(x)e−iωx, ω ∈[−π, π]. Then \ (f∗g)(ω) = b f(ω)bg(ω). 8. Convolution operators are Fourier multipliers If Tis a convolution operator Tf =f∗p, then c Tf(ω) = bp(ω)b f(ω). Thus the Fourier multiplier of Tis m(ω) = bp(ω). 9. Convolution as a Banach algebra In the absolutely summable space L1(Z), ∥f∗g∥L1≤ ∥f∥L1∥g∥L1. Thus L1(Z) is closed under convolution and forms a Banach algebra. 49 5.4 Fundamental Solution Let Lbe a linear shift operator, with (t, x)∈N0×Ω⊆Zn+1, and let u: N0×Ω→Rsolve the homogeneous linear evolution equation L(u)=0, with initial condition u(0,x) = f(x), together with suitable boundary conditions. If the initial data is the Kronecker delta f(x) = δ(x), then the corresponding solution G(t, x) is called the fundamental solution. In this paper, we also refer to Gas a Green’s function, or specifically, CauchyGreen’s Function, in the sense of initial value problems. For general initial data f, the solution is given by the discrete convolution u(t, x) = X s∈Ω f(s)G(t, x−s). 5.5 The Discrete Green’s Function Definition 5.4 (General Discrete Green’s Function).Let Ω⊆Zn, and consider a linear shift operator L:F(Z≥0×Ω) →F(Z≥0×Ω), acting on functions u:Z≥0×Ω→C. The discrete Green’s function is a kernel G: (Z≥0×Ω) ×(Z≥0×Ω) →C, defined by the property LG(t, x;τ, s)=δ(t−τ)δ(x−s),(t, x)∈Z≥0×Ω, subject to the causality condition G(t, x;τ, s)=0whenever t<τ, together with the prescribed boundary conditions on Ω. Proposition 5.5 (Representation of the Solution).Let f:Z≥0×Ω→Cbe a forcing term. Then the unique solution of L(u) = f, u(0,x)=0, is given by u(t, x) = t X τ=0 X s∈Ω G(t, x;τ, s)f(τ, s). 50 Corollary 5.6 (Shift-Invariant Case).If Lis shift-invariant in both time and space, then G(t, x;τ, s) = G(t−τ, x−s), so that the solution reduces to the discrete space–time convolution u(t, x)=(G∗f)(t, x) = t X τ=0 X s∈Ω G(t−τ, x−s)f(τ, s). Proposition 5.7 (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 with constant coefficient. Let G(t, x)be the discrete Green’s function, i.e., the solution of L(G) = δ(t)δ(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. 51 6 Discrete Fourier Transform 6.1 Motivation and Background In this section we develop the Fourier analytic framework necessary for the study of partial difference equations. The discrete Fourier series (DFS) and the discrete Fourier transform (DFT) provide natural tools for representing periodic discrete functions in terms of orthogonal exponential bases. These representations not only yield compact formulas for the coefficients and reconstruction of discrete functions, but also establish powerful theorems such as orthogonality relations, Parseval’s identity, and convergence results. Moreover, Fourier methods play a crucial role in solving linear partial difference equations by diagonalizing shift operators and constructing explicit solutions to initial–boundary value problems. The exposition here is inspired by classical results in Fourier analysis and adapted to the discrete function setting (cf. [16,28,29]). 6.2 Discrete Fourier Series We now introduce the discrete Fourier series (DFS), which provides an expansion of periodic discrete functions in terms of exponential basis functions. Let f: Z→Cbe a discrete function of period N, i.e. f(x+N) = f(x) for all x∈Z. Then fadmits the representation f(x) = N−1 X k=0 F(k)ei2π Nkx, where the Fourier coefficients are given explicitly by F(k) = 1 N N−1 X x=0 f(x)e−i2π Nkx, k = 0,1, . . . , N −1. This formula shows that every N-periodic discrete function can be decomposed into a finite linear combination of orthogonal exponential functions. In the following, we shall establish the orthogonality relations of the exponential basis, derive Parseval’s identity, and discuss the convergence properties of the discrete Fourier series. 6.3 Inner Product and Orthogonality We equip the space of N-periodic discrete functions with the normalized inner product ⟨f, g⟩:= 1 N N−1 X x=0 f(x)g(x). This structure makes the function space a finite-dimensional Hilbert space, naturally identified with CN. 52 Proposition 6.1 (Orthonormality of the exponential basis).Let ϕk(x) := ei2π Nkx, k = 0,1, . . . , N −1, x ∈Z. Then {ϕk}N−1 k=0 forms an orthonormal basis of L2(Ω) with Ω = {0,1, . . . , N −1}. In particular, ⟨ϕk, ϕm⟩=δkm,0≤k, m ≤N−1. Proposition 6.2 (Orthonormality of the exponential basis).Let ϕk(x) := ei2π Nkx, k = 0,1, . . . , N −1, x ∈Z. Then {ϕk}N−1 k=0 forms an orthonormal basis of L2(Ω) with Ω = {0,1, . . . , N −1}. In particular, ⟨ϕk, ϕm⟩=δkm,0≤k, m ≤N−1. Proof. By definition of the inner product, ⟨ϕk, ϕm⟩=1 N N−1 X x=0 ei2π N(k−m)x. If k=m, each term equals 1 and hence the sum equals 1. If k=m, the summand is a finite geometric progression with ratio ei2π N(k−m)= 1, and the sum vanishes. Thus ⟨ϕk, ϕm⟩=δkm. 6.4 Discrete Fourier Transform We now turn to the discrete Fourier transform (DFT), which is the natural Fourier analytic tool for functions defined on the entire integer lattice. Definition 6.3 (DFT and its inverse).Let f:Z→Cbe an infinite discrete function. The discrete Fourier transform of fis defined by F(ω) = X x∈Z f(x)e−iωx, ω ∈[−π, π]. The corresponding inverse transform is given by f(x) = 1 2πZπ −π F(ω)eiωx dω, x ∈Z. Remark 6.4 (Key properties). •The frequency variable ωis continuous and varies over the fundamental interval [−π, π]. •The spectrum F(ω)is periodic with period 2π, i.e. F(ω+ 2π) = F(ω). •This transform (commonly referred to as DTFT in engineering) is particularly useful in analyzing infinite discrete sequences, e.g. in digital signal processing, stability analysis, and discrete-time dynamical systems. 53 Convention. In this paper we refer to the Fourier transform of sequences defined on the full integer lattice Zas the Discrete Fourier Transform (DFT). This differs from the engineering convention, where “DFT” usually denotes the transform of finite sequences. The periodic finite case will be referred to here as the Discrete Fourier Series (DFS). 6.5 Parseval’s Identity Theorem 6.5 (Parseval’s identity: periodic discrete case).Let f:{0,1, . . . , N− 1} → Cand define the inner product ⟨f, g⟩:= 1 N N−1 X x=0 f(x)g(x). Let the exponential basis ϕk(x) = ei2π Nkx,k= 0, . . . , N −1, which is orthonormal under ⟨·,·⟩. Define Fourier coefficients F(k) = ⟨f, ϕk⟩=1 N N−1 X x=0 f(x)e−i2π Nkx. Then 1 N N−1 X x=0 |f(x)|2= N−1 X k=0 |F(k)|2. Equivalently, N−1 X x=0 |f(x)|2=N N−1 X k=0 |F(k)|2. Proof. Since {ϕk}is an orthonormal basis of the N-dimensional Hilbert space, the orthogonal expansion holds: f=PN−1 k=0 F(k)ϕkwith F(k) = ⟨f, ϕk⟩. Applying ∥f∥2=⟨f, f⟩and orthonormality, ⟨f, f⟩=DX k F(k)ϕk,X m F(m)ϕmE=X k,m F(k)F(m)⟨ϕk, ϕm⟩=X k|F(k)|2, i.e. 1 NPx|f(x)|2=Pk|F(k)|2. Theorem 6.6 (Parseval/Plancherel: infinite lattice (DFT in this paper)).Let f∈ L2(Z)and define its discrete Fourier transform (frequency continuous on [−π, π]) F(ω) = X x∈Z f(x)e−iωx, ω ∈[−π, π], with inverse f(x) = 1 2πZπ −π F(ω)eiωx dω. Then the Plancherel identity holds: X x∈Z|f(x)|2=1 2πZπ −π|F(ω)|2dω. 54 Proof sketch. Consider the isometric isomorphism F:L2(Z)→L2([−π, π]) given by f7→ Fabove. Using the orthogonality 1 2πRπ −πei(ω)(x−y)dω =δxy and Fubini/Tonelli, compute 1 2πZπ −π|F(ω)|2dω =1 2πZX x,y f(x)f(y)e−iω(x−y)dω =X x,y f(x)f(y)δxy =X x|f(x)|2. 6.6 Fourier Multiplier Operators We introduce one of the most fundamental notions in harmonic analysis. Definition 6.7. A linear operator Ton functions u:Zn→C(or u:Rn→C) is called a Fourier multiplier operator if its action in the Fourier domain is given by multiplication: c Tu(ξ) = m(ξ)bu(ξ), for some function m(ξ), called the multiplier or symbol of T. Example 6.1 (Differential operator).For the continuous derivative, ∂xeiξx =iξ eiξx, so m(ξ) = iξ. Example 6.2 (Shift operator).For the discrete shift, Exeiξx =eiξ eiξx, so m(ξ) = eiξ. Example 6.3 (Difference operator).Since ∆xu=Exu−u, we have ∆xeiξx = (eiξ −1) eiξx, so m(ξ) = eiξ −1. 6.7 Fourier Ansatz Consider the constant coefficient linear partial difference equation Etu(t, x) = X k∈S akEku(t, x), x ∈Zn,(1) where •S⊂Znis a finite index set, •ak∈Care constant coefficients, 55 •Ek=Ek1 x1···Ekn xndenotes the shift operator. We impose the initial condition u(0, x) = f(x), x ∈Zn, and assume no boundary conditions. We apply the discrete Fourier transform with respect to the spatial variable x: bu(t, ξ) = X x∈Zn u(t, x)e−iξ·x, ξ ∈[−π, π]n. Using the identities d Etu=Etbu, d Eku=eik·ξbu, equation (1) transforms into Etbu(t, ξ) = X k∈S akeik·ξ!bu(t, ξ). Definition 6.8 (Symbol).The symbol of the partial difference operator in (1) is defined by σ(ξ) := X k∈S akeik·ξ, ξ ∈[−π, π]n. For each fixed ξ, we obtain an ordinary difference equation in time: bu(t+ 1, ξ) = σ(ξ)bu(t, ξ). Solving this equation yields bu(t, ξ) = σ(ξ)tbu(0, ξ) = σ(ξ)tb f(ξ). Applying the inverse Fourier transform, we obtain the solution formula u(t, x) = 1 (2π)nZ[−π,π]nb f(ξ)σ(ξ)teiξ·xdξ. This expression provides the explicit solution of the constant coefficient linear partial difference equation (1). 56 7 Symbolic Analysis and Classification 7.1 Introduction In this chapter, we introduce the Cauchy problem for partial difference equations and investigate its well-posedness in the sense of Hadamard. In particular, we study the conditions under which a Cauchy problem admits existence, uniqueness, and continuous dependence on the initial data, as well as the mechanisms leading to ill-posedness. We then focus on linear partial difference equations and introduce the notion of the principal symbol. Using a low-frequency approximation of the symbol, we establish a classification theory analogous to that of partial differential equations. In particular, linear partial difference equations are classified into elliptic, parabolic, and hyperbolic types according to the spectral properties of the principal symbol. This symbolic classification provides a unified framework for understanding the qualitative behavior of solutions, including propagation, diffusion, and stability phenomena, in discrete dynamical systems. 7.2 Cauchy Problem Consider a partial difference equation of finite order F({Eαu}α∈S, x) = 0, x ∈Zn, where S⊂Znis a finite index set satisfying |α|1≤m, and Eα=Eα1 x1···Eαn xndenotes the shift operator. Let ϕ:Zn→Zbe a discrete scalar function. We define an (n−1)- dimensional discrete hypersurface Γ = {x∈Zn:ϕ(x)=0}. Let ∇cdenote the central discrete gradient. The (unit) normal direction to Γ is defined by v=∇cϕ |∇cϕ|, whenever ∇cϕ= 0. Let ∇denote the discrete gradient operator. We define the directional difference operator in the normal direction by ∆v=v·∇, and its iterates by ∆k vu= (v·∇)ku, k ∈N. 57 To analyze the large-scale behavior, we consider the low-frequency expansion around (ξ, η) = (0,0). Using the Taylor expansions sinη 2=η 2+O(|η|3),sinξ 2=ξ 2+O(|ξ|3), we obtain the quadratic approximation σ(ξ, η) = 1 4η2−c2ξ2+O(|(ξ, η)|4).(19) Up to an irrelevant constant factor, the principal symbol is therefore p(ξ, η) = η2−c2ξ2.(20) This quadratic form can be written in matrix form as p(ξ, η) = ξ ηAξ η, A =−c20 0 1. If c= 0 is real, then c2>0, and the matrix Ahas one positive and one negative eigenvalue. Therefore, the discrete wave equation is classified as a hyperbolic equation. 7.5 Trigonometric Polynomial and Quadric Surface This subsection serves as a visual and geometric guide explaining why low– frequency approximations and quadratic Taylor expansions are sufficient for the symbolic classification of partial difference equations. Rather than aiming at an exact global classification of trigonometric symbols, we focus on their local behavior near the origin in frequency space, which governs the large–scale dynamics of solutions. We begin by recalling several standard polynomial objects. Multivariable Polynomials A multivariable polynomial of degree mis a function of the form P(x1, . . . , xn) = X |α|≤m aαxα1 1xα2 2···xαn n, where α= (α1, . . . , αn)∈Nnand |α|=α1+α2+···+αn. Examples. 64 1. Quadratic function (one variable): P(x) = ax2+bx +c. 2. Conic section (two variables): P(x, y) = Ax2+Bxy +Cy2+Dx +Ey +F. 3. Quadric surface (three variables): P(x, y, z) = Ax2+By2+Cz2+Dxy +Eyz +Fzx +Gx +Hy +Iz +J. Such quadratic polynomials define conic sections in two dimensions and quadric surfaces in three dimensions. Trigonometric Polynomials Atrigonometric polynomial is a polynomial expression in sine and cosine functions, defined by Psin x1,...,sin xn,cos x1,...,cos xn, where Pis an ordinary multivariable polynomial. Example. P(sin x, cos x, sin y, cos y) = sin2x+ cos2y. Trigonometric polynomials arise naturally as symbols of partial difference equations, since shift operators act as phase factors eiξjunder Fourier transformation. Geometric Interpretation via Graphs Although trigonometric polynomials are globally periodic and bounded, their local behavior near low frequencies ξ= 0 is well approximated by quadratic polynomials through Taylor expansion: sin ξ≈ξ, sin2(ξ/2) ≈ξ2 4. This observation allows us to interpret the local geometry of a trigonometric symbol in terms of classical quadric surfaces. Example 7.5 (Hyperbolic Surface).Consider the quadratic surface z=x2−y2. Now consider the trigonometric polynomial z= 4 sin2(x/2) −4 sin2(y/2). 65 Figure 2: Graph of the quadric surface z=x2−y2. Figure 3: Graph of z= 4 sin2(x/2) −4 sin2(y/2). Near (x, y) = (0,0), we have the approximation 4 sin2(x/2) −4 sin2(y/2) ≈x2−y2. Hence, in a neighborhood of the origin, the trigonometric polynomial exhibits the same hyperbolic geometry as the corresponding quadric surface. 66 Figure 4: Graph of the elliptic paraboloid z=x2+y2. Figure 5: Graph of the trigonometric polynomial z= 4 sin2(x/2) + 4 sin2(y/2). Example 7.6 (Elliptic Paraboloid).The surface z=x2+y2 is an elliptic paraboloid, which serves as the canonical geometric model associated with elliptic quadratic forms. On the other hand, the trigonometric polynomial z= 4 sin2(x/2) + 4 sin2(y/2) 67 arises naturally as the symbol of the discrete Laplace operator. Using the Taylor expansion sin(x/2) = x 2+O(x3), we obtain, near (x, y) = (0,0), 4 sin2(x/2) + 4 sin2(y/2) = x2+y2+O(|(x, y)|4). Therefore, in a neighborhood of the origin, the trigonometric polynomial admits a quadratic approximation whose graph is locally indistinguishable from that of an elliptic paraboloid. This illustrates why low-frequency (long-wavelength) behavior of partial difference equations can be classified via quadratic forms and Hessian matrices, in complete analogy with the continuous PDE setting. This geometric correspondence explains why the classification of partial difference equations can be carried out using the quadratic approximation of their symbols: the low–frequency behavior captures the essential elliptic, parabolic, or hyperbolic nature of the equation, while higher–order terms only affect fine– scale oscillations. 68 7.6 Hadamard’s Well-Posedness Cauchy Problem for a Partial Difference Equation Let Xbe a space of initial data and Ya space of solutions. Consider a partial difference equation together with prescribed Cauchy data. (i) Existence. The Cauchy problem is said to satisfy existence if, for every initial datum u0∈X, there exists at least one solution u∈Ythat satisfies the partial difference equation together with the prescribed initial condition. (ii) Uniqueness. The Cauchy problem is said to satisfy uniqueness if, for every initial datum u0∈X, there exists at most one solution u∈Ysatisfying the equation and the initial condition. (iii) Solution Operator. If both existence and uniqueness hold, one may define the solution operator S:X→Y, S(u0) = u, where udenotes the unique solution corresponding to the initial datum u0. (iv) Continuous Dependence (Hadamard’s Third Condition). The Cauchy problem is said to satisfy Hadamard’s third condition if the solution operator S:X→Yis continuous with respect to the topologies of Xand Y. (v) Hadamard’s Well-Posedness. The Cauchy problem is said to be Hadamard’s well-posed if it satisfies existence, uniqueness, and continuous dependence on the initial data. Otherwise, the problem is said to be ill-posed. 7.7 Stability of Partial Difference Equations Definition 7.4 (Instability).Let Xbe a space of initial data and let u(t) = Stu0 denote the solution of a Cauchy problem for a partial difference equation, where St:X→Xis the associated solution operator. Definition 1 (Operator-Norm Definition). The Cauchy problem is said to be unstable if sup t≥0∥St∥X→X=∞. In other words, the solution operator is not uniformly bounded in time. 69 Equivalent Formulation. The system is unstable if and only if for every constant C > 0, there exist a time t≥0 and some initial datum u0= 0 such that ∥Stu0∥X> C ∥u0∥X. Exponential Instability. The system is said to be exponentially unstable if there exist constants c > 0 and u0= 0 such that ∥Stu0∥X∼ect ∥u0∥Xas t→ ∞. Equivalently, there exists c > 0 such that ∥St∥X→X≥ect for sufficiently large t. Remark (Distinction Between Instability and Ill-Posedness). Instability should be distinguished from ill-posedness in the sense of Hadamard: •The problem is unstable if ∥St∥X→X<∞for each fixed t≥0,but sup t≥0∥St∥X→X=∞. •The problem is ill-posed if there exists some t0>0 such that ∥St0∥X→X=∞. Thus, the essential distinction lies in whether the loss of control occurs at a fixed finite time or only asymptotically as time evolves. Summary. A Cauchy problem is said to be unstable if the associated solution operator is not uniformly bounded in time, that is, sup t≥0∥St∥X→X=∞. In particular, exponential growth of solutions implies instability, even though the problem may still be Hadamard well-posed. Example 7.7 (Discrete Backward Heat Equation).Consider the one–dimensional discrete backward heat equation ∆tu=−δx∆xu, with initial condition u(0, x) = f(x), and no boundary condition imposed. We seek solutions of the form u(t, x) = eλteiξx. 70 Substituting this ansatz into the equation yields eλ−1 = 4 sin2(ξ/2), and hence eλ= 1 + 4 sin2(ξ/2). Therefore, the solution can be written in Fourier representation as u(t, x) = 1 2πZπ −πb f(ξ)1 + 4 sin2(ξ/2)teiξx dξ. The amplification factor is given by G(ξ) := 1 + 4 sin2(ξ/2). Since sin2(ξ/2) ∈[0,1], it follows that G(ξ)∈[1,5], and in particular G(ξ)≥1for all ξ∈[−π, π]. Moreover, •G(ξ)=1if and only if ξ= 0 (the zero frequency mode), •G(ξ)>1for all ξ= 0. Time evolution. The solution satisfies u(t, x) = 1 2πZπ −πb f(ξ)G(ξ)teiξx dξ. Since G(ξ)t=1 + 4 sin2(ξ/2)t, every nonzero Fourier mode obeys G(ξ)t∼ec(ξ)t, c(ξ)>0. Hence, all nonzero frequency modes grow exponentially in time. Stability analysis. Let Stdenote the solution operator. Then sup t≥0∥St∥=∞, and in fact the operator norm grows exponentially in time due to the unbounded amplification of high-frequency modes. 71 Consequently, the discrete backward heat equation is exponentially unstable. Remark on well-posedness. For each fixed integer time t, the amplification factor satisfies |G(ξ)|t<∞, so the solution operator Stis well-defined and unique. However, sup t≥0∥St∥=∞, and therefore the problem fails to be stable in the sense of Hadamard. In the discrete-time setting, the backward heat equation is thus well-defined but severely unstable. 72 8 First Order Equations in Time 8.1 Introduction In this section, we study linear partial difference equations whose structural order in the time direction is one. Such systems evolve step by step in a manner where the next state depends only on the present state, and therefore they may naturally be referred to as Markov systems. These equations serve as the discrete-time analogue of first-order evolution equations in the continuous setting, and they provide a fundamental framework for analyzing propagation, transport, and probabilistic models in discrete space-time lattices. Definition 8.1 (Markov System).Let u:Z1+n→R(or C) be a function of discrete time t∈Zand spatial variables (x1, . . . , xn)∈Zn. A Markov system is a first–order time evolution equation of the form Etu(t, x) = F{Ek1 x1···Ekn xnu(t, x)}(k1,...,kn)∈S, t, x, where S⊂Znis a finite index set. This means that the state at time t+ 1 depends only on the configuration at time t, but not on any earlier times. 8.2 1D Pascal Evolution Equation We now consider the linear partial difference equation inspired by the recurrence relation in combinatorics [24] Etu=u+E−1 xu, u :Z2→R, with initial condition u(0, x) = f(x)∈ L2(Z), and no boundary conditions. Fourier Ansatz. We apply the Fourier ansatz u(t, x) = eλteikx. Substitution into the equation gives the dispersion relation eλ= 1 + e−ik. Integral Representation of the Solution. The general solution can be written as u(t, x) = Zπ −π A(k) (1 + e−ik)teikx dk. Using the initial condition t= 0, we obtain u(0, x) = f(x) = Zπ −π A(k)eikx dk. 73 In general, by induction, X(x) = X(0) (λ−1)(λ−2) ···(λ−x), x ≥1. In falling factorial notation, X(x) = X(0) (λ−1)x . Characteristic Solution. Thus the separated solution is uλ(t, x) = λt (λ−1)x . General Solution. By linear superposition, the general solution can be expressed as u(t, x) = X λ A(λ)λt (λ−1)x . Initial Condition. The initial condition u(0, x) = δ(x) imposes the constraint u(0, x) = X λ A(λ) (λ−1)x . - For x= 0: u(0,0) = 1 = X λ A(λ). - For x≥1: u(0, x) = 0 = X λ A(λ) (λ−1)x . This is precisely the structure of a binomial inversion formula. Solving for A(λ) gives A(λ) = (−1)x−λ x!x λ,0≤λ≤x. Substituting back, we obtain u(t, x) = 1 x! x X λ=0 (−1)x−λx λλt. Conclusion. The closed-form expression for u(t, x) is exactly the Stirling number of the second kind S(t, x). Hence we conclude that the Stirling numbers arise naturally as the Green’s function of the non-autonomous partial difference equation Etu=E−1 xu+xu. 80 9 Second Order Equations in Time 9.1 Introduction In this section, we study linear partial difference equations with structural order 2 in the time direction. Such systems depend not only on the present state but also on the previous time step, analogous to second-order evolution equations in the continuous setting. They naturally capture wave-like and oscillatory behaviour in discrete space-time lattices. Examples. •One-dimensional equation: Etu=aE−1 xu+bu +cExu+pE−1 tE−1 xu+qE−1 tu+rE−1 tExu. •Two-dimensional equation: Etu=u+E−1 xu+E−1 yu+E−1 tu. 9.2 Second–Order Pascal Evolution Equation We consider the second–order linear partial difference equation Etu=E−1 xu+Exu+E−1 tu, u :Z2→R. This equation extends the classical Pascal Evolution Equation by including a memory term through E−1 tu. Consequently, the Green’s function no longer produces the standard Pascal triangle, but instead generates a novel combinatorial structure that exhibits duplication and oscillatory patterns. We refer to this system as the Second–Order Pascal Evolution Equation. Initial Condition. u(0, x) = f(x), u(1, x) = g(x) f, g ∈ L2(Z) Fourier Ansatz. We seek solutions of the form u(t, x) = eλteikx. Substitution yields eλ=e−ik +eik +e−λ. Equivalently, with z=eλ, this gives the quadratic z2−2 cos(k)z−1=0, 81 whose roots are z±=z±(k) = cos(k)±p1 + cos2(k). General Solution. By superposition, the solution admits the Fourier integral representation u(t, x) = 1 2πZπ −πA(k)zt ++B(k)zt −eikx dk, where A(k), B(k) are determined from the initial conditions. Determination of Coefficients. At t= 0, u(0, x) = f(x) = 1 2πZπ −πA(k) + B(k)eikx dk, so that A(k) + B(k) = b f(k),b f(k) = X x∈Z f(x)e−ikx. At t= 1, u(1, x) = g(x) = 1 2πZπ −πA(k)z+(k) + B(k)z−(k)eikx dk, so that A(k)z+(k) + B(k)z−(k) = bg(k),bg(k) = X x∈Z g(x)e−ikx. Thus, A(k) = bg(k)−b f(k)z−(k) z+(k)−z−(k), B(k) = b f(k)z+(k)−bg(k) z+(k)−z−(k). Final Representation. The general solution is therefore u(t, x) = 1 2πZπ −π"bg(k)−b f(k)z−(k) z+(k)−z−(k)z+(k)t+b f(k)z+(k)−bg(k) z+(k)−z−(k)z−(k)t#eikx dk. Special Case: Delta Initial Data. For the choice f(x) = δ(x), g(x) = δ(x), the Fourier transforms are b f(k) = bg(k) = 1. Thus, u(t, x) = 1 2πZπ −π (1 −z−)zt ++ (z+−1)zt − z+−z− eikx dk. This Fourier integral represents a special solution associated with the second– order Pascal Evolution Equation and generates the novel combinatorial structure observed in the numerical pattern. 82 Figure 6: Integer array generated by the Second–Order Pascal Evolution Equation with two–delta initial condition. Notice the duplicated values along the diagonals and oscillatory structures emerging inside the triangle. 83 9.3 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. Solution We consider the discrete wave equation δt∆tu=c2∇2u, which in shift–operator form reads Etu+E−1 tu−2u=c2Exu+E−1 xu−2u, u :Z2→R. Initial and boundary conditions. u(0, x) = f(x), u(1, x) = g(x) f, g ∈ L2(Ω),Ω = {x∈Z: 1 ≤x≤L−1} u(t, 0) = u(t, L)=0. Separation of variables. We assume u(t, x) = T(t)X(x). Substitution gives EtT−2T+E−1 tT T=c2ExX−2X+E−1 xX X=−λ, for some separation constant λ. 84 Spatial problem (Discrete Sturm–Liouville). ExX−2X+E−1 xX=−λ c2X, X(0) = X(L) = 0. The eigenfunctions and eigenvalues are Xn(x) = sinnπ Lx, λn= 4c2sin2nπ 2L, n = 1,2, . . . , L −1. Temporal problem. EtT−2T+E−1 tT=−λT. Its characteristic polynomial is z2−(2 −λ)z+ 1 = 0, with roots z±(λ) = 2−λ±p(2 −λ)2−4 2. Hence Tn(t) = Anz+(λn)t+Bnz−(λn)t. Mode solutions. The n-th mode is un(t, x) = Anz+(λn)t+Bnz−(λn)tsinnπ Lx. General solution. u(t, x) = L−1 X n=1 Anz+(λn)t+Bnz−(λn)tsinnπ Lx. Determination of coefficients. From u(0, x) = f(x) and u(1, x) = g(x), we expand: f(x) = L−1 X n=1 (An+Bn) sinnπ Lx, g(x) = L−1 X n=1 (Anz+(λn) + Bnz−(λn)) sinnπ Lx. By orthogonality of the sine basis, An+Bn=2 L L−1 X x=1 f(x) sinnπ Lx, 85 Anz+(λn) + Bnz−(λn) = 2 L L−1 X x=1 g(x) sinnπ Lx. Explicit formulas. Solving this 2 ×2 linear system gives An=1 z+(λn)−z−(λn) 2 L L−1 X x=1 g(x)−z−(λn)f(x)sinnπ Lx!, Bn=1 z−(λn)−z+(λn) 2 L L−1 X x=1 g(x)−z+(λn)f(x)sinnπ Lx!. Thus the full solution is completely determined. Remark 9.1. The temporal part of the solution is expressed in terms of the characteristic roots Tn(t) = Anzt n,++Bnzt n,−, zn,±=2−λn±p(2 −λn)2−4 2. When the discriminant is negative, the roots zn,±are complex conjugates and the solution can be rewritten in trigonometric form using cos(ωnt)and sin(ωnt). Otherwise, when the roots are real, it is more natural to keep the exponential representation. Both forms are mathematically equivalent and depend only on the spectral parameter λn. 86 10 Steady State Problems 10.1 Introduction In this chapter, we study the stationary or steady-state problem for discrete field equations. Our goal is to develop a discrete analogue of the classical elliptic theory arising in continuous partial differential equations. We begin by introducing the discrete counterparts of the Laplace and Poisson equations, defined on integer lattices Znor on finite discrete domains with prescribed boundary conditions. In addition to the standard nearest-neighbour Laplacian, we also propose the Moore Laplacian and the corresponding Moore Poisson equation, obtained by extending the stencil to the full Moore neighbourhood. These operators retain many of the structural properties of their continuous elliptic counterparts, such as self-adjointness, positivity, and discrete maximum principles, while encoding richer geometric or combinatorial interactions. The discrete elliptic equations considered here take the general form Lu=f, where Lis one of the Laplace-type operators introduced above. We will examine the solvability, properties of solutions, and the relationship between these discrete formulations and the classical elliptic PDEs. In particular, we show that the discrete equations exhibit behavior closely analogous to continuous elliptic problems, including smoothing effects, uniqueness of solutions, and harmonicity on discrete domains. 87 10.2 2D Discrete Laplace Equation We now consider the two-dimensional discrete Laplace equation δx∆xu+δy∆yu= 0, u =u(x, y), subject to the Dirichlet boundary conditions u(0, y) = u(L, y) = u(x, M)=0, u(x, 0) = f(x), f(0) = f(L) = 0. Separation of Variables. We seek a separable solution of the form u(x, y) = X(x)Y(y). Substituting into the equation gives Y(y)δx∆xX+X(x)δy∆yY= 0, which can be rearranged as δx∆xX X(x)+δy∆yY Y(y)= 0. Thus, each term must equal a constant −λ, giving two ordinary difference equations: δx∆xX=−λX, δy∆yY=λY. Shift Operator Form. Using the shift operators Exand Ey, these can be written as ExX+E−1 xX−2X=−λX, (21) EyY+E−1 yY−2Y=λY. (22) Solution for X(x).Expanding Eq. (21) gives ExX+E−1 xX= (2 −λ)X. Assume a trial solution X(x) = rx, leading to the characteristic equation r+1 r= 2 −λ, or equivalently, r2−(2 −λ)r+ 1 = 0. The roots are r±=e±iθ,with cos θ= 1 −λ 2. Hence the general solution for X(x) is X(x) = Acos(θx) + Bsin(θx). Applying the Dirichlet conditions X(0) = X(L) = 0 yields A= 0, θn=nπ L, n = 1,2, . . . , L −1, so that Xn(x) = sinnπx L, λn= 4 sin2nπ 2L. 88 Solution for Y(y).Substituting λninto Eq. (22) gives EyY+E−1 yY= (2 + λn)Y. Assume Y(y) = ry, giving r+1 r= 2 + λn, r±=e±µn,where cosh µn= 1 + λn 2. Thus the general solution is Yn(y) = Cneµny+Dne−µny. Applying the boundary condition Y(M) = 0 gives Yn(y) = sinhµn(M−y), up to normalization, and we set Yn(0) = 1 for convenience: Yn(y) = sinhµn(M−y) sinh(µnM). Complete Solution. By the superposition principle, the full solution satisfying u(x, 0) = f(x) is u(x, y) = L−1 X n=1 Ansinnπx Lsinhµn(M−y) sinh(µnM),where cosh µn= 1+2 sin2nπ 2L. The coefficients Anare determined from the boundary data: f(x) = L−1 X n=1 Ansinnπx L. Hence, An=2 L L−1 X x=1 f(x) sinnπx L. Final Form. The discrete Laplace equation with Dirichlet boundary conditions admits the formal solution u(x, y) = L−1 X n=1 "2 L L−1 X x′=1 f(x′) sinnπx′ L#sinnπx Lsinhµn(M−y) sinh(µnM), where λn= 4 sin2nπ 2L,cosh µn= 1 + λn 2. 89 2D case. δt∆tu=c2δx∆xu+δy∆yu. 3D case. δt∆tu=c2δx∆xu+δy∆yu+δz∆zu. 4. Moore Heat Equation The Moore Laplacian in two dimensions is defined as ∇2 Mu:= δx∆xu+δy∆yu+δx∆xδy∆yu. The Moore heat equation is then ∆tu=α∇2 Mu. 2D case. ∆tu=αδx∆xu+δy∆yu+δx∆xδy∆yu. 3D case. The natural 3D Moore Laplacian includes all two–coordinate couplings: ∇2 Mu=∇2u+δx∆xδy∆yu+δy∆yδz∆zu+δz∆zδx∆xu +δx∆xδy∆yδz∆zu. Hence the 3D Moore heat equation becomes ∆tu=α∇2 Mu. 5. Moore Wave Equation The Moore wave equation is defined as δt∆tu=c2∇2 Mu. 2D case. δt∆tu=c2δx∆xu+δy∆yu+δx∆xδy∆yu. 3D case. δt∆tu=c2∇2 Mu, with ∇2 Mgiven by the full Moore coupling structure above. 96 10.6 List of Steady State Problems Steady state problems arise naturally from discrete evolution equations by imposing the condition ∆tu= 0. In this section we summarize several fundamental steady state equations in the theory of partial difference equations. 1. Discrete Laplace Equation The discrete Laplace equation is defined by ∇2u= 0. Explicitly: •1D: δx∆xu= 0. •2D: δx∆xu+δy∆yu= 0. •3D: δx∆xu+δy∆yu+δz∆zu= 0. 2. Discrete Poisson Equation The discrete Poisson equation is the discrete analogue of ∇2u=−f: ∇2u=−f. Explicit forms: •1D: δx∆xu=−f(x). •2D: δx∆xu+δy∆yu=−f(x, y). •3D: δx∆xu+δy∆yu+δz∆zu=−f(x, y, z). 97 3. Moore Laplace Equation Using the Moore neighbourhood, the Moore Laplacian satisfies ∇2 Mu= 0. Explicitly: •2D: δx∆xu+δy∆yu+δx∆xδy∆yu= 0. •3D: ∇2u+δx∆xδy∆yu+δy∆yδz∆zu+δz∆zδx∆xu+δx∆xδy∆yδz∆zu= 0. 4. Moore Poisson Equation The Moore Poisson equation is defined by ∇2 Mu=−f. Explicit forms: •2D: δx∆xu+δy∆yu+δx∆xδy∆yu=−f(x, y). •3D: ∇2u+δx∆xδy∆yu+δy∆yδz∆zu+δz∆zδx∆xu+δx∆xδy∆yδz∆zu=−f(x, y, z). 5. Discrete Biharmonic Equation The discrete biharmonic equation is defined by ∇4u= 0,∇4:= (∇2)2. Explicitly: •1D: ∇4u= (∇2)2u=δ2 x∆2 xu= 0. •2D: δ2 x∆2 xu+δ2 y∆2 yu+ 2 δx∆xδy∆yu= 0. •3D: ∇4u= (δx∆x+δy∆y+δz∆z)2u= 0. These steady state problems may be solved using separation of variables, eigenvalue methods, and discrete Fourier series expansions. 98 11 Discrete Evolution Equations 11.1 Introduction In this chapter, we develop the theoretical framework for discrete spatiotemporal dynamical systems and introduce the concept of discrete evolution equations. Our aim is to establish a rigorous parallel between continuous evolution equations, which arise in the study of partial differential equations (PDEs), and their fully discrete counterparts defined on the lattice Zn. We begin by defining discrete dynamical systems and formulating discrete evolution equations in operator form. The connection to semigroup theory is then discussed, providing the algebraic foundation for the time evolution of such systems. Subsequently, we classify well-posed problems into three main types: initial value problems,boundary value problems, and initial-boundary value problems. These settings serve as the natural discrete analogues of the classical PDE theory. Finally, we introduce the distinction between autonomous and non-autonomous systems, depending on whether the update operator depends explicitly on the independent variables. Together, these elements form a comprehensive theory of discrete evolution equations, laying the groundwork for subsequent analysis and applications. 11.2 Definitions Definition 11.1 (Discrete Dynamical System).Let Xbe a set (typically a metric space, topological space, or Banach space). A discrete dynamical system is a pair (X, φ)where φ:X→X is a mapping (often continuous if Xhas a topology). The dynamics is defined by iteration: un+1 =φ(un), n ∈Z≥0, u0∈X. Equivalently, for each n∈Z≥0, we define the n-th iterate φn(x) := φ◦φ◦···◦φ | {z } ntimes (x), so that the trajectory of x∈Xis given by {φn(x) : n∈Z≥0}. [9] Definition 11.2 (Discrete Spatiotemporal Dynamical System).Let u:N0×Zn→C denote the state of the system, where 99 •t∈N0is the discrete time, •x∈Znis the discrete spatial coordinate. Adiscrete spatiotemporal dynamical system is governed by an equation of the form Etu(t, x) = F{Em tEk xu(t, x)}(m,k)∈S, t, x, where •Em tu(t, x) = u(t+m, x)is the time–shift operator, •Ek xu(t, x) = u(t, x+k)is the spatial shift operator, •S⊂Zn+1 is a finite index set specifying which time–space shifts appear, •Fis a prescribed function, possibly nonlinear. Definition 11.3 (Discrete Evolution Equation).Let Xbe a Banach space, and let A:X→Xbe a linear operator. A discrete evolution equation is an iterative relation of the form Etu=Au +F(u, t), that is, u(t+ 1) = Au(t) + F(u(t), t), t ∈Z≥0, where •u:Z≥0→Xis the unknown sequence (or discrete trajectory), •A:X→Xis a linear operator, representing the linear part of the dynamics, •F:X×Z≥0→Xis a linear or nonlinear mapping, representing the forcing or nonlinear interaction term. An initial condition u(0) = u0∈X is prescribed, and the discrete evolution equation determines the trajectory {u(t)}t≥0. Example 11.1 (Linear Evolution Equation: Fibonacci Equation).Consider the recurrence relation Etx=x+E−1 tx, where x=x(t)and x:Z→R. This is a linear ordinary difference equation of order 2, commonly known as the Fibonacci equation. Example 11.2 (Nonlinear Evolution Equation: Logistic Map).Consider the nonlinear recurrence Etx=rx(1 −x), where x=x(t)and x:Z→R. This is a nonlinear ordinary difference equation, famously known as the logistic map, which plays a central role in the study of chaos theory. 100 [27] Example 11.3 (Discrete Evolution Equation: Vector Equation).Consider the discrete scalar field u=u(t, x)with t∈N0,x∈Z. The evolution is governed by the equation Etu=u+E−1 xu+Exu+E−1 tu+E−2 tu. Define auxiliary variables v(t, x) := E−1 tu(t, x) = u(t−1, x), w(t, x) := E−2 tu(t, x) = u(t−2, x). Then the system can be rewritten as Etu=u+E−1 xu+Exu+v+w, Etv=u, Etw=v. Introducing the vector-valued function u=u(t, x) :=  u(t, x) v(t, x) w(t, x) , the system takes the compact operator form Etu=Au, where Ais a linear operator acting on udefined by A u v w = u+E−1 xu+Exu+v+w u v . Example 11.4 (Nonlinear Evolution Equation: Rule 110).Consider the onedimensional cellular automaton Rule 110, which can be expressed as a nonlinear partial difference equation: Etu= mod2u+Exu+u Exu+u Exu E−1 xu, where u:N0×Z→Z2, u =u(t, x). This system is well-known for its Turing completeness. It can be regarded as a nonlinear evolution equation of the form Etu=Au +F(u, t), with A= 0, F(u, t) = mod2u+Exu+u Exu+u Exu E−1 xu. 101 [34] Example 11.5 (Langton’s Ant as a System of Coupled Difference Equations). Consider the following system: Etu=u+ (1 −2u)·δ(x−X)δ(y−Y), Etd= mod4d+ (1 −2u(t, X, Y )), EtX=X+ cosπ 2·Etd, EtY=Y−sinπ 2·Etd, where: •u=u(t, x, y)∈ {0,1}is the state of the lattice at time tand position (x, y), i.e. u:Z3→ {0,1}. •d=d(t)∈ {0,1,2,3}is the ant’s direction at time t, i.e. d:Z→Z4. •(X(t), Y (t)) ∈Z2is the position of the ant, with X, Y :Z→Z. •The modulo operator is defined as mod4(x) = x−4·x 4. This yields a coupled nonlinear evolution system consisting of: •One partial difference equation for the lattice state u(t, x, y). •Three ordinary difference equations for the ant’s internal state: the direction d(t), and the position (X(t), Y (t)). [7] 102 11.3 Semigroup Theory Definition 11.4 (Semigroup).Let Xbe a set. A semigroup is a pair (X, ·) where ·:X×X→Xis a binary operation satisfying the associativity property: (x·y)·z=x·(y·z),∀x, y, z ∈X. If there exists an element e∈Xsuch that e·x=x·e=x, ∀x∈X, then (X, ·)is called a monoid, and eis called the identity element. [31] Definition 11.5 (Discrete Operator Semigroup).Let Xbe a Banach space and T:X→Xa bounded linear operator. The family {Tn}n∈N0defined by Tn:= T◦T◦···◦T | {z } ntimes , T0:= I, is called a discrete operator semigroup generated by T. Theorem 11.6 (Solution of Linear Discrete Evolution Equation).Let Xbe a Banach space, T:X→Xa bounded linear operator, and u0∈X. Then the solution to Etu=u(t+ 1) = Tu(t), u(0) = u0, is given explicitly by u(t) = Ttu0, t ∈N0. Example 11.6 (Stability of the Discrete Heat Equation Semigroup).Consider the discrete heat equation Etu=u+α∇2u, u(0, x) = f(x)∈ L2(Z), where ∇2u=u(t, x + 1) −2u(t, x) + u(t, x −1) is the discrete Laplacian. Define the evolution operator T=I+α∇2:L2(Z)→ L2(Z). Then: 1. Tis a bounded linear operator. Therefore {Tt}t∈N0forms a discrete operator semigroup, and the solution always exists, given by u(t) = Ttf, t ∈N0. 2. The Fourier symbol of Tis b T(k)=1−4αsin2k 2, k ∈[−π, π]. 3. The system is stable (i.e. ∥Ttf∥L2≤ ∥f∥L2for all t) if and only if 0< α ≤1 2. If α > 1 2, then there exist modes that grow exponentially, and the solution explodes. 103 11.4 Initial Value Problems Definition 11.7 (Initial Value Problem for Discrete Evolution Equations).Let u=u(t, x)denote the state of the system, where (t, x)∈N0×Zn, with t∈N0representing discrete time and x∈Znrepresenting discrete space. Adiscrete initial value problem is given by Etu=Au+F(u, t, x), subject to the initial condition u(0,x) = f(x),x∈Zn, with no boundary conditions imposed. Here Ais a linear operator and Fis a (possibly nonlinear) function. Example 11.7 (Initial Value Problem: Rule 90 Cellular Automaton).Consider the Rule 90 cellular automaton, expressed as a partial difference equation: Etu= mod2E−1 xu+Exu, where u:Z2→Z2. The initial condition is given by u(0, x) = f(x)∈B(Z), with no boundary condition imposed. If we take f(x) = δ(x), where δ(x)is the Kronecker delta, then the explicit solution is u(t, x) = mod2C(2t, x +t), where C(n, k) =    n! k!(n−k)!,0≤k≤n, 0,otherwise. This solution generates the well-known Sierpinski triangle. [34] 104 11.5 Boundary Value Problems Definition 11.8 (Discrete Boundary via Level Set Function).Let φ:Zn→R be a discrete level set function defined on the integer lattice. We define the corresponding discrete domain Ω⊆Znby Ω = {x∈Zn:φ(x)<0}, and its complement by Ωc={x∈Zn:φ(x)>0}. Then the discrete boundary of Ωis defined as the set of lattice points where the level set function vanishes: ∂Ω = {x∈Zn:φ(x)=0}. [12] Definition 11.9 (Abstract Discrete Boundary Condition).Let Ω⊆Znbe a discrete domain with boundary ∂Ωdefined via a level set function φ:Zn→ R. A discrete boundary condition for a partial difference equation on Ωis expressed abstractly as B(u(t, x)) = 0,x∈∂Ω, where Bis a boundary operator acting on the discrete function u:N0×Zn→C. Types of Boundary Conditions Definition 11.10 (Discrete Dirichlet Boundary Condition).Let Ω⊆Znbe a discrete domain with boundary ∂Ω. A discrete Dirichlet boundary condition prescribes the value of the solution uon the boundary: u(t, x) = g(x),∀x∈∂Ω, t ∈N0, where g:∂Ω→Cis a given boundary function. Definition 11.11 (Discrete Neumann Boundary Condition).Let Ω⊆Znbe a discrete domain with boundary ∂Ωdefined via a level set function φ:Zn→R. Adiscrete Neumann boundary condition prescribes the normal discrete flux of uon the boundary: ∆nu(t, x) = g(t, x),∀x∈∂Ω, t ∈N0, where ∆nu:= ∇cu·n, n =∇cφ ∥∇cφ∥, with ∇cdenoting the central gradient. 105 Integral Representation of the Solution. Hence the general solution can be expressed in Fourier form as u(t, x, y, z) = 1 (2π)3Zπ −πZπ −πZπ −πb f(k, m, n) (1+e−ik+e−im+e−in)tei(kx+my+nz)dk dm dn, where b f(k, m, n) denotes the discrete Fourier transform of the initial data: b f(k, m, n) = X x,y,z∈Z f(x, y, z)e−i(kx+my+nz). Green’s Function. For the delta initial condition f(x, y, z) = δ(x)δ(y)δ(z), we obtain the Green’s function G(t, x, y, z) = 1 (2π)3Zπ −πZπ −πZπ −π (1+e−ik+e−im+e−in)tei(kx+my+nz)dk dm dn. Expanding the multinomial term and evaluating the integral, we find G(t, x, y, z) = C(t, x, y, z) where C(t, x, y, z) =      t! (t−x−y−z)! x!y!z!,if 0 ≤x, y, z ≤tand x+y+z≤t, 0,otherwise. which is the multinomial coefficient, provided all entries are nonnegative, and 0 otherwise. Convolution Form of the Solution. Therefore, the solution can also be written in convolution form: u(t, x, y, z) = X s,p,q∈Z f(s, p, q)Ct, x −s, y −p, z −q. Fourier–Multinomial Identity. This yields the identity 1 (2π)3Zπ −πZπ −πZπ −π (1+e−ik+e−im+e−in)tei(kx+my+nz)dk dm dn =t t−x−y−z, x, y, z . In other words, the Green’s function of the 3D Pascal Evolution Equation is exactly given by multinomial coefficients. 112 12.4 n-Dimensional Pascal Evolution Equation We now propose a general linear partial difference equation in nspatial dimensions, extending the binomial and trinomial cases. Equation. Etu=u+E−1 x1u+E−1 x2u+···+E−1 xnu, u :Zn+1 →R. Initial condition: u(0,x) = f(x)∈ L2(Zn), with no boundary conditions. Fourier Ansatz. We apply the Fourier ansatz u(t, x) = eλtei(k1x1+···+knxn). Substitution gives the dispersion relation eλ= 1 + e−ik1+e−ik2+···+e−ikn. Integral Representation. Thus the general solution has the Fourier integral form u(t, x) = 1 (2π)nZ[−π,π]nb f(k)1 + e−ik1+···+e−ikntei(k1x1+···+knxn)dk, where b f(k) is the discrete Fourier transform of f. Green’s Function. For the delta initial condition f(x) = δ(x), we obtain G(t, x) = 1 (2π)nZ[−π,π]n1 + e−ik1+···+e−ikntei(k1x1+···+knxn)dk. This integral can be identified with the multinomial coefficient in discrete form: G(t, x) = C(t, x) =      t! (t−Pn j=1 xj)! x1!···xn!,if Pn j=1 xj≤tand xj≥0, 0,otherwise. Convolution Form of the Solution. Therefore the general solution can also be written in convolution form: u(t, x) = X s∈Zn f(s)Ct, x−s. This establishes the n-dimensional Pascal Evolution Equation, where the Green’s function is given by the multinomial coefficient, generalizing the binomial and trinomial cases. 113 Theorem 12.1 (Fourier–Multinomial Identity).Let k= (k1, . . . , kn)∈[−π, π]n and x= (x1, . . . , xn)∈Zn. For t∈Nwe have 1 (2π)nZ[−π,π]n 1 + n X j=1 e−ikj  t ei(k·x)dk=C(t, x), where C(t, x)is the multinomial coefficient C(t, x) =      t! (t−Pn j=1 xj)! x1!···xn!,if Pn j=1 xj≤tand xj≥0, 0,otherwise. 114 13 Nonlinear Equations with mod n Nonlinearity 13.1 Introduction In this section, we focus exclusively on nonlinear partial difference equations with mod nnonlinearity. In particular, we demonstrate that several classical self–similar fractals can be realized as exact solutions of such evolutionary systems. We propose and analyze partial difference equations corresponding to three well–known fractals: •the Sierpinski triangle, •the Sierpinski carpet, •the Sierpinski pyramid. For each case, we derive the analytic form of the solution by applying the linear Green’s function representation, followed by reduction modulo n. This approach reveals that these fractals are not merely geometric constructs defined by iterated function systems (IFS), but can also be understood as solutions to explicitly defined nonlinear evolution equations. 13.2 Right–Side Sierpinski Triangle Equation We now consider the nonlinear partial difference equation Etu= mod2u+E−1 xu, u :Z2→R. Initial Condition. We impose u(0, x) = f(x)∈ L2(Z) with no boundary conditions. Relation to the Linear Case. From the previous section, the linear equation Etu=u+E−1 xu admits the convolution solution u(t, x) = X s∈Z f(s)C(t, x −s), where C(t, y) =   t y,0≤y≤t, 0,otherwise. 115 Nonlinear Modular System. For the modular system, the solution becomes u(t, x) = mod2 X s∈Z f(s)C(t, x −s)!. Delta Initial Condition. If f(x) = δ(x), then u(t, x) = mod2C(t, x). This evolution produces a tilted Sierpinski triangle pattern. Spatiotemporal Patterns. •For f(x) = δ(x), the evolution yields a perfect self–similar fractal structure (Sierpinski triangle). •For random initial conditions with f(x)∈ {0,1}, the dynamics still exhibit fractal–like triangular patterns. •For real–valued random initial conditions f(x)∈R, the system displays spatiotemporal chaos. 116 Figure 7: Spatiotemporal plot for δinitial condition. 117 Figure 8: Random initial condition with f(x)∈ {0,1}. 118 Figure 9: Random initial condition with f(x)∈R. 13.3 Sierpinski Carpet Equation A collaborator of mine, Wu Han (China), proposed the following nonlinear partial difference equation with modular reduction: Etu= mod3E−1 xu+Exu+E−1 tu, u :Z2→R. Initial Condition. We impose two delta initial conditions u(0, x) = u(1, x) = δ(x), with no boundary conditions. Numerical Observation. 119 Figure 10: The spatiotemporal plot of this system reveals a striking fractal pattern: the evolution generates a Sierpinski carpet. The precise mathematical reason for this emergence is not yet understood, but the numerical evidence is compelling. Relation to the Linear Case. From the previous section, the corresponding linear system without the modular reduction has the Green’s function K(t, x) = 1 2πZπ −π (1 −z−)zt ++ (z+−1) zt − z+−z− eikx dk, with characteristic roots z±=z±(k) = cos(k)±p1 + cos2(k). Modular Reduction. For the nonlinear modular system, the solution with delta initial conditions is given by u(t, x) = mod3K(t, x). Fractal Emergence. Remarkably, this modular reduction transforms the oscillatory integer structure of K(t, x) into a perfect self–similar fractal, the Sierpinski carpet. 13.4 Sierpinski Pyramid Equation We propose the nonlinear partial difference equation Etu= mod2u+E−1 xu+E−1 yu, u :Z3→R. Initial Condition. u(0, x, y) = f(x, y)∈ L2(Z2), 120 with no boundary conditions. Solution. From the previous section we already derived the solution for the linear version. For the nonlinear equation, the solution is u(t, x, y) = mod2 X (s,p)∈Z2 f(s, p)C(t, x −s, y −p) , where C(t, x, y) =      t! (t−x−y)! x!y!,0≤x, y, x +y≤t, 0,otherwise. In the special case f(x, y) = δ(x)δ(y), the solution reduces to u(t, x, y) = mod2C(t, x, y), which generates the well–known Sierpinski Pyramid structure. 121 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. 128 13.8 Fractals as Solutions to Evolution Equations Traditionally, many well–known fractals such as the Sierpinski triangle, the Sierpinski carpet, and the Sierpinski pyramid are generated via Iterated Function Systems (IFS) [14, 15] or purely geometric constructions. In this work, however, we have observed that these fractals can also emerge as exact solutions to suitable evolution equations. This provides a novel perspective, revealing deep connections between difference equations and fractal geometry. In particular, we noted that the Green’s functions of several linear partial difference equations coincide with classical combinatorial numbers (e.g. binomial, trinomial, and multinomial coefficients). When these coefficients are reduced modulo a small integer, the resulting solutions display exact fractal patterns. Hence, partial difference equations naturally unify Fourier analysis, combinatorics, fractal geometry, and chaos theory. Furthermore, the linear equations without modular reduction produce solutions that grow indefinitely, corresponding to a “stretching” mechanism. The modular reduction acts as a “folding” step that forces the solution values back into a bounded range. This interplay of stretching and folding is precisely the central mechanism responsible for chaotic dynamics. [2] We also observed that for many nonlinear or chaotic systems, when the initial condition or forcing term is a delta function, the solution is a perfect self–similar fractal. On the other hand, if the initial condition is random, the solution evolves into spatiotemporal chaos. From the convolution representations, it follows that the nonlinear solutions can be expressed as superpositions of evolution kernels, and we conjecture that some spatiotemporal chaos can be understood as the nonlinear superposition of many fractal kernels. 129 14 Oscillation Theory 14.1 Motivation In the theory of partial differential equations, a central theme is regularity theory. This emphasis arises from several fundamental features of PDEs: solutions are required to be continuous (or even smooth) functions, and differential operators are typically unbounded. As a result, a large portion of PDE analysis is devoted to questions of existence, uniqueness, and regularity of solutions in appropriate function spaces. The situation is markedly different in the setting of partial difference equations. In this discrete framework, the underlying functions are defined on lattices and are inherently discontinuous. There is no notion of smoothness in the classical sense, and difference operators are bounded operators on natural discrete function spaces. Consequently, many of the central difficulties that motivate regularity theory in PDEs are absent for partial difference equations. For this reason, we argue that regularity theory should not play the same foundational role in the study of partial difference equations and discrete dynamical systems. Instead, the primary phenomena of interest are the oscillatory behaviors of solutions. In the classical theory of one-dimensional difference equations, as developed for example by Elaydi, oscillation theory typically refers to the study of signchanging behavior of solutions. While this approach is effective for scalar or lowdimensional recurrence relations, it is insufficient for partial difference equations, whose behavior is more naturally interpreted as that of infinite-dimensional dynamical systems. Partial difference equations exhibit not only temporal oscillations, but also spatial oscillations and complex spatiotemporal patterns. These phenomena cannot be adequately captured by sign-based oscillation criteria alone. To address this richness of behavior, we propose to adopt ideas from harmonic analysis, viewing solutions through their frequency content, multiscale decompositions, and spectral representations. By decomposing discrete fields into frequencies, wave packets, or wavelets, one can analyze oscillations across different spatial and temporal scales. This perspective provides a natural and powerful framework for understanding the qualitative dynamics of partial difference equations and discrete field theories, beyond the scope of classical regularity-based methods. 14.2 Brillouin Zone and Nyquist Frequency Let us consider functions defined on the integer lattice, x∈Z. For a discrete Fourier mode of the form u(x) = eikx, 130 we observe that ei(k+2π)x=ei2πxeikx =eikx, x ∈Z, since ei2πx = 1 for all integers x. This shows that frequencies differing by integer multiples of 2πare indistinguishable on the integer lattice. In other words, the frequency spectrum is periodic with period 2π, and frequencies outside a fundamental interval are folded back (aliasing). For this reason, the interval [−π, π] is chosen as the fundamental frequency domain, and is called the Brillouin zone. A fundamental consequence of discretization is that the frequency spectrum is bounded. On the integer lattice, the largest possible magnitude of frequency is π. This corresponds to the fastest possible oscillation that a discrete function can exhibit. As a simple example, consider the function f(x) = cos(πx)=(−1)x, x ∈Z. This function alternates sign at every lattice point, representing the maximum possible oscillation rate in space. Written in trigonometric form, its frequency is exactly π. This maximal frequency is known as the Nyquist frequency. No discrete function defined on Zcan oscillate faster than this frequency without being aliased to a lower one. 14.3 Fourier Analysis 14.4 Wavelet Analysis 14.5 Discrete Conservation Law 131 15 System of Partial Difference Equations 15.1 Introduction In the preceding chapters, we focused primarily on scalar partial difference equations, that is, equations governing a single discrete field u:N0×Zn→R. Such equations suffice for modeling isolated diffusion processes, simple transport dynamics, or single–species population models. However, many natural and physical systems cannot be described by a single scalar equation. Phenomena in fluid mechanics, electromagnetism, reaction– diffusion systems, population interactions, and complex systems often require the simultaneous evolution of multiple interacting components. These interactions naturally lead to coupled equations, or equivalently, to the evolution of a vector–valued discrete field u= (u1, . . . , us) : N0×Zn→Rs. This chapter is devoted to the systematic formulation and analysis of systems of partial difference equations (P∆E systems). We introduce the general operator framework, define linear, semilinear, quasilinear, and fully nonlinear systems, and establish a unified notation for discrete vector fields and discrete differential operators. These systems form the natural discrete analogue of vector PDEs in the continuous setting and provide a general language for describing a wide range of coupled dynamical phenomena. 15.2 Classification of P∆E Systems In this section we introduce a systematic classification scheme for systems of partial difference equations (P∆E). Let x= (x1, x2, . . . , xn)∈Zn, k = (k1, k2, . . . , kn)∈Zn, and define the multi-indexed shift operator Ek:= Ek1 x1Ek2 x2···Ekn xn, together with the L1–norm |k|1:= |k1|+|k2|+···+|kn|. We consider vector-valued discrete fields u= (u1, u2, . . . , us) with u:N0×Zn→Rs, and write all systems in evolution form, Etu=Au +F(u, t, x). 132 Definition 15.1 (Linear System).A system is called linear if it is of the form Etu=A(t, x)u+F(t, x), where A(t, x)is an s×smatrix of shift operators with entries Aij(t, x) = X |k|1≤m aijk(t, x)Ek,1≤i, j ≤s. No nonlinear dependence on uoccurs. Definition 15.2 (Semilinear System).A system is called semilinear if Etu=A(t, x)u+F(u, t, x), where A(t, x)is as in the linear case, and the nonlinear term F(u, t, x)satisfies: no shift operators may appear in F. Thus Fmay be nonlinear in u, but must depend only on the pointwise value u(t, x)and not on its shifted versions. Definition 15.3 (Quasilinear System).A system is called quasilinear if Etu=A(t, x, u)u+F(u, t, x), where the entries of Aare given by Aij(t, x, u) = X |k|1≤m aijk(t, x, u)Ek, with the restriction: A(t, x, u)may depend nonlinearly on u, but no shift operators may appear inside A. Equivalently, the only admissible nonlinear transport-type terms are of the form f(u)Exu. Definition 15.4 (Fully Nonlinear System).A system is called fully nonlinear if Etu=F(Au, u, t, x), where Fmay depend nonlinearly on shifted versions of u. Typical examples include Exu Eyu, (Exu)2,pu+Exu, etc. In this case, shift operators appear inside the nonlinear term in genuinely nonlinear combinations. 133 Remark. This classification is fully parallel to the classical PDE classification. Indeed, replacing each shift operator Ekwith the differential operator ∂k=∂k1 x1∂k2 x2···∂kn xn, yields the corresponding notions of linear, semilinear, quasilinear, and fully nonlinear partial differential equation systems. Example 15.1 (Reformulation of a Second-Order P∆E into a First-Order System).Consider the three-dimensional discrete wave equation δt∆tu=c2∇2u, u :N0×Z3→R. Recall that δt∆tu=Etu+E−1 tu−2u. Therefore the equation can be rewritten in pure shift-operator form as Etu=c2∇2u+ 2u−E−1 tu. Introduce a new auxiliary variable v:= E−1 tu. Then the second-order equation becomes the coupled first-order system (Etu=c2∇2u+ 2u−v, Etv=u. This shows that any second-order P∆E in time can be equivalently formulated as a first-order system in an extended state vector. Example 15.2 (Quasilinear System: Discrete Navier–Stokes Equations).Consider the discrete analogue of the Navier–Stokes equations, which we refer to as the Discrete Navier–Stokes Equations (DNSE): ρ∆tu+u·∇cu=−∇cp+µ∇2u+f,∇c·u= 0, where the discrete vector field u= (u, v, w) = u(t, x, y, z),u:N0×Z3→R3, and the operators ∆t,∇c, and ∇2denote respectively the forward time difference, the central discrete gradient, and the discrete Laplacian on the cubic lattice Z3. Rewriting the equation in evolution form using Etu=u+ ∆tu, we obtain Etu=u−1 ρ∇cp+µ ρ∇2u−u·∇cu+1 ρf. 134 This system is quasilinear because the discrete gradient ∇cuappearing in the convective term is multiplied by the solution uitself: u·∇cu=u∆c xu+v∆c yu+w∆c zu, which is linear in the discrete derivatives of ubut nonlinear in u. Thus the DNSE fits precisely into the framework of quasilinear partial difference equation systems. Example 15.3 (Fully Nonlinear System: Forest Fire Model).Consider the two-dimensional Forest Fire Model, which can be written as Etu= mod3(2 δ(u−1) θ(S−1) + G(t, x, y)) , where the neighbourhood interaction term is S=X (i,j)∈M δEi xEj yu−2, and u:N0×Z2→Z3,Z3={0,1,2}. The states are interpreted as: 0 = empty,1 = tree,2 = burning tree. The forcing term G:N0×Z2→ {0,1} models external phenomena such as spontaneous tree growth or lightning strikes. This system is fully nonlinear because the highest-order spatial operator δEi xEj yu−2 depends nonlinearly on the shifted values of u, and therefore the nonlinearity appears inside the operator itself. 135 15.3 Linear Systems In this section, we study the class of linear partial difference equation systems of the form Etu=A(t, x)u+F(t, x), where u:N0×Zn→Rsis a vector-valued function, and A(t, x) is an s×s matrix whose entries are finite linear combinations of shift operators: Aij(t, x) = X ∥k∥1≤m aijk(t, x)Ek, k = (k1, . . . , kn)∈Zn. Here Ek:= Ek1 x1Ek2 x2···Ekn xn,∥k∥1=|k1|+···+|kn|, and u= (u1, u2, . . . , us)T is a discrete vector field. Reduction to Ordinary Difference Equations. If each entry Aij(t, x) contains no shift operator, i.e. Aij(t, x) = aij(t, x), then the system reduces to a system of ordinary difference equations (ODEs). Constant-Coefficient Linear Systems Consider the autonomous system Etu=Au, A ∈Rs×sconstant. The solution is well known and can be expressed using the discrete semigroup {At}t≥0: u(t) = Atu0, u0=u(0). Case 1: Ais diagonalizable. Suppose A=PΛP−1, where Λ = diag(λ1, . . . , λs). Then At=PΛtP−1, and hence the solution is u(t) = P   λt 1... λt s   P−1u0. 136 Case 2: Ais not diagonalizable. Let Ahave Jordan form A=PJP−1, J = diag (J1, J2, . . . , Jr), where each Jordan block is of the form Jℓ=       λℓ1 0 ··· 0 0λℓ1··· 0 . . ........ . . 0··· 0λℓ1 0··· 0 0 λℓ        . Then At=PJtP−1, where each block satisfies Jt ℓ=λt ℓ        1tt 2··· t kℓ−1 0 1 t··· t kℓ−2 . . ........ . . 0··· 0 1 t 0··· 0 0 1        . Thus the full solution is u(t) = PJtP−1u0, with polynomial prefactors t jarising from the nilpotent part of each Jordan block. System of Independent Equations If the coefficient matrix Ais diagonal, then the system Etu=Au contains no coupling between the components of u. In other words, each component evolves independently, and the system may be solved by treating the equations one by one. Example 15.4 (Discrete Vector Heat Equation).Consider the discrete vector heat equation ∆tu=α∇2u,u= (u, v, w), where u:N0×Z3→R3. Using the identity ∆t=Et−I, the system can be rewritten in evolution form: Etu=u+α∇2u, Etv=v+α∇2v, Etw=w+α∇2w. 137 By induction, we obtain the discrete Duhamel formula: u(t, x) = Atu0(x) + t−1 X s=0 At−1−sF(u(s), s, x). This representation is an implicit formulation of the semilinear system and will be fundamental in the study of existence, uniqueness, stability, and long– time dynamics. 144 16 Atlas of Nonlinear Partial Difference Equations 16.1 Introduction This chapter presents an atlas of nonlinear partial difference equations (P∆Es) that exhibit rich dynamical behavior, including pattern formation, spatiotemporal chaos, and a wide variety of fractal structures. The purpose of this chapter is twofold: 1. to illustrate how classical discrete models—such as cellular automata, sandpile dynamics, and lattice-based evolution rules—fit naturally into the unified framework of nonlinear P∆Es; and 2. to showcase several novel models introduced by the author, demonstrating the expressive power of P∆Es in generating complex and emergent behavior from simple algebraic rules. The examples in this atlas range from elementary one-dimensional evolution rules to high-dimensional systems with nonlinear couplings and mod-ninteractions. Despite their diversity, all models can be written compactly in the general form Etu=F(u, x, t), where Fmay involve nonlinear combinations of shift operators, local interactions, threshold or mod-nnonlinearities, or discrete approximations of classical differential operators. Many of the systems presented here produce striking phenomena: self-similar patterns, recursively generated fractals, intermittency, spatiotemporal turbulence, and long-range correlations. Through these examples, we highlight the central theme of this monograph: nonlinear partial difference equations provide a natural language for describing complexity in discrete spacetime. This atlas serves both as a reference and as a source of inspiration for future research on nonlinear discrete dynamics, offering a broad view of how simple algebraic rules can give rise to highly nontrivial emergent structures. 145 16.2 Coupled Map Lattice 1D Coupled Map Lattice The one-dimensional Coupled Map Lattice (CML), originally proposed by Kaneko [35], 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. 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. 146 Explanation Acoupled map lattice (CML) is a discrete-time, discrete-space dynamical system in which each lattice site evolves according to a prescribed local map and interacts with other sites via a coupling rule. The local map encodes the nonlinear dynamics at each site, while the coupling represents spatial interactions such as diffusion, transport, or synchronization. CMLs are capable of producing spatiotemporal chaos, where complex, irregular patterns emerge and evolve across both space and time. They occupy an intermediate position between continuous partial differential equations and fully discrete cellular automata: like PDE, they describe the evolution of a field; like CA, they are inherently discrete in space and time. In the context of this work, CMLs can be regarded as a particular subclass of partial difference equations in which the update operator contains a nonlinear local component together with a discrete coupling term. This perspective allows the mathematical machinery of discrete functional analysis to be applied directly to the study of CML dynamics. 147 16.3 Elementary Cellular Automata The study of one-dimensional cellular automata (1D CA) was significantly advanced by Stephen Wolfram [34], 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. 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. 148 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. 149 Example 1 This is the spatiotemporal plot of Rule 90 with the initial condition u(0, x) = δ(x) and no boundary conditions. Figure 16: 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. 150 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 17: 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. 151 Example 3 This is the spatiotemporal plot of Rule 90 with random initial condition and boundary conditions u(t, −200) = u(t, 200) = 0. Figure 18: Spatiotemporal plot of Rule 90 with random initial condition and boundary conditions. As shown, the system exhibits spatiotemporal chaos. 152 Rule 30 The evolution equation for Rule 30, derived using Boolean algebra, is given by: Etu= mod2E−1 xu+Exu+u+u·Exu, where Etis the time shift operator, Exis the spatial right shift operator. Example We consider the following simulation 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 shows the spatiotemporal evolution of the system under Rule 30: It is evident that the system exhibits spatiotemporal chaos, characterized by aperiodic and unpredictable patterns across both space and time. 153