EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 3, 2024, 1565-1584 ISSN 1307-5543 – ejpam.com Published by New York Business Global An efficient convergent approach for difference delayed reaction-diffusion equations M. Heidari1, M. Ghovatmand1,∗, M. H. Noori Skandari1, D. Baleanu2 1 Faculty of Mathematical Sciences, Shahrood University of Technology, Shahrood, Iran 2 Department of Computer Science and Mathematics, Lebanese American University, Beirut, Lebanon Abstract. It is usually not possible to solve partial differential equations, especially the delay type, with analytical methods. Therefore, in this article, we present an efficient method for solving differential equations of the difference delayed reaction-diffusion type, which can be generalized to other delayed partial differential equations. In the proposed approach, we first convert the delayed equation into an equivalent non-delayed equation by inserting the corresponding delay function with an effective technique. Then, using a pseudo-spectral method, we discretize the obtained equation in the Legendre-Gauss-Lobatto collocation points and present an algebraic system with an equal number of equations and unknowns which can be solved by quasi-Newton methods such as Levenderg-Marquardt algorithm. The approximate solutions can be obtained with exponential accuracy. The convergence analysis of the method is fully discussed and four examples are presented to evaluate the results and compare with one of the conventional methods used to solve partial differential equations, that is, the compact finite difference method. 2020 Mathematics Subject Classifications: 35K57, 65M70, 65N35 Key Words and Phrases: Difference delayed reaction-diffusion equations, Lagrange interpolat- ing polynomials, Legendre-Gauss-Lobatto points, Convergence Analysis, Modulus of continuity 1. Introduction Reaction-diffusion (RD) equation is one of the partial differential equations (PDEs) which appears in physical and chemical models. Many chemical systems are described by this equation, where the diffusion of matter competes with production of some kind of chemical reaction [3, 11, 13, 26, 32, 33]. Recently, some researchers have attempted to numerically solve this equation. Sharifi and Rashidian [28] have applied the collocation method that is a special case of the spectral methods for this equation. Also, in [2], the solvability of optimal control problems on both weak and strong solutions of a boundary ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v17i3.5197 Email addresses: Marziehheidari22@yahoo.com (M. Heidari), Mghovat@gmail.com (M. Ghovatmand), Math.Noori@yahoo.com (M. H. Noori Skandari ), Dumitru.Baleanu@lau.edu.lb (D. Baleanu) https://www.ejpam.com 1565 © 2024 EJPAM All rights reserved. M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1566 value problem for the nonlinear reaction–diffusion–convection (RDC) equation with vari- able coefficients is investigated and the requirements for smoothness of the multiplicative control are reduced. Some other methods for RD equation exist in Liu et al. [19], Koto [12], Priyanka et al. [25], Visuvasam et al. [31] and their references. The (difference) delayed RD (DRD) equation is a type of extended RD equations that are widely used in various sciences. One of the interesting applications is in biological sciences and for modeling the spread of bacteriophage infection [7]. The other applica- tion is the modeling of range of prey-predator systems [1]. Moreover, DRD equations are utilized for the modeling of virus infection of Hepatit C [34]. Further, DRD equations are widely proposed as models for population ecology and cell biology. For example, the diffusive Mackey–Glass equation [18] is used to describe a single-species population with age-structure and diffusion. Also, the diffusive Hematopoiesis model [35] is presented to investigate the dynamics of blood cell production. Most of the results in literature indi- cate that a delay term could change the dynamic properties of a system such as stability, oscillation, bifurcations, chaos etc. That is the reason that such equations became the focus of researchers in numerical analysis and simulation. So far many researchers have paid attention to the DRD equation and suggested many methods for solving this equation. Zhao and Ge [38] checked out the existing traveling wavefront solutions for DRD equations and qualitative behavior of these solutions. Wu and Lan [36] applied a method for RD with several delays by considering compact functions in a continuous space and applying the convergence theorem Lebseque’s dominated (Some other related works exist in [6, 14, 17, 21]). Moreover, Meleshko and Moyo [22] suggested the Lie group for solving this equation. Li et al. [16] proposed a Galerkin method with local discontinuity for solving this equation. They used a space of discontinuous piecewise polynomials to approximate the solution and estimated the delay term by an interpolation polynomial. This method is stable, conservative and has high accuracy but for delayed systems is weak. Polyanin and Zhurove [23] applied another method for DRD equation where the solution is considered as product of two functions that defined in the form additional functional limitations and delayed PDE. In this method, all of the solutions have some free parameters that are suitable for the main equations. Chen [4] proposed the linear compact Θ method for DRD equation. [30] utilized a method of reduction lines and finite differences to convert the DRD equation to a system of ordinary differential equations. By this method, error is reduced when all exact and numerical solutions tend to zero. L. Liu and J. Nieto [20] investigated the asymptotic behavior of fractional DRD equation. They proved the existence and uniqueness of the solution in the sense of weak solution with using the classic Galerkin approximation method and comparison principal. According to the knowledge that we have the existing methods for delayed partial differ- ential equations, especially the DRD equations, it is still felt to solve these equations with an efficient method with a suitable convergence rate and high accuracy. On the other hand, spectral and pseudo-spectral methods [8–10, 24, 37] are one of the most efficient and useful methods for solving continuous-time problems such as delayed/non-delayed or- dinary and partial differential equations. These methods have a good accuracy and high speed convergence compared with those of other methods. Hence, in this paper, we solve M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1567 DRD equation by using a new convergent pseudo-spectral collocation method. The DRD equation is first converted into two parts: after delay time and before delay time. An equivalent system of equations is then gained by inserting the delay time-function in DRD equation. The achieved system is discretized and converted into an algebraic system of equations using two-dimensional Lagrange polynomial and based on the Legendre-Gauss- Lobatto (LGL) points. Solving the last system, tends to the approximate solutions for the DRD equation. We analyze the convergence of our proposed method based on the modulus of continuity functions and its corresponding space and norm. Also, we show its efficiency with several comparative numerical examples. This article contains the following sections. In Section 2 and 3, we introduce a general form of DRD equation and implement our method for solving it, respectively. In Section 4, we discuss the convergence of method and in Section 5 we illustrate its validity by solving some numerical comparative examples. 2. Difference delayed reaction-diffusion equation Consider the following form of DRD equations ηζ(ζ, ι) = Z(ζ, ι, η(ζ, ι), η(ζ − µ, ι), ηιι), (ζ, ι) ∈ [0,M ]× [0, P ], (1) with following initial and boundary conditions η(ζ, ι) = φ(ζ, ι), (ζ, ι) ∈ [−µ, 0]× [0, P ], (2) η(ζ, 0) = g1(ζ), η(ζ, P ) = g2(ζ), ζ ∈ [0,M ] (3) where µ is the delay constant and Z , φ, g1, g2 are given smooth functions. The specific form of equation (1) by conditions (2) and (3) can be stated as the following system which has been studied by many researchers in recent years [15, 16] ηζ(ζ, ι) = ηιι(ζ, ι) + η(ζ − µ, ι), (ζ, ι) ∈ [0,M ]× [0, P ] η(ζ, ι) = φ(ζ, ι), (ζ, ι) ∈ [−µ, 0]× [0, P ], η(ζ, 0) = g1(ζ), η(ζ, P ) = g2(ζ), ζ ∈ [0,M ]. In this work, we implement our method for general form of DRD equation. Since the manner of system is different before and after delay time ζ = µ, we rewrite DRD equation (1) with conditions (2) and (3) as follows ηζ(ζ, ι) = { Z(ζ, ι, η(ζ, ι), φ(ζ − µ, ι), ηιι(ζ, ι)), (ζ, ι) ∈ [0, µ]× [0, P ], Z(ζ, ι, η(ζ, ι), η(ζ − µ, ι), ηιι(ζ, ι)), (ζ, ι) ∈ [µ,M ]× [0, P ], η(0, ι) = φ(0, ι), ι ∈ [0, P ], η(ζ, 0) = g1(ζ), η(ζ, P ) = g2(ζ), ζ ∈ [0,M ]. (4) M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1568 3. Implementing the method We explain the approximation solution of set (4) in form η(ζ, ι) ≃ ηQ(ζ, ι) = Q∑ i=0 Q∑ j=0 η̄Qijhi(ζ)hj(ι), (5) where η̄Qij are unknown coefficients and hν(.) is the Lagrange basis polynomial and defined in form hi(ζ) = Q∏ k=0,i ̸=k ζ − ζk ζi − ζk , hj(ι) = Q∏ l=0,j ̸=l ι− ιl ιj − ιl , (6) where {ιl}Ql=0 and {ζk}Qk=0 are the shifted LGL points in [0, P ] and [0,M ], respectively. These points are described as follows{ ιl = P 2 (vl + 1), l = 0, 1, · · · , Q, ζk = M 2 (vk + 1), k = 0, 1, · · · , Q, (7) where {vl}Ql=0 are the roots of following polynomial R(z) = (1− z2)ṘQ(z); z ∈ [−1, 1], where Legendre polynomial RQ(z) is described by recurrent formula{ RQ+1(z) = 2Q+1 Q+1 zRQ(z)− Q Q+1RQ−1(z), R0(z) = 1, R1(z) = z. (8) By using the approximation (5) we achieve{ ηζ(ζ, ι) ≃ ∑Q i=0 ∑Q j=0 η̄ Q ijh ′ i(ζ)hj(ι), ηιι(ζ, ι) ≃ ∑Q i=0 ∑Q j=0 η̄ Q ijhi(ζ)h ′′ j (ι). (9) according to following Lagrange polynomials property hi(zj) = { 1, i = j, 0, i ̸= j, (10) where {zj}Qj=0 are the LGL points. We can approximate the functions η(ζk, ιl), ηζ(ζk, ιl) and ηιι(ζk, ιl) in the interpolating points {ιl}Ql=0 and {ζk}Qk=0 as η(ζk, ιl) ≃ η̄Qkl, k, l = 0, 1, · · · , Q, ηζ(ζk, ιl) ≃ Q∑ i=0 Q∑ j=0 η̄Qijh ′ i(ζk)hj(ιl) = Q∑ i=0 η̄ilDki, l = 0, 1, · · · , Q, k = 0, 1, · · · , Q, (11) M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1569 ηιι(ζk, ιl) ≃ Q∑ i=0 Q∑ j=0 η̄Qijhi(ζk)h ′′ j (ιl) = Q∑ j=0 η̄Qkjh ′′ j (ιl) = Q∑ j=0 η̄QkjD (2) lj ; l = 0, 1, · · · , Q, k = 0, 1, · · · , Q, (12) where Dki and D2 lj are defined as follows Dki = h′i(ζk) =  − 2 M Q(Q+1) 4 , k = i = 0, RQ(ζk) RQ(ζi) 1 ζk−ζi , k ̸= i, 0, 1 ≤ k = i ≤ Q− 1, 2 M Q(Q+1) 4 , k = i = Q, (13) and D (2) lj = Q∑ p=0 DlpDpj . (14) Suppose index sµ is such that 0 = ζ0 < ζ1 ≤ ζ2 < · · · < ζsµ ≤ µ < ζsµ+1 < · · · < ζQ = M. Then system (1)-(3) can be approximated as the following discrete algebraic system ∑Q i=0 η̄ Q ilDki = { Z(ζk, ιl, η̄ Q kl, φ(ζk − µ, ιl), ∑Q j=0 η̄ Q kjD 2 lj), k = 1, · · · , sµ, l = 1, · · · , Q− 1, Z(ζk, ιl, η̄ Q kl, ∑Q i=0 η̄ Q il hi(ζk − µ), ∑Q j=0 η̄ Q kjD 2 lj), k = sµ+1, · · · , Q, l = 1, · · · , Q− 1, η̄Q0l = φ(0, ιl), l = 1, · · · , Q− 1, η̄Qk0 = g1(ζk), η̄ Q kQ = g2(ζk), k = 0, · · ·Q, (15) with unknowns η̄Qij for i, j = 0, 1, · · · , Q. Notice that the number of this equations is (Q+1)2 which is equal to the number of variables. By solving the system (15), we get the approximate solution η̄Qij (i, j = 0 · · · , Q) and continuous approximate solution ηQ(ζ, ι) presented with (5). 4. Convergence study Now we examine the convergence of the method. We define Γ̃ = [0,M ] × [0, P ] and show the space of all continuously differentiable functions from order k on Γ̃ by Ck(Γ̃). Definition 1. Function Ψ : R+ → R+ with the following characters [27] is called a modulus of continuity • Ψ is continuous and increasing, • Ψ(a) → 0 when a → 0, • Ψ(a1 + a2) ≤ Ψ(a1) + Ψ(a2) for any a1, a2 ∈ R, M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1570 • a ≤ bΨ(a) for some b > 0 and 0 < a ≤ 2 . One of the cases for modulus of continuity is Ψ(a) = aα, 0 < α ≤ 1. Show the unit circle in R2 by B2. The continuous function e on Γ̃, admits Ψ(.) as modulus of continuity if |e(., .)|Ψ = sup{ |e(ζ ′, ι′)− e(ζ ′′, ι′′)| Ψ(∥(ζ ′, ι′)− (ζ ′′, ι′′)∥∞) : (ζ ′, ι′), (ζ ′′, ι′′) ∈ Γ̃, (ζ ′, ι′) ̸= (ζ ′′, ι′′)}, (16) is finite and ∥(ζ ′, ι′)− (ζ ′′, ι′′)∥∞ = max{|ζ ′ − ζ ′′|, |ι′ − ι′′|}. Let C1 Ψ(B 2) includes functions with first continuous differentiation on the B2 that is equipped with the following norm ∥e(., .)∥1,Ψ = ∥e(., .)∥∞ + ∥eζ′(., .)∥∞ + ∥eι′(., .)∥∞ + |eζ′(., .)|Ψ + |eι′(., .)|Ψ. (17) Now we define C1 Ψ(Γ̃) = {e(., .) ∈ C1(Γ̃) : ∀(ζ ′′, ι′′) ∈ Γ̃,∃ map Θ : B2 → Γ̃ s.t (ζ ′′, ι′′) ∈ int(Θ(B2)), e ◦Θ(., .) ∈ C1 Ψ(B 2)}. (18) So if for some Θ1, · · · ,Θn Γ̃ = n⋃ i=1 int(Θi(B 2)), then e(., .) ∈ C1 Ψ(Γ̃) if and only if e ◦Θi(., .) ∈ C1 Ψ(B 2) for each i = 1, · · · , Q. Furthermore C1 Ψ(Γ̃) with the following norm is a Banach space ∥e(., .)∥1,Ψ = Q∑ i=1 ∥e ◦Θi∥1,Ψ. (19) We show pol(Q,Q, Γ̃) as the space of all polynomials, pol(Q,Q, Γ̃) = {Θ(ζ, ι) = Q∑ i=0 Q∑ j=0 c̄ijζ iιj ; (ζ, ι) ∈ Γ̃, c̄ij ∈ R}. (20) Lemma 1. For each e(., .) ∈ C1 Ψ(Γ̃), there exists a Θ(., .) ∈ pol(Q,Q, Γ̃) such that ∥e(., .)−Θ(., .)∥∞ ≤ c0c1 2Q Ψ( 1 2Q ), (21) where c1 = ∥e(., .)∥1,Ψ and constant c0 is independent of Q. M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1571 Proof. For proof see [27]. To guarantee the existence of solution for system (15), we change it in form | ∑Q i=0 η̄ Q ilDki − Z ( ζk, ιl, η̄ Q kl, φ(ζk − µ, ιl), ∑Q j=0 η̄ Q kjD (2) lj ) | ≤ √ Q 2Q−1Ψ( 1 2Q−1), k = 1, · · · sµ, l = 1, 2, · · · , Q− 1, | ∑Q i=0 η̄ Q ilDki − Z ( ζk, ιl, η̄ Q kl, ∑Q i=0 η̄ Q il hi(ζk − µ), ∑Q j=0 η̄ Q kjD (2) lj ) | ≤ √ Q 2Q−1Ψ( 1 2Q−1), k = sµ + 1, · · · , Q, l = 1, 2, · · · , Q− 1, |η̄Q0l − φ(0, ιl)| ≤ √ Q 2Q−1Ψ( 1 2Q−1), l = 1, 2, · · · , Q− 1, |η̄Qk0 − g1(ζk)| ≤ √ Q 2Q−1Ψ( 1 2Q−1), k = 0, 1, · · · , Q, |η̄QkQ − g2(ζk)| ≤ √ Q 2Q−1Ψ( 1 2Q−1), k = 0, 1, · · · , Q, (22) where Ψ(.) is a given modulus of continuity and Q is sufficiently big. Whereas limQ→∞ √ Q 2Q−1Ψ( 1 2Q−1) = 0, any solution η̄Qkl (k = 0, 1, · · · , Q, l = 0, 1, · · · , Q) of (22) is a solution for (15) when Q → ∞. In the next theorem, we prove that the system (22) has a solution. Theorem 1. Let η(., .) is a solution for set (1)-(3) where η(., .) ∈ C1 Ψ(Γ̃). Then there is a positive integer Q1 such that for Q ≥ Q1 the set (22) has a solution as η̄Q = (η̄Qkl; k = 0, 1, · · · , Q, l = 0, 1, · · · , Q), that satisfies relationship |η(ζk, ιl)− η̄Qkl| ≤ c 2Q− 1 Ψ( 1 2Q− 1 ), k = 0, 1, · · · , Q, l = 0, 1, · · · , Q, (23) where c > 0 is a constant, independent of Q. Proof. Let Θ(., .) ∈ pol(Q − 1, Q, Γ̃) be the best approximation for ηζ(., .). By using the Lemma 4.1 we have ∥ηζ(ζ, ι)−Θ(ζ, ι)∥∞ ≤ ξ 2Q− 1 Ψ( 1 2Q− 1 ), (ζ, ι) ∈ Γ̃, (24) where ξ > 0 is not dependent on Q. Now, we describe the function η̄(., .) in form η̄(ζ, ι) = η(0, ι) + ∫ ζ 0 Θ(τ, ι)dτ, (ζ, ι) ∈ Γ̃, (25) and η̄Qkl = η̄(ζk, ιl); k = 0, 1, · · · , Q, l = 0, 1, · · · , Q. (26) M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1572 We will confirm that η̄Q = (η̄Qkl; k = 0, 1, · · · , Q, l = 0, 1, · · · , Q) applies to system (22). With (24), (25) and (26) we have |η(ζ, ι)− η̄(ζ, ι)| = | ∫ ζ 0 (ηζ(τ, ι)−Θ(τ, ι))dτ | ≤ ∫ ζ 0 |ηζ(τ, ι)−Θ(τ, ι)|dτ ≤ ξ 2Q− 1 Ψ( 1 2Q− 1 ) ∫ ζ 0 dτ ≤ ξM 2Q− 1 Ψ( 1 2Q− 1 ). (27) Now, according to (25), η̄(., ι), ι ∈ [0, P ] is a polynomial of degree at most Q on [0,M ]. Thus, Q∑ i=0 η̄QilDki = η̄ζ(ζk, ιl); k = 1, · · · , Q, l = 1, · · · , Q. (28) Therefor, with (26), (27) and (28) we have for k = 1, · · · , sµ and l = 1, · · · , Q, | Q∑ i=0 η̄QilDki − Z ( ζk, ιl, η̄ Q kl, φ(ζk − µ, ιl), η̄ιι(ζk, ιl) ) | ≤ |η̄ζ(ζk, ιl)− ηζ(ζk, ιl)|+ |ηζ(ζk, ιl)− Z ( ζk, ιl, η̄ Q kl, φ(ζk − µ, ιl), η̄ιι(ζk, ιl) ) | ≤ |Θ(ζk, ιl)−ηζ(ζk, ιl)|+|Z ( ζk, ιl, η(ζk, ιl), φ(ζk−µ, ιl), ηιι(ζk, ιl) ) −Z ( ζk, ιl, η̄ Q kl, φ(ζk−µ, ιl), η̄ιι(ζk, ιl) ) | ≤ ξ 2Q− 1 Ψ( 1 2Q− 1 ) + L ( |η(ζk, ιl)− η̄Qkl|+ |φ(ζk − µ, ιl)− φ(ζk − µ, ιl)| ) ≤ ξ 2Q− 1 Ψ( 1 2Q− 1 ) + L( ξM 2Q− 1 )Ψ( 1 2Q− 1 ) = ξ(1 +MP ) 1 2Q− 1 Ψ( 1 2Q− 1 ). (29) Here, the function Z has L > 0 as the Lipschitz constant with respect to its third and fourth components. And, for k = sµ + 1, · · · , Q− 1 and l = 1, · · · , Q− 1, | Q∑ i=0 η̄QilDki − Z ( ζk, ιl, η̄ Q kl, η̄(ζk − µ, ιl), η̄ιι(ζk, ιl) ) | ≤ |η̄ζ(ζk, ιl)− ηζ(ζk, ιl)|+ |ηζ(ζk, ιl)− Z ( ζk, ιl, η̄ Q kl, η̄(ζk − µ, ιl), η̄ιι(ζk, ιl) ) | ≤ |Θ(ζk, ιl)−ηζ(ζk, ιl)|+|Z ( ζk, ιl, η(ζk, ιl), η(ζk−µ, ιl), ηιι(ζk, ιl) ) −Z ( ζk, ιl, η̄ Q kl, η̄(ζk−µ, ιl), η̄ιι(ζk, ιl) ) | ≤ ξ 2Q− 1 Ψ( 1 2Q− 1 ) + L ( |η(ζk, ιl)− η̄Qkl|+ |η(ζk − µ, ιl)− η̄(ζk − µ, ιl)| ) M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1573 ≤ ξ 2Q− 1 Ψ( 1 2Q− 1 ) + L( 2ξM 2Q− 1 )Ψ( 1 2Q− 1 ) = ξ(1 + 2MP ) 1 2Q− 1 Ψ( 1 2Q− 1 ). (30) Moreover, for k = 0, 1, · · · , Q we have |η̄Qk0 − g1(ζk)| ≤ |η̄Qk0 − u(ζk, 0)|+ |η(ζk, 0)− g1(ζk)| ≤ Mξ 2Q− 1 Ψ( 1 2Q− 1 ), (31) and |η̄kP − g2(ζk)| ≤ |η̄QkP − η(ζk, P )|+ |η(ζk, P )− g2(ζk)| ≤ Mξ 2Q− 1 Ψ( 1 2Q− 1 ). (32) Further, for l = 1, · · · , Q− 1 |η̄0l − φ(0, ιl)| ≤ |η̄Q0l − η(0, ιl)|+ |η(0, ιl)− φ(0, ιl)| ≤ Mξ 2Q− 1 Ψ( 1 2Q− 1 ). (33) If we select Q1 ∈ N as for all Q ≥ Q1, max{Mξ, ξ(1 + 2MP )} ≤ √ Q, then η̄Q = (η̄Qkl; k, l = 0, 1, · · · , Q, ) for Q ≥ Q1 satisfies system (22) by using (29)-(33). In next theorem, we show convergence theorem of solutions. Theorem 2. Let { η̄Qkl; k, l = 0, 1, · · · , Q }∞ Q=Q1 be the sequence of solutions of system (22) and { ηQ(., .) }∞ Q=Q1 be the sequence of polynomials presented in (5). Consider for any ι ∈ [0, P ] the sequence { (ηQ(0, ι), ηQζ (., .)) }∞ Q=Q1 has a subsequence { (ηQi(0, ι), ηQi ζ (., .)) }∞ i=1 that converges to { χ∞(ι), q(., .) } uniformly, where q(., .) ∈ C2(Γ̃), χ∞(.) ∈ C2([0, P ]) and limi→∞Qi = ∞. Then η̄(ζ, ι) = lim i→∞ ηQi(ζ, ι), (34) is a solution of the system (1). Proof. Define η̄(ζ, ι) = χ∞(ι) + ∫ ζ 0 q(τ, ι)dτ. (35) At first, we prove that η̄(., .) satisfies (1)-(3). By contradiction, we assume that there is index 1 ≤ l ≤ Q− 1 and τ ∈ [0,M ] such that η̄ζ(τ, ιl)− Z ( ζk, ιl, η̄(τ, ιl), φ(τ − µ, ιl), η̄(τ, ιl) ) ̸= 0, k = 1, · · · , sµ, (36) M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1574 or η̄ζ(τ, ιl)− Z ( ζk, ιl, η̄(τ, ιl), η̄(τ − µ, ιl), η̄ιι(τ, ιl) ) ̸= 0, k = sµ+1, · · · , Q. (37) Assume that we have (36). Since the shifted LGL points are dense in [0,M ] when Q → ∞, redthere are subsequences {ζkQi }∞i=1 that 0 < kQi < Qi and limi→∞ ζkQi = τ . So we have lim i→∞ ( η̄ζ(ζkQi , ιl)− Z ( ζk, ιl, η̄(ζkQi , ιl), φ(ζkQi , ιl), η̄ιι(ζkQi , ιl) )) = η̄ζ(τ, ιl)− Z ( ζk, ιl, η̄(τ, ιl), φ(τ − µ, ιl), η̄ιι(τ, ιl) ) ̸= 0. On the other hand, since limi→∞ √ Qi 2Qi−1Y ( 1 2Qi−1) = 0, from (22) we have lim i→∞ ( η̄ζ(ζkQi , ιl)− Z ( ζk, ιl, η̄(ζkQi , ιl), φ(ζkQi − µ, ιl), η̄ιι(ζkQi , ιl) )) = 0 and this contradicts the relation (36). Similarly, if we have (37), lim i→∞ ( η̄ζ(ζkQi , ιl)− Z ( ζk, ιl, η̄(ζkQi , ιl), η̄(ζkQi , ιl), η̄ιι(ζkQi , ιl) )) = η̄ζ(τ, ιl)− Z ( ζk, ιl, η̄(τ, ιl), η̄(τ − µ, ιl), η̄ιι(τ, ιl) ) ̸= 0. Since limi→∞ √ Qi 2Qi−1Y ( 1 2Qi−1) = 0, so lim i→∞ ( η̄ζ(ζkQi , ιl)− Z ( ζk, ιl, η̄(ζkQi , ιl), η̄(ζkQi − µ, ιl), η̄ιι(ζkQi , ιl) )) = 0, which contradicts the relation (37). It is also easy to see that η̄(ζ, ι) for ι = ιl (l = 1, 2, .., Q − 1) and ζ ∈ [0,M ] satisfies conditions (2) and (3), by using the assumptions of theorem and the last three equations of system (22). Hence, η̄(ζ, ι) for ι = ιl (l = 1, 2, .., Q−1) and ζ ∈ [0,M ] satisfies in all equations of system (22). We know that shifted LGL points {ιk}Qk=1 are dense in [0, P ] when Q tends to infinity, thus η̄(ζ, ι) for all [ζ, ι] ∈ Γ̃ satisfies system (1)-(3). 5. Illustrative examples To show the effectiveness of suggested approach, we solve four numerical examples. The system (15) is solved using MATLAB software with Levenberg-Marqurdt algorithm. The absolute error of obtained approximation solution ηQ(., .) for exact solution η(., .) is described by ERQ(ζ, ι) = |η(ζ, ι)− ηQ(ζ, ι)|, (ζ, ι) = [0,M ]× [0, P ]. Also we define the ERQ 2 and ERQ ∞ errors of approximations in forms ERQ 2 = ( Q∑ r=0 Q∑ s=0 |η(ζr, ιs)− ηQ(ζr, ιs)|2 ) 1 2 , M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1575 ERQ max = max{|η(ζr, ιs)− ηQ(ζr, ιs)| : r, s = 0, 1, · · · , Q}. Moreover, for a given time ζ = ζ̄, we define errors ĒR Q 2 = ( Q∑ s=0 |η(ζ̄, ιs)− ηQ(ζ̄, ιs)|2) 1 2 , ĒR Q max = max{|η(ζ̄, ιs)− ηQ(ζ̄, ιs)| : s = 0, 1, · · · , Q}. Example 1. Let the following DRD equation ηζ(ζ, ι) = ηιι(ζ, ι) + η2(ζ, ι) + η(ζ − µ, ι) + f(ζ, ι), (ζ, ι) ∈ [0, 1]× [0, 1], (38) with the conditions  η(ζ, ι) = ζ2e2ι, (ζ, ι) ∈ [−µ, 0]× [0, 1] η(ζ, 0) = ζ2, ζ ∈ [0, 1], η(ζ, 1) = ζ2e2, ζ ∈ [0, 1], (39) where µ = 0.5 and f(ζ, ι) = 2ζe2ι − 4ζ2e2ι − ζ4e4ι + (ζ − µ)2e2ι. The exact solution of this equation is η(ζ, ι) = ζ2e2ι. We estimate the solution of this system for Q = 8 by the presented method. Figure 1 shows the results. Figure 2 illustrates the estimate errors. It can be seen that the errors goes to zero when Q increases. Example 2. Let the following DRD equation ηt(ζ, ι) = ηιι(ζ, ι) + η(ζ, ι)− e−µη(ζ − µ), (ζ, ι) ∈ [0, 2]× [0, π], (40) with the conditions  η(ζ, ι) = e−ζsinι, (ζ, ι) ∈ [−µ, 0]× [0, π], η(ζ, 0) = 0, ζ ∈ [0, 2], η(ζ, π) = 0, ζ ∈ [0, 2], (41) where µ = 1. The exact solution is η(ζ, ι) = e−ζsinι. We present the gained results for Q = 10 in Figure 3. The logarithm of ERQ 2 and ERQ max for Q = 10 are presented in Figure 4. We can see that the errors converge to zero by increasing Q. Also Table 1 shows the compared results and superiority of our method with respect to the compact finite difference method [15]. M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1576 Example 3. Let the following DRD equation ηζ(ζ, ι) = ηιι(ζ, ι) + 2η(ζ, ι) + η(ζ − µ, ι) 1 + η2(ζ − µ, ι) + f(ζ, ι), (ζ, ι) ∈ [0, 0.5]× [0, 2π], (42) with the conditions  η(ζ, ι) = e−ζsinι, (ζ, ι) ∈ [−µ, 0]× [0, 2π], η(ζ, 0) = 0, ζ ∈ [0, 0.5], η(ζ, 2π) = 0, ζ ∈ [0, 0.5], (43) where µ = 0.1 and f(ζ, ι) = 2e−ζsinι− e−(ζ−µ)sinι 1 + e−2(ζ−µ)sin2ι . Here, function η(ζ, ι) = e−ζsinι is the exact solution. We approximate the solution for Q = 10 by the presented method and show the estimated solution and its absolute error in Figure 5. Figure 6 illustrates the estimate errors. Also, we compare the errors ERQ max and ERQ 2 at ζ = 0.5 for presented method with those of compact finite difference method [16], that are told in Table 2. The result in Table 2 presents that the error of presented method is less than that of the method [16]. Example 4. Assume the following DRD equation ηζ(ζ, ι) = ηιι(ζ, ι) + π2η + η(ζ − µ, ι) + f(ζ, ι), (44) with the conditions  η(ζ, ι) = e−ζ2sin(πι), (ζ, ι) ∈ [−µ, 0]× [0, 1], η(ζ, 0) = 0, ζ ∈ [0, 1], η(ζ, 1) = 0, ζ ∈ [0, 1], (45) where µ = 0.5 and f(ζ, ι) = −2ζe−ζ2sin(πι)− e−(ζ−µ)2sin(πι), Function η(ζ, ι) = e−ζ2sin(πι) is the exact solution. We show the results for Q = 10 in Figure 7. Moreover, the logarithm of errors ERQ 2 and ERQ max of estimate solutions are presented in Figure 8, which tends to zero by increasing Q. M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1577 0 0.5 1 0 0.5 1 −2 0 2 4 6 8 ζ � ι� η Q (. ,. ) 0 0.5 1 0 0.5 1 −13 −12 −11 −10 −9 −8 ζι lo g 1 0 (E R Q (. ,. )) Figure 1: The ηQ(., .) and logarithm of ERQ(., .) for Q = 8 in Example 1. Table 1: Comparison of error ERQ max at time ζ = 2 for Example 2. Number of points Method [15] Number of points Presented method 4× 200 1.52× 10−6 7× 7 3.58× 10−6 8× 400 3.80× 10−7 8× 8 2.57× 10−7 16× 800 9.50× 10−8 9× 9 1.68× 10−8 32× 1600 2.37× 10−8 10× 10 9.02× 10−10 3 4 5 6 7 8 −9 −8 −7 −6 −5 −4 −3 −2 Q lo g 1 0 (E R Q m a x ) 3 4 5 6 7 8 −9 −8 −7 −6 −5 −4 −3 −2 Q lo g 1 0 (E R Q 2 ) Figure 2: The logarithm of ERQ max and ERQ 2 for Example 1. Table 2: The comparison of errors at ζ = 0.5 in Example 3. Q = 10 ĒR Q 2 ĒR Q max Our approach 4.16× 10−9 7.63× 10−8 Compact finite difference method [16] 3.15× 10−3 5.98× 10−2 M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1578 0 1 2 0 2 4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 ζι η Q (. ,. ) 0 1 2 0 2 4 −11 −10.5 −10 −9.5 −9 ζι lo g 1 0 (E R Q (. ,. )) Figure 3: The ηQ(., .) and logarithm of ERQ(., .) for Q = 8 in Example 2. 2 4 6 8 10 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 q lo g 1 0 (E R Q m a x ) 2 4 6 8 10 −9 −8 −7 −6 −5 −4 −3 −2 −1 Q lo g 1 0 (E R Q 2 ) Figure 4: The logarithm of errors ERQ max and ERQ 2 for Q = 10 in Example 2. M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1579 0 0.5 0 5 10 −1 −0.5 0 0.5 1 ζι η Q (. ,. ) 0 0.5 0 5 10 −12 −11 −10 −9 −8 −7 ζι lo g 1 0 (E R Q (. ,. )) Figure 5: The ηQ(., .) and logarithm of ERQ(., .) for Q = 10 in Example 3. 2 4 6 8 10 −8 −7 −6 −5 −4 −3 −2 −1 Q lo g 1 0 (E R Q m a x ) 2 4 6 8 10 −7 −6 −5 −4 −3 −2 −1 0 Q lo g 1 0 (E R Q 2 ) Figure 6: The logarithm of errors ERQ max and ERQ 2 in Example 3. M. Ghovatmand et al. / Eur. J. Pure Appl. Math, 17 (3) (2024), 1565-1584 1580 0 0.5 1 0 0.5 1 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 ζι η Q (. ,. ) 0 0.5 1 0 0.5 1 −10.5 −10 −9.5 −9 −8.5 −8 ζι lo g 1 0 (E R Q (. ,. )) Figure 7: The ηQ(., .) and logarithm of ERQ(., .) for Q = 10 in Example 4. 2 4 6 8 10 −9 −8 −7 −6 −5 −4 −3 −2 −1 Q lo g 1 0 (E R Q m a x ) 2 4 6 8 10 −8 −7 −6 −5 −4 −3 −2 −1 0 Q lo g 1 0 (E R Q 2 ) Figure 8: The logarithm of errors ERQ max and ERQ 2 in Example 4. REFERENCES 1581 6. Conclusions In this article, we showed that pseudo-spectral methods with two-dimensional Lagrange- bases at the collocation points of Legendre-Gauss-Lobatto can be used for reaction-diffusion equations of the difference delayed type. The proposed technique was to insert the relevant delay function in the diffusion-reaction equation and converting it into a non-delayed equa- tion. This technique can be used for other delayed differential equations. The convergence analysis of the method showed that the approximate solutions obtained by considering mild conditions converge to the exact solution of the equation. Also, the use of modulus of continuity functions and its relevant space of continuously differentiable functions for convergence analysis of the proposed method can be extended to other pseudo-spectral methods and other time-continuous problems. For future research, we will use the method and its convergence analysis for fractional delayed/non-delayed one and two-dimensional diffusion-reaction equations. Conflicts of interest: This work does not have any conflicts of interest. Funding: There are no funders to report for this submission. AI statement: The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article. References [1] Ali, I., Rasool, G., & Alrashed, S. (2018). Numerical simulations of reaction diffusion equations modeling prey–predator interaction with delay. International Journal of Biomathematics, 11(04), 1850054. [2] Baranovskii, E. S., Brizitskii, R. V., Saritskaia, Z. V., 2024, Optimal control problems for the reaction–diffusion–convection equation with variable coefficients, Nonlinear Analysis: Real World Applications, 75, 103979. [3] Cardin, F., Favretti, M., Lovison, A., & Masci, L. (2019). Stochastic and geometric aspects of reduced reaction–diffusion dynamics. Ricerche di Matematica, 68(1), 103- 118. [4] Chen, L., Wei, H., & Wen, M. (2017). An interface-fitted mesh generator and virtual element methods for elliptic interface problems. Journal of Computational Physics, 334, 327-348. [5] Doha, E. H., & Bhrawy, A. H. (2008). Efficient spectral-Galerkin algorithms for di- rect solution of fourth-order differential equations using Jacobi polynomials. Applied Numerical Mathematics, 58(8), 1224-1244. REFERENCES 1582 [6] Faria, T., & Trofimchuk, S. (2006). Nonmonotone traveling waves in a single species reaction–diffusion equation with delay. Journal of Differential Equations, 228(1), 357- 376. [7] Gourley, S. A., & Kuang, Y. (2004). A delay reaction-diffusion model of the spread of bacteriophage infection. SIAM Journal on Applied Mathematics, 65(2), 550-566. [8] Huang, Y., Mohammadi Zadeh, F., Noori Skandari, M.H., Ahsani Tehrani, H., Tohidi E. (2021), Space–time Chebyshev spectral collocation method for nonlinear time- fractional Burgers equations based on efficient basis functions, 44, 5, 4117-4136. [9] Heidari, M., Ghovatmand, M., Noori Skandari, M.H., Baleanu, D. (2022) Numerical Solution of Reaction–Diffusion Equations with Convergence Analysis. J Nonlinear Math Phys, 30, 384–399. [10] Jafari, H., Mahmoudi, M., Noori Skandari, M.H., (2021), A new numerical method to solve pantograph delay differential equations with convergence analysis. Adv Differ Equ 2021, 129, https://doi.org/10.1186/s13662-021-03293-0. [11] Kirthiga, M., Balamurugan, S., Visuvasam, J., Rajendran, L., and Marwan, A., (2021). Theoretical Analysis of Single-Stage and Multi-Stage Monod Model of Landfill Degradation Through Mathematical Modelling, 7( 1), 48-62. [12] Koto, T. (2008). IMEX Runge–Kutta schemes for reaction-diffusion equations. Jour- nal of Computational and Applied Mathematics, 215(1), 182-195. [13] Kuttler, C. (2017). Reaction diffusion equations and their application on bacterial communication. In Handbook of Statistics, 37, 55-91, Elsevier. [14] Lan, K. Q., & Wu, J. H. (2003). Traveling wavefronts of scalar reaction-diffusion equations with and without delays. Nonlinear analysis: real world applications, 4(1), 173-188. [15] Li, D., Zhang, C., & Wen, J. (2015). A note on compact finite difference method for reaction–diffusion equations with delay. Applied Mathematical Modelling, 39(5-6), 1749-1754. [16] Li, D., Zhang, C., & Qin, H. (2011). LDG method for reaction-diffusion dynamical systems with time delay. Applied Mathematics and Computation, 217(22), 9173-9181. [17] Lin, G., & Li, W. T. (2009). Traveling wavefronts of Belousov-Zhabotinskii system with diffusion and delay. Applied mathematics letters, 22(3), 341-346. [18] Ling, Z., Lin, Z., (2010). Traveling wavefront in a hematopoiesis model with time delay, Appl. Math. Lett. 23, 426–431. REFERENCES 1583 [19] Liu, B., & Zhang, C., (2015). A spectral Galerkin method for nonlinear delay convection-diffusion-reaction equations. Computers & Mathematics with Applica- tions, 69(8), 709-724. [20] Liu, L., & Nieto, J. J., (2023). Dynamics of Fractional Delayed Reaction-Diffusion Equations. Entropy, 25(6), 950. [21] Lv, G., & Wang, M., (2012). Nonlinear stability of traveling wave fronts for nonlocal delayed reaction–diffusion equations. Journal of Mathematical Analysis and Applica- tions, 385(2), 1094-1106. [22] Meleshko, S. V., & Moyo, S., (2008). On the complete group classification of the reaction-diffusion equation with a delay. Journal of mathematical analysis and appli- cations, 338(1), 448-466. [23] Polyanin, A. D., & Zhurov, A. I., (2014). Functional constraints method for construct- ing exact solutions to delay reaction-diffusion equations and more complex nonlinear equations. Communications in Nonlinear Science and Numerical Simulation, 19(3), 417-430. [24] Peykrayegan, N., Ghovatmand, M., Noori Skandari, M.H., (2021), An efficient method for linear fractional delay integro-differential equations. Comp. Appl. Math. 40, 249. https://doi.org/10.1007/s40314-021-01640-1. [25] Priyanka, P., Arora, S.,, Mebrek-Oudina, F., Sahani, S., (2023). Super conver- gence analysis of fully discrete Hermite splines to simulate wave behaviour of Ku- ramoto–Sivashinsky equation, Wave Motion, 121, 103187. [26] Priyanka, P., Mebarek-Oudina, F., Sahani S., Arora, S., (2024). Travelling wave solution of fourth order reaction diffusion equation using hybrid quintic hermite splines collocation technique, Arabian Journal of Mathematics , https://doi.org /10.1007/s40065-024-00459-y. [27] Ragozin, D. L., (1970). Polynomial approximation on compact manifolds and homo- geneous spaces. Transactions of the American Mathematical Society, 150(1), 41-53. [28] Sharifi, S., & Rashidinia, J., (2019). Collocation method for Convection-Reaction- Diffusion equation. Journal of King Saud University-Science, 31(4), 1115-1121. [29] Shen, J., Tang, T., & Wang, L. L., (2011). Spectral methods: algorithms, analysis and applications (Vol. 41). Springer Science & Business Media. [30] Sorokin, V. G., & Vyazmin, A. V., (2022). Nonlinear Reaction–Diffusion Equations with Delay: Partial Survey, Exact Solutions, Test Problems, and Numerical Integra- tion. Mathematics, 10(11), 1886. [31] Visuvasam J., and Alotaibi H., Analysis of Von Karman swirling flows due to a porous rotating disk electrode, Micromashines, 14, 582. REFERENCES 1584 [32] Visuvasam J., Molina, A., Laborda, and E., Rajendran, L., (2018). Mathematical Models of the Infinite Porous Rotating Disk Electrode, Int. J. Electrochem. Sci., 13, 9999-10022. [33] Visuvasam, J., Meena, A., and Rajendran, L., (2020). New analytical method for solving nonlinear equation in rotating disk electrodes for second-order ECE reactions, 869, 15, 114106. [34] Wang, W., & Ma, W. (2018). Hepatitis C virus infection is blocked by HMGB1: A new nonlocal and time-delayed reaction–diffusion model. Applied Mathematics and Computation, 320, 633-653. [35] Wang, X., Li, Z., (2007). Dynamics for a type of general reaction–diffusion model, Nonlinear Anal. 67, 2699–2711. [36] Wu, F., Wang, Q., Cheng, X., & Chen, X. (2018). Linear Θ-Method and Compact Θ- Method for Generalized Reaction-Diffusion Equation with Delay. International Jour- nal of Differential Equations, 2018. [37] Xiaobing, P., Yang, X., Noori Skandari, M. H., Tohidi, E., Shateyi, S. (2021). A new high accurate approximate approach to solve optimal control problems of fractional order via efficient basis functions, Alexandria Engineering Journal, 61 (8), 5805-5818. [38] Zhao, Z., & Ge, W. (2011). Traveling wavefront solutions for reaction-diffusion equa- tion with small delay. Funkcialaj Ekvacioj, 54(2), 225-236.