EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 3, 2024, 1908-1936 ISSN 1307-5543 – ejpam.com Published by New York Business Global Modeling Rugose Spiraling Whitefly Infestation on Coconut Trees Using Delay Differential Equations: Analysis via HPM B. Dhivyadharshini1, R. Senthamarai1,* 1 Department of Mathematics, College of Engineering and Technology, SRM Institute of Science and Technology, Kattankulathur−603203, Tamilnadu, India Abstract. In this paper, a delay - induced pest control model is proposed. We have introduced a time delay in healthy trees and whitefly population in the infected tree density of the proposed system of equations to reduce the probability of healthy trees becoming infected, as well as the level of infection. We have analyzed the impact of time delay on the stability of the equilibrium and establish requirements to verify its asymptotic stability over all delays. The solutions of this system of non-linear ordinary differential equations(ODEs) and delay differential equations(DDEs) are presented by using homotopy perturbation method(HPM). Numerical simulation is also obtained for the same model of both ODE and DDE by using MATLAB software. The proposed method works very well and is easy to use, as shown by these findings. Our aim is to control the spread of whitefly such that harvest is not affected. 2020 Mathematics Subject Classifications: 65L05, 34A34, 92D45, 37N25, 92D25 Key Words and Phrases: Pest control model, Non linear differential equation, Delay differential equation, Homotopy perturbation method, Numerical simulation. 1. Introduction Rugose spiralling whitefly (RSW) is a highly polyphagous and recently introduced which is thought to have first appeared in Central America. Its presence is confined to Belize, Mexico, Guatemala, and Florida in Central and North America, and it has extended to neighbouring coun- tries in the oriental region which rise coconuts. Martin originally introduced RSW in 2004 using samples gathered from the leaves of coconut palms in Belize [11]. In India RSW has established its presence in Kerala (Palakkad, Malappuram, Thrissur, Idukki, Kozhikode, Kannur, Ernakulam, Kasaragod, Pathanamthitta, Alappuzha, Kollam, and Thiruvananthapuram districts), Tamil Nadu (Pollachi and Pattukottai), Karnataka (Udupi), and secluded areas of Andhra Pradesh (Kadiyam) have seen its expansion during a six−month period (August 2016 − January 2017) [6]. RSW was initially discovered on coconut palm leaves in the Coimbatore District of Tamil Nadu [22]−[26]. It was then appeared on several plants, namely monocots and dicots, in Karnataka, Kerala, Andhra Pradesh, and Assam. One of the most significant palm crops in tropical, subtropical, and warm ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v17i3.5249 Email addresses: senthamr@srmist.edu.in (R. Senthamarai), db0558@srmist.edu.in (B. Dhivyadharshini) https://www.ejpam.com 1908 © 2024 EJPAM All rights reserved. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1909 temperate regions is the coconut (Cocos nucifera). It is commonly known as the “tree of heaven” or “kalpavriksha” due to the fact that it offers a wider range of beneficial products to the public. More than 90 countries, mostly in Asia, grow coconuts. Covering an area of 12 million hectares, the Pacific Islands and South America are predicted to produce 70 billion nuts annually. Approx- imately 50 % of its yearly production is used in India. Many pests present and create different kinds of problems when growing these coconut trees. In August of 2016, reports of the dangerous invasive pest on coconuts (Cocos nucifera L.) have been detected in Pollachi, Tamil Nadu, India [22]−[26]. In south India’s coconut and oil palm growing regions, the pest’s entrance has been taken severely due to its polyphagous nature and ability to consume over 200 host plants. Because of the way this whitefly feeds, the leaves lose water and nutrients. This has an impact on the host tree. Additionally, it departure behind sooty mould, which envelops the leaf surface and may in- hibit photosynthesis, so affecting the production and growth [12]. It has been explained how mathematical models are developed and used in epidemiology [2]. In agriculture, pest control is essential for producing healthy, high-yield crops [18]. The ecology of pest control is complicated, involving complex interactions between plants, pests, natural enemies, other living things, and the surrounding environment. For the African cassava mosaic virus disease (ACMD), the dynamics of plant and vector populations within a region have been investigated. A class of models that connects virus epidemiology with vector dynamics for ACMD has not yet been utilised, and it is based on a differential equation system [16]. The effect of incubation delay in plant-vector interaction has been examined [19]. Numerous authors have included various forms of time delays in biological models [21]−[4]. Compared to ordinary differential equations, delay-differential equations typically display for more complex dynamics especially when the complexity of the delay terms increase. For in- stance, in biology many changes don’t occur instaneously, that is, think about the birth of human. The birth rate varies based on available resources at the period. Such kind of situations were cap- tured by DDE’s [17] and [1]. A DDE is an ODE in which derivative depends on past values of the state at time t, the evolution of system depends on the current time t, current state of system and at some time τ > 0 in the past. In this paper we have considered a constant DDE since the time delay parameter τ is constant. A generalized delay-induced epidemic model has been analyzed, taking into account a nonlinear incidence rate, latency, and relapse [29]. Delay-differential equa- tion modeling has been investigated as a means of explaining HIV infection of CD4+ T−cells [7]. As far as we are aware, no delay differential equations system exists that simulates RSWs that impact coconut trees. Inspired by its concept of plant-vector interaction [19], we have investigated the disease dynamics with the time delay using the mathematical modeling concept and examined the linearized system’s transcendental characteristic equation. We evaluated the parameter values primarily in light of the Pollachi tract in Tamil Nadu. Its stability was noted at the equilibrium points. A next-generation matrix was used to determine the reproduction number. For the sensi- tivity analysis, we investigated the factors influencing the system. Multiple techniques exist for resolving this type of non−linear issue. Various analytical techniques, such as the variational it- eration method, Adomian decomposition method, homotopy analysis method and others, can be employed to get approximate solutions for nonlinear differential equations [28]−[23]. The ho- motopy perturbation method (HPM), was introduced by Dr. Ji Huan He (1999b) has successfully been applied to solve numerous types of linear and nonlinear functional equations [14] and [15]. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1910 Using the homotopy perturbation method, an approximate analytical solution has been found for the proposed system that proved it’s suitability for getting solutions that are valid for all param- eters appeared on the system [8]−[3]. It has a substantial benefit in that provides an analytical approximation solution to a large range of nonlinear problems in applied sciences [20]−[24]. The parameter values were taken from [25] based on the facts and methods. Numerical simulation is also acquired by using MATLAB software. It is determined that the analytical and numerical simulation outcomes are in perfect agreement. 2. Mathematical formulation In order to examine the effects of RSWs on coconut trees, a plant−vector interaction model [16] was promptly created by taking into account the whitefly population and coconut plants. The population of trees was splitted into two groups: healthy trees (X) and infected trees (Y ). C stands for the number of whiteflies density per square meter. 2.1. The system of ordinary differential equations The mathematical model is given as follows [8]: dX dt = bX ( 1− X +Y v ) −aXC, (1) dY dt = aXC−µY, (2) dC dt = wY − sC. (3) with the initial conditions X(0) = l > 0, Y (0) = m > 0, C(0) = n > 0. (4) 2.2. The system of delay differential equations The Mathematical model with delay is given as follows: dX(t) dt = bX(t) ( 1− X(t)+Y (t) v ) −aX(t)C(t), (5) dY (t) dt = aX(t − τ)C(t − τ)−µY (t), (6) dC(t) dt = wY (t)− sC(t). (7) with the initial values X(θ) = l, Y (θ) = m, C(θ) = n θ ∈ [−τ,0] (8) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1911 3. Analysis of the model 3.1. Positivity The following concise notation can be used to express the system of Eqs.(5)−(7): dU dt = ψ (U) . Here, U (θ) = (X(θ),Y (θ),C(θ))T ∈ Z where Z = Z ( [−τ,0] ,R3 + ) represents the Banach space of continuous functions, and ψ = (ψ1,ψ2,ψ3) T are the right sides of the system of Eqs.(5)−(7). According to the fundamental theory of functional differential equations [13], there exists a unique solution (X(t),Y (t),C(t)) to the system of Eqs.(5)−(7) with initial conditions given by Eq.(8). Theorem 1. All the solutions of Eqs.(5)−(7) with initial conditions Eq.(8) are positive. Proof. The principle mentioned above has been proven using the approach proposed by Bodnar, as described in [5] and [30]. It is simple to verify in system of Eqs.(5)−(7) that when selecting X(θ) ∈ R+ that if X = 0,Y = 0,C = 0, then ψi(U)|Ui=0,U∈R3 + ≥ 0. By applying Lemma 2 in [30] and Theorem 1.1 in [5], we can conclude that any solution x(t) = x(t,x(θ)) of Eqs.(5)−(7) with x(θ) ∈ Z satisfies U(t) ∈ R3 + for all t ≥ 0. Therefore, the solution to the system (5)−(7) occurs within the region R3 +, and all solutions remain non-negative for all t > 0. Thus, the positive cone R3 + is an invariant region. 3.2. Boundedness The total tree population N = X +Y satisfies dX dt + dY dt = bX ( 1− X +Y v ) −aXC+aX(t − τ)C(t − τ)−µY, d (X +Y ) dt +ηX +ηY = bX ( 1− X +Y v ) −aXC+aX(t − τ)C(t − τ)−µY +ηX +ηY, dN dt +ηN ≤−bX ( N v ) −aXC+aX(t − τ)C(t − τ)+(b+η)X +(η −µ)Y, dN dt +ηN ≤−bX2 v +(b+η)X , Here η < µ . It seems clear that −bX2 v +(b+η)X is quadratic in X and its maximum value is (b+η)2v 4b . dN dt +ηN ≤ l, B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1912 where l = (r+η)2v 4b , 0 ≤ N(t)≤ e−ηt ( N(0)− l η ) + l η (9) AS t−→∞, N(t)−→ l η since supt−→∞ C(t)= 1 wη as a result of using the bound of Y. Therefore, the following positive invariant set represents the biologically viable range of the sys- tem represented by Eqs.(5) – (7). Ω = { (X ,Y,C) ∈ R3 +|0 ≤ X ,Y ≤ l η ,C ≤ l wη } 4. Reproduction number If τ = 0, disease class Eqs.(6), (7) describe the population inputs that are related to the size of the population and the number of initial infections. 4.1. Disease class F = [ aXC 0 ] V = [ −µY wY − sC ] Constructing F matrix Let f (Y1c) = aXC, g(Y1c) = 0 F = [ ∂ f ∂Y ∂ f ∂C ∂g ∂Y ∂g ∂C ] = [ 0 aX 0 0 ] ∂ f ∂Y = 0 ∂ f ∂C = aX ∂g ∂Y = 0 ∂g ∂C = 0 Constructing V matrix Let f (Y1c) = −µY , g(Y1c) = wY − sC V = [ ∂ f ∂Y ∂ f ∂C ∂g ∂Y ∂g ∂C ] = [ −µ 0 w −s ] ∂ f ∂Y =−µ ∂ f ∂C = 0 ∂g ∂Y = w ∂g ∂C =−s |FV−1 −λ I|= 0 (10) By substituting and solving Eq.(10) using next-generation matrix invented by Diekmann et al. [10]. By using it, we can find the dominant eigen value which is called the basic reproduction number (R0), since avw µs is the dominant eigen value. R0 = avw µs (11) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1913 5. Analysis of the system’s stability with delay We express linearized system of Eqs.(5)−(7) in matrix form as follows: d dt X(t) Y (t) C(t) = A1 X(t) Y (t) C(t) +A2 X(t − τ) Y (t − τ) C(t − τ)  , where A1 and A2 are, A1 = b ( 1− X+Y v ) − bX v −aC −bX v −aX 0 −µ 0 0 w −s  , A2 =  0 0 0 aC 0 aX 0 0 0  , (12) The system’s characteristic equation can be expressed as △(λ ) = |λ I −A1 − e−λτA2|= 0, (13) that is, the system given by Eqs.(5)−(7) possess three equilibria: (i) Trivial equilibrium: E0 = (0,0,0) (ii) Pest - free equilibrium: E1 = (v,0,0) (iii) Coexistence equilibrium: Ē = (X∗,Y ∗,C∗) whereas, X∗ = µs aw Y ∗ = sb(avw−µs) aw(sb+avw) C∗ = b(avw−µs) a(sb+avw) (14) 5.1. Trivial Equilibrium: E0 = (0,0,0) Lemma 1. The system of non−linear delay differential equations Eqs.(5)−(7) at the point E0 = (0,0,0) is always unstable. 5.2. Pest - free equilibrium: E1 = (v,0,0) By substituting the equilibrium points in the above matrixces Eq.(12), we can get A1 and A2: A1 = −b −b −av 0 −µ 0 0 w −s  , A2 = 0 0 0 0 0 av 0 0 0  , (15) Substituting A1 and A2 in Eq.(13), the transcendental equation of the matrix is ψ(λ ,τ) = λ 3 +D1λ 2 +D2λ +D3 + e−λτ [B1λ +B2] = 0, (16) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1914 with coefficients: D1 = µ + s, D2 = µs+b, D3 = bµs, B1 =−avw, B2 =−abvw. Theorem 2. The system given by Eqs.(5)−(7) around the pest-free equilibrium E1 is locally asymptotically stable (LAS), provided R0 < 1. Proof. At E1 then τ = 0 in Eqn.(16) becomes, λ 3 +D1λ 2 +D2λ +D3 +B1λ +B2 = 0 which possesses the form λ 3 +D1λ 2 +(D2 +B1)λ +D3 +B2 = 0 The system is LAS according to the Routh-Hurwitz (R-H) criterion if D1 > 0, D3 +B2 > 0 and D1 (D2 +B1)>D3+B2. Hence Pest - free equilibrium E1 is locally asymptotically stable if avw< µs i.e., R0 < 1. 5.3. Coexistence equilibrium: Ē = (X∗,Y ∗,C∗) By Substituting the coexistence equilibrium points in the matrices (12). We can get A1 and A2: A1 = b ( 1− X∗+Y ∗ v ) − bX∗ v −aC∗ −bX∗ v −aX∗ 0 −µ 0 0 w −s  , A2 =  0 0 0 aC∗ 0 aX∗ 0 0 0  , (17) Substituting A1 and A2 in (13), the transcendental equation of the matrix is ψ(λ ,τ) = λ 3 +D1λ 2 +D2λ +D3 + e−λτ [B1λ +B2] = 0, (18) with coefficients: D1 =C∗a+ 2X∗b v + X∗b v −b+µ + s, D2 =C∗aµ +C∗as+ 2X∗bµ v + 2X∗bs v + Y ∗bs v −bs+µs, D3 =C∗aµs+ 2X∗bµs v −bµs. B1 = C∗X∗ab v −X∗aw, B2 =C∗X∗a2w−C∗X∗a2w+ C∗X∗abs v − 2(X∗)2 abw v − X∗Y ∗abw v +X∗abw. It is very difficult to determine the sign of roots when τ is present. If there are negative real parts in each of the roots of the related characteristic equation, then Ē is known to be locally asymptotically stable. If a root has a positive real part, it is unstable. If there are only purely imaginary roots, B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1915 stability flips. There are infinitely many complex roots for the transcendental equation Eq.(18). As a result, we start our analysis with the delay set to zero and subsequently determine the stability requirements when τ > 0. Case I: when τ = 0 Eq.(18) becomes, ψ(λ ,τ) = λ 3 +D1λ 2 +(D2 +B1)λ +(D3 +B2) = 0, (19) According to the Routh−Hurwitz criterion, all eigenvalues of Eq.(19) have negative real parts if only if D1 > 0, D3 +B2 > 0, D1 (D2 +B1)− (D3 +B2)> 0. (20) All the criteria in Eq.(20) are satisfied, and the infected steady state Ē is asymptotically stable. Case II: when τ > 0 Eq.(18) has an infinitely many roots. The roots of the characteristic Eq.(18) with negative real parts are necessary for stability. If there are only purely imaginary solutions to the characteristic Eq.(18), then Ē experiences stability changes. Assume that a root of Eq.(18) is iω . From it, we obtain −iω3 −D1ω 2 + iD2ω +B1ω (sinωτ + icosωτ)+B2 (cosωτ − isinωτ)+D3 = 0. (21) When we separate the real and imaginary parts, we get D1ω 2 −D3 = B1ω sinωτ +B2 cosωτ, (22) ω 3 −D2ω = B1ω cosωτ −B2 sinωτ. (23) The above two equations can be squared and added to get ω 6 + ( D2 1 −2D2 ) ω 4 + ( D2 2 −2D1D3 −B2 1 ) ω 2 + ( D2 3 −B2 2 ) = 0. (24) Let z = ω 2, ζ1 = D2 1 −2D2, ζ2 = D2 2 −2D1D3 −B2 1, ζ3 = D2 3 −B2 2. Then Eq.(24) becomes, h(z) = z3 +ζ1z2 +ζ2z+ζ3 = 0. (25) Given that ζ3 = D2 3 −B2 2 > 0 for the parameter values shown in Table 1. We assume that ζ3 ≥ 0 and make the following claim. Claim 1. If ζ3 ≥ 0 (26) and ζ2 > 0, (27) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1916 Therefore, there are no positive real roots in Eq. (25). In fact, observe that dh(z) dt = 3z2 +2ζ1z+ζ2. Set 3z2 +2ζ1z+ζ2 = 0. (28) Then the roots of Eq.(29) can be represented as z1,2 = −ζ1 ± √ ζ 2 1 −3ζ2 3 . (29) If ζ2 > 0, then ζ 2 1 −3ζ2 < ζ 2 1 ; that is, √ ζ 2 1 −3ζ2 < ζ1. Hence, neither z1 nor z2 is positive. Thus, Eq.(25) has no positive roots. As a result, Claim 1 suggests that there is no ω such that iω is an eigenvalue of the characteristic Eq.(18). Consequently, for all delay of τ ≥ 0, the real parts of all eigenvalue of Eq.(18) are nega- tive. Thus, the analysis presented above leads to the following Theorem. Theorem 3. Suppose that (i) D1 > 0, B1 +B2 > 0, D1 (D2 +B2)− (B1 +D3)> 0; (ii) ζ3 ≥ 0 and ζ2 > 0. Then the infected steady state Ē of the delay model Eqs.(5)−(7) is perfectly stable; that is, Ē is asymptotically stable for all τ ≥ 0. Proof. Observe that all the conditions in Theorem 3 is satisfied for the specified parameter values in Table 1. For all τ ≥ 0, the infected steady state Ē is thus asymptotically stable. Remark. Theorem 3 states that if the parameters satisfy the conditions (i) and (ii), then the steady state of the delay model Eqs.(5)−(7) is asymptotically stable for all delay values; that is, independent of the delay. However, it is important to note that if the conditions (condition (ii)) in Theorem 3 are not satisfied, then the stability of the steady state depends on the delay value and the delay may even cause oscillations. An example, if (a) ζ3 < 0, then from Eq.(25) we have h(0) < 0 and limt→∞ h(z) = ∞. Thus, Eq.(25) has at least one positive root, say z0. Consequently, Eq.(24) has at least one positive root, denoted by ω0. If (b) ζ2 < 0, then √ ζ 2 1 −3ζ2 > ζ1. By Eq.(29), z1 = 1 3 ( −ζ1 + √ ζ 2 1 −3ζ2 ) > 0. Accordingly, Eq.(25), hence Eq.(24), has a positive root ω0. This suggests that there are pair of purely imaginary roots to the characteristic Eq.(18) ±iω0. Let λ (τ) = η (τ)+ iω (τ) represent the eigenvalue of Eq.(18) so that η (τ0) = 0,ω (τ0) = ω0. By means of Eqs.(22) and (23) we have cosωτ = 1 △ ∣∣∣∣D1ω2 −D3 B2ω ω3 −D2ω −B1 ∣∣∣∣ B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1917 = ( − 1 △ )( B2ω 4 +(D1 −D2B2)ω 2 −D3B1 ) , where △= ∣∣∣∣−B1 B2ω B2ω B1 ∣∣∣∣ = B2 1 +B2 2ω 2 > 0. sinωτ = 1 △ ∣∣∣∣B2ω D1ω2 −D3 −B1 ω3 −D2ω ∣∣∣∣ = B2ω 4 +(B1D1 −D2B2)ω 2 −B1D3. τ j = 1 ω0 ( cos−1 ( B2ω4 +(B1D1 −D2B2)ω2 −B1D3 B2 1 +B2 2ω2 ) +2 jπ ) , j = 0,1,2,3,4, ..., and τ0 = 1 ω0 cos−1 ( B2ω4 +(B1D1 −D2B2)ω2 −B1D3 B2 1 +B2 2ω2 ) , j = 0. The set of ordered pair is (ω0,τ0). Furthermore, we are able to confirm that the subsequent transversality requirement: d dt Reλ (τ)| τ=τ0 = d dt η (τ) ∣∣∣∣ τ=τ0 (30) holds. By continuity, when τ > τ0, the real part of λ (τ) becomes positive, making the steady state unstable. This leads to the bifurcation at τ = τ0, as shown by Eqs.(5)−(7). 6. Analysis of sensitivity parameters This section provides information on how changing parameter values affect the functional value of the reproduction number R0. The crucial parameter needs to be identified since it may serve as a significant threshold for the treatment of diseases. The following are the algebraic representations of the R0 sensitivity index to the parameters v,s,b,a,w, µ: ∂R0 ∂v = µsaw (µs)2 , ∂R0 ∂ s = −avwµ (µs)2 , ∂R0 ∂b = 0, ∂R0 ∂a = µsvw (µs)2 , ∂R0 ∂w = µsav (µs)2 , ∂R0 ∂ µ = −avws (µs)2 . The conclusion is that there exist positive partial derivatives, and that the basic reproductive num- ber R0 grows as any of the above positive value parameters v,a,w increase. The proportional response to the proportional perturbation is used to estimate the elasticity. We’ve Ev = v R0 ∂R0 ∂v = ( vµs avw )( µsaw (µs)2 ) = 1, Ea = a R0 ∂R0 ∂a = (aµs avw )( µsvw (µs)2 ) = 1, B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1918 Ew = w R0 ∂R0 ∂w = (wµs avw )( µsav (µs)2 ) = 1. It can be seen from the following expressions that Ev,Ea, and Ew are all positive. This suggests that increasing the values of the parameters v, a, and w increases the value of the basic reproduction number R0. A large variation in the fundamental reproduction number might result from even the smallest adjustment in these parameters. It is important to carefully calculate an extremely sensitive parameter since even a small variation might cause significant quantitative changes in the system. Figure 1: Nature and extent of infestation in the coconut plant 7. Mathematical analysis 7.1. Approximate analytical solution of ODE Eqs.(1)−(4) using homotopy perturba- tion method Several authors using the homotopy perturbation method [8] and [9] to produce approximate analytical solution for different non-linear problems established the effectiveness of the homotopy perturbation method for solving various engineering, physical, chemical and biological problems. This method plays a vital role in bio-mathematical sciences. We obtain the following analytical solutions of the of Eqs.(1)−(4) by using homotopy perturbation method. X(t) = ( l − l(lebt−mbe−µt µ − vane−µt s ) v + l(l − mb µ − van s ) v ) ebt B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1919 + ( 1 v2µs(µ +b− s)(−s+µ) ( l(µ +b− s) ( bs(−s+µ) ( lµe−µt b − vm2be−sµt 2µ )) − 1 µ (( −ms(−s+µ)(vm+1)b3 − (((1+(−l +m)v)s)+ v2an)µ − s2(vm+1))(−s+µ)mb2 −µv(m(van− ls)µ2)−2 ( − lsm+an ( vm− l 2 )) sµ +(−ls2m+an(vm− l)s+ vamw)s ) b− v2amµws(−s+µ)ebt ) + tmb2s3 + tlµb2s2 + tmµb3s+ tmµ 2b2s+ tavnµ 2b2 + tavnµ 3b+ tavnµbs2 + valnµbs(−s+µ)e(b−s)t b− s + v2amnµb(−s+µ)(µ +b− s)e−(µ+s)t −µ − s − tmb3s2 − tlµbs3 +2tlµ2bs2 − tlµ2b2s− tlµ3bs−2tmµb2s2 − tavnµb2s−2tavnµ 2bs− vlmµbs(−s+µ)(µ +b− s)e−(µ−b)t −µ +b )) − 1 v2µs(µ +b− s)(−s+µ) ( l ( (µ +b− s) ( bs(−s+µ)( lµ b − vm2b 2µ ) − µ(−n(−s+µ)b+ vmws)av s ) − 1 µ ( −ms(−s+µ)(vm+1)b3 − (((1+(−l +m)v)s+ v2an)µ − s2(vm+1))(−s+µ)mb2 −µv ( m(van− ls)µ2 −2 ( − lsm+an ( vm− 1 2 )) sµ +(−ls2m +an(vm− l)s+ vamw)s)b− v2amµws(µ − s) ) + valnµbs(µ − s) b− s + v2amnµb(−s+µ)(µ +b− s) −µ − s − vlmµbs(−s+µ)(µ +b− s) b−µ ))) ebt (31) Y (t) = ( m+ a lnet(µ+b−s) µ +b− s − lna µ +b− s ) e−µt + ( 1 (µ +b− s)vµs ( lna ( ((van− ls)µ +mbs)e(µ+b−s)t + lµs(va+µ +b− s)e(µ+2b−s)t µ +2b− s − s ( mb(µ +b− s)e(b−s)t b− s + valµebt b ) − vnµa(µ +b− s)e(µ+2b−s)t µ +b−2s )) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1920 − lna ( (van− ls)µ +mbs+ lµs(va+µ+b−s) µ+2b−s − s ( mb(µ+b−s) b−s + valµ b ) − vnµa(µ+b−s) µ+b−2s ) (µ +b− s)vµs ) e−µt (32) C(t)= ( n− wme−t(µ−s) µ − s + wm µ − s ) e−st + −q2m ( −t + e−t(p−s) p−s ) p− s + q2m (p− s)(s− p) e−st (33) 7.2. Approximate analytical solution of DDE Eqs.(5)−(8) using homotopy perturba- tion method By applying homotopy perturbation method [27], we are able to derive the following analytical solutions for the delay system of Eqs.(5)−(8). X(t) = ( l − l(lebt − mbe−µt µ − vane−µt s ) v + l(l − mb µ − van s ) v ) ebt + ( 1 µsv2(µ +b− s)(−s+µ) ( l ( tavnµ 3b+tavnµ 2b2+tavnµbs2 + tmµ 2b2s+ tmb2s3 + avnµs(−s+µ)e(−b+s)τ−t(µ+b) −µ −b +(−µ −b + s) ( − m(e−µt µb2s+ vb2s2(−e−µt − e−µt)+ e−µtav2nµ2b+ e−µtvmµb2s − e−µts2b2 − e−µtvmb2(s2 − vµb2s(−e−µt − e−µt)− e−µtav2nµbs µ ) −b2s(−s+µ) ( µt2 2 − vm2e−2µt 2µ ) + vµa(−nbµ + s(vmw+nb))e−st s − av2mnµb(−s+µ)e−t(µ+s) −µ − s ) − tmb3s2 −2tmµb2s2 − tavnµb2s −2tavnµ 2bs− avlnµbs(−s+µ)e(t−τ)(b−s) b− s )) − 1 µsv2(µ +b− s)(−s+µ) ( l ( avnµs(−s+µ)e(−b+s)τ −µ − s +(−µ −b+ s) ( − m(av2nµ2b−av2nµbs+av2µws+ vmµb2s − vmb2s2 + vµb2s− vb2s2 +µb2s−b2s2) µ + b2s(−s+µ)vm2 2µ + vµa(−nbµ + s(vmw+nb)) s − av2mnµb(−s+µ) −µ − s ) − avlnµbs(−s+µ)e−τ(b−s) b− s ))) (34) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1921 Y (t) = me−µt + ( a lnebt−bτ−st+sτ µ +b− s − e−µta lne−bτ+sτ µ +b− s ) + ( 1 (µ +b− s)(−s+µ)µsv ( al ( (µ +b− s) ( avn2µ(−s+µ)e(µ+b−2s)t+2sτ µ +b−2s + mnbs(−s+µ)e(µ+s)τ+t(b−s) b− s − n((avn− ls)µ +mbs)(−s+µ)e(µ+b−s)t+sτ µ +b− s − lnµs(−s+µ)e(µ+2b−s)t+sτ µ +2b− s − vmµwse(µ−b)τ+bt b ) + µ(−nµ2 +(mw+ sn)µ +mw(b− s))vse(µ+b−s)t−τ(b−s) µ +b− s + te−τ(b−s)vnµ 3s− te−τ(b−s)vnµ 2s2 )) − ( 1 (µ +b− s)(−s+µ)µsv ( al ( (µ +b− s) ( avn2µ(−s+µ)e2sτ µ +b−2s + mnbs(−s+µ)e(µ+s)τ b− s − n((avn− ls)µ +mbs)(−s+µ)esτ µ +b− s − lnµs(−s+µ)esτ µ +2b− s − vmµwse(µ−b)τ b ) + µ(−nµ2 +(mw+ sn)µ +mw(b− s))vse(µ+b−s)t−τ(b−s) µ +b− s ))) e−µt (35) C(t) = ( n+ −wme−t(−s+µ) −s+µ + wm −s+µ ) e−st + a ln ( webt−bτ+sτ b ) − ( e−µt−bτ+st+sτ −µ+s ) µ +b− s − a ln ( we−bτ+sτ b ) − ( e−bτ+sτ −µ+s ) µ +b− s e−st (36) 8. Numerical Example For the study of our mathematical results, we numerically extract the disease caused by the rugose spiraling whitefly on coconut trees. The system of both ordinary differential equations (ODE) and delay differential equations (DDE) were numerically solved using MATLAB software for a range of parameter values in order to test the accuracy of the homotopy perturbation method results under specified sets of parameters. Addtionally, we contrasted the numerical and analytical results. We have seen that there is an effective level of consensus for all parameter values. B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1922 Table 1: The parameters and their values are given in the following table [25]. Symbol Meaning Values taken for the analysis Range Unit v Tree density 0.0138 0.0138 - 0.0178 m−1 s Death rate for RSWs 0.06 0 - 0.01 day−1 b Replanting rate 0.0005 0 - 0.008 day−1 a Contact rate 0.0002 0 - 0.002 pest−1day−1 w Whitefly birth rate 0.2 0.1-0.3 day−1 µ Mortality rate for trees 0.0001 0 - 0.008 day−1 8.1. Numerical Results and Discussion The Eqs.(31) − (33) denote the approximate analytical expressions for the non linear differen- tial equations by homotopy perturbation method. and Eqs.(34) − (36) represent the approximate analytical expressions for the delay differential equations by homotopy perturbation method. Us- ing MATLAB coding, we conducted a comparative analysis between the numerical simulation and our approximate analytical results for the proposed model. The graph displayed in figures (2a) − (5b) represent the both numerical and approximate analytical results of the system of non-linear ordinary differential equations Eqs.(1) − (4) and figures (6a) − (8b) for a delay differential equa- tions Eqs.(5)−(8) From figure (2a), it is evident that the whitefly death rate s increases when healthy tree density varies from 0.006 to 0.01. From figure (2b), the contact rate a decreases when healthy tree density varies from 0.00015 to 0.002. From figure (3a), it is found that the contact rate a increases the infected tree density varies from 0.006 to 0.01. In figure (3b), the death rate of RSW s decreases when infected tree density from 0.00015 to 0.002. The mortality rate for trees µ decreases the infected tree density varies from this 0.00001 to 0.001 and it is indicated in figure (4). From figure (5a), the death rate for RSWs s decreases the whitefly density from 0.006 to 0.01. From figure (5b), the whitefly birth rate w increases the whitefly density from 0.1 to 2.5. In figure (6a), when the time delay τ = 0.5 and the healthy tree density varies from 0.00015 to 0.002, the contact rate a decreases whereas the whitefly death rate s increases when the healthy tree density varies from 0.006 to 0.01 with the time delay τ = 0.5 and it is shown in figure (6b). In figure (7a), the infected tree density increases when the time delay τ increases from 0.5 to 20. By knowing the range of the infection rate, we can control the spread of the whitefly population. In figure (7b), the mortality rate of trees µ decreases when the infected tree density with the time delay τ = 0.5 ranging from 0.00001 to 0.001 whereas the contact rate a increases when the infected tree density with the time delay τ = 0.5 ranging from 0.00015 to 0.002. and it is shown from (7c), and the death rate of RSWs s decreases when the infected tree density with the time delay τ = 0.5 in the ranging from 0.006 to 0.01 is shown in figure (7d). From figure (8a), the whitefly birth rate w increases when the whitefly density with the time delay τ = 0.5 varies from 0.1 to 2.5. From figure (8b), the death rate of RSW s decreases when the whitefly density with the time delay τ = 0.5 between 0.006 and 0.01. It is found from all of these figures that our analytical results closely match with the numer- ical results. Table 1 contains the parametric values that were used in the analysis. The homotopy perturbation method result and the numerical simulation result of Eq.(1) are find to derive the er- B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1923 ror approximation in Table 2. The error estimation achieved by finding the numerical simulation result and the homotopy perturbation method result of Eq.(5) is displayed with the time delay τ = 0.5 in Table 3. Thus, the results obtained by homotopy perturbation method seems fine. 0 20 40 60 80 100 120 Time (t) (days) 4.98 4.99 5 5.01 5.02 5.03 H ea lt h y t re e d en si ty X (t ) (m -2 ) 10 -3 *** Numerical Simulation Analytical Approximation s = 0.006, 0.007, 0.008, 0.009, 0.01 Increases b = 0.0005 = 0.0001 w = 0.2 a = 0.0002 v = 0.0138 (a) 0 20 40 60 80 100 120 Time (t) (days) 5 5.02 5.04 5.06 5.08 5.1 5.12 5.14 5.16 5.18 H ea lt h y t re e d en si ty X (t ) (m -2 ) 10 -3 *** Numerical Simulation Analytical Approximation Decreases b = 0.0005 = 0.0001 w = 0.2 s = 0.006 v = 0.0138 a = 0.00001, 0.00002, 0.00003, 0.00005, 0.0001 (b) Figure 2: Graphs of the healthy tree density X over time t, determined by numerical simulations and the homotopy perturbation method solution. The numerical simulations are indicated by the curves (—) using Eq.(1) and the analyitcal solution of homotopy perturbation method is plotted by (***) using Eq.(31). 0 20 40 60 80 100 120 Time (t) (days) 7 7.5 8 8.5 9 9.5 In fe ct ed t re e d en si ty Y (t ) (m -2 ) 10 -4 Increases *** Numerical Simulation Analytical Approximation b = 0.0005 = 0.0001 w = 0.2 s = 0.006 v = 0.0138 a = 0.00015, 0.00011, 0.00010, 0.003, 0.002 (a) 0 20 40 60 80 100 120 Time (t) (days) 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 In fe ct ed t re e d en si ty Y (t ) (m -2 ) 10 -4 Numerical Simulation Analytical Approximation s = 0.006, 0.007, 0.008, 0.009, 0.01 *** Decreases b = 0.0005 = 0.0001 w = 0.2 a = 0.0002 v = 0.0138 (b) Figure 3: Graphs of the infected tree density Y over time t, determined by numerical simulations and the homotopy perturbation method solution. The numerical simulations are showed by the curves (—) using Eq.(2) and the analytical solution of homotopy perturbation method is plotted by (***) using Eq.(32). B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1924 0 20 40 60 80 100 120 Time (t) (days) 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 10 -4 In fe ct ed t re e d en si ty Y (t ) (m -2 ) ** * ** ** * * * ** ** **** * ** * ** ** ** Numerical Simulation Analytical Approximation*** Decreases b = 0.0005 s = 0.006 w = 0.2 a = 0.0002 v = 0.0138 = 0.00001, 0.0002, 0.0005, 0.0008, 0.001 Figure 4: Graphs of the infected tree density Y over time t, determined by numerical simulations and the homotopy perturbation method solution. The numerical simulations are indicated by the curves (—) using Eq.(2) and the analyitcal solution of homotopy perturbation method is plotted by (***) using Eq.(32). 0 20 40 60 80 100 120 Time (t) (days) 0.6 0.8 1 1.2 1.4 1.6 1.8 2 W h it ef ly d en si ty C (t ) (m -2 ) Decreases s = 0.006, 0.007, 0.008, 0.009, 0.01 *** Numerical Simulation Analytical Approximation b = 0.0005 = 0.0001 w = 0.2 a = 0.0002 v = 0.0138 (a) 0 20 40 60 80 100 120 Time (t) (days) 0.8 1 1.2 1.4 1.6 1.8 2 W h it ef ly d en si ty C (t ) (m -2 ) Numerical Simulation Analytical Approximation *** * * * * * * * ** * * * w = 0.1,0.5,1.0,2.0,2.5 b = 0.0005 = 0.0001 a = 0.0002 v = 0.0138 s = 0.006 Increases (b) Figure 5: Graphs of the whitefly density C over time t, determined by numerical simulations and the homotopy perturbation method solution. The numerical simulations are showed by the curves (—) using Eq.(3) and the analyitcal solution of homotopy perturbation method is plotted by (***) using Eq.(33). B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1925 0 20 40 60 80 100 120 Time (t) (days) 5 5.02 5.04 5.06 5.08 5.1 5.12 5.14 5.16 5.18 H ea lt h y t re e d en si ty X (t ) (m -2 ) 10 -3 *** Numerical Simulation Analytical Approximation Decreases b = 0.0005 = 0.0001 w = 0.2 s = 0.006 v = 0.0138 = 0.5 a = 0.00001, 0.00002, 0.00003, 0.00005, 0.0001 (a) 0 20 40 60 80 100 120 Time (t) (days) 4.98 4.99 5 5.01 5.02 5.03 H ea lt h y t re e d en si ty X (t ) (m -2 ) 10 -3 *** Numerical Simulation Analytical Approximation Increases b = 0.0005 = 0.0001 w = 0.2 a = 0.0002 v = 0.0138 = 0.5 s = 0.006, 0.007, 0.008, 0.009, 0.01 (b) Figure 6: Graphs of the healthy tree density X over time t with the time delay, determined by numerical simulations and the homotopy perturbation method solution. The numerical simulations are indicated by the curves (—) using Eq.(5) and the analyitcal solution of homotopy perturbation method is plotted by (***) using Eq.(34). B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1926 0 20 40 60 80 100 120 Time (t) (days) 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 9 10 -4 In fe ct ed t re e d en si ty Y (t ) (m -2 ) Numerical Simulation Analytical Approximation*** = 0.5, 5, 10, 15, 20 b = 0.0005 = 0.0001 w = 0.2 a = 0.0002 v = 0.0138 s=0.006 Increases (a) 0 20 40 60 80 100 120 Time (t) (days) 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 10 -4 In fe ct ed t re e d en si ty Y (t ) (m -2 ) ** * ** ** * * * ** ** **** * ** * ** ** ** Decreases Numerical Simulation Analytical Approximation*** b = 0.0005 s = 0.006 w = 0.2 a = 0.0002 v = 0.0138 = 0.5 = 0.00001, 0.0002, 0.0005, 0.0008, 0.001 (b) 0 20 40 60 80 100 120 Time (t) (days) 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 In fe ct ed t re e d en si ty Y (t ) (m -2 ) 10 -3 *** Increases Numerical Simulation Analytical Approximation a = 0.00015, 0.00011, 0.00010, 0.003, 0.002 b = 0.0005 = 0.0001 w = 0.2 s = 0.006 v = 0.0138 = 0.5 (c) 0 20 40 60 80 100 120 Time (t) (days) 7 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 10 -4 In fe ct ed t re e d en si ty Y (t ) (m -2 ) s = 0.006, 0.007, 0.008, 0.009, 0.01 *** Numerical Simulation Analytical Approximation Decreases b = 0.0005 = 0.0001 w = 0.2 a = 0.0002 v = 0.0138 = 0.5 (d) Figure 7: Graphs of the infected tree density Y over time t with the time delay, determined by numerical simulations and the homotopy perturbation method solution. The numerical simulations are showed by the curves (—) using Eq.(6) and the analyitcal solution of homotopy perturbation method is plotted by (***) using Eq.(35). B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1927 0 20 40 60 80 100 120 Time (t) (days) 0.8 1 1.2 1.4 1.6 1.8 2 W h it ef ly d en si ty C (t ) (m -2 ) Numerical Simulation Analytical Approximation *** * * * * * * * ** * * * b = 0.0005 = 0.0001 s = 0.006 a = 0.0002 v = 0.0138 = 0.5 w = 0.1,0.5,1.0,2.0,2.5 Increases (a) 0 20 40 60 80 100 120 Time (t) (days) 0.6 0.8 1 1.2 1.4 1.6 1.8 2 W h it ef ly d en si ty C (t ) (m -2 ) s = 0.006, 0.007, 0.008, 0.009, 0.01 *** Numerical Simulation Analytical Approximation Decreases b = 0.0005 = 0.0001 w = 0.2 a = 0.0002 v = 0.0138 = 0.5 (b) Figure 8: Graphs of the whitefly density C over time t with the time delay, determined by nu- merical simulations and the homotopy perturbation method solution. The numerical simulations are indicated by the curves (—) using Eq.(7) and the analyitcal solution of homotopy perturbation method is plotted by (***) using Eq.(36). 9. Approximate analytical solution of ODE by using homotopy perturbation method The analytical solutions of the system of ODE Eqs.(1) − (3) using Homotopy perturbation method. Linear and non-linear terms in the given equation can be separated by using homotopy perturbation method theory. [1−P] [ dX dt −bX ] +P [ dX dt −bX + 1 v ( bX2 +bXY ) +aXC ] = 0 (37) [1−P] [ dY dt +µY ] +P [ dY dt +µY −aXC ] = 0 (38) [1−P] [ dC dt + sC ] +P [ dC dt + sC−wY ] = 0 (39) And the initial conditions are as follows, X0(0) = l, Y0(0) = m, C0(0) = n. (40) X = X0 + pX1 + p2X2 + ... (41) Y = Y0 + pY1 + p2Y2 + ... (42) C =C0 + pC1 + p2C2 + ... (43) p0 : dX0 dt −bX0 = 0, (44) p1 : dX1 dt −bX1 + bX0 v + bX0Y0 v +aX0C0 = 0, (45) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1928 p2 : dX2 dt −bX2 + bX1 v + bX0Y1 v +bX1Y0 +aX0C1 = 0 (46) p0 : dY0 dt + pY0 = 0, (47) p1 : dY1 dt + pY1 −aX0C0, (48) p2 : dY2 dt − pY2 −aX0C1 −aX1C0 = 0, p0 : dC0 dt + sC0 = 0, (49) p1 : dC1 dt + sC1 −wY0 = 0, (50) p2 : dC2 dt + sC2 −wY1. (51) X0 = lebt , Y0 = me−µt , C0 = ne−st . (52) X1(t) = −l ( lebt − mbe−µt µ − kane−µt s ) k + l ( l − mb µ − kan s ) k ebt , (53) Y1(t) = ( alnet(µ+b−s) µ +b− s − lna µ +b− s ) e−µt , (54) C1(t) = ( −wme−t(−s+µ) −s+µ + wm −s+µ ) e−st . (55) X2(t) = ( 1 v2µs(µ +b− s)(−s+µ) ( l(µ +b− s) ( bs(−s +µ) ( lµe−µt b − vm2be−sµt 2µ )) − 1 µ (( −ms(−s+µ)(vm+1)b3 − (((1+(−l +m)v)s)+ v2an)µ − s2(vm+1))(−s+µ)mb2 −µv(m(van− ls)µ2)−2 ( − lsm+an ( vm− l 2 )) sµ +(−ls2m+an(vm− l)s+ vamw)s ) b− v2amµws(−s+µ)ebt ) + tmb2s3 + tlµb2s2 + tmµb3s+ tmµ 2b2s+ tavnµ 2b2 + tavnµ 3b + tavnµbs2 + valnµbs(−s+µ)e(b−s)t b− s + v2amnµb(−s+µ)(µ +b− s)e−(µ+s)t −µ − s − tmb3s2 B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1929 − tlµbs3 +2tlµ2bs2 − tlµ2b2s− tlµ3bs−2tmµb2s2 − tavnµb2s −2tavnµ 2bs− vlmµbs(−s+µ)(µ +b− s)e−(µ−b)t −µ +b )) − 1 v2µs(µ +b− s)(−s+µ) ( l ( (µ +b− s) ( bs(−s+µ)( lµ b − vm2b 2µ ) − µ(−n(−s+µ)b+ vmws)av s ) − 1 µ ( −ms(−s+µ)(vm+1)b3 − (((1+(−l +m)v)s+ v2an)µ − s2(vm+1))(−s+µ)mb2 −µv ( m(van− ls)µ2 −2 ( − lsm+an ( vm− 1 2 )) sµ +(−ls2m+an(vm− l)s+ vamw)s)b− v2amµws(−s+µ) ) + valnµbs(−s+µ) b− s + v2amnµb(−s+µ)(µ +b− s) −µ − s − vlmµbs(−s+µ)(µ +b− s) −µ +b ))) ebt , (56) Y2(t) = ( 1 (µ +b− s)vµs ( lna ( ((van− ls)µ +mbs)e(µ+b−s)t + lµs(va+µ +b− s)eµ+2b−st µ +2b− s − s ( mb(µ +b− s)e(b−s)t b− s + valµebt b ) − vnµa(µ +b− s)e(µ+2b−s)t µ +b−2s )) − lna ( (van− ls)µ +mbs+ lµs(va+µ+b−s) µ+2b−s − s ( mb(µ+b−s) b−s + valµ b ) − vnµa(µ+b−s) µ+b−2s ) (µ +b− s)vµs ) e−µt , (57) C2(t) = −w2m ( −t + e−t(−s+µ) −s+µ ) µ − s + w2m (µ − s)(s−µ) e−st . (58) X(t) = X0 +X1 +X2 (59) Y (t) = Y0 +Y1 +Y2 (60) C(t) =C0 +C1 +C2 (61) X(t) = lim p→1 X(t) = X0 +X1 +X2 (62) Y (t) = lim p→1 Y (t) = Y0 +Y1 +Y2 (63) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1930 C(t) = lim p→1 C(t) =C0 +C1 +C2 (64) Final solution is given in the text as Eqs. (31) − (33). 10. Approximate analytical solution of DDE by using homotopy perturbation method The analytical solutions of the system of DDE Eqs.(5) − (7) using Homotopy perturbation method. Linear and non-linear terms in the given equation can be separated by using homotopy perturbation method theory with the initial conditions given in Eq.(8). The time delay is introduced in Healthy trees and Whitefly population in Eq.(6). [1−P] [ dY dt +µY ] +P [ dY dt + pY −αX(t − τ)C(t − τ) ] = 0 (65) with the initial conditions given in Eq.(8) X(t − τ) = X0(t − τ)+ pX1(t − τ)+ p2X2(t − τ)+ ... (66) C(t − τ) =C0(t − τ)+ pC1(t − τ)+ p2C2(t − τ)+ ... (67) p0 : dX0(t − τ) dt −bX0(t − τ) = 0, (68) p1 : dX1(t − τ) dt −bX1(t − τ)+ bX0(t − τ) v + bX0(t − τ)Y0 v +aX0(t − τ)C0(t − τ) = 0, (69) p2 : dX2(t − τ) dt −bX2(t − τ)+ bX1(t − τ) k + bX0(t − τ)Y1 v +bX1(t − τ)Y0 +aX0(t − τ)C1(t − τ) = 0 (70) p0 : dY0 dt +µY0 = 0, (71) p1 : dY1 dt + pY1 −αX0(t − τ)C0(t − τ) = 0, (72) p2 : dY2 dt − pY2 −αX0(t − τ)C1(t − τ)−αX1(t − τ)C0(t − τ) = 0, p0 : dC0 dt + sC0 = 0, (73) p1 : dC1(t − τ) dt + sC1(t − τ)−wY0 = 0, (74) p2 : dC2(t − τ) dt + sC2(t − τ)−wY1 = 0. (75) B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1931 X0 = lebt , Y0 = me−µt , C0 = ne−st . (76) X1(t) = ( − l(lebt − mbe−µt µ − vane−µt s ) v + l(l − mb µ − van s ) v ) ebt (77) Y1(t) = ( a lnebt−bτ−st+sτ µ +b− s − e−µta lne−bτ+sτ µ +b− s ) (78) C1(t) = ( −wme−t(−s+µ) −s+µ + wm −s+µ ) e−st (79) X2(t) = ( a lnebt−bτ−st+sτ µ +b− s − e−µta lne−bτ+sτ µ +b− s ) + ( 1 (µ +b− s)(−s+µ)µsv ( al ( (µ +b − s) ( avn2µ(−s+µ)e(µ+b−2s)t+2sτ µ +b−2s + mnbs(−s+µ)e(µ+s)τ+t(b−s) b− s − n((avn− ls)µ +mbs)(−s+µ)e(µ+b−s)t+sτ µ +b− s − lnµs(−s+µ)e(µ+2b−s)t+sτ µ +2b− s − vmµwse(µ−b)τ+bt b ) + µ(−nµ2 +(mw+ sn)µ +mw(b− s))vse(µ+b−s)t−τ(b−s) µ +b− s + te−τ(b−s)vnµ 3s− te−τ(b−s)vnµ 2s2 )) − ( 1 (µ +b− s)(−s+µ)µsv ( al ( (µ +b− s) ( avn2µ(−s+µ)e2sτ µ +b−2s + mnbs(−s+µ)e(µ+s)τ b− s − n((avn− ls)µ +mbs)(−s+µ)esτ µ +b− s − lnµs(−s+µ)esτ µ +2b− s − vmµwse(µ−b)τ b ) + µ(−nµ2 +(mw+ sn)µ +mw(b− s))vse(µ+b−s)t−τ(b−s) µ +b− s ))) e−µt (80) Y2(t) = ( 1 (µ +b− s)(−s+µ)µsv ( al( ( µ +b− s) ( avn2µ(−s+µ)e(µ+b−2s)t+2sτ µ +b−2s + mnbs(−s+µ)e(µ+s)τ+t(b−s) b− s − n((avn− ls)µ +mbs)(−s+µ)e(µ+b−s)t+sτ µ +b− s B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1932 − lnµs(−s+µ)e(µ+2b−s)t+sτ µ +2b− s − vmµwse(µ−b)τ+bt b ) + µ(−nµ2 +(mw+ sn)µ +mw(b− s))vse(µ+b−s)t−τ(b−s) µ +b− s + te−τ(b−s)vnµ 3s− te−τ(b−s)vnµ 2s2 )) − ( 1 (µ +b− s)(−s+µ)µsv ( al ( (µ +b− s) ( avn2µ(−s+µ)e2sτ µ +b−2s + mnbs(−s+µ)e(µ+s)τ b− s − n((avn− ls)µ +mbs)(−s+µ)esτ µ +b− s − lnµs(−s+µ)esτ µ +2b− s − vmµwse(µ−b)τ b ) + µ(−nµ2 +(mw+ sn)µ +mw(b− s))vse(µ+b−s)t−τ(b−s) µ +b− s ))) e−µt (81) C2(t) = a ln ( webt−bτ+sτ b ) − ( e−µt−bτ+st+sτ −µ+s ) µ +b− s − a ln ( we−bτ+sτ b ) − ( e−bτ+sτ −µ+s ) µ +b− s e−st (82) Following same as Eqs.(59) − (64). Final solution is given in the text as Eqs.(34) − (36). B. Dhivyadharshini, R. Senthamarai / Eur. J. Pure Appl. Math, 17 (3) (2024), 1908-1936 1933 Ta bl e 2: O bs er va tio n of nu m er ic al re su lts an d th e ho m ot op y pe rt ur ba tio n m et ho d so lu tio n of E q. (1 ) fo r va ri ou s va lu es of a as in Fi gu re .( 2b ) a = 0. 00 00 5 a = 0. 00 08 a = 0. 00 1 t N um er ic al H PM E rr or % N um er ic al H PM E rr or % N um er ic al H PM E rr or % 0 0. 00 50 0 0. 00 50 0 0. 01 79 3 0. 00 50 0 0. 00 50 0 0. 01 79 6 0. 00 50 0 0. 00 50 0 0. 01 79 8 20 0. 00 50 2 0. 00 50 0 0. 82 84 7 0. 00 50 1 0. 00 50 0 0. 01 79 6 0. 00 50 1 0. 00 50 0 0. 22 69 9 40 0. 00 50 4 0. 00 50 0 1. 25 69 9 0. 00 50 3 0. 00 50 0 0. 61 64 0 0. 00 50 2 0. 00 50 0 0. 47 47 9 60 0. 00 50 6 0. 00 50 0 1. 69 80 9 0. 00 50 5 0. 00 50 0 0. 95 73 4 0. 00 50 4 0. 00 5 0. 75 71 80 0. 00 50 9 0. 00 50 0 2. 14 96 4 0. 00 50 7 0. 00 50 0 1. 32 14 8 0. 00 50 5 0. 00 50 0 1. 06 96 8 10 0 0. 00 51 1 0. 00 50 0 2. 14 96 4 0. 00 50 9 0. 00 50 0 1. 70 54 4 0. 00 50 7 0. 00 50 0 1. 40 83 0 A ve ra ge er ro r% 1. 06 09 7 0. 82 01 2 0. 65 91 4 Ta bl e 3: O bs er va tio n of nu m er ic al re su lts an d th e ho m ot op y pe rt ur ba tio n m et ho d so lu tio n of E q. (5 ) fo r va ri ou s va lu es of a w ith th e tim e de la y τ = 0. 5 as in Fi gu re .( 6a ) a = 0. 00 00 5 a = 0. 00 08 a = 0. 00 1 t N um er ic al H PM E rr or % N um er ic al H PM E rr or % N um er ic al H PM E rr or % 0 0. 00 50 0 0. 00 50 0 0. 00 01 0 0. 00 50 1 0. 00 50 1 0. 00 31 6 0. 00 07 3 0. 00 50 0 0. 61 42 20 0. 00 50 2 0. 00 50 0 0. 39 66 0 0. 00 50 1 0. 00 50 1 0. 00 31 6 0. 00 07 3 0. 00 50 0 0. 58 05 0 40 0. 00 50 4 0. 00 50 0 0. 00 50 3 0. 01 21 3 0. 00 50 3 0. 00 07 7 0. 00 07 7 0. 00 05 0 0. 53 14 6 60 0. 00 50 6 0. 00 50 0 1. 23 88 0 0. 00 50 5 0. 00 50 5 0. 02 66 2 0. 00 07 9 0. 00 00 5 0. 53 14 6 80 0. 00 50 9 0. 00 50 0 1. 67 97 0 0. 00 50 7 0. 00 50 6 0. 04 63 4 0. 00 08 1 0. 00 05 0 0. 51 35 0 10 0 0. 00 51 1 0. 00 50 0 0. 13 11 0 0. 00 50 9 0. 00 50 8 0. 07 09 8 0. 00 08 1 0. 00 50 0 0. 49 86 9 A ve ra ge er ro r% 1. 04 27 0 0. 02 65 4 0. 02 34 1 REFERENCES 1934 11. Conclusion This paper examines how the behaviour of interaction among species populations and param- eters affect the time delay system. We concentrated on how RSWs and coconut palms interacted in pollachi with the time delay. Similarly, wherever we need a pest control we can extend our research in the same way. The equilibrium points and the conditions to be LAS and sensitiv- ity parameters have been studied. We have described the approximate analytical and numerical solution of ordinary differential equations and delay differential equations using He’s homotopy perturbation method. The homotopy perturbation method is extensibly enhanced to derive explicit solutions of the ODEs and DDEs using symbolic computation software such as Matlab. We dis- covered strong agreement between the analytical results and the numerical simulation results. The study cited above offers convincing evidence that one important influencing element is the time delay τ . Therefore, by finding the range of days of trees getting infected with the help of time delay, we can introduce control measures like pesticides, insecticides and by introducing a preda- tor in the particular range of tree becoming infected so that we can enormously assist farmers in controlling the disease. Even homotopy perturbation method does not has a controlling parameter h̃ it has given as the convergent series solution for our model since it has less non-linearity. Acknowledgements We sincerely thank the reviewers for their valuable comments, which were of great help in revising the manuscript. The authors are pleased to acknowledge the financial support of the Se- lective Excellence Research Initiative (SRMIST/R/AR(A)/SERI2023/174/31). It is our pleasure to thank the College of Engineering and Technology, SRM IST for its valuable support and constant encouragement. References [1] M. Abolhasani, H. Ghaneai, and M. Heydari. Modified homotopy perturbation method for solving delay differential equations. Appl. Sci. Rep., 16(2):89–92, 2016. [2] L. J. Allen, F. Brauer, P. Van den Driessche, and J. Wu. Mathematical Epidemiology. Springer, Berlin, 2019. [3] Naveed Anjum and Ji-Huan He. Homotopy perturbation method for n/mems oscillators. Mathematical Methods in the Applied Sciences, pages 1–15, 2020. [4] Fahad Al Basir, Ezio Venturino, Santanu Ray, and Priti Kumar Roy. Impact of farming awareness and delay on the dynamics of mosaic disease in jatropha curcas plantations. Com- putational and Applied Mathematics, 37:6108–6131, 2018. [5] M. Bodnar. The nonnegativity of solutions of delay differential equations. Appl Math Lett, 13(6):91–95, 2000. REFERENCES 1935 [6] P. Chowdappa. Invasive rugose spiralling whitefly on coconut. ICAR-CPCRI Technical Bulletin, 117, 2017. [7] R.V. Culshaw and S. Ruan. A delay-differential equation model of hiv infection of cd4+ t-cells. Mathematical Biosciences, 165:27–39, 2000. [8] B. Dhivyadharshini and R. Senthamarai. Mathematical analysis of a non linear prey predator system: Analytical approach by hpm. In AIP Conference Proceedings, volume 2516, 2022. [9] B. Dhivyadharshini and R. Senthamarai. A non-linear prey-predator model: Detailed analy- sis by hpm and ham. Research Trends in Mathematics and Statistics, 25:31–57, 2023. [10] O. Diekmann, Heesterbeek JAP, and Metz JAJ. On the definition and the computation of the basic reproduction ratio r0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology, 28(4):365–382, 1990. [11] K. Elango, S. J. Nelson, and A. Aravind. Rugose spiralling whitefly, aleurodicus rugiop- erculatus martin (hemiptera, aleyrodidae): An invasive foes of coconut. J. Entomol. Res, 44:261–266, 2020. [12] K. Elango, S. J. Nelson, S. Sridharan, V. Paranidharan, and S. Balakrishnan. Biology, dis- tribution and host range of new invasive pest of india coconut rugose spiralling whitefly aleurodicus rugioperculatus martin in tamil nadu and the status of its natural enemies. Int. J. Agricul. Sci, 11:8423–8426, 2019. [13] J. Hale. Theory of Functional Differential Equations. Springer, Heidelberg, 1977. [14] Ji-Huan He. On the frequency-amplitude formulation for nonlinear oscillators with general initial conditions. Int. J. Appl. Comput. Math, 111(7), 2021. [15] Ji-Huan He, G. M. Moatimid, and D. R. Mostapha. Nonlinear instability of two streaming- superposed magnetic reiner-rivlin fluids by he-laplace method. Journal of Electroanalytical Chemistry, 895, 2021. [16] Holt, M. J. Jeger, J. M. Thresh, and G. W. Otim-Nape. An epidemiological model incorpo- rating vector population dynamics applied to african cassava mosaic virus disease. J. Appl. Ecol, 34:793–806, 1997. [17] Y. Kuang. Delay Differential Equations with Applications in Population Dynamics. Aca- demic Press, Inc., New York, 1993. [18] S. Pathak and A. Maiti. Pest control using virus as control agent: A mathematical model. Nonlinear Analysis: Modelling and Control, 17:67–90, 2012. [19] S. Ray and F. A. Basir. Impact of incubation delay in plant-vector interaction. Math. Comput. Simul, 170:16–31, 2020. REFERENCES 1936 [20] F. A. Rihan, D. H. Abdel Rahman, S. Lakshmanan, and A. S. Alkhajeh. A time delay model of tumor-immune system interactions: Global dynamics, parameter estimation, sensitivity analysis. Applied Mathematics and Computation, 232:606–623, 2014. [21] F. A. Rihan, C. Tunc, S.H. Saker, S. Lakshmanan, and R. Rakkiyappan. Application of delay differential equations in biological systems. Complexity, 2018. [22] S. Shanas, J. Job, T. Joseph, and G. Anju Krishnan. First report of the invasive rugose spiraling whitefly, aleurodicus rugioperculatus martin (hemiptera: Aleyrodidae) from the old world. Entomon, 41:365–368, 2016. [23] M. Sivakumar, M. Mallikarjuna, and R. Senthamarai. A kinetic non-steady state analysis of immobilized enzyme systems with external mass transfer resistance. AIMS Mathematics, 9(7):18083–18102, 2024. [24] M. Sivakumar and R. Senthamarai. Mathematical model of epidemics: Analytical approach to sirw model using homotopy perturbation method. In AIP Conference Proceedings, volume 2277, 2020. [25] G. Suganya and R. Senthamarai. Mathematical modeling and analysis of the effect of the rugose spiraling whitefly on coconut trees. AIMS Mathematics, 7:13053–13073, 2022. [26] R. Sundararaj and K. Selvaraj. Invasion of rugose spiraling whitefly, aleurodicus rugioper- culatus martin (hemiptera: Aleyrodidae): A potential threat to coconut in india. Phytopara- sitica, 45:71–74, 2017. [27] Nural Atiqah Talib, Normah Maan, and Aminu Barde. Analytical approximation solution for logistic delay differential. Malaysian Journal of Fundamental and Applied Sciences, 16(3):368–373, 2020. [28] T. Vijayalakshmi and R. Senthamarai. An analytical approach to the density dependent prey- predator system with reddington-deangelies functional response. In AIP Conference Pro- ceedings, volume 2112, 2019. [29] S. Wang, Z. Ma, X. Li, and T. Qi. A generalized delay-induced sirs epidemic model with relapse. AIMS Math, 7:6600–6618, 2022. [30] X. Yang, L. Chen, and J. Chen. Permanence and positive periodic solution for the single species nonautonomus delay diffusive model. Comput Math Appl, 32:109–116, 1996.