EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 1, Article Number 5793 ISSN 1307-5543 – ejpam.com Published by New York Business Global Exploring Fractional-Order Models in Computational Finance via an Efficient Hybrid Approach Imtiaz Ahmad1,∗, Rashid Jan2,3, Normy Norfiza Abdul Razak2, Aziz Khan4, Thabet Abdeljawad4,5,6,7,∗ 1 Institute of Informatics and Computing in Energy (IICE), Universiti Tenaga Nasional (UNITEN), Kajang 43000, Selangor, Malaysia. 2 Institute of Energy Infrastructure (IEI), Department of Civil Engineering, College of Engineering, Universiti Tenaga Nasional (UNITEN), Putrajaya Campus, Jalan IKRAM- UNITEN, 43000 Kajang, Selangor, Malaysia. 3 Mathematics Research Center, Near East University TRNC, Mersin 10, Nicosia 99138, Turkey. 4 Department of Mathematics and Sciences, Prince Sultan University, P.O. Box 66833, 11586 Riyadh, Saudi Arabia. 5 Department of Mathematics, Saveetha School of Engineering, Saveetha Institute of Medical and Technical Sciences, Saveetha University, Chennai 602105, Tamil Nadu, India. 6 Center for Applied Mathematics and Bioinformatics (CAMB), Gulf University for Science and Technology, Hawally, 32093, Kuwait. 7 Department of Mathematics and Applied Mathematics, Sefako Makgatho Health Sciences University, Garankuwa, Medusa 0204, South Africa. Abstract. In this study, a hybrid numerical method is applied to solve the time-fractional Black- Scholes model for various options, including traditional (European and American) as well as non- standard options (such as butterfly spread, double barrier, and digital options). The method combines the fractional Liouville-Caputo scheme for time derivatives with the Strang splitting al- gorithm, while a meshless approach based on Lucas and Fibonacci polynomials is used for spatial derivatives. Numerical experiments are conducted using the L∞ error norm to evaluate the ac- curacy and effectiveness of the method, with the double mesh technique employed for validation when exact solutions are unavailable. In addition, the model performance is further evaluated through the computation of key sensitivity measures, specifically the Greeks (delta and gamma). The accuracy and robustness of the proposed solution are validated by benchmarking its results against those obtained using alternative methods reported in the literature. 2020 Mathematics Subject Classifications: 65M70, 91G60, 33C45, 35R11, 65D99 Key Words and Phrases: Hybrid numerical method, Fractional PDEs, Vanilla options, Exotic options, Strang splitting algorithm ∗Corresponding author. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5793 Email addresses: imtiazkakakhil@gmail.com (I. Ahmad), tabdeljawad@psu.edu.sa (T. Abdeljawad) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 2 of 27 1. Introduction Mathematical modeling and simulation of financial derivatives have become increas- ingly popular in recent years. When used correctly, financial derivatives help investors optimize gains while reducing losses. These instruments serve two primary purposes: they allow one party to hedge against risk, while offering the other party the potential for substantial returns [44]. Derivatives are used to manage significant risks associated with various financial assets, including fluctuations in the prices of commodities, stocks, bonds, indices, and changes in exchange rates. Common examples of financial derivatives in- clude European, American, and Asian options [44]. This study focuses on vanilla options, including European and American options, as well as exotic options such as butterfly spreads, double barriers, and digital options. Exotic options are frequently used in finan- cial contracts due to their distinct characteristics and, unlike standardized options, are customized to meet the specific requirements of investors. A European option, which can only be exercised on its expiration date, is commonly referred to as a vanilla option. In contrast, an American option offers more flexibility to investors, as it can be exercised at any point before or on the expiration date. However, this flexibility adds a computational challenge, as it requires determining both the option’s value and the optimal exercise point at each time step. The American option model is reformulated as a moving boundary problem. To manage the variable boundary inherent in the American put option model, a penalty source term approach is employed. This method introduces a minor, continuous penalty term [40] into the Black-Scholes equa- tion [26], thereby removing the free boundary and reformulating the problem as a fixed boundary value problem. The Black-Scholes (B-S) model represents an innovative method, categorized as a convection-diffusion partial differential equation (PDE). This classifica- tion arises from the assumption within the risk-neutral probability framework that stock prices follow a Brownian motion. While some alternative models offer closed-form solu- tions, many realistic scenarios require numerical methods due to their complex dynamics. Thus, developing efficient and accurate numerical techniques is crucial for effectively eval- uating these PDE models in option pricing. Various numerical techniques have been documented in recent publications, including finite element methods [51] and finite difference methods [11]. These approaches have demonstrated considerable success in computational finance, typically generating results with a high degree of accuracy. A recent study [11] highlighted an effective extrapolation method for one-dimensional European and digital options that achieved improved preci- sion. Nonetheless, radial basis function (RBF) based numerical algorithms present strong alternatives due to their spectral accuracy, straightforward coding, adaptability for both uniform and scattered nodes, and ease of implementation in higher-dimensional spaces. Various researchers have explored meshless methods to accurately evaluate PDE models in computational finance [1, 12, 29]. For instance, in [42], the authors introduced the radial basis point interpolation technique for solving the B-S model related to American and European options with a single asset, proposing various numerical methods. Similarly, other meshless methods [2, 25, 45, 48] have been developed to efficiently address American Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 3 of 27 and other option models. Solving partial differential equations (PDEs) can be quite challenging due to their complexity, which often renders traditional methods ineffective. As a result, there is a growing demand for efficient computational techniques that provide accurate approxima- tions. Among various techniques, meshless methods have demonstrated effectiveness and accuracy, making them suitable for a diverse array of both fractional and non-fractional PDEs. Several meshless approaches are available, including the moving least squares method, the natural element method, the reproducing kernel particle method, the mesh- less RBF methods, and meshless methods based on Lucas and Fibonacci polynomials. The primary objective of this study is to implement a hybrid method that utilizes the connection between Fibonacci and Lucas polynomials to numerically solve the proposed fractional model. The method integrates the Strang splitting technique with the Liouville- Caputo fractional derivative framework to address the temporal aspect. This approach fa- cilitates the application of higher-order derivatives, improves accuracy even with a coarser number of collocation points, and ultimately reduces computational costs. Additionally, these polynomials are known for their beneficial applications in differential equations (DEs). Previous research has shown that effectively tackling boundary value problems necessitates an understanding of the relationship between Lucas and Chebyshev polyno- mials. For example, the Lucas polynomials have proven useful in solving higher-order DEs [19]. Additionally, studies indicate that Fibonacci polynomials can be employed to solve Volterra-Fredholm integral DEs effectively [39]. Moreover, a hybrid method that integrates Taylor and Lucas polynomials for solving delay difference equations has also been proposed [13]. The integration of hybrid Fibonacci and Lucas polynomial algorithms has enabled innovative solutions for time-dependent PDEs [41]. Furthermore, the integra- tion of Lucas polynomials with finite difference techniques has yielded efficient numerical schemes for diverse PDE models [4, 6]. 1.1. Motivation This study focuses on developing a relatively new numerical method for solving the time-fractional Black-Scholes model (TFBSM), applicable to various financial derivatives, including digital options, butterfly spreads, double barriers, and traditional options (both European and American). The proposed solution integrates the Strang splitting technique with the Liouville-Caputo fractional derivative framework to handle the temporal aspect. For spatial derivatives, a meshless method using Fibonacci and Lucas polynomials is em- ployed. Given the complexity of fractional PDE models, accurate and efficient numerical methods are essential, as they often hinder the availability of closed-form analytical so- lutions. This study develops a robust numerical framework for managing the TFBSM applied to complex option structures. Additionally, the non-orthogonality of Lucas and Fibonacci polynomials, unlike orthogonal polynomials such as Chebyshev, removes the necessity for interval transformations. This simplification enhances the approximation of higher-order derivatives for unknown functions. As a result, the method boosts com- putational accuracy, even with coarser discretizations, and offers a stable approach to Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 4 of 27 addressing the complexities inherent in fractional PDEs. The structure of this paper is organized as follows: Section 2 introduces essential concepts and foundational principles in fractional calculus and polynomial theory. Section 3 presents an overview and literature review of the Black-Scholes model. Section 4 details the proposed numerical methodology. In Section 5, the numerical results are presented and analyzed. Lastly, Section 6 provides concluding remarks. 2. Basic Concepts in Polynomial Theory and Fractional Calculus Fractional derivatives are essential operators in fractional calculus [9, 36, 37], with applications spanning various fields such as physics, biology, magnetohydrodynamics, en- gineering, and finance. This mathematical framework enables the effective modeling of intricate processes exhibiting fractal characteristics, memory effects, and non-local inter- actions. Below are some widely utilized definitions of fractional derivatives. Definition 1: The Caputo’s fractional derivative [15]: ∂ηP(A, T ) ∂T η = 1 Γ (1− η) T∫ 0 ∂P(A, ζ) ∂ζ (T − ζ)−η dζ, 0 < η < 1. (1) Definition 2: The Riemann-Liouville derivative [32]: ∂ηP(A, T ) ∂T η = 1 Γ (1− η) d dT T∫ T (P(A, ϑ)− P(A, T )) (ϑ− T )η dϑ, 0 < η < 1. (2) 2.1. Lucas and Fibonacci Polynomial Theory This section is dedicated to defining and applying Lucas and Fibonacci polynomials for the approximation of unknown functions and their derivatives. Lucas polynomials [10]: Lucas polynomials can be expressed using a three-term re- currence relation: Lk(A) = kLk−1(A) + Lk−2(A), k ≥ 2, (3) with the initial values L0(A) = 2 and L1(A) = A. By allowing Eq. (3) to generate a sequence of Lucas numbers A = 1. Fibonacci polynomials [10]: Fibonacci polynomials, which are an extension of Fibonacci numbers, are defined using a three-term recurrence relation: Fk(A) = kFk−1(A) + Fk−2(A), k ≥ 2. (4) Starting with F0(A) = 0 and F1(A) = 1, the sequence proceeds. Equation (4) yields the well-known Fibonacci number sequence for A = 1. Lemma[10]: The derivative of order m of the kth Lucas polynomial, represented by Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 5 of 27 Lk(A), can be formulated using the kth Fibonacci polynomial, Fk(A), in the following way: L (m) k (A) = kFk(A)Dm−1, Dm−1 = D ×D ×D · · · D︸ ︷︷ ︸ (m−1)time . (5) In this context, D denotes the (M + 1)× (M + 1) matrix, which is written as follows: D =  0 0 . . . 0 0 ... d 0  , where the computation of d follows [10]: dij = { i sin (j−i)π 2 , if j > i, 0, otherwise. 2.2. Approximation of Function Let P(A) be a continuous function, and suppose that P ∈ L2(R). In this case, P can be represented as the following linear combination of the kth Lucas polynomials: P(A) = ∞∑ k=0 ΛkLk(A). (6) Here, Λk refers to the unknown coefficients, while Lk(A) signifies the Lucas polynomials. Likewise, a linear combination of the kth Fibonacci polynomials allows for the expansion of P(A) under the same conditions, as shown below: P(A) = ∞∑ k=0 ΛkFk(A). Here, the Fibonacci polynomials are represented as Fk(A), and the unknown coefficients are indicated by Λk. The first-order derivative of the function P(A) can be obtained by expanding it using the Lucas polynomial series. P ′(A) = ∞∑ k=0 ΛkL ′ k(A). (7) For the function P(A), the associated mth order derivative is as follows: Pm(A) = ∞∑ k=0 ΛkL (m) k (A), (8) Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 6 of 27 where Pm(A) = dmP(A) dAm , L (m) k (A) = dmLk(A) dAm . The Eqs. (7) and (8) can be expressed using the relationship defined in (5): P ′(A) = ∞∑ k=0 ΛkkFk(A), (9) In the same way, the mth derivative can be found using the formula below: P(m)(A) = ∞∑ k=0 ΛkkFk(A)Dm−1, (10) in which D and Dm−1 are previously defined. Remark 1: Truncated Lacus and Fibonacci polynomial series are commonly utilized in numerical computations to express P(A) and its mth derivative. Specifically, we examine the following: P(A) ≃ M∑ k=0 ΛkLk(A), M ∈ N, and P(m)(A) ≃ M∑ k=0 ΛkL (m) k (A) = M∑ k=0 ΛkkFk(A)Dm−1, M ∈ N. (11) 3. Fractional Black-Scholes Models Fractional derivatives serve as an effective tool for capturing the memory and hered- itary characteristics of various materials [8, 30, 31]. These derivatives are particularly relevant when examining differential equations within the framework of fractal geometry and dynamics [3, 5, 7, 14, 43]. The discovery of fractal structures in finance has led to a growing interest in fractional models. As a result, fractional-order PDEs have been incorporated into financial modeling to better represent complex financial dynamics by replacing classical Brownian motion with fractional Brownian motion in stochastic pro- cesses. To effectively depict memory effects, one must comprehend the non-local character of integrals and fractional derivatives [24]. In their work, Jumarie [33] introduced time- and space-fractional B-S models utilizing a fractional Taylor series expansion approach. Farhadi et al. [24] employed a TFBSM to account for the influence of memory effects on trend dynamics, while Liang et al. [38] developed bifractional B-S models tailored to op- tion pricing applications. Further research in this area remains essential and valuable, as fractional B-S partial differential equations continue to pose increasing challenges in terms of solution complexity. Various numerical methods have been applied to price European options using fractional B-S models. Yavuz et al. [49] utilized the Laplace homotopy Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 7 of 27 analysis method for pricing options. In their work, Cen et al. [17] applied an integral ap- proach for time discretization alongside central differencing on a piecewise uniform mesh for spatial discretization. In a finite-moment log-stable framework, Chen [22] solved fractional B-S model in term of the American option. Using a power penalty technique in conjunction with a finite- difference scheme. For both hazardous and risk-free assets, the conventional B-S model has proven successful in pricing options in a full market with no transaction costs. But when transaction costs are present, this strategy falls short because complete hedging is unachievable [23]. Edeki et al. [46] used modified differential transform method for the solution of the nonlinear TFBSM. Furthermore, Chen et al. [20] implemented the penalty method to address fractional-order American put options. A predictor-corrector technique specifically designed for American options was presented by Chen et al. [22]. Zhou et al. [52] used a combination of the finite difference method and the Laplace transform technique to solve American options. Cartea et al. [16] addressed exotic option pricing using the shifted Grűnwald-Letnikov scheme, while Chen et al. [21] provided a closed-form analytical solution for double barrier options, highlighting their significance, applications, and advantages. These options are typically less expensive than standard vanilla options but still allow holders to benefit from fluctuations in the spot price, whether upward or downward. The payout to the holder depends on the underlying asset breaching two barriers, Bl and Br. For pricing double barrier options, Zhang et al. [50] introduced an implicit discrete scheme, Golbabai et al. [27] proposed the moving least-squares technique, and [28] presented a meshless radial basis function method. In model (12), it is assumed that the underlying asset adheres to a conventional ge- ometric Brownian motion, while the fluctuations in option prices are modeled based on fractal transmission patterns [21, 50]. The boundary and terminal conditions for the time-fractional Black-Scholes model used in option pricing are defined as follows [21, 50]: ∂ηP ∂tη = −1 2 ς2A2 ∂ 2P ∂A2 − ϵA∂P ∂A + rP, 0 < η ≤ 1, 0 ≤ t ≤ T,A ≥ 0, P(A, T ) = F(A), P(0, t) = F1(t), P(+∞, t) = F2(t). (12) The payoff function serves as the initial condition for an initial value problem that reformu- lates model (12) through a linear transformation T = T − t of the variable. Consequently, this leads to a redefinition of the boundary and initial conditions for the fractional PDE model. ∂ηP ∂T η = 1 2 ς2A2 ∂ 2P ∂A2 + ϵA∂P ∂A − rP, P(A, 0) = F(A), P(0, T ) = F1(T ), P(+∞, T ) = F2(T ). (13) When η = 1, Eq. (13) simplifies to the standard Black-Scholes model. In this scenario, we implement the proposed method to derive the numerical solution by examining three Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 8 of 27 models. For European options, we utilize the TFBSM [18]: ∂ηP ∂T η − 1 2 ς2A2 ∂ 2P ∂A2 − ϵA∂P ∂A + ϵP = 0, (A, T ) ∈ (0,Amax)× (0, T ), P(A, 0) = F(A), P(0, T ) = F1(T ), P(Amax, T ) = F2(T ). (14) For American put options with penalty term approach [35] is given as follow using TFBSM. This formulation extends the model presented in [25] from a fixed domain to incorporate time-fractional orders. ∂ηP ∂T η = 1 2 ς2A2 ∂ 2P ∂A2 + ϵA∂P ∂A − ϵP + µC P + µ− q(A) , (A, T ) ∈ (0,Amax)× (0, T ), P(A, 0) = max(E − A, 0), P(0, T ) = E , lim A→∞ P(A, T ) = 0. (15) The initial stock price is denoted as A0, while the current stock price is represented by A, with the condition that 0 < η ≤ 1. The volatility of returns is indicated by ς(≥ 0), and ϵ signifies the risk-free interest rate where E is exercise price, T refers to the time until expiration, and ρ represents market liquidity. The option price is expressed as P(A, T ), with the additional conditions that C ≥ ϵE , 0 < µ ≪ 1, and q(A) = E − A. 4. Proposed Methodology The hybrid meshless method proposed in this section aims to approximate the two- term time-fractional model Eqs. (14)-(15). To facilitate the discussion, we provide the following notation: Pn+1(A) = P(A, T n+1), Pn+1 i = P(Ai, T n+1), here the collocation points are denoted by Ai and the time-step size is represented by ∆T , and T n = n×∆T : Ai = a+ ihA (i = 1, 2, · · · ,M, M ∈ N). In the spatial domain hA = (b−a)/M , the spatial step sizes are denoted by hA = (b−a)/M . 4.1. Temporal Discretization An approximation error of order O(∆T 2−η) is obtained for the discrete formulation at the (n+ 1)th time level by employing the established L1 formula, with the condition that 0 < η ≤ 1 [4]. ∂ηP(A, T n+1) ∂T η = 1 Γ(1− η) ∫ T n+1 0 ∂P(A, ζ) ∂ζ (T n+1 − ζ)−ηdζ, Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 9 of 27 = 1 Γ(1− η) n∑ k=0 ∫ (k+1)×∆T k×∆T ∂P(A, ζ) ∂ζ (T n+1 − ζ)−ηdζ, = 1 Γ(1− η) n∑ k=0 [ Pk+1A−PkA ∆T +O(∆T ) ] ∫ (k+1)×∆T k×∆T ((k + 1)∆T − ζ)−ηdζ. Following integration, it gives: ∂ηP(A, T n+1) ∂T η =  Aη ∑n k=0Kη(k) [ Pn−k+1A−Pn−kA ] +O(∆T 2−η), 0 < η < 1, Pn+1A−PnA ∆T +O(∆T ), η = 1, (16) where Aη = ∆T −η Γ(2−η) and Kη(k) = (k + 1)1−η − (k)1−η. Thus, by omitting the error term, we can write: ∂ηP(A, T n+1) ∂T η = Aη [ Pn+1A−PnA ] +Aη n∑ k=1 Kη(k) [ Pn−k+1A−Pn−kA ] , (17) with Kη(k) = 1, k = 0. Let’s rewrite Eq. (14) prior to implementing the θ-weighted scheme: ∂ηP(A, T ) ∂T η = LP(A, T ), A ∈ Ω ⊂ R, 0 < η ≤ 1, T > 0, (18) with the conditions P(A, 0) = P0A, P(A, T ) = f1(A, T ), A ∈ ∂Ω, (19) where L = 1 2 ς 2A2 ∂2P ∂A2 + ϵA∂P ∂A − ϵP. We now apply the θ-weighted rule, which converts Eq. (18) to: ∂ηP(A, T ) ∂T η = θLPn+1(A, T ) + (1− θ)LPn(A, T ). (20) Using Eq. (16), we get: AηPn+1(A)− θLPn+1(A) = AηPn(A) + (1− θ)LPn(A)−Qn η (A), (21) where Qn η (A) = Aη n∑ k=1 Kη(k) [ Pn−k+1(A)− Pn−k(A) ] . Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 10 of 27 Thus far, Eq. (21) serves as the time-discrete representation of Eq. (18). Furthermore, the Strang splitting algorithm, a well-established technique for decom- position, is applied in the time dimension and is formulated as follows [44, 47]: ∂P1 ∂T = L1P1(T ) with T ∈ [T n, T n+1/2] and P1(T n) = Pn sp, (22) ∂P2 ∂T = L2P2(T ) with T ∈ [T n, T n+1] and P2(T n) = P1(T n+1/2), (23) ∂P3 ∂T = L1P3(T ) with T ∈ [T n+1/2, T n+1] and P3(T n+1/2) = P2(T n+1), (24) where L1P(T ) = 1 2 ς 2A2 ∂2P ∂A2 , L2P(T ) = rA∂P ∂A − rP, and Pn sp = P0 is the initial solution defined for the Eq. (18). 4.2. Space Discretization An alternative way to express the function P(A) is through a truncated series expan- sion involving Lucas polynomials. P(A) ≃ M∑ k=0 ΛkLk(A). (25) Equation (25) can be write as: P(A) ≃ LT(A)Λ, (26) with Λ = [Λ0,Λ1,Λ2, · · · ,ΛM ]T and LT(A) = [L0(A), L1(A), L2(A), · · · , LM (A)]. Partial derivatives of P(A) can be calculates as: ∂P(A) ∂A ≃ M∑ k=0 Λk dLk(A) dA = LT A(A)Λ, and ∂2P(A) ∂A2 ≃ M∑ k=0 Λk d2Lk(A) dA2 = LT AA(A)Λ, where LT A(A) = { ∂ ∂ALk(A)}Mk=0 = {dLk(A) dA }Mk=0 and LT AA(A) = { ∂2 ∂A2Lkm(A)}Mk=0 = {d2Lk(A) dA2 }Mk=0. One can refined further as: ∂P ∂A ≃ ( LT A(A) ) Λ. (27) ∂2P ∂A2 ≃ ( LT AA(A) ) Λ. (28) Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 11 of 27 The function P(A, T ) and its partial derivatives can be expressed as follows: P(A, T ) ≈ LT(A)Λ(T ), ∂P(A, T ) ∂A ≈ LT A(A)Λ(T ), ∂2P(A, T ) ∂A2 ≈ LT AA(A)Λ(T ). (29) Let Λ(T ) represent the vector of unknown coefficients that vary with time. By combining the results from Eq. (29) and Eq. (21), the following expression is obtained: (Aη)L(A)Λn+1 − θ ( 1 2 ς2A2LT AA(A) + ϵALT A(A)− ϵLT(A) ) Λn+1 (30) = (Aη)L(A)Λn + (1− θ) ( 1 2 ς2A2LT AA(A) + ϵALT A(A)− ϵLT(A) ) Λn −Qn η (A). The boundary conditions: BL(A)Λn+1 ≈ f1(A, T n+1) = fn+1 1 (A), ∀A ∈ ∂Ω, (31) where Λn = Λ(T n). 4.3. Full Discretization We utilize collocation on Eqs. (30) and (31) at the discrete mesh points Ai associated with the corresponding time levels to formulate the discrete representation of the model given in Eq. (18). As a result, we obtain the following: (Aη)L(Ai)Λ n+1 − θ ( 1 2 ς2A2 iL T AA(Ai) + ϵALT A(Ai)− ϵLT(Ai) ) Λn+1 = (Aη)L(Ai)Λ n + (1− θ) ( 1 2 ς2A2 iL T AA(Ai) + ϵALT A(Ai)− ϵLT(Ai) ) Λn −Qn η (Ai), and BL(Ai)Λ n+1 ≈ f1(Ai, T n+1) = fn+1 1 (Ai), ∀Ai ∈ ∂Ω, where i = 0, 1, 2, · · · ,M . In matrix-vector form: GΛn+1 = HΛn +Qn+1. (32) The entries of the matrix in this case are defined as follows: (G)i = { (Aη)L(Ai)− θ ( 1 2 ς 2A2 iL T AA(Ai) + ϵALT A(Ai)− ϵLT(Ai) ) , Ai ∈ Ω, BL(Ai), Ai ∈ ∂Ω, (33) (H)i = { (Aη)L(Ai) + (1− θ) ( 1 2 ς 2A2 iL T AA(Ai) + ϵALT A(Ai)− ϵLT(Ai) ) , Ai ∈ Ω, 0, Ai ∈ ∂Ω, (34) Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 12 of 27 and (Q)i = { −Qn η (Ai), Ai ∈ Ω, f1(Ai, T n+1), Ai ∈ ∂Ω. (35) The methodological framework for Eq. (18) has been completed up to this stage. Accord- ing to [41], it is assured that the linear system defined by Eq. (32) possesses a solution. At each time step, we can determine the vector of unknown coefficients by solving Eq. (32). The procedure is initiated using the first condition as follows: LT(Ai)Λ 0 = P(Ai, 0) = P0 =⇒ AΛ0 = P0. Upon solving, we obtain Λ0. Equation (32) demonstrates that this iterative procedure continues until the desired time level is achieved. Subsequently, the coefficient vector is computed to determine the solution at the specified time step. P(Ai, t n) = LT(Ai)Λ n, ∀Ai ∈ Ω, (n ≥ 0). 5. Numerical Analysis of Computational Simulations This study assesses the proposed method in terms of accuracy, convergence order, and applicability. The evaluation includes both vanilla and exotic options. We employed the maximum error (L∞) norm to measure accuracy. L∞ = max|Pex − P|, where Pex and P denote the precise and numerical values respectively. Consider the following model equation [50]: ∂ηP ∂T η = 1 2 ς2A2 ∂ 2P ∂A2 + ϵA∂P ∂A − ϵP + f(A, T ), (A, T ) ∈ (0,Amax)× (0, T ), P(0, T ) = 0, P(Amax, T ) = 0, P(A, 0) = A2(1−A). (36) We define f as follows: f = ( 2T 2−η Γ(3− η) + 2T 1−η Γ(2− η) ) A2(1−A)− (T +1)2[x(2−6A)+y(2A−3A2)− zA2(1−A)]. This expression is selected so that P = (T + 1)2A2(1 − A) serves as the exact solution to Eq. (36). The parameters for this scenario are defined as follows: ς = 0.25, ϵ = 0.05, Amax = 1, x = 1 2 ς 2, y = ϵ− a, z = ϵ, and T = 1. The convergence rates and slopes of the Problem 5 are illustrated in Figure 1, where they are compared to the implicit discrete scheme (IDS) [50] and the local RBF collocation Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 13 of 27 method (LRBFCM) [1]. The figure demonstrates that the results obtained from the pro- posed method exhibit greater accuracy than those from the other methods. Additionally, the numerical solution for nodal points M = 64, fractional-order η = 0.7, and η = 1 is presented in Figure 2. The numerical outcomes obtained from the proposed methodology exhibit a high degree of concordance with the exact solutions. 10 −2 10 −1 10 0 10 −6 10 −5 10 −4 10 −3 10 −2 Problem 1 log(h) lo g( L ∞ ) LRBFCM IDS Proposed Method Figure 1: Comparative analysis of the proposed method with the LRBFCM [1] and IDS [50] in terms of L∞ for various spatial steps h in Problem 5. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 A P (A , T ) Exact solution Numerical solution 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 A P (A , T ) Exact solution Numerical solution Figure 2: Numerical and exact solution for η = 0.7 (left) and η = 1 (right) for Problem 5. Initially, examine the time-fractional model described by Eq. (14), focusing on the butterfly spread call option. This option is defined by three strike prices: S1, S2 and S3 where S1 < S2 < S3 and S2 = S1+S2 2 . ∂ηP ∂T η = 1 2 ς2A2 ∂ 2P ∂A2 + ϵA∂P ∂A − ϵP, (A, T ) ∈ (0, Br)× (0, T ), P(0, T ) = 0, P(Br, T ) = 0, P(A, 0) = max(A− S1, 0)− 2max(A− S2, 0) + max(A− S3, 0), (37) Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 14 of 27 Table 1: Comparative results of the proposed approach with the LMCM [29] are presented using L∞ norms for Problem 5. η Method M = 61 M = 121 M = 241 M = 481 M = 961 0.3 LMCM [29] 1.3954e−04 3.5386e−05 8.8776e−06 2.3228e−06 1.0259e−06 Proposed method 7.8320e−05 9.1639e−06 5.4617e−06 9.5813e−07 8.4187e−07 0.5 LMCM [29] 1.2142e−04 3.0580e−05 7.6584e−06 2.0221e−06 9.7377e−07 Proposed method 7.8355e−05 9.6832e−06 4.8274e−06 9.2025e−07 6.4002e−07 0.7 LMCM [29] 8.7415e−05 2.1697e−05 5.4133e−06 1.4676e−06 8.7138e−07 Proposed method 4.5277e−05 8.1539e−06 1.7841e−06 8.4630e−07 6.2771e−07 0.9 LMCM [29] 3.0188e−05 7.5350e−06 1.8837e−06 5.6087e−07 7.0604e−07 Proposed method 1.0605e−05 4.5635e−06 8.2638e−07 4.0256e−07 5.8730e−07 where strike prices S1 = 0.4, S2 = 0.5, S3 = 0.6, ϵ = 0.1, ς = 0.5 and T = 1 (year) are the parameters values according to [44]. Using a double mesh approach for Problem 5, the findings shown in Table 1 contrast the suggested approach with the LMCM [29], emphasizing the maximum error at different η values. The accuracy of the proposed method improves as the number of nodal points increases, producing results that surpass those reported in [29]. Figure 3 illustrates the outcomes of the proposed approach for different fractional order values of η. It indicates that when the asset price A is near the strike prices S1 and S3, the model tends to yield lower option prices; conversely, higher prices are observed when A is close to the strike price S2. This suggests that the time-fractional Black-Scholes process provides a more realistic portrayal of the characteristics associated with significant price movements compared to the traditional integer-order Black-Scholes model. Additionally, Figure 4 displays the solution of the butterfly spread option in term of asset (A), time T versus option pricing value for fractional orders η = 0.3 and η = 0.7. Furthermore, we compute the gamma of the options portfolio, which reflects how the delta of the portfolio changes in response to fluctuations in the underlying asset price, along with the delta of an option itself, which measures the sensitivity of the option’s value to shifts in the asset price. The numerical values of delta and gamma for the butterfly spread option across a range of η values are shown in Figure 5. Examine the model represented by Eq. (14), which involves a digital call or cash-or- nothing call option. ∂ηP ∂T η = 1 2 ς2A2 ∂ 2P ∂A2 + ϵA∂P ∂A − ϵP, (A, T ) ∈ (0, Bϵ)× (0, T ), (38) with payoff function boundary conditions: P(A, T ) =  0 for A < S, 0.5 for A = S, 1 for A > S. (39) Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 15 of 27 0 0.2 0.4 0.6 0.8 1 0 0.005 0.01 0.015 0.02 0.025 A O p ti o n v a lu e η = 1 η = 0.7 η = 0.5 η = 0.3 Figure 3: Numerical solution of the fractional-order η butterfly spread option for Problem 5. 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.02 0.04 0.06 0.08 AT O p ti o n p ri ce 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.02 0.04 0.06 0.08 AT O p ti o n p ri ce Figure 4: Numerical solution of the butterfly spread option for fractional orders η = 0.3 (left) and η = 0.7 (right) in Problem 5. V (A, 0) = { 0 for A = 0, e−r(T−T ) for A → ∞, (40) where the parameters ϵ = 0.05, S = 0.5, ς = 0.2, and T = 1 (year) [44]. The results obtained from applying the double mesh method to Problem 5 and comparing them with those from the LMCM [29] are presented in Table 2 for different fractional orders η. The data indicate that the proposed method outperforms the LMCM [29], as reflected in the table. Furthermore, the accuracy improves with an increasing number of nodes. Figure 6 displays the option prices at various η values. Additionally, Figure 7 illustrates the three-dimensional profile of the digital call option for fractional orders η = 0.3 and η = 0.7. Lastly, Figure 8 presents the delta and gamma of the digital call option across different η values. Examine the equation presented in Eq. (14) pertaining Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 16 of 27 0 0.2 0.4 0.6 0.8 1 −0.04 −0.02 0 0.02 0.04 0.06 A D el ta (∆ ) η = 0.5 η = 0.7 η = 0.9 η = 1 0 0.2 0.4 0.6 0.8 1 −3 −2 −1 0 1 A G a m m a (Γ ) η = 1 η = 0.9 η = 0.7 η = 0.5 Figure 5: Numerical solution of ∆ and Γ for the butterfly spread option with fractional orders η in Problem 5. Table 2: A comparison of the proposed approach with the LMCM [29] is made using L∞ norms for Problem 5. η Method M = 61 M = 121 M = 241 M = 481 M = 961 0.3 LMCM [29] 1.0159e−03 2.5886e−04 6.5008e−05 1.6627e−05 7.1306e−06 Proposed method 7.6382e−04 9.0435e−05 4.0917e−05 7.3491e−06 4.5375e−06 0.5 LMCM [29] 8.4988e−04 2.1518e−04 5.3940e−05 1.3876e−05 6.7375e−06 Proposed method 3.5562e−04 8.5784e−05 2.5383e−05 6.3471e−06 2.3566e−06 0.7 LMCM [29] 6.1047e−04 1.5205e−04 3.7943e−05 9.9106e−06 6.2241e−06 Proposed method 3.4860e−04 6.9042e−05 1.1684e−05 7.1095e−06 2.4261e−06 0.9 LMCM [29] 3.9657e−04 1.0060e−04 2.5239e−05 6.3455e−06 5.4104e−06 Proposed method 9.9648e−05 6.7810e−05 1.1105e−05 3.5633e−06 1.7319e−06 to the double barrier knock-out call option model [50] ∂ηP ∂T η = 1 2 ς2A2 ∂ 2P ∂A2 + (ϵ−D)A∂P ∂A − ϵP, (A, T ) ∈ (Bl, Bϵ)× (0, T ), P(Bl, T ) = p1(T ), P(Bϵ, T ) = p2(T ), P(A, 0) = q(A), (41) where q(A) = max(A−S, 0), p1(T ) = p2(T ) = 0 and other parameter values are ς = 0.45, D = 0.01, ϵ = 0.03, Br = 15, Bl = 3, T = 1 (year) and S = 10 as given in [21, 50]. The outcomes are depicted in Figure 9 for various values of η. The model demonstrates that option prices decline when the asset price A is below the strike price S. In contrast, for in-the-money options, where A > S the prices are elevated. Furthermore, Figure 10 presents the numerical solution for the double barrier knock-out option at η = 0.3 and η = 0.7. The European call options is consider in the following equation using TFBSM Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 17 of 27 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 A O p ti o n v a lu e η = 1 η = 0.7 η = 0.5 η = 0.3 Figure 6: Numerical solution of the fractional-order η digital option for Test problem 5. 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 AT O p ti o n p ri ce 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 AT O p ti o n p ri ce Figure 7: Numerical solution of the digital option for fractional orders η = 0.3 (left) and η = 0.7 (right) in Problem 5. [18]: ∂ηP ∂T η = 1 2 ς2A2 ∂ 2P ∂A2 + ϵA∂P ∂A − ϵP, (A, T ) ∈ (0,Amax)× (0, T ), P(0, T ) = 0, P(Amax, T ) = Amax − Ee−rT , T ∈ (0, T ), P(A, 0) = max(A− E , 0). (42) with ς = 0.3, ϵ = 0.04, E = 10, T = 1, and Amax = 40. The outcomes produced by the proposed method are compared with those from the CDM [18] and LRBFCM [1], and they are evaluated using the L∞ error norm, as presented in Table 3. Given that there is no exact solution for Problem 5, the double mesh technique was employed to approximate the error norm. The results indicate that the proposed approach demonstrates superior accuracy compared to both CDM [18] and LRBFCM [1], as detailed in Table 3. Additionally, the numerical solutions obtained from the proposed Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 18 of 27 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 A D el ta (∆ ) η = 1 η = 0.9 η = 0.7 η = 0.5 0 0.2 0.4 0.6 0.8 1 −80 −60 −40 −20 0 20 40 A G a m m a (Γ ) η = 1 η = 0.9 η = 0.7 η=0.5 Figure 8: Numerical solution of ∆ and Γ for the digital option with fractional orders η in Problem 5. method for various fractional orders are illustrated in Figure 11. Specifically, Figure 12 depicts the numerical solution for a European call option with η = 0.7 and M = 41. Table 3: A comparison of the proposed approach with the methods in reported in [1, 18] is made using L∞ norms for Problem 5. M = 64 M = 128 η Proposed Method CDM[18] LRBFCM [1] Proposed Method CDM[18] LRBFCM [1] 0.1 9.7328e−4 1.2029e−3 3.7775e−3 1.8301e−4 2.9366e−4 9.4820e−4 0.3 9.4581e−4 1.1477e−3 3.6010e−3 1.5016e−4 2.7890e−4 9.0293e−4 0.5 8.5899e−4 1.0952e−3 3.4144e−3 1.0164e−4 2.6816e−4 8.5492e−4 0.7 8.1536e−4 1.0831e−3 3.2257e−3 9.6493e−5 2.6511e−4 8.0467e−4 0.9 7.9455e−4 1.1384e−3 3.0512e−3 8.8433e−5 2.7875e−4 7.7231e−4 Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 19 of 27 0 5 10 15 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 A O p ti o n P ri ce η=0.5 η=0.7 η=0.9 η=1 Figure 9: Numerical solution of the fractional-order η the double barrier option for Test problem 5. 0 5 10 15 0 0.5 1 0 2 4 6 8 10 AT O p ti o n p ri ce 0 5 10 15 0 0.5 1 0 2 4 6 8 10 AT O p ti o n p ri ce Figure 10: Numerical solution of the double barrier option for fractional orders η = 0.3 (left) and η = 0.7 (right) in Problem 5. 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 0 10 20 30 40 TA O p ti o n p ri ce Figure 12: Numerical solution of the European option for M = 41 and η = 0.7 in Problem 5. Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 20 of 27 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 A O p ti o n p ri ce η=0.3 η=0.5 η=0.7 η=0.9 η=1 16.15 16.16 16.17 16.18 16.19 16.2 16.21 16.22 6.65 6.66 6.67 6.68 6.69 6.7 6.71 6.72 A O p ti o n p ri ce η=0.3 η=0.5 η=0.7 η=0.9 η=1 Figure 11: Numerical solution for various values of η for Problem 5. The American put option is taken into account using the following TFBSM, utilizing the penalty term method: ∂ηP ∂T η = 1 2 ς2A2 ∂ 2P ∂A2 + ϵA∂P ∂A − ϵP + µC P + µ− q(A) , (A, T ) ∈ (0,Amax)× (0, T ), P(0, T ) = E , P(Amax, T ) = 0, T ∈ (0, T ), P(A, 0) = max(E − A, 0). (43) with parameters E = 1, ϵ = 0.1, ς = 0.2, C = rE , µ = 0.01, q(A) = E − A, Amax = 2 and T = 1. The model equation presented in Eq. (43) does not have an exact solution. To evaluate accuracy, we applied the double mesh procedure. Table 4 presents a comparison for η = 1 between the proposed method and the reference price [25], along with other current methodologies [1, 25, 34]. Additionally, the numerical solutions for various values ofM and η are provided in Table 5. Figure 13 displays the matching slopes and convergence rates. Moreover, Figure 14 shows a surface plot for an American put option with parameters η = 0.5 and M = 21. Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 21 of 27 Table 4: A comparison of the proposed approach with a reference solution [25], obtained using very fine nodes, as well as with other methods, is presented for T = 1 and η = 1 values in Problem 5. Reference Proposed LRBFCM A price [25] method [1] [25] [34] 0.6 0.4000037 0.4000039 0.4000043 0.4000176 0.4000185 0.7 0.3001161 0.3001180 0.3001216 0.3001007 0.3002333 0.8 0.2020397 0.2020307 0.2020324 0.2019901 0.2022428 0.9 0.1169591 0.1169463 0.1168272 0.1165422 0.1154885 1.0 0.0602833 0.0602845 0.0601380 0.0597033 0.0580422 1.1 0.0293272 0.0293238 0.0292004 0.0287648 0.0276763 1.2 0.0140864 0.0140863 0.0140160 0.0136840 0.0132349 1.3 0.0070408 0.0070417 0.0070086 0.0068192 0.0066983 1.4 0.0038609 0.0038609 0.0038484 0.0037485 0.0037539 Table 5: Comparative results of the proposed approach with the LMCM [29] are presented using L∞ norms for Problem 5. η Method M = 21 M = 41 M = 81 M = 161 0.9 LRBFCM[1] 5.7315e−3 8.3826e−4 1.6169e−4 4.7204e−5 Proposed Method 8.7368e−4 2.3298e−4 6.5822e−5 1.0377e−5 0.7 LRBFCM[1] 5.8620e−3 1.1410e−3 2.3723e−4 2.0000e−4 Proposed Method 8.9310e−4 5.2644e−4 7.5681e−5 7.4367e−5 0.5 LRBFCM[1] 6.0068e−3 1.3576e−3 7.5219e−4 5.6710e−4 Proposed Method 9.0376e−4 5.5783e−4 8.2655e−5 9.6883e−5 0.3 LRBFCM[1] 2.3432e−3 2.5194e−3 1.5109e−3 6.7021e−4 Proposed Method 9.7833e−4 6.8201e−4 5.9836e−4 1.4356e−4 Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 22 of 27 10 −1 10 −6 10 −5 10 −4 10 −3 10 −2 log(h) lo g( L ∞ ) η=0.3 η=0.5 η=0.7 η=0.9 Figure 13: Numerical solution the proposed method in terms of L∞ for various spatial steps h in Problem 5. Figure 14: Numerical solution of the American option for M = 21 and η = 0.5 in Problem 5. 6. Conclusion This study presents an effective numerical method for solving the time-fractional Black- Scholes model for various options, including both vanilla options (like European and Amer- ican options) and exotic options (such as butterfly spreads, double barriers, and digital op- tions). The proposed hybrid meshless technique integrates the Strang splitting algorithm Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 23 of 27 with Lucas and Fibonacci polynomials to address fractional partial differential equation challenges. The method’s effectiveness is showcased through detailed tables and figures that provide an objective assessment of its performance across various criteria. Significant insights were gained during the numerical experiments. Tests were conducted with varying fractional parameter values η to highlight the benefits of the time-fractional Black-Scholes model in comparison to the traditional Black-Scholes model for both call and put op- tions. Furthermore, the proposed method demonstrated greater accuracy and efficiency than existing techniques reported in the literature. Consequently, this approach shows substantial potential for application in various option pricing models, as well as in other areas involving partial differential equations in finance and fields that utilize fractional derivatives. Declarations Competing interests: Authors declare no conflict of interest. Availability of data and materials: Data will be provided on request to the corresponding author. Acknowledgments and Funding: Aziz Khan and Thabet Abdeljawad would like to thank Prince Sultan University for pay- ing the APC and support through TAS research lab. References [1] Hijaz Ahmad, Muhammad Nawaz Khan, Imtiaz Ahmad, Mohamed Omri, and Maged F Alotaibi. A meshless method for numerical solutions of linear and non- linear time-fractional Black-Scholes models. AIMS Math, 8(8):19677–19698, 2023. [2] Imtiaz Ahmad, Hijaz Ahmad, Phatiphat Thounthong, Yu-Ming Chu, and Clemente Cesarano. Solution of multi-term time-fractional PDE models arising in mathematical biology and physics by local meshless method. Symmetry, 12(7):1195, 2020. [3] Imtiaz Ahmad, Ihteram Ali, Rashid Jan, Sahar Ahmed Idris, and Mohamed Mousa. Solutions of a three-dimensional multi-term fractional anomalous solute transport model for contamination in groundwater. Plos one, 18(12):e0294348, 2023. [4] Imtiaz Ahmad, Abdulrahman Obaid Alshammari, Rashid Jan, Normy Norfiza Ab- dul Razak, and Sahar Ahmed Idris. An efficient numerical solution of a multi- dimensional two-term fractional order PDE via a hybrid methodology: The Caputo– Lucas–Fibonacci approach with strang splitting. Fractal and Fractional, 8(6):364, 2024. [5] Imtiaz Ahmad, Asmidar Abu Bakar, Hijaz Ahmad, Aziz Khan, and Thabet Abdel- jawad. Investigating virus spread analysis in computer networks with Atangana– Baleanu fractional derivative models. Fractals, 2440043:17, 2024. Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 24 of 27 [6] Imtiaz Ahmad, Asmidar Abu Bakar, Ihteram Ali, Sirajul Haq, Salman Yussof, and Ali Hasan Ali. Computational analysis of time-fractional models in energy infrastruc- ture applications. Alexandria Engineering Journal, 82:426–436, 2023. [7] Imtiaz Ahmad, Ibrahim Mekawy, Muhammad Nawaz Khan, Rashid Jan, and Salah Boulaaras. Modeling anomalous transport in fractal porous media: A study of fractional diffusion PDEs using numerical method. Nonlinear Engineering, 13(1):20220366, 2024. [8] Imtiaz Ahmad, Sakhi Zaman, et al. Local meshless differential quadrature collocation method for time-fractional PDEs. Discrete & Continuous Dynamical Systems-Series S, 13(10), 2020. [9] Saim Ahmed, Ahmad Taher Azar, Mahmoud Abdel-Aty, Hasib Khan, and Jehad Alz- abut. A nonlinear system of hybrid fractional differential equations with application to fixed time sliding mode control for Leukemia therapy. Ain Shams Engineering Journal, 15(4):102566, 2024. [10] Ihteram Ali, Sirajul Haq, Kottakkaran Sooppy Nisar, and Dumitru Baleanu. An effi- cient numerical scheme based on Lucas polynomials for the study of multidimensional Burgers-type equations. Advances in Difference Equations, 2021(1):1–24, 2021. [11] L. V. Ballestra. Repeated spatial extrapolation: An extraordinarily efficient approach for option pricing. Journal of Computational and Applied Mathematics, 256:83–91, 2014. [12] L. V. Ballestra and G. Pacelli. Pricing European and American options with two stochastic factors: A highly efficient radial basis function approach. Journal of Eco- nomics Dynamics and Control, 37:1142–1167, 2013. [13] N Baykuş-Savaşaneril and Mehmet Sezer. Hybrid Taylor-Lucas collocation method for numerical solution of high-order Pantograph type delay differential equations with variables delays. Applied Mathematics and Information Sciences, 11(6):1795–1801, 2017. [14] Salah Boulaaras, Rashid Jan, Amin Khan, Ali Allahem, Imtiaz Ahmad, and Salma Bahramand. Modeling the dynamical behavior of the interaction of T-cells and human immunodeficiency virus with saturated incidence. Communications in Theoretical Physics, 76(3):035001, 2024. [15] Michele Caputo. Linear models of dissipation whose Q is almost frequency independent-II. Geophysical Journal International, 13(5):529—-539, 1967. [16] Alvaro Cartea and Diego del Castillo-Negrete. Fractional diffusion models of option prices in markets with jumps. Physica A: Statistical Mechanics and its Applications, 374(2):749–763, 2007. [17] Zhongdi Cen, Jian Huang, Aimin Xu, and Anbo Le. Numerical approximation of a time fractional Black-Scholes equation. Computers & Mathematics with Applications, 75(8):2874–2887, 2018. [18] Cen Zhongdi, Jian Huang, Aimin Xu, and Anbo Le. Numerical approximation of a time-fractional Black–Scholes equation. Computers & Mathematics with Applications, 75:2874–2887, 2018. [19] Muhammed Çetin, Mehmet Sezer, and Coşkun Güler. Lucas polynomial approach Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 25 of 27 for system of high-order linear differential equations and residual error estimation. Mathematical problems in engineering, 2015, 2015. [20] Wen Chen and Song Wang. A penalty method for a fractional order parabolic varia- tional inequality governing American put option valuation. Computers & Mathematics with Applications, 67(1):77–90, 2014. [21] Wenting Chen, Xiang Xu, and Song-Ping Zhu. Analytically pricing double barrier options based on a time-fractional Black-Scholes equation. Computers & Mathematics with Applications, 69(12):1407–1419, 2015. [22] Wenting Chen, Xiang Xu, and Song-Ping Zhu. A predictor–corrector approach for pricing American options under the finite moment log-stable model. Applied Numer- ical Mathematics, 97:15–29, 2015. [23] Rafael Company, Lucas Jódar, and José-Ramón Pintos. A numerical method for European option pricing with transaction costs nonlinear equation. Mathematical and Computer Modelling, 50(5-6):910–920, 2009. [24] A Farhadi, M Salehi, and GH Erjaee. A new version of Black-Scholes equation pre- sented by time-fractional derivative. Iranian Journal of Science and Technology, Transactions A: Science, 42(4):2159–2166, 2018. [25] Gregory Eric Fasshauer, Abdul QayyumMasud Khaliq, and David Albert Voss. Using meshfree approximation for multi-asset American options. Journal of the Chinese Institute of Engineers, 27(4):563–571, 2004. [26] Black Fischer and Myron Scholes. The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81:637–654, 1973. [27] Ahmad Golbabai and Omid Nikan. A computational method based on the moving least-squares approach for pricing double barrier options in a time-fractional Black– Scholes model. Computational Economics, 55(1):119–141, 2020. [28] Ahmad Golbabai, Omid Nikan, and Touraj Nikazad. Numerical analysis of time fractional Black-Scholes european option pricing model arising in financial market. Computational and Applied Mathematics, 38(4):173, 2019. [29] Mustafa Inc, Muhammad Nawaz Khan, Imtiaz Ahmad, Shao-Wen Yao, Hijaz Ahmad, and Phatiphat Thounthong. Analysing time-fractional exotic options via efficient local meshless method. Results in Physics, 19:103385, 2020. [30] Rashid Jan, Salah Boulaaras, Mohammad Alnegga, and Farah Aini Abdullah. Fractional-calculus analysis of the dynamics of typhoid fever with the effect of vaccina- tion and carriers. International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, 37(2):e3184, 2024. [31] Rashid Jan, Normy Norfiza Abdul Razak, Salah Boulaaras, and Ziad Ur Rehman. Fractional insights into Zika virus transmission: Exploring preventive measures from a dynamical perspective. Nonlinear Engineering, 12(1):20220352, 2023. [32] Guy Jumarie. Stock exchange fractional dynamics defined as fractional exponential growth driven by (usual) gaussian white noise. application to fractional Black-Scholes equations. Insurance: Mathematics and Economics, 42(1):271–287, 2008. [33] Guy Jumarie. Derivation and solutions of some fractional Black-Scholes equations in coarse-grained space and time. Application to Merton’s optimal portfolio. Computers Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 26 of 27 & mathematics with applications, 59(3):1142–1164, 2010. [34] Mohan K Kadalbajoo, Alpesh Kumar, and Lok Pati Tripathi. Application of the local radial basis function-based finite difference method for pricing American options. International Journal of Computer Mathematics, 92(8):1608–1624, 2015. [35] AQM Khaliq, DA Voss, and SHK Kazmi. A linearly implicit predictor–corrector scheme for pricing American options using a penalty method approach. Journal of Banking & Finance, 30(2):489–502, 2006. [36] Hasib Khan, Jehad Alzabut, JF Gómez-Aguilar, and Abdulwasea Alkhazan. Essential criteria for existence of solution of a modified-ABC fractional order smoking model. Ain Shams Engineering Journal, 15(5):102646, 2024. [37] Hasib Khan, Altaf Hussain Rajpar, Jehad Alzabut, Muhammad Aslam, Sina Etemad, and Shahram Rezapour. On a fractal–fractional-based modeling for influenza and its analytical results. Qualitative Theory of Dynamical Systems, 23(2):70, 2024. [38] Jin-Rong Liang, Jun Wang, Wen-Jun Zhang, Wei-Yuan Qiu, and Fu-Yao Ren. The solution to a bi-fractional Black-Scholes-Merton differential equation. International Journal of Pure and Applied Mathematics, 58(1):99–112, 2010. [39] Farshid Mirzaee and Seyede Fatemeh Hoseini. Application of Fibonacci collocation method for solving Volterra–Fredholm integral equations. Applied Mathematics and Computation, 273:637–644, 2016. [40] Bjørn Fredrik Nielsen, Ola Skavhaug, and Aslak Tveito. Penalty and front-fixing methods for the numerical solution of American option problems. Journal of Com- putational Finance, 5(4):69–98, 2002. [41] Ömer Oruç. A new algorithm based on Lucas polynomials for approximate solu- tion of 1D and 2D nonlinear generalized Benjamin–Bona–Mahony–Burgers equation. Computers & Mathematics with Applications, 74(12):3042–3057, 2017. [42] J. A. Rad, K. Parand, and L. V. Ballestra. Pricing European and American options by radial basis point interpolation. Applied Mathematics and Computation, 251:363–377, 2015. [43] Ziad Ur Rehman, Salah Boulaaras, Rashid Jan, Imtiaz Ahmad, and Salma Bahra- mand. Computational analysis of financial system through non-integer derivative. Journal of Computational Science, 75:102204, 2024. [44] Siraj-ul-Islam and Imtiaz Ahmad. A comparative analysis of local meshless formu- lation for multi-asset option models. Engineering Analysis with Boundary Elements, 65:159–176, 2016. [45] Siraj-ul-Islam and Imtiaz Ahmad. Local meshless method for PDEs arising from models of wound healing. Applied Mathematical Modelling, 48:688–710, 2017. [46] S.O.Edeki, O.O Ugbebor, and E.A Owoloko. Analytical solutions of a time-fractional nonlinear transaction-cost model for stock option valuation in an illiquid market set- ting driven by a relaxed Black-Scholes assumption. Cogent Mathematics, 4:1–14, 2017. [47] G. Strang. On the construction and comparision of difference schemes. SIAM Journal on Numerical Analysis, 5:506–517, 1968. [48] Fuzhang Wang, Juan Zhang, Imtiaz Ahmad, Aamir Farooq, and Hijaz Ahmad. A Imtiaz Ahmad et al. / Eur. J. Pure Appl. Math, 18 (1) (2025), 5793 27 of 27 novel meshfree strategy for a viscous wave equation with variable coefficients. Fron- tiers in Physics, 9:701512, 2021. [49] Mehmet Yavuz and Necati Özdemir. European vanilla option pricing model of frac- tional order without singular kernel. Fractal and Fractional, 2(1):3, 2018. [50] H Zhang, F Liu, I Turner, and Q Yang. Numerical solution of the time fractional Black Scholes model governing European options. Computers & Mathematics with Applications, 71:1772–1783, 2016. [51] R. Zhang, Q. Zhang, and H. Song. An effcient finite element method for pricing Amer- ican multi-asset put options. Communications in Nonlinear Science and Numerical Simulation, 29:25–36, 2015. [52] Zhiqiang Zhou and Xuemei Gao. Numerical methods for pricing American options with time-fractional PDE models. Mathematical Problems in Engineering, 2016, 2016.