Jurnal Matematika MANTIK Vol. 7, No. 1, May 2021, pp. 9-19 CONTACT: Sadeq Taha Abdulazeez, sadiq.taha@uod.ac Department of Mathematics, College of Basic Education, University of Duhok, Duhok, Iraq The article can be accessed here. https://doi.org/10.15642/mantik.2021.7.1.9-19 A New Computational Method Based on Integral Transform for Solving Linear and Nonlinear Fractional Systems Diyar Hashim Malo1, Rogash Younis Masiha2, Muhammad Amin Sadiq Murad3, Sadeq Taha Abdulazeez4* 1,3Department of Mathematics, College of Sciences, University of Duhok, Duhok, Iraq 2Department of Statistics, Van Yüzüncü Yıl Üniversitesi, Van, Turkey 4Department of Mathematics, College of Basic Education, University of Duhok, Duhok, Iraq Article history: Received Mar 16, 2021 Revised, Apr 30, 2021 Accepted, May 24, 2021 Kata Kunci: Sistem fraksional stiff, Persamaan diferensial, EHPM, Metode Kernel Hilbert space, Solusi aproksimasi Abstrak. Pada artikel ini, Metode Elzaki homotopy perturbation (EHPM) diterapkan untuk menyelesaikan sistem fraksional stiff. Metode EHPM adalah kombinasi dari modifikasi transformasi integral Laplace yang disebut transformasi Elzaki dan metode gangguan homotopi. Metode yang diusulkan telah diterapkan pada beberapa contoh sistem fraksional stiff linier dan nonlinier. Hasil yang diperoleh dengan metode ini dibandingkan dengan hasil yang diperoleh dengan metode Kernel Hilbert space (KHSM) menunjukkan bahwa metode EHPM merupakan metode yang efektif dan akurat untuk menyelesaikan sistem persamaan diferensial orde fraksional. Keywords: Fractional stiff systems, Differential equations, Elzaki homotopy perturbation methods, Kernel Hilbert space method, Approximate solutions Abstract. In this article, the Elzaki homotopy perturbation method is applied to solve fractional stiff systems. The Elzaki homotopy perturbation method (EHPM) is a combination of a modified Laplace integral transform called the Elzaki transform and the homotopy perturbation method. The proposed method is applied for some examples of linear and nonlinear fractional stiff systems. The results obtained by the current method were compared with the results obtained by the kernel Hilbert space KHSM method. The obtained result reveals that the Elzaki homotopy perturbation method is an effective and accurate technique to solve the systems of differential equations of fractional order. How to cite: D. H. Malo, R. Y. Masiha, M. A. S. Murad, and S. T. Abdulazeez, “A New Computational Method Based Integral Transform for Solving Linear and Nonlinear Fractional Systems”, J. Mat. Mantik, vol. 7, no. 1, pp. 9-19, May. 2021. Jurnal Matematika MANTIK Vol. 7, No. 1, May 2021, pp. 9-19 ISSN: 2527-3159 (print) 2527-3167 (online) mailto:sadiq.taha@uod.ac https://doi.org/10.15642/mantik.2021.7.1.9-19 https://doi.org/10.15642/mantik.2021.7.1.9-19 http://u.lipi.go.id/1458103791 Jurnal Matematika MANTIK Vol. 7, No. 1, May 2021, pp. 9-19 10 1. Introduction Many scientific disciplines can be modeled by initial value problems of fractional order which lead to a better understanding in characterizing the natural phenomena in various fields of science such as engineering, biology, economy, computer science, and physics. The exact and analytical solutions of fractional order differential equations are difficult to be found, thus the numerical methods are used to investigate the solutions of differential equations of fractional orders. The proposed system at first was highlighted by Hirschfelder and Curtiss [1]. Then stiff systems were studied by many researchers in different branches [2]–[4]. Several numerical methods have been developed to obtain the analytical and approximate solutions of stiff systems of integer and fractional orders, such as Adomian decomposition method [5], homotopy perturbation method [6], homotopy analysis method [7], variation iteration method [8], rational homotopy perturbation method[9], Block method [10], Laplace adomian decomposition method and modified decomposition method [11], multistage Bernstein method[12], and fractional power series method [13]. Homotopy perturbation method (HPM) is an important semi-analytic technique for solving differential equations [14]–[16] and an efficient technique to study the different types of nonlinear functional equations. Volterra-Fredholm nonlinear systems were solved by HPM [17], hyperbolic PDEs [18], Zakharov-Kuznetsov[19], and system of nonlinear differential equations[20]. So far, many modified techniques have been described to solve linear and nonlinear differential equations of integer and fractional orders. Those of which are Laplace homotopy perturbation method, Laplace Adomian decomposition method, and recently, Elzaki homotopy transformation perturbation method is used to solve a range of problems such as a family of Fisher’s equation [21], spatial diffusion of Biological population [22], nonlinear oscillators [23], system of linear and nonlinear PDEs of fractional orders [24], time-fractional Navier–Stokes equations [25], and An algorithm for solving the Burgers– Huxley equation using the Elzaki transform [26] Elzaki integral transform is a modification of the Laplace and Sumudu transforms which was invented by Tariq [27], Elzaki transformation is an efficient and powerful technique that has found the exact solutions to several differential equations which cannot be solved by Sumudu transform [28]. Elzaki integral equation is a powerful and efficient technique that has been used to solve many differential equations of integer and fractional orders [25][29]–[33]. The objective of this paper is to expand the applications of EHPM and illustrate the efficiency of the proposed method, thus we consider the stiff systems fractional ordinary differential equations: 𝐷𝜇𝑖 𝑧𝑖(𝑡) + 𝐹𝑖(𝑡, 𝑧1(𝑡), 𝑧2(𝑡), … , 𝑧𝑛 (𝑡)) = 𝑓𝑖 (𝑡), 𝑧𝑖 (𝑡0) = 𝑎𝑖,0, 𝑚 − 1 < 𝜇𝑖 < 𝑚 , 𝑚 ∈ 𝑁. (1) 2. Preliminaries In this section, we introduce some definitions and properties of fractional calculus and Elzaki transform which are used in this article. Definition 2.1. [4] A real valued function 𝑔(𝑦), 𝑦 > 0 is belongs to the space ∁𝜎, 𝜎 ∈ 𝑅 if there exists at least a real number 𝑑 > 𝜎, such that 𝑔(𝑦) = 𝑦 𝑑 𝑔1(𝑦) where 𝑔1(𝑦) ∈ ∁(0, ∞), and it is said to be in the space ∁𝜎 𝑛 if 𝑔𝑛 ∈ 𝑅𝜎 , 𝑛 ∈ 𝑁. Definition 2.2. [34] The function 𝑓(𝑢) is called R-L fractional integral of order ∝≥ 0 if it is defined as: 𝐽𝛼 𝑓(𝑢) = 1 𝛤(𝛼) ∫ 𝑢 0 (𝑢 − 𝑡)𝛼−1𝑓(𝑡)𝑑𝑡, 𝛼 > 0, 𝑡 > 0. In particular 𝐽0𝑓(𝑢) = 𝑓(𝑢). For 𝜃 ≥ 0 and 𝜗 ≥ −1, some properties of the operator 𝐽𝛼 Diyar H. Malo, Rogas Y. Masiha, Muhammad A. S. Murad, and Sadeq T. Abdulazeez A New Computational Method Based Integral Transform for Solving Linear and Nonlinear Fractional Systems 11 1. 𝐽𝛼 𝐽𝜃 𝑓(𝑢) = 𝐽𝛼+𝜃 𝑓(𝑢), 2. 𝐽𝛼 𝐽𝜃 𝑓(𝑢) = 𝐽𝜃 𝐽𝛼 𝑓(𝑢), 3. 𝐽𝛼 𝑥𝜗 = 𝛤(𝜗+1) 𝛤(𝛼+𝜗+1) 𝑥𝛼+𝜗. Definition 2.3. [34] The function 𝑓 ∈ 𝐶−1 𝑛 , 𝑛 ∈ 𝑁 is called Caputo fractional derivative if it is defined as 𝐷𝛼 𝑓(𝑢) = 1 𝛤(𝑛−𝛼) ∫ 𝑢 0 (𝑢 − 𝑡)𝑛−𝛼−1𝑓 𝑛(𝑡)𝑑𝑡, 𝑛 − 1 < 𝛼 ≤ 𝑛. Definition 2.4. [27] The Elzaki-transform of the function 𝑓(𝑢) is defined as: 𝐸[𝑓(𝑢)] = 𝑇(𝑣) = 𝑣 ∫ ∞ 0 𝑓(𝑢)𝑒 −𝑢 𝑣 𝑑𝑢 𝑢 > 0. Suppose that 𝑓 is piecewise continuous, then we can calculate 𝐸 [ 𝜕𝑓 𝜕𝑥 ] as follows: 𝐸 [ 𝜕𝑓(𝑥,𝑢) 𝜕𝑥 ] = ∫ ∞ 0 𝑣𝑒 −𝑢 𝑣 𝜕𝑓(𝑥,𝑢) 𝜕𝑥 𝑑𝑢 = 𝜕 𝜕𝑥 ∫ ∞ 0 𝑣𝑒 −𝑢 𝑣 𝑓(𝑥, 𝑢)𝑑𝑢 = 𝜕 𝜕𝑥 𝑇(𝑥, 𝑣), Similarly, we can have: 𝐸 [ 𝜕2𝑓(𝑥, 𝑢) 𝜕𝑥2 ] = 𝑑2𝑇(𝑥, 𝑢) 𝑑𝑥2 . Assume that 𝜕𝑓 𝜕𝑢 = ℎ, then we have: 𝐸 [ 𝜕2𝑓(𝑥, 𝑢) 𝜕𝑢2 ] = 𝐸 [ 𝜕ℎ(𝑥, 𝑢) 𝜕𝑢 ] = 1 𝑣 𝐸[ℎ(𝑥, 𝑢)] − 𝑣ℎ(𝑥, 0) 𝐸 [ 𝜕2𝑓(𝑥, 𝑢) 𝜕𝑢2 ] = 𝑇(𝑥, 𝑢) 𝑣 2 − 𝑓(𝑥, 0) − 𝑣 𝜕𝑓 𝜕𝑢 (𝑥, 0) By mathematical induction one can extend this result to the 𝑛𝑡ℎpartial derivative as follows: 𝐸 [ 𝜕𝑛𝑓(𝑥,𝑢) 𝜕𝑢𝑛 ] = 𝑇(𝑥,𝑡) 𝑣𝑛 − ∑𝑛−1𝑖=0 𝑣 2−𝑛+𝑖 𝜕 𝑖𝑓(𝑥,0) 𝜕𝑢𝑖 . (2) 3. Elzaki homotopy perturbation method analysis Assume the following system of nonlinear differential equations 𝛭𝑖(𝑧𝑖) + 𝛮𝑖 (𝑧𝑖) = 𝑔𝑖 (𝑡), 𝑡 ∈ 𝛬, 𝑖 = 1,2, … , 𝑛 (3) where 𝑧𝑖 and 𝑔𝑖(𝑡) are sought and known functions respectively, 𝛭𝑖 and 𝛮𝑖 are linear and nonlinear operators respectively, where: 𝛮𝑖(𝑧𝑖) = ∑ ∞ 𝑘=0 𝐴𝑖𝑘 , where 𝐴𝑖𝑘 represents the Adomian polynomial as follows: 𝐴𝑖𝑘 = 1 𝑚! 𝑑𝑚 𝑑𝜏𝑚 (𝛮𝑖 (∑ ∞ 𝑘=0 𝑧1𝑘𝜏 𝑚, … , ∑ ∞ 𝑘=0 𝑧𝑛𝑘𝜏 𝑚 ) , Now, we assume the following homotopy И𝑖 (ϒ𝑖 , 𝑝) = (1 − 𝑝)(𝛭𝑖(ϒ𝑖 ) − 𝑧𝑖,0) + 𝑝(𝛭𝑖(ϒ𝑖 ) + 𝛮𝑖(ϒ𝑖 ) − 𝑔𝑖), or equivalently, И𝑖 (ϒ𝑖 , 𝑝) = 𝛭𝑖(ϒ𝑖 ) − 𝑧𝑖,0 + 𝑝𝑧𝑖,0 + 𝑝(𝛮𝑖(ϒ𝑖 ) − 𝑔𝑖), (4) Since 𝛭𝑖(ϒ𝑖 ) − 𝑧𝑖,0 = 0 and 𝛭𝑖(ϒ𝑖 ) + 𝛮𝑖(ϒ𝑖 ) − 𝑔𝑖 = 0 which are equivalent to the operator equations И𝑖(𝑧, 0) = 0 and И𝑖(𝑧, 1) = 0 respectively, also 𝑝 ∈ [0,1] is a homotopy parameter, 𝑧𝑖,0 are initial approximations and ϒ𝑖 : (𝑡, 𝑝): 𝛬 × [0,1] → 𝑅. Now, we apply the Elzaki transform on (4), we obtain 𝐸[𝛭𝑖(ϒ𝑖 ) − 𝑧𝑖,0 + 𝑝𝑧𝑖,0 + 𝑝(𝛮𝑖 (ϒ𝑖 ) − 𝑔𝑖 )] = 0 By differential property of of Elzaki transform, we get 1 𝑣𝜇 𝐸[ϒ𝑖 ] − 𝑣 2−𝜇 ϒ𝑖,0 − 𝑣 3−𝜇 ϒ′ 𝑖,0 − ⋯ − 𝑣 𝑛+1−𝜇 ϒ𝑖,0 (𝑛−1) = 𝐸[𝑧𝑖,0 − 𝑝𝑧𝑖,0 + 𝑝(𝛮𝑖(ϒ𝑖 ) − 𝑔𝑖)] Jurnal Matematika MANTIK Vol. 7, No. 1, May 2021, pp. 9-19 12 Applying the inverse Elzaki transform, we obtain ϒ𝑖 = 𝐸 −1 [𝑣 2ϒ𝑖,0 + 𝑣 3ϒ′𝑖,0 + ⋯ + 𝑣 𝑛+1ϒ𝑖,0 (𝑛−1) ] + 𝐸−1 [𝑣 𝜇 𝐸[𝑧𝑖,0 − 𝑝𝑧𝑖,0 + 𝑝(𝛮𝑖 (ϒ𝑖) − 𝑔𝑖)]]. (5) By HPM, the following series represent the solution of equation (5), where 𝑝 ∈ [0,1]. ϒ𝑖 (𝑡) = ∑ ∞ 𝑗=0 𝑝 𝑗 ϒ𝑖,𝑗 , 𝑖 = 1,2, … , 𝑛 (6) Substituting (6) in (5), we obtain ∑ ∞ 𝑛=0 𝑝𝑛ϒ𝑖,𝑛 = 𝐸 −1 [𝑣 2ϒ𝑖,0 + 𝑣 3ϒ′𝑖,0 + ⋯ + 𝑣 𝑛+1ϒ𝑖,0 (𝑛−1) ] + 𝐸−1 [𝑣 𝜇 𝐸 [𝑧𝑖,0 − 𝑝𝑧𝑖,0 + 𝑝 (𝛮𝑖 (∑ ∞ 𝑛=0 𝑝𝑛ϒ𝑖,𝑛) − 𝑔𝑖)]] Now, by comparing the coefficients of the same powers of the embedded parameter 𝑝 , leads to 𝑝0 ∶ ϒ𝑖,0 = 𝐸 −1 [𝑣 2ϒ𝑖,0 + 𝑣 3ϒ′ 𝑖,0 + ⋯ + 𝑣 𝑛+1ϒ𝑖,0 (𝑛−1) + 𝐸[𝑧0]] 𝑝1 ∶ ϒ𝑖,1 = 𝐸 −1 [𝑣 𝜇 𝐸[𝛮𝑖 (ϒ𝑖,0) − 𝑧𝑖,0 − 𝑔𝑖]] 𝑝2 ∶ ϒ𝑖,2 = 𝐸 −1 [𝑣 𝜇 𝐸[𝛮𝑖(ϒ𝑖,0, ϒ𝑖,1)]] 𝑝3 ∶ ϒ𝑖,3 = 𝐸 −1 [𝑣 𝜇 𝐸[𝛮𝑖 (ϒ𝑖,0, ϒ𝑖,1, ϒ𝑖,2)]] ⋮ 𝑝 𝑗 ∶ ϒ𝑖,𝑗 = 𝐸 −1 [𝑣 𝜇 𝐸[𝛮𝑖(ϒ𝑖,0, ϒ𝑖,1, ϒ𝑖,2, … , ϒ𝑖,𝑗−1 )]] ⋮ Assume that ϒ𝑖,0 = 𝛽𝑖,0, ϒ ′ 𝑖,0 = 𝛽𝑖,1, … , ϒ𝑖,0 (𝑛−1) = 𝛽𝑖,𝑛−1 are the initial approximations, where 𝑖 = 1, 2, … , 𝑛. Finally, the solution of (3) is obtained by the following series: 𝑧𝑖(𝑡) = ϒ𝑖 = ϒ𝑖,0 + ϒ𝑖,1 + ϒ𝑖,2 + ⋯ +. (7) The convergence and uniqueness of the series (7) is discussed in [21]. 4. Applications of Elzaki homotopy perturbation method The numerical patterns are employed to confirm the efficiency and accuracy of the Elzaki homotopy transform perturbation technique for solving the fractional stiff systems. For this purpose, at the first example, the comparison of error analysis with the kernel Hilbert space technique [13] [35] is made, at the second example, the exact solution of the proposed problem is achieved where 𝛽 = 1, and at the third example the comparison between the proposed method and (ADM and MLDM [11]) is made. Example 4.1 Consider the following linear stiff system of fractional order 0 < 𝛽 ≤ 1 𝐷𝑡 𝛽 𝑧(𝑡) = −𝑧(𝑡) + 95𝑤(𝑡), 𝐷𝑡 𝛽 𝑤(𝑡) = −𝑧(𝑡) − 97𝑤(𝑡), (8) with the initial conditions 𝑧(0) = 1, 𝑤(0) = 1. The exact solution of the system when 𝛽 = 1 is 𝑧(𝑡) = 95 47 𝑒−2𝑡 − 48 47 𝑒−96𝑡 , 𝑤(𝑡) = 48 47 𝑒−96𝑡 − 1 47 𝑒−2𝑡 . To solve system (8) by the EHPM, we use the following homotopy 𝐷𝑡 𝛽 𝑍(𝑡) − 𝑧0(𝑡) + 𝑝(𝑧0(𝑡) + 𝑍(𝑡) − 95𝑊 (𝑡)) = 0 𝐷𝑡 𝛽 𝑊(𝑡) − 𝑤0(𝑡) + 𝑝(𝑤0(𝑡) + 𝑍(𝑡) + 97𝑊 (𝑡)) = 0. Applying the Elzaki transform, we have Diyar H. Malo, Rogas Y. Masiha, Muhammad A. S. Murad, and Sadeq T. Abdulazeez A New Computational Method Based Integral Transform for Solving Linear and Nonlinear Fractional Systems 13 𝐸 [𝐷𝑡 𝛽 𝑍(𝑡) − 𝑧0(𝑡) + 𝑝(𝑧0(𝑡) + 𝑍(𝑡) − 95𝑊(𝑡))] = 0 𝐸 [𝐷𝑡 𝛽 𝑊(𝑡) − 𝑤0(𝑡) + 𝑝(𝑤0(𝑡) + 𝑍(𝑡) + 97𝑊 (𝑡))] = 0. Using the differential property of Elzaki transform, we obtain 1 𝑣 𝛽 𝐸[𝑍(𝑡)] = 𝑣 2−𝛽 𝑍0(𝑡) + 𝐸[𝑧0(𝑡) − 𝑝(𝑧0(𝑡) + 𝑍(𝑡) − 95𝑊(𝑡))] 1 𝑣𝛽 𝐸[𝑊(𝑡)] = 𝑣 2−𝛽𝑊0(𝑡) + 𝐸[𝑤0(𝑡) − 𝑝(𝑧0(𝑡) + 𝑍(𝑡) + 97𝑊 (𝑡))]. Applying inverse Elzaki transform, we have 𝑍(𝑡) = 𝐸−1 [𝑣 2𝑍0(𝑡) + 𝑣 𝛽𝐸[𝑧0(𝑡) − 𝑝(𝑧0(𝑡) + 𝑍(𝑡) − 95𝑊(𝑡))]] 𝑊(𝑡) = 𝐸−1 [𝑣 2𝑊0(𝑡) + 𝑣 𝛽 𝐸[𝑤0(𝑡) − 𝑝(𝑧0(𝑡) + 𝑍(𝑡) + 97𝑊(𝑡))]]. (9) Here, the following form is the solution of equation (9): 𝑍(𝑡) = 𝑍0(𝑡) + 𝑝𝑍1(𝑡) + 𝑝 2𝑍2(𝑡) + ⋯, 𝑊(𝑡) = 𝑊0(𝑡) + 𝑝𝑊1(𝑡) + 𝑝 2𝑊2(𝑡) + ⋯, (10) Substituting (10) in to (9) and collecting the coefficients of equivalent powers of embedded parameter 𝑝 , we obtain the following: 𝑝0 ∶ {𝑍0(𝑡) = 𝐸 −1 [𝑣 2𝑍0(0) + 𝑣 𝛽𝐸[𝑧0(𝑡)]] , 𝑊0(𝑡) = 𝐸 −1 [𝑣 2𝑊0(0) + 𝑣 𝛽𝐸[𝑤0(𝑡)]], 𝑝1 ∶ {𝑍1(𝑡) = 𝐸 −1 [−𝑣 𝛽𝐸[(𝑧0(𝑡) + 𝑍0(𝑡) − 95𝑊0(𝑡))]] , 𝑊1(𝑡) = 𝐸−1 [−𝑣 𝛽𝐸[(𝑧0(𝑡) + 𝑍0(𝑡) + 97𝑊0(𝑡))]], 𝑝2 ∶ {𝑍2(𝑡) = 𝐸 −1 [−𝑣 𝛽𝐸[(𝑍1(𝑡) − 95𝑊1(𝑡))]] , 𝑊2(𝑡) = 𝐸 −1 [−𝑣 𝛽 𝐸[(𝑍1(𝑡) + 97𝑊1(𝑡))]], 𝑝3 ∶ {𝑍3(𝑡) = 𝐸 −1 [−𝑣 𝛽𝐸[(𝑍2(𝑡) − 95𝑊2(𝑡))]] , 𝑊3(𝑡) = 𝐸 −1 [−𝑣 𝛽𝐸[(𝑍2(𝑡) + 97𝑊2(𝑡))]], ⋮ Here, we set 𝑍0(0) = 𝑧0(𝑡) = 1 and 𝑊0(0) = 𝑤0(𝑡) = 1. Thus, the above equations lead to the following results: 𝑍0(𝑡) = 1 + 𝑡𝛽 𝛤(𝛽+1) , 𝑊0(𝑡) = 1 + 𝑡𝛽 𝛤(𝛽+1) , 𝑍1(𝑡) = 93 𝛤(𝛽+1) 𝑡𝛽 + 94 𝛤(2𝛽+1) 𝑡 2𝛽 , 𝑊1(𝑡) = − 99 𝛤(𝛽+1) 𝑡𝛽 − 98 𝛤(2𝛽+1) 𝑡 2𝛽 , 𝑍2(𝑡) = − 9498 𝛤(𝛽+1) 𝑡 2𝛽 − 9404 𝛤(3𝛽+1) 𝑡 3𝛽 , 𝑊2(𝑡) = 9510 𝛤(2𝛽+1) 𝑡 2𝛽 + 9412 𝛤(3𝛽+1) 𝑡 3𝛽 𝑍3(𝑡) = 9129486 𝛤(3𝛽+1) 𝑡 3𝛽 + 9035446 𝛤(4𝛽+1) 𝑡 4𝛽 , 𝑊3(𝑡) = − 912972 𝛤(3𝛽+1) 𝑡 3𝛽 − 903560 𝛤(4𝛽+1) 𝑡 4𝛽 , Thus, the solution of system (8) using the series (7) is 𝑧(𝑡) = 1 + 94 𝛤(𝛽 + 1) 𝑡𝛽 − 9404 𝛤(2𝛽 + 1) 𝑡 2𝛽 + 903544 𝛤(3𝛽 + 1) 𝑡 3𝛽 …, 𝑤(𝑡) = 1 − 98 𝛤(𝛽 + 1) 𝑡𝛽 + 9412 𝛤(2𝛽 + 1) 𝑡 2𝛽 − 903560 𝛤(3𝛽 + 1) 𝑡 3𝛽 …, Jurnal Matematika MANTIK Vol. 7, No. 1, May 2021, pp. 9-19 14 Figure 1. The behavior of EHPM solution of system (8) for different values of 𝛽 and exact solution. Figure 1: illustrates the behavior of the EHPM solution of system (8) for different values of 𝛽 and exact solution. It can be seen that there is a good agreement between the exact solution and the approximate solution using EHPM, especially when 𝛽 = 1. Table 1. The Absolute Error of kernel Hilbert space method KHSM and EHPM of (8) 𝒕 𝒛(𝒕)EHPM 𝒛(𝒕)KHSM 𝒘(𝒕)EHPM 𝒘(𝒕)KHSM 0.00 0 0 0 0 0.5 3.0163443 × 10−15 4.76706570 × 10−4 3.0163483 × 10−15 4.44439990 × 10−4 0.1 5.0178941 × 10−14 4.11985309 × 10−5 3.4785609 × 10−13 1.37132721 × 10−5 0.15 6.4189477 × 10−11 2.37714435 × 10−5 1.1162905 × 10−10 4.49704615 × 10−7 0.2 7.3684661 × 10−9 1.99138385 × 10−5 1.2169510 × 10−9 2.1210257 × 10−7 0.25 3.4096780 × 10−7 1.67357026 × 10−5 3.6738962 × 10−7 1.76192995 × 10−7 Table 1: shows the comparison between the absolute errors of the approximate solution of the proposed technique at different values of 𝑡 when 𝛽 = 1 for 20 order approximation and the approximate solution of KHSM. The results illustrate that the proposed technique is superior to KHSM and an efficient method for investigating the solution of proposed fractional systems. Diyar H. Malo, Rogas Y. Masiha, Muhammad A. S. Murad, and Sadeq T. Abdulazeez A New Computational Method Based Integral Transform for Solving Linear and Nonlinear Fractional Systems 15 Example 4.2: Consider the following non-linear stiff system of fractional order 𝐷𝑡 𝛽 𝑧(𝑡) = −(𝜌−1 + 2)𝑧(𝑡) + 𝜌−1𝑤 2(𝑡), 0 < 𝛽 ≤ 1, 𝐷𝑡 𝛽 𝑤(𝑡) = 𝑧(𝑡) − 𝑤(𝑡) − 𝑤 2(𝑡), 𝑡 ∈ [0,2] (11) with the initial conditions 𝑧(0) = 1, 𝑤(0) = 1. The exact solution of the system when 𝛽 = 1 is 𝑧(𝑡) = 𝑒−2𝑡 , 𝑤(𝑡) = 𝑒𝑡 . To solve system (11) by the EHPM, we use the following homotopy 𝐷𝑡 𝛽 𝑍(𝑡) − 𝑧0(𝑡) + 𝑝(𝑧0(𝑡) + (𝜌 −1 + 2)𝑍(𝑡) − 𝜌−1𝑊2(𝑡)) = 0 𝐷𝑡 𝛽 𝑊(𝑡) − 𝑤0(𝑡) + 𝑝(𝑤0(𝑡) − 𝑍(𝑡) + 𝑊(𝑡) + 𝑊 2 (𝑡)) = 0. Applying the Elzaki transform, we have 𝐸 [𝐷𝑡 𝛽 𝑍(𝑡) − 𝑧0(𝑡) + 𝑝(𝑧0(𝑡) + (𝜌 −1 + 2)𝑍(𝑡) − 𝜌−1𝑊2(𝑡))] = 0 𝐸 [𝐷𝑡 𝛽 𝑊(𝑡) − 𝑤0(𝑡) + 𝑝(𝑤0(𝑡) − 𝑍(𝑡) + 𝑊(𝑡) + 𝑊 2(𝑡))] = 0. Using the differential property of Elzaki transform, we obtain 1 𝑣 𝛽 𝐸[𝑍(𝑡)] = 𝑣 2−𝛽 𝑍0(𝑡) + 𝐸[𝑧0(𝑡) − 𝑝(𝑧0(𝑡) + (𝜌 −1 + 2)𝑍(𝑡) − 𝜌−1𝑊2(𝑡))] 1 𝑣𝛽 𝐸[𝑊(𝑡)] = 𝑣 2−𝛽𝑊0(𝑡) + 𝐸[𝑤0(𝑡) − 𝑝(𝑤0(𝑡) − 𝑍(𝑡) + 𝑊(𝑡) + 𝑊 2 (𝑡))]. Applying inverse Elzaki transform, we have 𝑍(𝑡) = 𝐸−1 [𝑣 2𝑍0(𝑡) + 𝑣 𝛽𝐸[𝑧0(𝑡) − 𝑝(𝑧0(𝑡) + (𝜌 −1 + 2)𝑍(𝑡) − 𝜌−1𝑊2(𝑡))]] 𝑊(𝑡) = 𝐸−1 [𝑣 2𝑊0(𝑡) + 𝑣 𝛽 𝐸[𝑤0(𝑡) − 𝑝(𝑤0(𝑡) − 𝑍(𝑡) + 𝑊(𝑡) + 𝑊 2(𝑡))]]. (12) Here, the following form is the solution of system (11): 𝑍(𝑡) = 𝑍0(𝑡) + 𝑝𝑍1(𝑡) + 𝑝 2𝑍2(𝑡) + ⋯, 𝑊(𝑡) = 𝑊0(𝑡) + 𝑝𝑊1(𝑡) + 𝑝 2𝑊2(𝑡) + ⋯. (13) Substituting (13) in to (12) and collecting the coefficients of equivalent powers of embedded parameter 𝑝 , we obtain the following: 𝑝0 ∶ {𝑍0(𝑡) = 𝐸 −1 [𝑣 2𝑍0(0) + 𝑣 𝛽𝐸[𝑧0(𝑡)]] , 𝑊0(𝑡) = 𝐸 −1 [𝑣 2𝑊0(0) + 𝑣 𝛽𝐸[𝑤0(𝑡)]], 𝑝1 ∶ {𝑍1(𝑡) = 𝐸 −1 [−𝑣 𝛽𝐸[𝑧0(𝑡) + (𝜌 −1 + 2)𝑍0(𝑡) − 𝜌 −1𝑊0 2(𝑡)]] , 𝑊1(𝑡) = 𝐸−1 [−𝑣 𝛽𝐸[𝑤0(𝑡) − 𝑍0(𝑡) + 𝑊0(𝑡) + 𝑊0 2(𝑡)]], 𝑝 𝑗 : { 𝑍2(𝑡) = 𝐸−1 [−𝑣 𝛽𝐸 [(𝜌−1 + 2)𝑍𝑗−1(𝑡) − 𝜌−1 ∑ 𝑗−1 𝑙=0 𝑊𝑗 (𝑡)𝑊𝑗−𝑙−1(𝑡)]] , 𝑊2(𝑡) = 𝐸 −1 [−𝑣 𝛽𝐸 [−𝑍𝑗−1(𝑡) + 𝑊𝑗−1(𝑡) + ∑ 𝑗−1 𝑙=0 𝑊𝑗 (𝑡)𝑊𝑗−𝑙−1(𝑡)]] , 𝑗 = 2,3, …, ⋮ Here, we set 𝑍0(0) = 𝑧0(𝑡) = 1 and 𝑊0(0) = 𝑤0(𝑡) = 1. Thus, the above equations lead the following results. 𝑍0(𝑡) = 1 + 𝑡𝛽 𝛤(𝛽 + 1) , Jurnal Matematika MANTIK Vol. 7, No. 1, May 2021, pp. 9-19 16 𝑊0(𝑡) = 1 + 𝑡𝛽 𝛤(𝛽 + 1) , 𝑍1(𝑡) = −3 𝑡𝛽 𝛤(𝛽 + 1) + (𝜌−1 − 2) 𝑡 2𝛽 𝛤(2𝛽 + 1) + 𝜌−1𝛤(2𝛽 + 1) 𝑡 3𝛽 𝛤 2(𝛽 + 1)𝛤(3𝛽 + 1) , 𝑊1(𝑡) = −2 𝑡𝛽 𝛤(𝛽 + 1) − 2 𝑡 2𝛽 𝛤(2𝛽 + 1) − 𝛤(2𝛽 + 1)𝑡 3𝛽 𝛤 2(𝛽 + 1)𝛤(3𝛽 + 1) , 𝑍2(𝑡) = (6 − 𝜌 −1) 𝑡 2𝛽 𝛤(2𝛽 + 1) − ( 1 𝜌2 + 4𝜌−1 + 4𝜌−1𝛤(2𝛽 + 1) 𝛤 2(𝛽 + 1) − 4) 𝑡 3𝛽 𝛤(3𝛽 + 1) − ( 1 𝜌2 + 4𝜌−1 + 4𝜌−1𝛤(3𝛽 + 1)𝛤(𝛽 + 1) 𝛤 2(2𝛽 + 1) ) 𝛤(2𝛽 + 1) 𝑡 4𝛽 𝛤 2(𝛽 + 1)𝛤(4𝛽 + 1) − 2𝜌−1𝛤(2𝛽 + 1)𝛤(4𝛽 + 1) 𝑡 5𝛽 𝛤 3(𝛽 + 1)𝛤(3𝛽 + 1)𝛤(5𝛽 + 1) , 𝑊2(𝑡) = 3 𝑡 2𝛽 𝛤(2𝛽 + 1) + (𝜌−1 + 4 + 4𝛤(2𝛽 + 1) 𝛤 2(𝛽 + 1) ) 𝑡 3𝛽 𝛤(3𝛽 + 1) + (2𝜌−1 + 6 + 4𝛤(3𝛽 + 1) 𝛤 2(𝛽 + 1)𝛤(2𝛽 + 1) ) 𝑡 4𝛽 𝛤(4𝛽 + 1) + 2𝛤(2𝛽 + 1)𝛤(4𝛽 + 1)𝑡 5𝛽 𝛤 3(𝛽 + 1)𝛤(3𝛽 + 1)𝛤(5𝛽 + 1) , ⋮ Thus, the solution of system (11) using the series (7) is 𝑧(𝑡) = 𝑍0(𝑡) + 𝑍1(𝑡) + 𝑍2(𝑡) + ⋯, 𝑤(𝑡) = 𝑊0(𝑡) + 𝑊1(𝑡) + 𝑊2(𝑡) + ⋯, 𝑧(𝑡) = 1 − 2𝑡𝛽 𝛤(𝛽 + 1) + 4𝑡 2𝛽 𝛤(2𝛽 + 1) + 𝜌−1𝛤(2𝛽 + 1) 𝑡 3𝛽 𝛤 2(𝛽 + 1)𝛤(3𝛽 + 1) − ( 1 𝜌2 + 4𝜌−1 + 4𝜌−1𝛤(2𝛽 + 1) 𝛤 2(𝛽 + 1) − 4) 𝑡 3𝛽 𝛤(3𝛽 + 1) − ( 1 𝜌2 + 4𝜌−1 + 4𝜌−1𝛤(3𝛽 + 1)𝛤(𝛽 + 1) 𝛤 2(2𝛽 + 1) ) 𝛤(2𝛽 + 1) 𝑡 4𝛽 𝛤 2(𝛽 + 1)𝛤(4𝛽 + 1) − 2𝜌−1𝛤(2𝛽 + 1)𝛤(4𝛽 + 1) 𝑡 5𝛽 𝛤 3(𝛽 + 1)𝛤(3𝛽 + 1)𝛤(5𝛽 + 1) + ⋯, 𝑤(𝑡) = 1 − 𝑡𝛽 𝛤(𝛽 + 1) + 𝑡 2𝛽 𝛤(2𝛽 + 1) − 𝛤(2𝛽 + 1)𝑡 3𝛽 𝛤 2(𝛽 + 1)𝛤(3𝛽 + 1) + (𝜌−1 + 4 + 4𝛤(2𝛽 + 1) 𝛤 2(𝛽 + 1) ) 𝑡 3𝛽 𝛤(3𝛽 + 1) + (2𝜌−1 + 6 + 4𝛤(3𝛽 + 1) 𝛤 2(𝛽 + 1)𝛤(2𝛽 + 1) ) 𝑡 4𝛽 𝛤(4𝛽 + 1) + 2𝛤(2𝛽 + 1)𝛤(4𝛽 + 1)𝑡 5𝛽 𝛤 3(𝛽 + 1)𝛤(3𝛽 + 1)𝛤(5𝛽 + 1) + ⋯. Here, setting 𝛽 = 1 and follow the above solution, the following results are obtained: 𝑍0(𝑡) = 1 + 𝑡, 𝑊0(𝑡) = 1 + 𝑡 𝑍1(𝑡) = −3𝑡 + (𝜌 −1 − 2) 𝑡2 2! + 2 𝜌−1 𝑡 3 3! , 𝑊1(𝑡) = −2𝑡 − 2 𝑡 2 − 2𝑡 3 3! , 𝑍2(𝑡) = (6 − 𝜌 −1) 𝑡 2 2 − ( 1 𝜌2 + 4𝜌−1 + 8𝜌−1 − 4) 𝑡 3 3! − ( 2 𝜌2 + 20𝜌−1) 𝑡 4 4! − 2𝜌−1 𝑡 5 15 , 𝑊2(𝑡) = 3 𝑡 2 2! + (𝜌−1 + 12) 𝑡 3 3! + (2𝜌−1 + 3 4 ) 𝑡 4𝛽 4! + 2 𝑡 5 15 , Diyar H. Malo, Rogas Y. Masiha, Muhammad A. S. Murad, and Sadeq T. Abdulazeez A New Computational Method Based Integral Transform for Solving Linear and Nonlinear Fractional Systems 17 One can express the above values in a series after finding the other terms of the solution 𝑧(𝑡) = 1 − 2(𝑡 − 𝑡 2) − 8 6 𝑡 3 − 4 6 𝑡 4 + ⋯ = ∑ ∞ 𝑘=0 (2𝑡)𝑘 (−1)𝑘 𝑘! = 𝑒 −2𝑡 𝑤(𝑡) = 1 − 𝑡 + 𝑡 2 2 − 𝑡 3 6 + 𝑡 4 24 + ⋯ = ∑ ∞ 𝑘=0 (𝑡)𝑘 (−1)𝑘 𝑘! = 𝑒−𝑡 , which is the exact solution of system (11). Table 2. The Absolute Error of kernel Hilbert space method KHSM and EHPM of (8) 𝑡 𝑧(𝑡)EHPM 𝑧(𝑡)KHSM 𝑤(𝑡)EHPM 𝑤(𝑡)KHSM 0.00 0 0 0 1.2 × 10−6 0.4 5.55 × 10−17 1.20 × 10−6 0 2.47 × 10−2 0.8 5.55 × 10−16 1.28 × 10−6 0 1.19 × 10−1 1.2 1.70 × 10−12 9.10 × 10−7 5.55 × 10−17 2.62 × 10−1 1.6 6.90 × 10−10 5.59 × 10−7 0 4.38 × 10−1 2.0 7.30 × 10−8 3.18 × 10−7 5.56 × 10−16 1.2 × 10−6 Table 2: illustrates the comparison between the absolute errors of the numerical solution of the proposed technique at different values of 𝑡 where 𝛽 = 1 for 20 terms iterations and the solution obtained by KHSM. It is clear that the EHPM technique is more accurate and converges faster than the KHSM technique. 5. Conclusions In this work, a new computational technique namely Elzaki homotopy perturbation method are constructed to solve linear and nonlinear systems of differential equation of fractional order. The construction scheme of the current method is discussed and implemented to show the performance of the computational method for this problem. The results obtained showed that the proposed method is a powerful and efficient method for solving linear and nonlinear stiff systems of fractional orders and compared the obtained results of the current method with the results obtained by the other methods. Table 1 and table 2 observed that the suggested technique is superior to kernel Hilbert space method KHSM. References [1] C. F. Curtiss and J. O. Hirschfelder, “Integration of Stiff Equations.,” Proc. Natl. Acad. Sci. U. S. A., vol. 38, no. 3, pp. 235–243, Mar. 1952, doi: 10.1073/pnas.38.3.235. [2] J. E. Flaherty and R. E. O’Malley, “The Numerical Solution of Boundary Value Problems for Stiff Differential Equations,” Math. Comput., vol. 31, no. 137, p. 66, 1977, doi: 10.2307/2005781. [3] D. Zwillinger, “Stiff Equations” Handb. Differ. Equations, vol. 38, no. 1950, pp. 690–694, 1992, doi: 10.1016/b978-0-12-784391-9.50181-0. [4] R. V Slonevskii and R. R. Stolyarchuk, “Rational-fractional methods for solving stiff systems of differential equations,” J. Math. Sci., vol. 150, no. 5, pp. 2434–2438, 2008, doi: 10.1007/s10958-008-0141-x. [5] H. Jafari and V. Daftardar-Gejji, “Solving a system of nonlinear fractional differential equations using Adomian decomposition,” J. Comput. Appl. Math., vol. Jurnal Matematika MANTIK Vol. 7, No. 1, May 2021, pp. 9-19 18 196, no. 2, pp. 644–651, 2006, doi: 10.1016/j.cam.2005.10.017. [6] H. Aminikhah and M. Hemmatnezhad, “An effective modification of the homotopy perturbation method for stiff systems of ordinary differential equations,” Appl. Math. Lett., vol. 24, no. 9, pp. 1502–1508, 2011, doi: 10.1016/j.aml.2011.03.032. [7] A. Rani, M. Saeed, Q. M. Ul-Hassan, M. Ashraf, M. Y. Khan, and K. Ayub, “Solving system of differential equations of fractional order by homotopy analysis method,” J. Sci. Arts, vol. 17, no. 3, pp. 457–468, 2017. [8] M. T. Darvishi, F. Khani, and A. A. Soliman, “The numerical simulation for stiff systems of ordinary differential equations,” Comput. Math. with Appl., vol. 54, no. 7–8, pp. 1055–1063, 2007. [9] J. Biazar, M. A. Asadi, and F. Salehi, “Rational homotopy perturbation method for solving stiff systems of ordinary differential equations,” Appl. Math. Model., vol. 39, no. 3–4, pp. 1291–1299, 2015, doi: 10.1016/j.apm.2014.09.003. [10] O. A. Akinfenwa, B. Akinnukawe, and S. B. Mudasiru, “A family of Continuous Third Derivative Block Methods for solving stiff systems of first order ordinary differential equations,” J. Niger. Math. Soc., vol. 34, no. 2, pp. 160–168, 2015, doi: 10.1016/j.jnnms.2015.06.002. [11] O. H. Mohammed and H. A. Salim, “Computational methods based laplace decomposition for solving nonlinear system of fractional order differential equations,” Alexandria Eng. J., vol. 57, no. 4, pp. 3549–3557, 2018, doi: 10.1016/j.aej.2017.11.020. [12] M. H. T. Alshbool and I. Hashim, “Multistage Bernstein polynomials for the solutions of the Fractional Order Stiff Systems,” J. King Saud Univ. - Sci., vol. 28, no. 4, pp. 280–285, 2016, doi: 10.1016/j.jksus.2015.06.001. [13] A. Freihet, S. Hasan, M. Al-Smadi, M. Gaith, and S. Momani, “Construction of fractional power series solutions to fractional stiff system using residual functions algorithm,” Adv. Differ. Equations, vol. 2019, no. 1, 2019, doi: 10.1186/s13662- 019-2042-3. [14] M. Kumar and A. S. Saxena, “New iterative method for solving higher order KDV equations,” pp. 246–257. [15] M. Javidi and B. Ahmad, “Numerical solution of fourth-order time-fractional partial differential equations with variable coefficients,” J. Appl. Anal. Comput., vol. 5, no. 1, pp. 52–63, 2015, doi: 10.11948/2015005. [16] D. H. Shou, “The homotopy perturbation method for nonlinear oscillators,” Comput. Math. with Appl., vol. 58, no. 11–12, pp. 2456–2459, 2009, doi: 10.1016/j.camwa.2009.03.034. [17] J. Biazar, B. Ghanbari, M. G. Porshokouhi, and M. G. Porshokouhi, “He’s homotopy perturbation method: A strongly promising method for solving non- linear systems of the mixed Volterra–Fredholm integral equations,” Comput. Math. with Appl., vol. 61, no. 4, pp. 1016–1023, 2011, doi: https://doi.org/10.1016/j.camwa.2010.12.051. [18] J. Biazar and H. Ghazvini, “Homotopy perturbation method for solving hyperbolic partial differential equations,” Comput. Math. with Appl., vol. 56, no. 2, pp. 453– 458, 2008, doi: https://doi.org/10.1016/j.camwa.2007.10.032. [19] J. Biazar, F. Badpeima, and F. Azimi, “Application of the homotopy perturbation method to Zakharov–Kuznetsov equations,” Comput. Math. with Appl., vol. 58, no. 11, pp. 2391–2394, 2009, doi: https://doi.org/10.1016/j.camwa.2009.03.102. Diyar H. Malo, Rogas Y. Masiha, Muhammad A. S. Murad, and Sadeq T. Abdulazeez A New Computational Method Based Integral Transform for Solving Linear and Nonlinear Fractional Systems 19 [20] T. M. Elzaki and J. Biazar, “Homotopy perturbation method and Elzaki transform for solving system of nonlinear partial differential equations,” World Appl. Sci. J., vol. 24, no. 7, pp. 944–948, 2013, doi: 10.5829/idosi.wasj.2013.24.07.1041. [21] A. C. Loyinmi and T. K. Akinfe, “Exact solutions to the family of Fisher’s reaction‐ diffusion equation using Elzaki homotopy transformation perturbation method,” Eng. Reports, vol. 2, no. 2, pp. 1–32, 2020, doi: 10.1002/eng2.12084. [22] J. Ul Rahman, D. Lu, M. Suleman, J. H. He, and M. Ramzan, “He-elzaki method for spatial diffusion of biological population,” Fractals, vol. 27, no. 5, 2019, doi: 10.1142/S0218348X19500695. [23] N. Anjum, M. Suleman, D. Lu, J. H. He, and M. Ramzan, “Numerical iteration for nonlinear oscillators by Elzaki transform,” J. Low Freq. Noise Vib. Act. Control, 2019, doi: 10.1177/1461348419873470. [24] D. Lu, M. Suleman, J. H. He, U. Farooq, S. Noeiaghdam, and F. A. Chandio, “Elzaki projected differential transform method for fractional order system of linear and nonlinear fractional partial differential equation,” Fractals, vol. 26, no. 3, 2018, doi: 10.1142/S0218348X1850041X. [25] R. M. Jena and S. Chakraverty, “Solving time-fractional Navier–Stokes equations using homotopy perturbation Elzaki transform,” SN Appl. Sci., vol. 1, no. 1, pp. 1– 13, 2019, doi: 10.1007/s42452-018-0016-9. [26] A. C. Loyinmi and T. K. Akinfe, “An algorithm for solving the Burgers–Huxley equation using the Elzaki transform,” SN Appl. Sci., vol. 2, no. 1, p. 7, 2020. [27] T. M. Elzaki, “The New Integral Transform ‘ELzaki Transform’,” vol. 7, no. 1, pp. 57–64, 2011. [28] E. M. A. Hilal, “Elzaki and Sumudu Transforms for Solving Some,” vol. 8, no. 2, pp. 167–173, 2012. [29] T. M. Elzaki, “Solution of Nonlinear Differential Equations Using Mixture of Elzaki Transform and Differential Transform Method,” vol. 7, no. 13, pp. 631–638, 2012. [30] D. Ziane and M. H. Cherif, “Resolution of Nonlinear Partial Differential Equations by Elzaki Transform Decomposition Method Laboratory of mathematics and its applications (LAMAP),” vol. 5, pp. 17–30, 2015. [31] O. E. Ige, R. A. Oderinu, and T. M. Elzaki, “Adomian polynomial and Elzaki transform method for solving sine-gordon equations,” IAENG Int. J. Appl. Math., vol. 49, no. 3, pp. 1–7, 2019. [32] Q. Branch, N. Branch, and M. Benchohra, “Applications of homotopy perturbation method and Elzaki transform for solving nonlinear partial differential equations of fractional order,” vol. 2015, no. 6, pp. 91–104, 2016. [33] N. Shawagfeh, “Decomposition method for fractional partial differential equations,” no. December, 2017, doi: 10.5829/idosi.wasj.2019.18.24. [34] A. Prakash and V. Verma, “Numerical method for fractional model of newell- whitehead-segel equation,” Front. Phys., vol. 7, no. FEB, pp. 1–10, 2019, doi: 10.3389/fphy.2019.00015. [35] O. Abu Arqub and M. Al-Smadi, “Numerical algorithm for solving time-fractional partial integrodifferential equations subject to initial and Dirichlet boundary conditions,” Numer. Methods Partial Differ. Equ., vol. 34, no. 5, pp. 1577–1597, 2018, doi: 10.1002/num.22209.