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 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 (Laurent 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. 1 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 6 1 Introduction to Partial Difference Equations 7 1.1 Motivation and Background . . . . . . . . . . . . . . . . . . . . . 7 1.2 DiscreteFunction........................... 8 1.3 Operators............................... 8 1.4 Ordinary Difference Equations . . . . . . . . . . . . . . . . . . . 9 1.5 Partial Difference Equations . . . . . . . . . . . . . . . . . . . . . 10 1.6 The Order of a Partial Difference Equation . . . . . . . . . . . . 10 1.7 Classification and Examples . . . . . . . . . . . . . . . . . . . . . 12 1.8 Notations ............................... 13 1.9 Notation Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2 Introduction to Linear Partial Difference Equations 15 2.1 Definition ............................... 15 2.2 Discrete 1D Transport Equation . . . . . . . . . . . . . . . . . . 16 2.3 3D Discrete Transport Equation . . . . . . . . . . . . . . . . . . 18 2.4 Existence and Uniqueness of Solutions to P∆E . . . . . . . . . . 19 3 Discrete Function Spaces and Operators 20 3.1 Introduction.............................. 20 3.2 Discrete Function Spaces . . . . . . . . . . . . . . . . . . . . . . . 20 3.3 TypesofOperators.......................... 25 3.4 OperatorTheory ........................... 30 4 Discrete Functionals and Convergence 39 4.1 Introduction.............................. 39 4.2 Discrete Functionals . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.3 Compactness of the Unit Ball in Discrete Function Spaces . . . . 40 4.4 Hahn–Banach Theorem . . . . . . . . . . . . . . . . . . . . . . . 41 4.5 Riesz Representation Theorem . . . . . . . . . . . . . . . . . . . 42 4.6 Banach Fixed Point Theorem . . . . . . . . . . . . . . . . . . . . 43 4.7 Types of Convergence . . . . . . . . . . . . . . . . . . . . . . . . 43 5 The Discrete Green’s Function 45 5.1 The Kronecker Delta Function . . . . . . . . . . . . . . . . . . . 45 5.2 Fundamental Solution . . . . . . . . . . . . . . . . . . . . . . . . 46 5.3 The Discrete Green’s Function . . . . . . . . . . . . . . . . . . . 47 6 Discrete Fourier Transform 49 6.1 Motivation and Background . . . . . . . . . . . . . . . . . . . . . 49 6.2 Discrete Fourier Series . . . . . . . . . . . . . . . . . . . . . . . . 49 6.3 Inner Product and Orthogonality . . . . . . . . . . . . . . . . . . 49 6.4 Discrete Fourier Transform . . . . . . . . . . . . . . . . . . . . . 50 6.5 Parseval’s Identity . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3 6.6 Justification of the Fourier Ansatz for Linear P∆E . . . . . . . . 52 6.7 Classification of Second Order Linear P∆E . . . . . . . . . . . . 52 6.8 Discrete Fourier Analysis and Hadamard’s Well-posedness . . . . 55 7 First Order Equations in Time 57 7.1 Introduction.............................. 57 7.2 1D Pascal Evolution Equation . . . . . . . . . . . . . . . . . . . . 57 7.3 1D Discrete Heat Equation . . . . . . . . . . . . . . . . . . . . . 60 7.4 Stirling Second Equation . . . . . . . . . . . . . . . . . . . . . . . 63 8 Second Order Equations in Time 65 8.1 Introduction.............................. 65 8.2 Second–Order Pascal Evolution Equation . . . . . . . . . . . . . 65 8.3 Discrete 1D Wave Equation . . . . . . . . . . . . . . . . . . . . . 68 9 Steady State Problems 71 9.1 Introduction.............................. 71 9.2 2D Discrete Laplace Equation . . . . . . . . . . . . . . . . . . . . 72 9.3 2D Discrete Poisson Equation . . . . . . . . . . . . . . . . . . . . 74 9.4 Two-Dimensional Moore Laplace Equation . . . . . . . . . . . . . 76 9.5 List of Discrete Evolution Equations . . . . . . . . . . . . . . . . 79 9.6 List of Steady State Problems . . . . . . . . . . . . . . . . . . . . 81 10 Discrete Evolution Equations 83 10.1Introduction.............................. 83 10.2Definitions............................... 83 10.3SemigroupTheory .......................... 88 10.4 Initial Value Problems . . . . . . . . . . . . . . . . . . . . . . . . 89 10.5 Boundary Value Problems . . . . . . . . . . . . . . . . . . . . . . 90 10.6 Initial-Boundary Value Problems . . . . . . . . . . . . . . . . . . 92 10.7 Autonomous and Non–Autonomous Systems . . . . . . . . . . . . 93 11 Higher Dimensional Problems 95 11.1Introduction.............................. 95 11.2 2D Pascal Evolution Equation . . . . . . . . . . . . . . . . . . . . 95 11.3 3D Pascal Evolution Equation . . . . . . . . . . . . . . . . . . . . 96 11.4 n-Dimensional Pascal Evolution Equation . . . . . . . . . . . . . 98 12 Nonlinear Equations with mod n Nonlinearity 100 12.1Introduction..............................100 12.2 Right–Side Sierpinski Triangle Equation . . . . . . . . . . . . . . 100 12.3 Sierpinski Carpet Equation . . . . . . . . . . . . . . . . . . . . . 104 12.4 Sierpinski Pyramid Equation . . . . . . . . . . . . . . . . . . . . 105 12.5 Sierpinski Fan Equation . . . . . . . . . . . . . . . . . . . . . . . 110 12.6 Mod 5 Sierpi´nski Equation . . . . . . . . . . . . . . . . . . . . . . 111 12.7 Mod 4 Sierpinski Fractal Equation . . . . . . . . . . . . . . . . . 112 4 12.8 Fractals as Solutions to Evolution Equations . . . . . . . . . . . 114 13 System of Partial Difference Equations 115 13.1Introduction..............................115 13.2 Classification of P∆E Systems . . . . . . . . . . . . . . . . . . . 115 13.3LinearSystems ............................119 13.4 Constant Coefficient Linear P∆E System . . . . . . . . . . . . . 122 13.5 Non-homogeneous Linear System . . . . . . . . . . . . . . . . . . 125 13.6 Semilinear Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 126 14 Atlas of Nonlinear Partial Difference Equations 128 14.1Introduction..............................128 14.2 Elementary Cellular Automata . . . . . . . . . . . . . . . . . . . 131 14.3 Abelian Sandpile Model . . . . . . . . . . . . . . . . . . . . . . . 140 14.4 Kuramoto Firefly Model . . . . . . . . . . . . . . . . . . . . . . . 143 14.5IsingModel..............................145 14.6 Discrete Logistic Diffusion Equation . . . . . . . . . . . . . . . . 146 14.7 Bacterial Colony Model . . . . . . . . . . . . . . . . . . . . . . . 155 14.8 Discrete Navier–Stokes Equations . . . . . . . . . . . . . . . . . . 161 15 Summary 165 15.1 From Numerical Methods to the Language of Complexity . . . . 165 15.2 Limitations and Future Work . . . . . . . . . . . . . . . . . . . . 166 Acknowledgements 167 References 168 5 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. 6 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,22,24], 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. 7 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 Discrete Function In this section, we introduce the concept of a discrete function, which serves as the fundamental object in our framework for discrete dynamical systems. Definition 1.1 (Discrete Function).Adiscrete function is a mapping f: Ω ⊆Zn→C, where Znis the discrete n-dimensional integer lattice, and Cis the codomain of the function (real, integer, or complex values depending on context). In this paper, we primarily focus on the case where f:Zn→R, i.e., realvalued discrete functions. Example 1.1 (Discrete Single-variable Function).Let f:Z→Rbe defined as f(x) = sin(x). This is a real-valued function defined on the integer lattice Z. Example 1.2 (Discrete Multivariable Function).Let u:Z2→Zbe defined as u(t, x) = ⌊x+ sin(t)⌋. Here, umaps discrete spacetime coordinates (t, x)to an integer value via the floor function applied to a real expression. Definition 1.2 (Discrete Vector Field).Adiscrete vector field on a discrete domain Ω⊆Znis a mapping F: Ω →Rm, that assigns to each discrete point x∈Ωa vector F(x) = (F1(x), F2(x), . . . , Fm(x)). 1.3 Operators Definition 1.3 (Shift Operator).Let x:Z→C. We define the shift operator Ekacting on a discrete function x(t)by Ekx:= x(t+k) for any integer k∈Z. 8 Examples. Ex =x(t+ 1) E2x=x(t+ 2) E−1x=x(t−1) [12] Definition 1.4 (Partial Shift Operator).Given a discrete scalar field u(x1, x2, . . . , xn), the partial shift operator Eki xiacts on uby shifting the i-th coordinate: Eki xiu:= u(x1, . . . , xi+ki, . . . , xn) for any integer ki∈Z. Examples. Etu=u(t+ 1, x, y, z) E2 xu=u(t, x + 2, y, z) Definition 1.5 (Difference Operator).Let x:Z→C. We define the (forward) difference operator ∆by ∆x:= x(t+ 1) −x(t) We may define the shift operator Eas Ex := x(t+ 1), so that the difference operator can be written compactly as ∆x=Ex −x Definition 1.6 (Partial Difference Operator).Let u:Zn→C, and denote u=u(x1, x2, . . . , xn). The partial difference operator with respect to the variable xiis defined as: ∆xiu:= Exiu−u where Exiis the partial shift operator acting on the i-th coordinate: Exiu:= u(x1, . . . , xi+ 1, . . . , xn) 1.4 Ordinary Difference Equations Definition 1.7 (Ordinary Difference Equation of Order k).An ordinary difference equation of order kis an equation involving a single-variable function x:Z→C, expressed in terms of shift operators: F(Enx, En−1x, . . . , Ex, x, E−1x, . . . , E−mx, t)=0, where Eix:= x(t+i),n, m ∈Z+,k=m+nis the total order, and Fis a (possibly nonlinear) function. 9 2.2 Discrete 1D Transport Equation We consider the Discrete Transport Equation: Etu=Ek xu where u:Z2→Ris defined on discrete spacetime. Initial condition: u(0, x) = f(x) We impose no boundary condition. Interpretation: Each value is transported rigidly along a discrete path. The effective velocity is −k, since data moves from x+kto x. Characteristic curve: x(t) = x0+kt Solution We consider the partial difference equation Etu=Ek xu, k ∈Z, with no boundary conditions. Let us define a new variable: ξ=x+kt, and define a new function v(t, ξ) := u(t, x). Then we compute the time-shift of v: Etv=v(t+1, ξ) = v(t+1, x+k(t+1)) = v(t+1, x+kt+k) = v(t+1, ξ+k) = EtEk ξv=EtEk xu Etv=EtEk xu v=Ek xu On the other hand, from the original equation: Etu=u(t+ 1, x) = Ek xu(t, x) = u(t, x +k). Substitute u(t, x) = v(t, ξ), we get: Etu=u(t+ 1, x) = v(t+ 1, ξ) = Etv. 16 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. 17 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). 18 2.4 Existence and Uniqueness of Solutions to P∆E For linear partial difference equations, the question of existence and uniqueness of solutions is considerably simpler than in the continuous case. Since the equations are algebraic recursions on a discrete lattice, solutions can be constructed step by step as long as the recursion is well-defined. In particular, a solution exists and is unique provided the following conditions hold: •The equation is well-defined, i.e. the shift operators involved do not lead to undefined values. •Sufficient initial conditions are given: if the equation has structural order kin the time direction, then one must prescribe u(0,·), u(1,·), . . . , u(k− 1,·) on the spatial domain. •Boundary conditions are specified so that whenever a shifted point falls outside the computational domain, the value of uis still uniquely determined. Under these assumptions, the solution can be constructed inductively in time, and at each step the new layer of values is uniquely determined by previously known layers and the prescribed boundary data. Therefore, unlike the case of partial differential equations, no additional regularity conditions are required to guarantee existence and uniqueness. 19 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,12,21,28]. 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 α: 20 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. [20] 21 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)|. [16] 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. 22 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). 23 Proposition 3.11. The discrete function space F0(Z) := {f:Z→R|supp(f)is finite} is not complete with respect to the L1norm. Proof. Define a sequence {fn}∞ n=1 ⊂F0(Z) by fn(x) := (e−x2,−n≤x≤n, 0,otherwise,x∈Z. Clearly each fnhas finite support, hence fn∈F0(Z). Consider the metric induced by the L1norm: d(f, g) = ∥f−g∥L1=X x∈Z|f(x)−g(x)|. For m>n, we compute ∥fm−fn∥L1=X x∈Z|fm(x)−fn(x)|=X n<|x|≤m e−x2. Since Px∈Ze−x2<∞, the tail sum Pn<|x|≤me−x2→0 as n→ ∞. Thus {fn} is a Cauchy sequence in (F0(Z),∥·∥L1). Its pointwise and L1limit is the function f(x) = e−x2, x ∈Z, which belongs to L1(Z) since Px∈Ze−x2<∞, but f /∈F0(Z) because supp(f) = Zis infinite. Therefore, F0(Z) is not complete. Its completion is L1(Z). Proposition 3.12. The space of finitely supported functions F0(Z) = {f:Z→R|supp(f)is finite} is dense in L1(Z). Proof. Let f∈ L1(Z), i.e. X x∈Z|f(x)|<∞. For each N∈N, define the truncated function fN(x) := (f(x),|x| ≤ N, 0,|x|> N. Clearly fN∈F0(Z), since its support is contained in {−N, −N+ 1, . . . , N}. 24 Now estimate the approximation error: ∥f−fN∥L1=X x∈Z|f(x)−fN(x)|=X |x|>N |f(x)|. Since f∈ L1(Z), the series Px∈Z|f(x)|converges, hence the tail sum P|x|>N |f(x)| → 0 as N→ ∞. Therefore, lim N→∞ ∥f−fN∥L1= 0, which proves that F0(Z) is dense in L1(Z). Definition 3.13 (Integer Valued Function Space).Let Ω⊆Zn. We define the integer valued function space as Q(Ω) = {f: Ω →Z}. Remark 3.14. The integer valued function space Q(Ω) is not a vector space, since scalar multiplication over R(or C) is not defined. Instead, it forms a module over the integer ring Z, where the scalars come from Z. Definition 3.15 (Mod nInteger Valued Function Space).Let Ω⊆Zn. We define the Mod nInteger Valued Function Space as Qn(Ω) = {f: Ω →Zn}, where Zn={0,1,2, ..., n −1}, n ∈Z}denotes the ring of integers modulo n. Definition 3.16 (Boolean Function Space).Let Ω⊆Zn. We define the Boolean Function Space as B(Ω) = Q2(Ω) = {f: Ω →Z2}, where Z2={0,1}denotes the Boolean ring with two elements. Any f∈B(Ω) is called a Boolean function. 3.3 Types of Operators Definition 3.17 (Discrete Operator).Let V, W be discrete function spaces, i.e., subspaces of F(Ω), where Ω⊆Zn. Adiscrete operator is a map T:V→W, which assigns to each function f∈Va function T f ∈W. Definition 3.18 (Shift Operator).Let x:Z→C. We define the shift operator Ekacting on a discrete function x(t)by Ekx:= x(t+k) for any integer k∈Z. 25 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,+). 32 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). 33 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(Ω). 34 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. 35 (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(Ω). 36 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. 37 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. 38 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,21] 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∥. 39 Definition 4.3 (Dual Space).Let (X, ∥·∥)be a normed vector space over Ror C. The dual space of X, denoted by X∗, is the collection of all linear functionals F:X→R(or C) that are bounded, i.e., there exists a constant C≥0such that |F(x)| ≤ C∥x∥,∀x∈X. The operator norm of F∈X∗is defined by ∥F∥:= sup x∈X, x=0 |F(x)| ∥x∥. Equipped with this norm, the dual space (X∗,∥·∥)is itself a Banach space, even if Xis not complete. 4.3 Compactness of the Unit Ball in Discrete Function Spaces We illustrate the phenomenon using the discrete Hilbert space L2(Ω). Theorem 4.4 (Compactness of the Unit Ball).Let Ω⊆Zn. Consider the unit ball B:= {f∈ L2(Ω) : ∥f∥L2≤1}. Then: 1. If Ωis finite (say |Ω|=N < ∞), then L2(Ω) ∼ =RN(or CN). By the Heine–Borel theorem, every closed and bounded set in RNis compact. Hence Bis compact. 2. If Ωis infinite, then L2(Ω) is infinite-dimensional. In this case, Bis not compact in the norm topology. Counterexample in the infinite case. Let Ω = Zand consider the standard orthonormal basis {en}n∈Z⊂ L2(Z), where en(m) := (1, m =n, 0, m =n. Clearly, ∥en∥L2= 1 for all n, so each en∈B. However, for n=m, we have ∥en−em∥L2=√2. Thus, the sequence {en}has no Cauchy subsequence, and therefore no convergent subsequence in norm. Hence, Bis not compact. Theorem 4.5 (Banach–Alaoglu Theorem).Let Xbe a normed vector space over Ror C, and let X∗denote its dual space. Consider the closed unit ball in the dual space, B∗:= {f∈X∗:∥f∥ ≤ 1}. Then B∗is compact in the weak-* topology on X∗, i.e. the topology of pointwise convergence on X. 40 4.4 Hahn–Banach Theorem Definition 4.6 (Sublinear functional).Let Xbe a real vector space. A mapping p:X→Ris called a sublinear functional if it satisfies: 1. Positive homogeneity: p(λx) = λp(x)for all x∈Xand λ > 0. 2. Subadditivity: p(x+y)≤p(x) + p(y)for all x, y ∈X. Theorem 4.7 (Hahn–Banach Theorem).Let Xbe a real vector space, p: X→Ra sublinear functional, and let Y⊂Xbe a linear subspace. Suppose f:Y→Ris a linear functional such that f(y)≤p(y),∀y∈Y. Then there exists a linear functional F:X→Rextending f(i.e. F|Y=f) such that F(x)≤p(x),∀x∈X. Example 4.2 (Hahn–Banach Extension on Z).Consider the Banach space X=L1(Z) = nf:Z→RX x∈Z|f(x)|<∞o. Subspace. Let F0(Z) = {f:Z→R|supp(f)is finite} be the space of finitely supported functions. Functional on the subspace. Define g(f) = f(0), f ∈F0(Z), that is, evaluation at the origin. Boundedness. For any f∈F0(Z), |g(f)|=|f(0)| ≤ X x∈Z|f(x)|=∥f∥L1. Hence gis a bounded linear functional. Hahn–Banach extension. By the Hahn–Banach theorem: 1. The functional gdefined on F0(Z)can be extended to a bounded linear functional Gon all of L1(Z). 2. The extension satisfies |G(f)|≤∥f∥L1,∀f∈ L1(Z). Conclusion. In fact, every such extension has the form G(f) = X x∈Z w(x)f(x), w ∈ L∞(Z), 41 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. 48 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. [15,27,28]). 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. 49 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. 50 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ω. 51 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 Justification of the Fourier Ansatz for Linear P∆E The use of the Fourier ansatz in solving linear partial difference equations is not an ad hoc guess, but rather a direct consequence of the spectral properties of shift and difference operators. Consider first the spatial shift operator (Exu)(x) = u(x+ 1). Its eigenfunctions are exponential functions of the form eikx, since Exeikx=eik(x+1) =eik eikx. Thus eikx is an eigenfunction of Exwith eigenvalue eik. Similarly, in the temporal direction, the forward difference operator ∆tu(t) = u(t+ 1) −u(t) acts on exponential functions eλt by ∆teλt =eλ(t+1) −eλt = (eλ−1) eλt. Hence eλt is also an eigenfunction, with eigenvalue eλ−1. More generally, any linear difference operator constructed from shifts and finite linear combinations thereof has exponential functions as its eigenfunctions. Since linear P∆Es are built from such operators, solutions may be represented as superpositions of separable modes of the form u(t, x) = eλteikx. This explains the standard Fourier ansatz: each exponential mode diagonalizes the operators involved, reducing the partial difference equation to an algebraic relation between λand k. By the principle of superposition, the general solution is then obtained by combining these modes according to the Fourier expansion of the initial data. 6.7 Classification of Second Order Linear P∆E For simplicity, let us consider a linear evolution equation with two independent variables (t, x), of structural order 2, and without mixed shift operators: aEtu+bE−1 tu+cExu+dE−1 xu+e u = 0, u =u(t, x). 52 Fourier–Laplace Ansatz. We employ the Fourier–Laplace ansatz u(t, x) = eλteikx, where λ∈Cand k∈Rdenote the temporal growth rate and spatial frequency, respectively. Substituting into the equation yields aeλ+be−λ+ceik +de−ik +e= 0. Let z=eλ, w =eik, so that the relation becomes the Laurent polynomial Q(z, w) = az +bz−1+cw +dw−1+e= 0. Classification. Solving for zin terms of wgives the dispersion relation z=eλ=f(k). The qualitative behaviour of the solution is then determined by the modulus of z: •If |z|<1, the solution decays in time, corresponding to a diffusive/parabolic type. •If |z|= 1, the solution oscillates without growth or decay, corresponding to a wave/hyperbolic type. •If |z|>1, the solution exhibits exponential growth in time, leading to instability or blow-up. Steady-state problems. For purely steady-state problems of the form aExu+bE−1 xu+cEyu+dE−1 yu+e u = 0, u =u(x, y), we instead use the Fourier ansatz u(x, y) = eikxeimy,(k, m)∈R2, which leads to the symbol Q(w, v) = aw +bw−1+cv +dv−1+e= 0, with w=eik,v=eim. The structure of the zero set {(w, v)∈C2:Q(w, v)=0} characterizes the admissible steady-state solutions. Example 6.1 (Discrete Diffusion Equation).Consider the discrete diffusion equation Etu=1 2Exu+E−1 xu. 53 Fourier–Laplace Ansatz. Take u(t, x) = eλteikx, k ∈R. Substituting into the equation gives eλ=1 2eik +e−ik= cos(k). Dispersion relation. Thus z=eλ= cos(k). Since |cos(k)|<1for almost all k∈(0, π), the temporal growth rate satisfies ℜ(λ)<0in general. Classification. Therefore the solution decays in time and smooths out oscillations. This behaviour corresponds to the parabolic/diffusive type, where the solution tends to flatten as t→ ∞. Example 6.2 (Discrete Laplace Equation).Consider the two–dimensional discrete Laplace equation ∇2u= 0, which in shift–operator form reads Exu+E−1 xu+Eyu+E−1 yu−4u= 0. Fourier Ansatz. Take u(x, y) = eikxeimy,(k, m)∈R2. Substituting into the equation yields eik +e−ik +eim +e−im −4=0, or equivalently 2 cos(k) + 2 cos(m)−4=0. Frequency condition. This reduces to cos(k) + cos(m)=2. Since cos(θ)≤1, the only solution is k≡0 (mod 2π), m ≡0 (mod 2π). Conclusion. Hence the only admissible Fourier modes are constant frequencies. The discrete Laplace equation thus admits only constant solutions (up to boundary conditions), consistent with the continuous case where ∆u= 0 admits only harmonic functions with trivial oscillatory modes. 54 6.8 Discrete Fourier Analysis and Hadamard’s Well-posedness Consider the two-dimensional discrete Laplace equation Exu+E−1 xu+Eyu+E−1 yu−4u= 0. Suppose we prescribe the initial conditions u(0, y) = f(y), u(1, y) = g(y), with no boundary conditions in the y-direction. Fourier Ansatz. We take u(x, y) = ekxeimy, m ∈[−π, π], and substitute into the equation. The resulting dispersion relation is cosh(k)=2−cos(m). Hence ek=e±cosh−1(2−cos m). Observation. For almost all frequencies m, we have 2 −cos(m)>1, so that cosh−1(2 −cos m)>0. Therefore |ek|>1 for generic m, with the only neutral case being m= 0 (mod 2π), where k= 0. This means that most Fourier modes lead to exponential growth in the x-direction. General Solution. By Fourier inversion, the full solution can be written as u(x, y) = 1 2πZπ −πA(m)excosh−1(2−cos m)+B(m)e−xcosh−1(2−cos m)eimy dm, where the coefficients A(m), B(m) are determined by the initial data f(y), g(y). Conclusion. For generic initial conditions, unless all coefficients A(m) vanish, the solution contains growing modes of the form excosh−1(2−cos m)and therefore exhibits unbounded growth (“blow-up”). This demonstrates that interpreting the discrete Laplace equation as an evolution equation leads to an ill-posed problem in the sense of Hadamard: the solution is unstable with respect to the initial data. Remark 6.7. The Fourier analysis shows that the only stable mode is the trivial zero-frequency mode (m= 0), while almost all other frequencies produce exponential growth. Thus the discrete Laplace equation, if reinterpreted as an evolution law, is fundamentally unstable and ill-posed. 55 Hadamard’s Well-posedness in Partial Difference Equations In analogy with the classical theory of partial differential equations (PDEs), we extend the concept of Hadamard well-posedness to partial difference equations (P∆Es). A problem is said to be well-posed if it satisfies the following three criteria: 1. Existence: For every admissible set of initial or boundary data, a solution exists. 2. Uniqueness: The solution is unique within the prescribed function space. 3. Stability (Continuous Dependence): The solution depends continuously on the initial and boundary data; in particular, small perturbations in the data lead to only small changes in the solution. If any of these three conditions fails, the problem is referred to as an ill-posed problem. This framework allows us to analyze the stability properties of P∆Es in direct analogy with PDEs. 56 7 First Order Equations in Time 7.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 7.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. 7.2 1D Pascal Evolution Equation We now consider the linear partial difference equation inspired by the recurrence relation in combinatorics [23] 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. 57 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. 64 8 Second Order Equations in Time 8.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. 8.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, 65 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. 66 Figure 1: 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. 67 8.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 λ. 68 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, 69 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 8.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. 70 9 Steady State Problems 9.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. 71 9.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, (1) EyY+E−1 yY−2Y=λY. (2) Solution for X(x).Expanding Eq. (1) 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. 72 Solution for Y(y).Substituting λninto Eq. (2) 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. 73 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. 80 9.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). 81 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. 82 10 Discrete Evolution Equations 10.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. 10.2 Definitions Definition 10.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 10.2 (Discrete Spatiotemporal Dynamical System).Let u:N0×Zn→C denote the state of the system, where 83 •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 10.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 10.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 10.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. 84 [26] Example 10.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 10.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. 85 [33] Example 10.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] Observations and Phenomena The following figure illustrates the evolution of Langton’s Ant over discrete time steps. 86 Figure 2: The evolution of Langton’s Ant over time. The black squares represent flipped cells (state 1), the white background denotes unvisited cells (state 0), and the red dot indicates the current position of the ant at the final time step. The system exhibits transient chaos, followed by the emergence of a repeating structure called the highway, which acts as a periodic attractor in the phase space. 87 10.3 Semigroup Theory Definition 10.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. [30] Definition 10.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 10.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 10.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. 88 10.4 Initial Value Problems Definition 10.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 10.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. [33] 89 Expanding the trinomial yields (1 + e−ik +e−im)t=X a+b+c=t a,b,c≥0 t! a!b!c!e−ikae−imb, so that G(t, x, y) = t t−x−y, x, y. Convolution Form of the Solution. Therefore, the general solution can be written in convolution form: u(t, x, y) = X (s,p)∈Z2 f(s, p)C(t, x −s, y −p), where C(t, x, y) =      t! (t−x−y)! x!y!,if 0 ≤x, y ≤tand x+y≤t, 0,otherwise. Fourier–Combinatorial Identity. This establishes the identity 1 (2π)2Zπ −πZπ −π1 + e−ik +e−imtei(kx+my)dk dm =t t−x−y, x, y, which shows that the trinomial coefficient arises naturally as the Green’s function of the 2D Pascal Evolution Equation. 11.3 3D Pascal Evolution Equation We now consider a simple linear partial difference equation in three spatial dimensions: Etu=u+E−1 xu+E−1 yu+E−1 zu, u :Z4→R. The initial condition is prescribed as u(0, x, y, z) = f(x, y, z)∈ L2(Z3), with no boundary conditions. Fourier Ansatz. We apply the Fourier ansatz u(t, x, y, z) = eλtei(kx+my+nz). Substitution into the equation gives the dispersion relation eλ= 1 + e−ik +e−im +e−in. 96 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. 97 11.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. 98 Theorem 11.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. 99 12 Nonlinear Equations with mod n Nonlinearity 12.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. 12.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. 100 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. 101 Figure 3: Spatiotemporal plot for δinitial condition. 102 Figure 4: Random initial condition with f(x)∈ {0,1}. 103 Figure 5: Random initial condition with f(x)∈R. 12.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. 104 Figure 6: 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. 12.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), 105 Figure 11: Spatiotemporal evolution of the Mod 5 Sierpi´nski Equation 12.7 Mod 4 Sierpinski Fractal Equation We consider the following nonlinear partial difference equation: Etu= mod4E−1 xu+u+Exu+δ(x) with initial condition: u(0, x) = δ(x) and no boundary condition. u=u(t, x)∈ {0,1,2,3} Here, δ(x) is the discrete delta function defined by: δ(x) = (1,if x= 0 0,otherwise The values of u(t, x) are restricted to {0,1,2,3}, and we visualize them using the following colour scheme: Value Colour 0 White 1 Green 2 Purple 3 Yellow The spatiotemporal plot of this equation is shown below: 112 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. 113 12.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) [13, 14] 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. 114 13 System of Partial Difference Equations 13.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. 13.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). 115 Definition 13.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 13.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 13.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 13.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. 116 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 13.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 13.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. 117 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 13.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. 118 13.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. 119 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 13.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. 120 Hence the matrix Ais diagonal: A= diag I+α∇2, I +α∇2, I +α∇2. It is now clear that the system decouples into three scalar equations, each of which may be solved independently. Thus the vector heat equation is simply three copies of the scalar discrete heat equation. 121 14 Atlas of Nonlinear Partial Difference Equations 14.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. 128 Coupled Map Lattice 1D Coupled Map Lattice The one-dimensional Coupled Map Lattice (CML), originally proposed by Kaneko [34], 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. 129 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. 130 14.2 Elementary Cellular Automata The study of one-dimensional cellular automata (1D CA) was significantly advanced by Stephen Wolfram [33], 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. 131 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. 132 Example 1 This is the spatiotemporal plot of Rule 90 with the initial condition u(0, x) = δ(x) and no boundary conditions. Figure 12: 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. 133 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 13: 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. 134 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 14: Spatiotemporal plot of Rule 90 with random initial condition and boundary conditions. As shown, the system exhibits spatiotemporal chaos. 135 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. 136 The Rule 153 We define the evolution equation as follows: Etu= mod21−E−1 xu 1 + Exu+Exu·E−1 xu+3E−1 xu−4Exu+u Example •Initial condition: u(0, x) is assigned a random binary value (0 or 1) at each spatial point. •Boundary condition: u(t, −200) = u(t, 200) = 0 for all t. •Colour scheme: 0 is shown as white, and 1 is shown as black. The following image illustrates the space-time evolution from t= 0 to t= 400, and x∈[−200,200], with a random initial condition: This system exhibits spatiotemporal chaos. 137 This model shows interesting connections to: •Phase synchronization phenomena •Topological vortex dynamics •Pattern formation and self-organization Figure 19: Spatiotemporal evolution of the Kuramoto model. Topological defects (vortices) emerge and move, and their annihilation leads to large-scale synchronization. 144 14.5 Ising Model Governing Equation We propose the following discrete evolution equation for the well-known Ising Model [29]: Etu= (1 −δ(sign(JSvu)−sign(u))) θ(F−ε) sign(JSvu) +δ(sign(JSvu)−sign(u)) θ(F−F0) (−sign(JSvu)) (3) where Svu=X (i,j)∈V Ei xEj yu, •Vis the von Neumann neighbourhood. •u=u(t, x, y)∈ {−1,+1}is the binary state at discrete time tand spatial coordinate (x, y). •Etu:= u(t+ 1, x, y) is the next time step value. •F=F(t, x, y)∼ N(µ, c2) is a random driving field at site (x, y), drawn from a normal distribution centered at µ=kT, where Tis the temperature and kis a constant. •εis a small positive number (baseline activation threshold). •F0is a larger threshold required for flipping a stable site. •Jis the coupling constant determining interaction strength between neighbouring sites. Physical Interpretation The function u(t, x, y) represents the spin (or state) of a particle or site at position (x, y) and time t. The variable Etudenotes the updated state at time t+1. The update is governed by local interactions and stochastic environmental drive F. •The first term activates when the site’s current state is different from the local neighbourhood majority (sign(JSvu)= sign(u)), allowing it to flip easily when F > ε. •The second term activates when the site’s current state matches the neighbourhood majority, but may still flip if the driving force exceeds a higher threshold F0, modeling thermal noise or instability. This model expresses the Ising-like behaviour in the language of partial difference equations (P∆E), linking local deterministic update rules with stochastic driving force under thermal control. 145 14.6 Discrete Logistic Diffusion Equation In this section, we propose a discrete-time, discrete-space population model called the Discrete Logistic Diffusion Equation, designed to capture the expansion and growth of a population over a spatial domain. This equation is a discrete analogue of the Fisher-KPP Equation. ∆tu=jα∇2u+ru 1−u Kk (4) where the discrete Laplacian is defined as ∇2u:=  X (i,j)∈V Ei xEj yu −4u with the von Neumann neighbourhood V={(1,0),(−1,0),(0,1),(0,−1)} Here, the state variable u=u(t, x, y)∈Z≥0represents the number of individuals at site (x, y) at time t, and ∆tu=u(t+ 1, x, y)−u(t, x, y). This is a Nonlinear Autonomous Partial Difference Equation Biological Interpretation We consider a two-dimensional discrete lattice where each site (x, y) corresponds to a unit habitat area. Time tis discrete and may represent one minute, one hour, or one day depending on the modeling scale. The variable u(t, x, y)∈Z≥0 denotes the population size at location (x, y) and time t. The evolution equation consists of the following components: •u: current population size at each site. •α∇2u: discrete diffusion term. Individuals move to neighboring sites in the von Neumann neighborhood. The coefficient αcontrols how strongly the population spreads in space. •ru 1−u K: logistic growth term. Population increases locally at growth rate r, but the growth slows as uapproaches the carrying capacity K.K is the maximum number of individuals allowed in each unit area.K∈Z •⌊·⌋: the floor function ensures that the updated population remains an integer, reflecting the discreteness of individuals in biological populations. •Etu=u(t+ 1, x, y): forward time evolution operator, representing the state of the system at the next time step. This model captures the competition between local growth and spatial dispersal in a biologically plausible manner. It also respects the fact that population counts are discrete and limited by spatial constraints. 146 Simulation Results Delta Initial Condition We consider the following parameter setup for the simulation: •Diffusion coefficient: α= 0.05 •Carrying capacity: K= 1000 •Initial condition: u(0, x, y) = δ(x)δ(y), meaning a single individual is placed at the center of the grid. •Boundary condition: (u(t, −100, y) = u(t, 100, y) = 0, u(t, x, −100) = u(t, x, 100) = 0, representing a bounded domain with zero population on the edges. We simulate the evolution of the system under the Discrete Logistic Diffusion Equation for various values of the growth rate r. Below are the population distributions u(t, x, y) at time t= 250, corresponding to different values of r: 147 Random Initial Condition In this simulation, we keep all parameters and boundary conditions unchanged: •α= 0.05 •K= 1000 148 •Boundary conditions: u(t, −100, y) = u(t, 100, y)=0, u(t, x, −100) = u(t, x, 100) = 0 However, we change the initial condition to a random configuration: u(0, x, y) = random integers (e.g., uniformly sampled from {0,1,2, ...}) Below are the evolution plots for different values of r: 149 From the simulation results under the delta initial condition, where the system is initialized with u(0, x, y) = δ(x)δ(y), we observe a nearly circular outward expansion due to the localized starting point. As the growth rate rincreases, the resulting spatial patterns exhibit a rich sequence of transitions: •For small values of r, the solution reaches a spatially homogeneous steady state. •As rincreases, the system begins to form checkerboard-like structures. •With further increase in r, the patterns evolve into concentric wave-like ripples. •At higher values of r, the symmetry breaks down and complex, irregular structures emerge. Under the random initial condition, where the system is initialized with spatially heterogeneous random values, a similar progression is observed. The system transitions from a steady state to increasingly disordered and chaotic spatial patterns as rincreases. These observations suggest that the model exhibits a type of infinitedimensional bifurcation behavior, where complex pattern dynamics emerge through successive instabilities driven by the control parameter r. Simple Model with Infinite Complexity If we remove the floor function, set the carrying capacity K= 1, and restrict the model to one spatial dimension, the discrete equation simplifies to: ∆tu=α∇2u+ru(1 −u),(5) 150 where u=u(t, x)∈R≥0represents the population density at site xand time t, and the discrete Laplacian is given by: ∇2u:= Exu+E−1 xu−2u. We call this Logistic Diffusion Equation. We fix the diffusion coefficient at: α= 0.05, and use the following settings: •Initial condition: u(0, x) is a randomly generated real-valued function on the spatial domain. •Boundary condition: u(t, −100) = u(t, 100) = 0 for all t. Below we present the spatiotemporal evolution of u(t, x) for different values of the growth rate r: 151 152 These simulations illustrate a rich spectrum of behaviours, including steady states, perfect periodicity, quasiperiodicity, spatiotemporal intermittency, and full spatiotemporal chaos, depending on the value of r. As the parameter rincreases, the system exhibits a wide variety of dynamical behaviours: •r∈[0,1.8]: the system converges to a steady state, resembling a fixed point. •r≈1.9: the solution begins to bifurcate. •r≈2.0: the system forms perfect checkerboard-like patterns. •r∈[2.2,2.3]: the system exhibits spatially periodic structures. •r≈2.4: the solution appears almost perfectly periodic, but upon closer inspection reveals quasiperiodic behaviour. •r≈2.5: the system begins to show localized irregularities — part quasiperiodic, part disordered — suggesting the onset of chaotic features. 153