EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 4, 2024, 3826-3846 ISSN 1307-5543 – ejpam.com Published by New York Business Global Employing the Sadik Residual Power Series Method to Analyze a System of Nonlinear Caputo Time-Fractional Partial Differential Equations Prapart Pue-on Mathematics and Applied Mathematics Research Unit, Department of Mathematics, Faculty of Science, Mahasarakham University, Maha Sarakham, 44150,Thailand Abstract. The present study proposes an approximate analytical solution to a nonlinear system of time-fractional partial differential equations. The Sadik residual power series method, which integrates the two-part Sadik integral transform with the residual power series technique, is em- ployed to solve the fractional differential equation in the Caputo sense. Nonlinear problems with known and unknown solutions are examined to demonstrate the capacity of the technique. Nu- merical simulations and 3D visualizations are conducted for various values of the fractional order to further understand the solution’s behavior. Additionally, the results are validated against exact solutions or existing methodologies to ensure their reliability and accuracy. A key advantage of the proposed method is its ability to generate results without the need for Adomian polynomials, perturbation techniques, discretization, or linearization, enabling a more efficient. 2020 Mathematics Subject Classifications: 26A33, 35R11, 35A22 Key Words and Phrases: Residual Power Series, Sadik transform, Nonlinear System, Caputo fractional derivatives, Approximate solution, method 1. Introduction A fractional differential equation is a generalization of an integer-order differential equation that appears in various scientific frameworks [18]. The study of this sort of equation has gained increasing interest among scholars because it can describe real-world physical phenomena in greater detail and accuracy than integer differential equations offer, especially in the physical processes involving memory and heritability properties. Ana- lyzing non-linear fractional differential equations has become an essential challenge for a century. Scientists have attempted to devise and enhance novel approaches, whether ana- lytical or numerical, for conquering such circumstances. Among these intriguing methods are the Adomian decomposition method [15], [24], [30], homotopy-perturbation method [1], the homotopy analysis method [16], the differential transform method [12], the Chebyshev DOI: https://doi.org/10.29020/nybg.ejpam.v17i4.5494 Email address: prapart.p@msu.ac.th (P. Pue-on) https://www.ejpam.com 3826 Copyright: © 2024 The Author(s). (CC BY-NC 4.0) P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3827 Wavelets method [19]. Although the procedures described above are reliable and verifi- able, they come with practical constraints and require a long time to generate results. As a result, establishing a strategy for effectively addressing issues, reducing determining time, and making it advantageous is challenging for researchers. Sadik integral transformation is one of the most vital tools used to address the issues dealing with fractional calculus in science and engineering [25], [28], [26]. The transform is considered to generalize other important integral transformations since this transformation can be reduced to others by assigning the parameter. Due to their straightforward use, the Sadik transformation has continued to gain popularity among academics since it was first introduced in 2018 [27]. The advancement of the transform in solving fractional differential equations is found in much literature and references therein, such as Alabdala et al. [2] used the Sadik transform to solve the Caputo delayed fractional differential equation; Derakhshan [9] employed the Sadik transform to obtain the solution for the equal width fractional differential equation; Pue-on [20] utilized the Sadik decomposition method to generate the solution for a system of nonlinear fractional Volterra Integro differential equations of convolution type; Pue-on [21], [22] also obtained the solution of the space and time fractional telegraph equation and well-known fractional partial differential equation by double Sadik transform. The residual power series approach is a highly effective semi-analytic procedure for handling scientific challenges. The aforementioned strategy is based on the generalized Taylor series expansion. The concept is to view the solution as a fractional power series with unknown coefficients and establish a residual function. A traditional residual power series method and an improved version have been widely utilized in literature to address a system of fractional differential equations, both linear and non-linear. Recently, some of these developments in non-linear system can be observed in [3], [5], [13], [14] and [23]. The inspiration for this work was to devise an analytical method that leveraged two widely recognized mathematical tools, the Sadik integral transform and the residual power series technique, to solve a nonlinear system of fractional order differential equations. The primary objective of the current study is to deploy the Sadik residual power series approach to solve the initial value problem for a system of nonlinear Caputo fractional partial differential equations with the form Dγ t ui(x, t) = Ni(uj(x, t),x, t) + fi(x, t), i = 1, . . . , n (1) ui(x, 0) = gi(x). (2) 0 < γ ⩽ 1. Here Dγ t (·) denotes the Caputo partial fractional derivative of order γ, Ni represents the general nonlinear differential operator and fi is the source term. This system is highly significant in the fields of science and technology. The structure of the manuscript is designed as follows: Section 2 presents fundamen- tal definitions and theorems of fractional calculus, the Sadik transform, and the residual power series method. These theoretical foundations support the methodology employed in this study. Section 3 describes the implementation of the Sadik residual power se- ries method to nonlinear systems of fractional partial differential equations. This section P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3828 also comprises illustrative examples that serve to demonstrate and validate the proposed technique. Section 4 shows a comprehensive discussion of the findings and presents the conclusions from this research. 2. Basic Definitions and Theorems of Fractional Calculus, Sadik Integral Transform and Residual Power Series Method This section reviews the basic notions and essential properties of fractional calculus. Furthermore, the definitions and theorems of the Sadik integral transform and residual power series technique are discussed. 2.1. Fractional Calculus Definition 1. [18] The fractional integral operator of order γ(γ ⩾ 0) of ϕ(t) of Riemann- Liouville type is defined as Iγϕ(t) =  1 Γ(γ) ∫ t 0 ϕ(τ) (t− τ)1−γ dτ, γ > 0, t > 0, ϕ(t), γ = 0 Definition 2. The Caputo fractional derivative operator Dγ of order γ, (n − 1 < γ ⩽ n, n ∈ N) is defined in the following form, Dγϕ(t) = 1 Γ(n− γ) ∫ t 0 (t− τ)−γ+n−1ϕ(n)(τ)dτ, (3) γ > 0, t > 0, where the function ϕ(t) has absolutely continuous derivatives up to order n− 1. Moreover, the following basic properties can be proved 1.) DγC = 0, C is a constant 2.) For p ∈ N0 = {0, 1, 2, . . .} and the ceiling function ⌈γ⌉ refers to the smallest integer greater than or equal to γ Dγxp =  0, for p < ⌈γ⌉ Γ(p+ 1) Γ(p+ 1− γ) xp−γ , for p ⩾ ⌈γ⌉ 3.) DγIγϕ(t) = ϕ(t), 4.) IγDγϕ(t) = ϕ(t)− n−1∑ k=1 ϕ(k)(0) k! tk. P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3829 Definition 3. [7] The fractional derivative of order γ > 0 in the Caputo sense is stated as Dγ t u(x, t) = ∂γu(x, t) ∂tγ =  1 Γ(n−γ) ∫ t 0 (t− τ)n−γ−1 ∂ nu(x,τ) ∂τn dτ, n− 1 < γ < n ∂nu(x,t) ∂tn , γ = n ∈ N. Definition 4. [17] The Mittag-Leffler function Ep(z) with p > 0 is defined as Ep(z) = ∞∑ n=0 zn Γ(pn+ 1) . 2.2. The Sadik Transform Definition 5. [27] If f(t) is piecewise continuous function on the interval 0 ⩽ t ⩽ A for any A > 0 and |f(t)| ⩽ KeBt when t ⩾ M, for any real constant B and some positive constant K and M . Then Sadik transform of f(t) is defined by S [f(t)] = 1 vβ ∫ ∞ 0 f(t)e−tvαdt = F (vα, β) (4) where v is complex variable, α is any non zero real number, and β is any real number. Here S is called the Sadik transform operator. Remark 1. The Sadik transform can be evolved to the other intrgral transform by assign- ing the parameter α and β as follows : α = 1, β = 0 (Laplace Transform ), α = 1, β = 1 (Aboodh Transform ), α = 1, β = −1 (Laplace–Carson transform ), α = −1, β = 0 ( Kamal TransformX, α = −1, β = 1 (Sumudu Transform), α = −1, β = −1 (Elzaki Transform), α = −2, β = 1 (Tarig Transform). Remark 2. The basic property of Sadik transform and its transform of elemetary functions are shown in [27]. Theorem 1. [25] Let F (vα, β) denote the Sadik transform of f(t) and f ′(t), f ′′(t), f ′′′(t), . . . , f (n−1)(t) are continuous on [0,∞.) Then the Sadik transform of f (n)(t) is S [ f (n)(t) ] = vnαF (vα, β)− n−1∑ k=0 vkα−βf (n−1−k)(0). (5) Theorem 2. [25] Let n − 1 < γ < n, (n = [γ] + 1) and f(t), f ′(t), . . . , f (n−1)(t) are continuous on [0,∞) and of exponential order, while Dγf(t) is piecewise continuous on [0,∞). Then Sadik transform of Caputo fractional derivative of order γ of function f is given by S[Dγf(t)] = vγαF (vα, β)− n−1∑ k=0 v(γ−n+k)α−βf (n−1−k)(0+). P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3830 2.3. The Residual Power Series Method Definition 6. [10, 11] If n ∈ N where n − 1 < γ ⩽ n. A power series expansion of the form ∞∑ n=0 an(t− t0) nγ = a0 + a1(t− t0) γ + a2(t− t0) 2γ + . . . , t ⩾ t0 is called fractional power series about t0. Theorem 3. [10, 11] Suppose that u(t) has a fractional power series representation at t = t0 of the form u(t) = ∞∑ n=0 an(t− t0) nγ , 0 ⩽ m− 1 < γ ⩽ m, t0 ⩽ t < t0 +R. If u(t) ∈ C[t0, t0 +R) and Dnγu(t) ∈ C(t0, t0 +R), n = 0, 1, 2, . . . , then an = Dnγu(t) Γ(nγ+1) . Definition 7. [10, 11] A power series of the form ∞∑ n=0 fn(x)(t− t0) nγ = f0(x) + f1(x)(t− t0) γ + f2(x)(t− t0) 2γ + . . . is called multiple fractional power series about t = t0, where t is a variable and fm’s are functions of x called the coefficients of the series. Theorem 4. [10, 11] Suppose that u(x, t) has a multiple fractional power series represen- tation at t = t0 of the form u(x, t) = ∞∑ n=0 un(x, t) = ∞∑ n=0 fn(x)(t− t0) nγ , 0 ⩽ m − 1 < γ ⩽ m, x ∈ I, t0 ⩽ t < t0 + R. If Dnγ t u(x, t), n = 0, 1, 2, . . . are continuous on I × (t0, t0 +R), then fn(x) = Dnγ t u(x, t) Γ(nγ + 1) , n = 0, 1, 2, . . . . Here Dnγ t (·) = ∂nγ ∂tnγ (·) = ∂γ ∂tγ ( ∂γ ∂tγ (. . . ( ∂γ ∂tγ (·)))) (n − times), and R = minc∈I Rc, in which Rc is the radius of convergence of the fractional power series ∑∞ k=0 fn(c)(t− t0) nγ . According to the convergence of the classic residual power series method, there is a real number 0 < λ < 1, such that ||un+1(x, t)|| ⩽ λ||un(x, t)||, t ∈ (t0, t0 +R). 3. Main Results This section explains the Sadik residual power series method and presents revealed examples that support its suggested concept. P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3831 3.1. Implementation of the Sadik residual power series method in a sys- tem of time-fractional PDEs Consider the nonlinear system of time-fractional partial differential equation (1) Dγ t ui(x, t) = Ni(uj(x, t),x, t) + fi(x, t), i = 1, . . . , n with initial conditions ui(x, 0) = gi(x). Applying the Sadik transform and along with the initial condition yields vγαUi(x, v α, β)− v(γ−1)α−βgi(x) = S [ Ni(uj(x, t),x, t) ] + Fi(x, v α, β). Rearranging the equation obtained Ui(x, v α, β) = 1 vα+β gi(x) + 1 vγα S [ Ni(uj(x, t),x, t) ] + 1 vγα Fi(x, v α, β). Performing the inverse Sadik transform results in ui(x, t) = gi(x) + S−1 [ 1 vγα S [ Ni(uj(x, t),x, t) ]] + S−1 [ 1 vγα Fi(x, v α, β) ] . Assume the solution to (1) is provided within the infinite series ui(x, t) = ∞∑ n=0 ai,n(x) tnγ Γ(nγ + 1) . (6) According to the residual power series strategy, we define the residual function as follows Resi(x, t) = ui(x, t)− gi(x)− S−1 [ 1 vγα S [ Ni(uj(x, t),x, t) ]] − S−1 [ 1 vγα Fi(x, v α, β) ] , and the kth truncated series of (6) is stated as uk(x, t) = k∑ n=0 ai,n(x) tnγ Γ(nγ + 1) . Since the solution satisfies the initial condition, ai,0(x) = gi(x). Hence, the k- order ap- proximated solution turns into ui,k(x, t) = gi(x) + k∑ n=1 ai,n(x) tnγ Γ(nγ + 1) . (7) Furthermore, the k- residual function is expressed as Resi,k(x, t) = ui,k(x, t)− gi(x)− S−1 [ 1 vγα S [ Ni(uj,k(x, t),x, t) ]] − S−1 [ 1 vγα Fi(x, v α, β) ] . P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3832 For the calculation of the coefficients ai,n(x), n = 1, 2, 3, . . . in equation (7), construct the relation for recurrence ui,0(x, t) = gi(x) Resi,k(x, t) = ui,k(x, t)− gi(x)− S−1 [ 1 vγα S [ Ni(uj,k(x, t),x, t) ]] − S−1 [ 1 vγα Fi(x, v α, β) ] (8) and and afterward follow the condition [29] t−kα · Resi,k(x, t) ∣∣∣ t=0 = 0. (9) Convergence criteria for Sadik residual power series method Theorem 5. [8] Let the Banach space B ≡ C(X × [0, T ]) is defined on X × [0, T ], where X = [a1, b1]× [a2, b2]× . . .× [an, bn], then Eq. (6) as u(x, t) = ∑k n=0 uk(x, t) is convergent series, if u0 ∈ B is bounded and ∀uk ∈ B, ∥uk+1∥ ⩽ λ∥uk∥, where 0 < λ < 1 and ∥u∥ = supx∈X,t∈[0,T ] |u(x, t)|. Theorem 6. [8] The maximum absolute error for (6) is defined by ∥u(x, t)− k∑ n=0 un(x, t)∥ ⩽ λk+1 1− λ ∥u0(x, t)∥. Remark 3. The procedure for proving the aforementioned theorems is the same as that used in the reference, which will not be examined further here. 3.2. Illustrative Examples Example 1. [4] Consider the nonlinear system of fractional partial differential equation Dγ t u(x, t) + w2(x, t)u2x(x, t) + u(x, t) = 1, Dγ t w(x, t) + u2(x, t)w2 x(x, t)− w(x, t) = 1, with initial conditions u(x, 0) = ex, w(x, 0) = e−x. The exact solution to the problem when γ = 1 is u(x, t) = ex−t, w(x, t) = e−x+t. Here N1(u,w) = −w2u2x − u,N2(u,w) = −u2w2 x + w, f1(x, t) = f2(x, t) = 1, g1(x) = ex, g2(x) = e−x. The kth− truncated term series solution for this problem is uk(x, t) = k∑ n=0 an(x) tnγ Γ(nγ + 1) , P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3833 wk(x, t) = k∑ n=0 bn(x) tnγ Γ(nγ + 1) . The recurrence relation to the above issue is u0(x, t) = ex, w0(x, t) = e−x, Res1,k(x, t) = uk(x, t)− ex − S−1 [ 1 vγα S [ − w2 k−1(uk−1) 2 x − uk−1 ]] − tγ Γ(γ + 1) , Res2,k(x, t) = wk(x, t)− e−x + S−1 [ 1 vγα S [ (uk−1) 2(wk−1) 2 x − wk−1 ]] − tγ Γ(γ + 1) , k ⩾ 1. For k = 1. The first-order approximate solution is assumed u1(x, t) = ex + a1(x) tγ Γ(γ + 1) , w1(x, t) = e−x + b1(x) tγ Γ(γ + 1) . Substituting into the above iteration yields Res1,1(x, t) = u1(x, t)− ex − S−1 [ 1 vγα S [ − w2 0(u0) 2 x − u0 ]] − tγ Γ(γ + 1) , = a1(x) tγ Γ(γ + 1) + ex tγ Γ(γ + 1) , Res2,1(x, t) = w1(x, t)− e−x + S−1 [ 1 vγα S [ (u0) 2(w0) 2 x − w0 ]] − tγ Γ(γ + 1) , = b1(x) tγ Γ(γ + 1) − ex tγ Γ(γ + 1) . Figure 1: Subplot (a) displays the approximate solution u5(x, t), and subplot (b) presents the approximate solution w5(x, t) for Example 1 with different γ values at x = 0.5. P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3834 Figure 2: The graphics represent the approximate solution u5(x, t) with different γ values: (a) γ = 0.6, (b) γ = 0.8, (c) γ = 1, and (d) the exact solution for Example 1. Figure 3: The graphics represent the approximate solution w5(x, t) with different γ values: (a) γ = 0.6, (b) γ = 0.8, (c) γ = 1, and (d) the exact solution for Example 1. P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3835 Using the condition (9) with k = 1, t−γRes1,1(x, t) ∣∣ t=0 = 0, t−γRes2,1(x, t) ∣∣ t=0 = 0, produces the coefficients that are a1(x) = −e−x and b1(x) = ex. As a result, the first- order truncated series solution is u1(x, t) = ex − ex tγ Γ(γ + 1) , w1(x, t) = e−x + e−x tγ Γ(γ + 1) . Table 1: Comparison of the approximate solution u5(x, t) for γ = 1 and error analysis in Example 1 using different t, x values. t x u5(x, t) Exact Soln. Abs.Error (SRPSM) Abs.Error (LRPSM) [4] 0.2 1.105170916 1.105170918 2× 10−9 2× 10−9 0.4 1.349858806 1.349858808 2× 10−9 2× 10−9 0.1 0.6 1.648721268 1.648721271 3× 10−9 3× 10−9 0.8 2.013752704 2.013752707 3× 10−9 3× 10−9 1.0 2.459603107 2.459603111 4× 10−9 4× 10−9 0.2 0.9512290245 0.9512294245 4× 10−7 4× 10−7 0.4 1.161833755 1.161834243 4.88× 10−7 4.88× 10−7 0.25 0.6 1.419066952 1.419067549 5.97× 10−7 5.97× 10−7 0.8 1.733252289 1.733253018 7.29× 10−7 7.29× 10−7 1.0 2.116999126 2.117000017 8.90× 10−7 8.90× 10−7 Table 2: Comparison of the approximate solution w5(x, t) for γ = 1 and error analysis in Example 1 using different x, t values. t x w5(x, t) Exact Soln. Abs.Error (SRPSM) Abs.Error (LRPSM) [4] 0.2 0.9048374172 0.9048374180 8× 10−10 8× 10−10 0.4 0.7408182199 0.7408182207 8× 10−10 8× 10−10 0.1 0.6 0.6065306591 0.6065306597 6× 10−10 6× 10−10 0.8 0.4965853033 0.4965853038 5× 10−10 5× 10−10 1.0 0.4065696594 0.4065696597 3× 10−10 3× 10−10 0.2 1.051270808 1.051271096 2.88× 10−7 2.88× 10−7 0.4 0.8607077406 0.8607079764 2.358× 10−7 2.358× 10−7 0.25 0.6 0.7046878967 0.7046880897 1.93× 10−7 1.93× 10−7 0.8 0.5769496523 0.5769498104 1.581× 10−7 1.581× 10−7 1.0 0.4723664234 0.4723665527 1.293× 10−7 1.293× 10−7 For k ⩾ 2. By employing the offered procedure and condition (9), the coefficients of the series solution are ak(x) = (−1)kex, bk(x) = e−x. Hence, the k-truncated series solution is uk(x, t) = k∑ n=0 (−1)nex tnγ Γ(nγ + 1) = ex k∑ n=0 (−tγ)n Γ(nγ + 1) , P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3836 wk(x, t) = k∑ n=0 e−x tnγ Γ(nγ + 1) = e−x k∑ n=0 tnγ Γ(nγ + 1) . It is important to note that when k → ∞, the approximate solution converges on the exact solution u(x, t) = ex ∞∑ n=0 (−tγ)n Γ(nγ + 1) = exE(−tγ), w(x, t) = e−x k∑ n=0 tnγ Γ(nγ + 1) = e−xE(tγ). If γ = 1, the solutions is the exact solution of first-order nonlinear PDEs system u(x, t) = ex−t, w(x, t) = e−x+t. The approximate solution for Example 1 are shown in Figure 1. Figures 2 and 3 compare the visualizations to the approximate and exact solutions. Figure 4: subplot (a) depicts the approximate solution u5(x, y, t) and subplot (b) presents the approximate solution w5(x, y, t) for Example 2 at different γ values for x = y = 0.2. Example 2. [6] Consider the Fractional Reaction-Diffusion Brusselator System Dγ t u(x, y, t) = −2u(x, y, t) + u2(x, y, t)w(x, y, t) + 1 4 [ uxx(x, y, t) + uyy(x, y, t) ] , Dγ t w(x, y, t) = u(x, y, t)− u2(x, y, t)w(x, y, t) + 1 4 [wxx(x, y, t) + wyy(x, y, t)], subject to the initial conditions u(x, y, 0) = e−x−y, w(x, y, 0) = ex+y. The exact solution when γ = 1 is u(x, y, t) = e−x−y− 1 2 t, w(x, y, t) = ex+y+ 1 2 t. P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3837 Figure 5: The graphics represent the approximate solution u5(x, y, t) when t = 1, 0 ⩽ x, y ⩽ 1, with different γ values: (a) γ = 0.6, (b) γ = 0.8, (c) γ = 1, and (d) the exact solution for Example 2. Note that N1(u,w) = −2u + u2w + 1 4(uxx + uyy), N2(u,w) = u − u2w2 x + 1 4(wxx + wyy), f1(x, t) = f2(x, t) = 0, g1(x, y) = e−x−y, g2(x, y) = ex+y. The kth− truncated term series solution to the problem is uk(x, y, t) = k∑ n=0 an(x, y) tnγ Γ(nγ + 1) and wk(x, y, t) = k∑ n=0 bn(x, y) tnγ Γ(nγ + 1) . The recursive relation for this problem is u0(x, y, t) = e−x−y, w0(x, y, t) = ex+y, Res1,k(x, y, t) = uk(x, y, t)− e−x−y + 2S−1 [ 1 vγα S [ uk−1 ]] − S−1 [ 1 vγα S [ u2k−1wk−1 ]] − 1 4 S−1 [ 1 vγα S [ (uk−1)xx + (uk−1)yy ]] , k ⩾ 1 Res2,k(x, y, t) = wk(x, y, t)− ex+y − S−1 [ 1 vγα S [ uk−1 ]] + S−1 [ 1 vγα S [ (uk−1) 2wk−1 ]] − 1 4 S−1 [ 1 vγα S [ (wk−1)xx + (wk−1)yy ]] , k ⩾ 1. Proceed to continue the recurrence to obtain the coefficients, a1(x, y) = −1 2 e−x−y, a2(x, y) = 1 4 e−x−y, . . . , ak(x, y) = 1 (−2)k e−x−y P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3838 Figure 6: The graphics represent the approximate solution w5(x, y, t) when t = 1, 0 ⩽ x, y ⩽ 1, with different γ values: (a) γ = 0.6, (b) γ = 0.8, (c) γ = 1, and (d) the exact solution for Example 2. b1(x, y) = 1 2 ex+y, b2(x, y) = 1 4 ex+y, . . . , bk(x, y) = 1 2k ex+y. Hence, the k-order approximate series solution is uk(x, y, t) = e−x−y [ 1− 1 2 tγ Γ(γ + 1) + 1 4 t2γ Γ(2γ + 1) + . . .+ 1 (−2)k tkγ Γ(kγ + 1) ] wk(x, y, t) = ex+y [ 1 + 1 2 tγ Γ(γ + 1) + 1 4 t2γ Γ(2γ + 1) + . . .+ 1 2k tkγ Γ(kγ + 1) ] . When k → ∞, the approximate solution converges to u(x, y, t) = e−x−y ∞∑ n=0 (− tγ 2 ) n Γ(nγ + 1) = e−x−yE( −tγ 2 ) w(x, y, t) = ex+y ∞∑ n=0 ( t γ 2 ) n Γ(nγ + 1) = ex+yE( tγ 2 ). If γ = 1, the result is reduced to the solution of the first-order equation u(x, y, t) = e−x−y− 1 2 t and w(x, y, t) = ex+y+ 1 2 t. The approximate solution for Example 2 are shown in Figure 4. Figures 5 and 6 compare the visualizations to the approximate and exact solutions. P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3839 Table 3: Comparison of the approximate solution u5(x, y, t) for x = y = 0.2 in Example 2 using different γ values. t γ = 0.2 γ = 0.4 γ = 0.6 γ = 0.8 γ = 1 Exact Soln. 0.2 0.477488627 0.512326854 0.547063183 0.579385607 0.606530658 0.606530659 0.4 0.456444828 0.474696013 0.496237125 0.521635809 0.548811578 0.548811636 0.6 0.443011054 0.449925766 0.460073110 0.475641940 0.496584653 0.496585303 0.8 0.432811651 0.430960053 0.431579433 0.437159809 0.449325358 0.449328964 1.0 0.424420175 0.415283298 0.407895963 0.404110210 0.406556090 0.406569659 Table 4: Comparison of the approximate solution w5(x, y, t)for x = y = 0.2 in Example 2 using different γ values. t γ = 0.2 γ = 0.4 γ = 0.6 γ = 0.8 γ = 1 Exact soln. 0.2 2.417846010 2.074372749 1.865586354 1.734071972 1.648721268 1.648721271 0.4 2.646725351 2.342078614 2.111643502 1.944446048 1.822118665 1.822118800 0.6 2.816085119 2.572130993 2.344840870 2.160190740 2.013751129 2.013752707 0.8 2.956794558 2.785737340 2.577721089 2.387602322 2.225531932 2.225540928 1.0 3.079989436 2.990722585 2.815404246 2.630080531 2.459568271 2.459603111 Table 5: Comparison of the absolute error for u2(x, y, t) when x = y = 0.2 and x = y = 0.5 in Example 2 using different t values. (x,y) t u2(x, y, t) Exact Soln. Abs. Error (SRPSM) Abs. Error (CRPS)[6] 0.2 0.6066396416 0.6065306597 1.089819× 10−4 1.0898× 10−4 (0.2, 0.2) 0.4 0.5496624377 0.5488116361 8.508016× 10−4 8.5080× 10−4 0.6 0.4993884343 0.4965853038 2.8031305× 10−3 2.8031× 10−3 0.8 0.4558176313 0.4493289641 6.4886672× 10−3 6.4886× 10−3 0.2 0.3329308943 0.3328710837 5.98106× 10−5 5.9810× 10−5 (0.5, 0.5) 0.4 0.3016611418 0.3011942119 4.669299× 10−4 4.6692× 10−4 0.6 0.2740701836 0.2725317930 1.5383906× 10−3 1.5383× 10−3 0.8 0.2501580200 0.2465969639 3.5610561× 10−3 3.5610× 10−3 Table 6: Comparison of the absolute error for w2(x, y, t) when x = y = 0.2 and x = y = 0.5 in Example 2 using different t values. (x, y) t w2(x, y, t) Exact Soln. Abs. Error (SRPSM) Abs. Error (CRPS)[6] 0.2 1.648466291 1.648721271 2.54980× 10−4 2.5497× 10−4 (0.2, 0.2) 0.4 1.820026132 1.822118800 2.092668× 10−3 2.0926× 10−3 0.6 2.006504218 2.013752707 7.248489× 10−3 7.2484× 10−3 0.8 2.207900553 2.225540928 1.7640375× 10−2 1.7640× 10−2 0.2 3.003701420 3.004166024 4.64604× 10−4 4.6460× 10−4 (0.5, 0.5) 0.4 3.316303831 3.320116923 3.813092× 10−3 3.8130× 10−3 0.6 3.656089058 3.669296668 1.3207610× 10−2 1.3207× 10−3 0.8 4.023057105 4.055199967 3.2142862× 10−2 3.2142× 10−2 P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3840 Figure 7: subplot (a) depicts the approximate solution u2(x, t) and subplot (b) presents the approximate solution w2(x, t) for Example 3 at different γ values for x = 0.5. Example 3. Consider the coupled system of reaction-diffusion equations Dγ t u(x, t) = u(x, t)− u2(x, t)− u(x, t)w(x, t) + uxx(x, t), Dγ t w(x, t) = wxx(x, t)− u(x, t)w(x, t), subject to the initial conditions u(x, 0) = epx (1 + e0.5px)2 , w(x, 0) = 1 1 + e0.5px where p is constant. Here N1(u,w) = u−u2−uw+ux, N2(u,w) = wxx−uw, f1(x, t) = f2(x, t) = 0, g1(x) = epx (1+e0.5px)2 , g2(x) = 1 1+e0.5px . The kth− truncated term series solution to the problem is uk(x, t) = k∑ n=0 an(x) tnγ Γ(nγ + 1) and wk(x, t) = k∑ n=0 bn(x) tnγ Γ(nγ + 1) . The iteration for this problem is u0(x, t) = epx (1 + e0.5px)2 , w0(x, t) = 1 1 + e0.5px , Res1,k(x, t) = uk(x, t)− epx (1 + e0.5px)2 − S−1 [ 1 vγα S [ uk−1 ]] + S−1 [ 1 vγα S [ u2k−1 ]] + S−1 [ 1 vγα S [ uk−1wk−1 ]] − S−1 [ 1 vγα S [ (uk−1)xx ]] , k ⩾ 1 Res2,k(x, t) = wk(x, t)− 1 1 + e0.5px − S−1 [ 1 vγα S [ (wk−1)xx ]] + S−1 [ 1 vγα S [ uk−1wk−1 ]] k ⩾ 1. After a few iterations, the coefficients are obtained: a1(x) = epx (1 + e0.5px)4 [ e0.5px − 0.5p2e0.5px + p2 ] , P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3841 b1(x) = 0.25 (1 + e0.5px)3 [ p2epx − 4epx − p2e0.5px ] , and a2(x) = 1 8(1 + e0.5px)6 [ − 32p2e2px + 16e2px + 28p2e1.5px − 33p4e1.5px + 4p2e2.5px + 18p4e2px − p4e2.5px + 8p4epx ] , b2(x) = 1 16(1 + e0.5px)5 [ p4e2px + 16e2px − 8p2e2px + 11p4epx − 11p4e1.5px + 40p2e1.5px − p4e0.5px − 32p2epx − 16e1.5px ] , Hence, the 2-nd order approximate solution is Table 7: Comparison of the approximate solution u2(x, t) in Example 3 using different γ values. t x γ = 0.2 γ = 0.4 γ = 0.6 γ = 0.8 γ = 1 0.2 0.3346368216 0.3075765318 0.2909673743 0.2809869555 0.2750395108 0.4 0.3545182352 0.3265194339 0.3093155911 0.2989687897 0.2927994540 0.1 0.6 0.3747087628 0.3458387005 0.3280801409 0.3173903381 0.3110127944 0.8 0.3951398996 0.3654736576 0.3472053559 0.3361990391 0.3296288598 1.0 0.4157415419 0.3853612031 0.3666326584 0.3553391449 0.3485936321 Table 8: Comparison of the approximate solution w2(x, t) in Example 3 using different γ values. t x γ = 0.2 γ = 0.4 γ = 0.6 γ = 0.8 γ = 1 0.2 0.3968241630 0.4264368572 0.4475212727 0.4616009993 0.4705466352 0.4 0.3798767149 0.4091653431 0.4303163383 0.4445647827 0.4536616511 0.1 0.6 0.3632965320 0.3921407562 0.4132861957 0.4276608484 0.4368839467 0.8 0.3471168287 0.3754019837 0.3964705401 0.4109280813 0.4202514134 1.0 0.3313660940 0.3589848110 0.3799068930 0.3944037158 0.4038005919 u2(x, t) = epx (1 + e0.5px)2 + epx(e0.5px − 0.5p2e0.5px + p2) (1 + e0.5px)4 · tγ Γ(γ + 1) + (18p4e2px − 32p2e2px + 16e2px + 28p2e1.5px − 33p4e1.5px + 4p2e2.5px − p4e2.5px + 8p4epx) 8(1 + e0.5px)6 × t2γ Γ(2γ + 1) , w2(x, t) = 1 1 + e0.5px + 0.25(p2epx − 4epx − p2e0.5px) (1 + e0.5px)3 · tγ Γ(γ + 1) + (p4e2px + 16e2px − 8p2e2px + 11p4epx − 11p4e1.5px + 40p2e1.5px − p4e0.5px − 32p2epx − 16e1.5px) 16(1 + e0.5px)5 × t2γ Γ(2γ + 1) . P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3842 The approximate solution for Example 3 where p = 2/3 are shown in Figure 7. Figures 8 and 9 provide the visualizations to the approximate solutions. Figure 8: The graphics represent the approximate solution u2(x, t) with different γ values: (a) γ = 0.4, (b) γ = 0.6, (c) γ = 0.8, and (d) γ = 1.0 for Example 3. 4. Discussion Based on the findings in Example 1, the suggested method provides a solution in the form of an infinite series that converges to exact solutions. Figure 1 compares the approximate solutions at different γ values with the exact solutions when x = 0.5. Figures 2 and 3 compare the visualizations to the approximate and exact solutions. These graphs show how the solution shapes will approach the precise solution as gamma approaches 1. Tables 1 and 2 deliver a comparison of numerical solutions derived with accurate solutions. Furthermore, an error analysis of the Sadik residual power series method (SRPSM) and the Laplace residual power series method (LRPSM) is carried out, demonstrating that both performed similarly. In Example 2, the approximate solution in the form of an infinite convergent series is found. Figure 4. provides the outcome at different values of γ when x = y = 0.2 and the 3D P. Pue-on / Eur. J. Pure Appl. Math, 17 (4) (2024), 3826-3846 3843 Figure 9: The graphics represent the approximate solution w2(x, t) with different γ values: (a) γ = 0.4, (b) γ = 0.6, (c) γ = 0.8, and (d) γ = 1.0 for Example 3. graphic of solutions are displayed in Figures 5 and 6. Moreover, the numerical simulation when x = y = 0.2 at different values of γ is shown in Tables 3 and 4. The Sadik residual power series method and conformable residual power series method (CRPS) [6] errors are analyzed at the corresponding points shown in Tables 5 and 6. This demonstrates that both have equivalent potential. In Example 3, the graph of approximate solutions at different values of γ when x = 0.5 are indicated in Figure 7 and the 3D graphics of approximate solutions are manifested in Figures 8 and 9. In addition, Tables 7 and 8 present the numerical solution at different values of γ. 5. Conclusion This investigation successfully addressed the system of nonlinear fractional partial differential equations by using the Sadik residual power series method. The employed methodology offers a solution to the nonlinear problem without necessitating the utiliza- tion of Adomian polynomials, linearization techniques, or perturbation processes. The ease of use of the technique due to condition (9) is an advantage compared to traditional power series approaches, which involve differentiation to obtain the series solution coef- ficients. Therefore, the recommended approach requires less time for computation. The REFERENCES 3844 obtained results are expressed as a power series that converges to the closed form via Mittag-Leffler functions, as shown in Example 1 and Example 2, while in Example 3, the outcome is written in the form of a truncated power series. Numerical simulation results and 3D graphics confirm that the method is an efficient, reliable, and accurate method for solving non-linear problems. Based on the information above, the Sadik residual power series approach is an efficient, dependable, and practical method for solving nonlinear fractional order partial differential equation systems. It is suitable for tackling other types of non-linear fractional calculus issues. Acknowledgements This research project was financially supported by Thailand Science Research and Innovation (TSRI). The author is grateful to the reviewers for their informative suggestions and feedback. References [1] O. Abdulaziz, I. Hashim, and S. Momani. Solving systems of fractional differential equations by homotopy-perturbation method. Physics Letters A, 372(4):451–459, 2008. [2] A. T. Alabdala, B. N. Abood, S. S. Redhwan, and S. Alkhatib. Caputo delayed fractional differential equations by Sadik transform. Nonlinear Functional Analysis and Applications, 28(2):439–448, 2023. [3] M. Alaroud. Application of Laplace residual power series method for approximate solutions of fractional IVP’s. Alexandria Engineering Journal, 61(2), 2022. [4] A. A. Alderremy, R. Shah, N. Iqbal, S. Aly, and K. Nonlaopon. Fractional series solu- tion construction for nonlinear fractional reaction-diffusion brusselator model utilizing Laplace residual power series. Symmetry, 14(1944), 2022. [5] A. S. Alshehry, R. Ullah, N. A. Shah, R. Shah, and K. Nonlaopon. Implementation of Yang residual power series method to solve fractional non-linear systems. AIMS Mathematics, 8(4):8294–8309, 2023. [6] A. Arafa. A different approach for conformable fractional biochemical reaction- diffusion models. Applied Mathematics-A Journal Chinese Universities., 35:452–467, 2020. [7] M. Caputo. Linear models of dissipation whose q is almost frequency independent. part ii. Geophysical Journal of the Royal Astronomical Society, 13:529–539, 1967. REFERENCES 3845 [8] J. P. Chauhan and S. R. Khirsariya. A semi-analytic method to solve nonlinear differential equations with arbitrary order. Results in Control and Optimization, 12(100267), 2023. [9] M. Derakhshan. Analytical solutions for the equal width equations containing gener- alized fractional derivative using the efficient combined method. International Journal of Differential Equations, 2021(7066398):14 pages, 2021. [10] A. El-Ajou, O. A. Arqub, and S. Momani. Approximate analytical solution of the nonlinear fractional KdV-Burgers equation: A new iterative algorithm. Journal of Computational Physics., 293:81–95, 2015. [11] A. El-Ajou, O. A. Arqub, Z. A. Zhour, and S. Momani. New results on fractional power series: theories and applications. Entropy, 15:5305–5323, 2013. [12] V. S. Ertürk and S. Momani. Solving systems of fractional differential equations using differential transform method. Journal of Computational and Applied Mathematics, 215(1):142–151, 2008. [13] A. R. Hadhoud, A. A. Rageh, and T. Radwan. Employing the Laplace residual power series method to solve (1+1)- and (2+1)-dimensional time-fractional nonlinear differential equations. Fractal and Fractional, 8(7), 2024. [14] E. Hesameddini and A. Rahimi. A novel iterative method for solving systems of fractional differential equations. Journal of Applied Mathematics, 2013(1), 2012. [15] H. Jafari and V. Daftardar-Gejji. Solving a system of nonlinear fractional differen- tial equations using Adomian decomposition. Journal of Computational and Applied Mathematics, 196(2):644–651, 2006. [16] H. Jafari and S. Seifi. Solving a system of nonlinear fractional partial differential equations using homotopy analysis method. Communications in Nonlinear Science and Numerical Simulation, 14(5):1962–1966, 2009. [17] G. Mittag-Leffler. “sur la nouvelle fonction Eα(x)”. Comptes Rendus de l’Academie des Sciences Paris, 137:554–558, 1903. [18] I. Podlubny. Fractional Differential Equations. Academic Press, New York, 1999. [19] S.N.T. Polat and A.T. Dincel. Solution method for systems of nonlinear fractional differential equations using third kind Chebyshev wavelets. Axioms, 12(546):12 pages, 2023. [20] P. Pue-on. The modified Sadik decomposition method to solve a system of non- linear fractional Volterra integrodifferential equations of convolution type. WSEAS TRANSACTIONS on MATHEMATICS, 20:335–343, 2021. REFERENCES 3846 [21] P. Pue-on. The exact solutions of the space and time fractional telegraph equations by the double Sadik transform method. Mathematics and Statistics, 10(5):995–1004, 2022. [22] P. Pue-on. Exploring the remarkable properties of the double Sadik transform and its applications to fractional Caputo partial differential equations. International Journal of Analysis and Applications., 21(118), 2023. [23] A. Qazza, A. Burqan, and R. Saadeh. Application of ARA-residual power series method in solving systems of fractional differential equations. Mathematical Problems in Engineering, 2022(1), 2021. [24] A. Raheem and A. Afreen. Study of a nonlinear system of fractional differential equations with deviated arguments via Adomian decomposition method. Interna- tional Journal of Applied and Computational Mathematics, 8(269):17 pages, 2022. [25] S. S. Redhwan, S. L. Shaikh, and M. S. Abdo. Some properties of Sadik transform and its applications of fractional-order dynamical systems in control theory. Advances in the Theory of Nonlinear Analysis and its Applications, 4:51–66, 2020. [26] S. S. Redhwan, S. L. Shaikh, M. S. Abdo, and S. Y. Al-Mayyahi. Sadik transform and some result in fractional calculus. Malaya Journal of Matematik, 8(2):536–543, 2020. [27] S. L. Shaikh. Introducing a new integral transform: Sadik transform. American In- ternational Journal of Research in Science, Technology, Engineering & Mathematics, 22:100–102, 2018. [28] K. Yuxiao M. Shuhua Z. Yonghong. Variable order fractional grey model and its application. Applied Mathematical Modelling, 97:619–635, 2021. [29] J. Zhang, X. Chen, L. Li, and C. Zhou. Elzaki transform residual power series method for the fractional population diffusion equations. Engineering Letters, 29(4):12 pages, 2021. [30] E. Ziada. Analytical solution of nonlinear system of fractional differential equations. Journal of Applied Mathematics and Physics, 9(10):2544–2557, 2021.