Fractional calculus is the non-integer or fractional-order differentiation and integration that has received significant attention over the last two decades for its application in describing real-world problems. Non-linear fractional-order differential equations (NL-FDEs) play a major role in applied sciences and social sciences, such as material science, signal processing, control theory, finance, and food supplements [1–3]. Recently, many advanced and reliable approaches have been discussed by the researchers. For example, Shen et al. [4] discussed the analytical and numerical solutions for a 2D multi-term time-fractional DWE. They used the approach of the variable of separation to drive analytical solutions, as well as the properties of Mittag–Leffler functions, and showed that stability and convergence analyses were studied. Ruzhansky et al. [5] considered the multi-term DWE and used the Caputo derivative, a non-local initial problem, and the Mittag–Leffler function in this study. Fan et al. [6] investigated the inverse problem to recognize the initial value for the space–time FDWE. They utilized the Landweber iterative regularization method for the time–space FDWE and calculated the error between validity and stability. The approximate solutions for 1D and 2D multi-term time-fractional sub-diffusion and DWEs with smooth and non-smooth solutions have been discussed by Rashidinia et al. [7]. They applied the Legendre collocation method and the Caputo derivative for solving the proposed equation. Moreover, they found the convergence analysis and concluded that this method has the benefit of using limited Legendre polynomials to obtain accurate and acceptable results. Feng et al. [8] researched a novel approach based on the finite element method (FEM) 2D diffusion type non-integer order equation. The utilized scheme is approximated for three different types of equations. The matrix form of the scheme is generated using the FEM, and the formulation of the equation is determined. The theoretical analysis, such as stability and convergence, is established. The suggested technique is versatile and resilient and may be used to solve various multi-term time-fractional diffusion problems directly. The Caputo derivative is in the temporal direction of the 2D time-dependent FDWE, as studied by Yang et al [9]. They used the finite difference method (FDM) to discretize the fractional-order derivative and provide a meshless approximation in spatial directions using moving least squares (MLS), which may be used to solve more complicated problem domains. The study also investigated the time-related convergence and stability properties of the semi-discretized scheme. In another study, Yang et al. [9] considered a 2D non-linear time-fractional DWE solved using a Crank–Nicolson Legendre spectral technique. The time-stepping is performed using the Crank–Nicolson difference method, while the spatial discretization is performed using the Legendre spectral approach. Lyu et al. [10] suggested a fast and linearized FDM for solving the non-linear multi-term FDWE. The suggested technique is based on a weighted approach, fast L 2 − 1 σ discretization, and multi-term L 2 − 1 σ -type discretization. The obtained truncation error in the suggested weighted discretization to the multi-term Caputo derivative rigorously proved unconditional convergence by demonstrating certain crucial aspects of the refined coefficients of full discretization. Salehi [11] investigated multi-term FDWS in 2D and solved it using a meshless collocation approach. The shape functions for spatial approximation are constructed using the moving least squares reproducing kernel particle approximation. They also created a semi-discrete technique by discretizing Caputo’s time derivatives using a finite difference approximation. The obtained difference schemes demonstrated unconditional stability and convergence. Ghafoor et al. [12] developed a technique based on the approximations of finite difference and Haar wavelets. The approach is used to numerically solve 1 + 1 -D and 1 + 2 -D time fractional PDEs. Zhuang and Liua [13] worked on a finite domain and the 2D temporal fractional diffusion equation (2D-TFDE). The 2D-TFDE is modeled using an implicit difference approximation. The mathematical induction approach is used to assess the stability and convergence analysis. Heydari et al. [14] investigated the approximate solution to the variable-order (VO) space–time fractional non-linear diffusion-wave problem. The established approach combines the collocation and tau methods with the Chebyshev cardinal functions and their operational matrix of VO-FDs. A resilient and reliable numerical approach for solving a 2D VO non-linear FDWE on arbitrary domains was suggested by Shekari et al. [15] using the MLS meshless procedure for the space domain and the FDS for the time domain . The method was devised in such a way that it is independent of the uniformity of the domain in consideration and the solution of the algebraic system of equations. The proposed method was tested using a variety of space domains, including quadrilateral, rounded, triangular, and polar domains, as well as different types of non-linearity. Kumar et al. [16]provided the local collocation method that depends on radial basis functions to investigate the solution of time-fractional non-linear DWEs. They demonstrated the numerical schemes that are unconditionally stable and converge in semi-discrete. Ding [17] researched the creation of a high-order numerical approach for the 2D time-space FDWEs. A new approach with order O   τ 2 + h 4 + h 4 is derived based on the fourth-order fractional-compact difference operator, where the temporal step size\tau and h 1 and h 2 are the spatial step sizes, respectively. The energy approach is used to analyze the algorithm’s stability and convergence, and a numerical experiment is conducted to confirm the numerical algorithm’s viability. Li et al. [18] analyzed a 2D non-linear FDWE in both time and space. The spatial component is discretized using the Galerkin FEM, while the temporal part is discretized using the new ADI method,” which also proved stability and convergence. They presented a 2D FDWE with the fractional derivative of order α   1 < α < 2 . Li et. al. [19] considered the ADI analysis based on the Crank–Nicolson method and the Galerkin FEM, both of which are analyzed. The ADI scheme is unconditionally stable, and L 2 norm convergence is rigorously illustrated. Datsko et al. [20] proposed that the TFDE with mass absorption in a sphere is observed under the harmonic impact on the surface of a sphere. The TFD of “Caputo” is implemented. The Mittag–Leffler function is also used to express the Laplace transform with respect to time and the finite sin-Fourier transform with respect to the spatial coordinates. Ren and Sun [21] discussed the fourth-order compact algorithm for solving the TFDWE with Neumann boundary conditions. For the time-fractional derivative, L 1 discretization is used, whereas the compact difference methodology is used for spatial discretization. The compact difference scheme is unconditionally stable, and global convergence is carefully shown. In addition, the Crank–Nicolson scheme with second-order spatial precision is described and perhaps an error estimate as well. The 2D case is solved using the compact ADI difference algorithm with Robin boundary conditions. Yang et al. [22] considered the initial value problem with mixed initial conditions of TFDWEs. They used the truncated regularization method and investigated conditional stability and error approximations, as well as some acceptable numerical examples, to establish the method’s validity. The convergence estimates in this publication are not saturated when compared with other authors’ works, and the non-homogeneous factors are correlated. Ali and Abdullah [23] discussed the FDWEs of arbitrary order.” They explain efficiency and fruitfulness through numerical test examples. There are many related studies available in the literature discussing different types of fractional-order models. For instance, Nawaz et al. [24] considered two numerical approaches, found the solution for two fractional-order models, and compared their solution with classical solutions. The fractional-order operator is in the Caputo sense, and we obtained the results numerically and graphically. They claimed that the methods are quickly convergent and yield encouraging results. Farid et al. [25] examined the Laplace transform with an iterative method for the space–time fractional-order models. They compared the obtained results with other existing schemes in the literature, which are more feasible and effective. Sayevand and Jafari [26] introduced the fractional-order KdV model in a fractal domain and applied the transformation to convert the fractional order into an ordinary-order derivative. They analyzed the theoretical analyses and presented many numerical examples to confirm the accuracy of the suggested approach. Li et al. [27] studied the significant properties of the Caputo fractional derivative in the real line and further developed it in the complex plane, which is used in signal processing. Guariglia [28] derived the functional equation for the fractional-order derivative Hurwitz function and proved the relation between the fractional-order zeta function and Bernoulli numbers. In another study, Guariglia [29] presented the functional equation with the Grunwald–Letnikov fractional derivative and discussed the link with the distribution of prime numbers. The other comprehensive literature can be studied in [30–37].

According to the previously described literature, fractional calculus is still a relatively new field that requires more accurate numerical approaches to examine the more practicable FDEs. The goal of this research is to develop a more accurate and reliable numerical technique for the FDWE. An attempt to discretize the R-L integral operator numerically and implement it in the Riemann–Liouville fractional-order derivative to approximate the FDWE has not yet been made. This approach reduces computational complexity and increases accuracy. Analysis, including aspects such as stability, consistency, and convergence, was also investigated based on the Fourier method and Taylor’s series expansion.

The remainder of the paper is arranged as follows: the associated preliminaries are described in Section 2. Section 3 explains the proposed method. The stability analysis is discussed in Section 4. Consistency and convergence are addressed in Section 5. Section 6 presents the numerical results, and the conclusion is provided in Section 7.

Here, we consider the FDWE as follows [23]: ∂ 2 w x , y , t ∂ t 2 = D t 2 − q ∂ 2 w x , t ∂ x 2 + ∂ 2 w y , t ∂ y 2 + F x , y , t ,   q ϵ 0,1 ,   x ϵ 0 , L ,   y ϵ 0 , L ,   t ϵ 0 , T , (1) where the conditions are w x , y , 0 = φ 1 x , y , w t x , y , 0 = φ 2 x , y , (2) w 0 , y , t = φ 3 y , t , w L , y , t = φ 4 y , t , w x , 0 , t = φ 5 x , t , w x , L , t = φ 6 x , t , (3) where F x , y , t represents the non-homogenous term and D t 2 − q denotes the R-L fractional-order derivative operator of order q lying between 0 and 1.

The R-L fractional differential operator is the most important extension of the classical differential operator. The R-L fractional-order derivative is defined as follows [38]: D t 2 − q w x , y , t = 1 Γ q d 2 d t 2   ∫ 0 t w x , y , z t − z 1 − q d z = d 2 d t 2   I 0 q w x , y , t . (4)

The R-L integral operator can be defined as follows: I 0 q w x , y , t = 1 Γ q   ∫ 0 t w x , y , z t − z 1 − q   d z , (5) where I 0 q represents the R-L fractional-order integral operator of order q lying between 0 and 1. Equation 5 can also be written as I 0 q w x , y , t n = 1 Γ q   ∫ 0 t n t n − z q − 1   w x , y , z d z . (6)

We use the Jumarie property in [25] as = 1 q Γ q ∫ 0 t n W x , y , z d z q .

By discretization of the aforementioned equation, we obtain the following equations: = 1 Γ 1 + q ∑ k = 0 n − 1 ∫ t k t k + 1 w x , y , z d z q , = 1 Γ 1 + q ∑ k = 0 n − 1 w x , y , t n − k ∫ t k t k + 1 z 0 d z q . (7)

Applying the Jumarie property ∫ 0 t z m d z n = Γ 1 + m Γ 1 + n Γ 1 + m + n t m + n , = τ q Γ 1 + q ∑ k = 0 n − 1 w ( x , y , t n − k ) k + 1 q − k q , = τ q Γ 1 + q ∑ k = 0 n − 1 C k q w x , y , t n − k , (8) where Δ t = τ is the time step and can be defined as t k + 1 = k + 1 Δ t = k + 1 τ and C k q = k + 1 q − k q , where k = 0 , 1 , 2 , … , n − 1 .

Lemma 1

The q 0 < q < 1 -order R-L fractional integral of w x , y , t in 0 , T is defined as I 0 q w x i , y j , t k = τ q Γ 1 + q ∑ v = 0 k − 1 C v q w x i , y j , t k − v . (9)

Lemma 2

The coefficient C k q k = 0,1,2 … satisfies the following properties [39]: • C 0 q = 1 , C k q > 0 ,   k = 0,1 , … . , • C k − 1 q > C k q ,   k = 1,2 , … . , • τ ≤ C 1 C k q τ q ,   w h e r e   C 1 > 0 , • ∑ k = 0 n C k q τ q = n + 1 q ≤ T q .

In this section, the implicit difference scheme (IDS) for the FDWE is constructed using Lemma 1 for the fractional-order part, and the space-derivative is reduced to the central difference approximation. The step for space is x i = i Δ x and y j = j Δ y , and the step for time is t k = k τ , where 1 ≤ i ≤ M − 1 , Δ x = L M , Δ y = V H , 0 ≤ k ≤ N , and  τ = T N , respectively. We substitute Eq. 2 in Eq. 1 at mesh point w x i , y j , t k as follows: ∂ 2 w x i , y j , t k ∂ t 2 = d 2 d t 2 I 0 q δ x 2 w x i , y j , t k Δ x 2 + d 2 d t 2 I 0 q δ y 2 w x i , y j , t k Δ y 2 + F x i , y j , t k . (10)

To eliminate the second-order time derivatives, we applied the backward difference approximation, then integrated from t k to t k + 1 , and also, used the trapezoidal rule for the forcing term. Finally, we obtained w i , j k + 1 − 2 w i , j k + w i , j k − 1 = I 0 q Δ x 2 δ x 2 w i , j k + 1 − 2 δ x 2 w i , j k + δ x 2 w i , j k − 1 + I 0 q Δ y 2 δ y 2 w i , j k + 1 − 2 δ y 2 w i , j k + δ y 2 w i , j k − 1 + τ 2 2 F i , j k + 1 + F i , j k . (11)

We substituted Lemma 1, and after simplification, we obtained the approximated scheme for the FDWE as follows: w i , j k + 1 − 2 w i , j k + w i , j k − 1 = [ S 1 w i + 1 , j k + 1 − 2 w i , j k + 1 + w i − 1 , j k + 1 − S 1 w i + 1 , j k − 2 w i , j k + w i − 1 , j k + ∑ v = 0 k − 1 C v + 1 q − C v q w i + 1 , j k − v − 2 w i , j k − v + w i − 1 , j k − v − S 1 ∑ v = 1 k − 1 C v q − C v − 1 q w i + 1 , j k − v − 2 w i , j k − v + w i − 1 , j k − v ] + [ S 2 w i , j + 1 k + 1 − 2 w i , j k + 1 + w i , j − 1 k + 1 − S 2 w i , j + 1 k − 2 w i , j k + w i , j − 1 k + S 2 ∑ v = 0 k − 1 C v + 1 q − C v q w i , j + 1 k − v − 2 w i , j k − v + w i , j − 1 k − v − S 2 ∑ v = 0 k − 1 C v q − C v − 1 q w i , j + 1 k − v − 2 w i , j k − v + w i , j − 1 k − v ] + τ 2 2 F i , j k + 1 + F i , j k . (12)

We know that w t x , y , 0 = φ 2 x , y . Therefore, w i , j − 1 = w i , j 1 + 2 τ φ 2 x i , y j , (13) w i , j 0 = φ 1 x i , y i , w 0 , j k = w M x , j k = w 0 , j k = w 0 , M y k = 0 , (14) where S 1 = τ q г q + 1 ∆ x 2 , S 2 = τ q г q + 1 ∆ y 2 , i = 1 , 2 , … . M x − 1 , j = 1 , 2 , … , M y − 1 , a n d   k = 1 , 2 , … , N − 1 .

The von Neumann method is used to determine the stability of the proposed scheme, and the approach in [40]is followed. Let w i , j k represent the exact solution for Eq. 12, and we obtain w i + 1 , j k + 1 − 2 w i , j k + 1 + w i − 1 , j k + 1 = [ S 1 w i + 1 , j k + 1 − 2 w i , j k + 1 + w i − 1 , j k + 1 − S 1 w i + 1 , j k − 2 w i , j k + w i − 1 , j k + S 1 ∑ v = 0 k − 1 C v + 1 q − C v q w i + 1 , j k − v − 2 w i , j k − v + w i − 1 , j k − v − S 1 ∑ v = 1 k − 1 C v q − C v − 1 q w i + 1 , j k − v − 2 w i , j k − v + w i − 1 , j k − v ] + [ S 2 w i , j + 1 k + 1 − 2 w i , j k + 1 + w i , j − 1 k + 1 − S 2 w i , j + 1 k − 2 w i , j k + w i , j − 1 k + S 2 ∑ v = 0 k − 1 C v + 1 q − C v q w i , j + 1 k − v − 2 w i , j k − v + w i , j − 1 k − v − S 2 ∑ v = 0 k − 1 C v q − C v − 1 q w i , j + 1 k − v − 2 w i , j k − v + w i , j − 1 k − v ] . (15)

The error is defined as E i , j k = w i , j k − W i , j k , and E i , j k satisfies Eq. 15 as follows: E i + 1 , j k + 1 − 2 E i , j k + 1 + E i − 1 , j k + 1 = [ S 1 E i + 1 , j k + 1 − 2 E i , j k + 1 + E i − 1 , j k + 1 − S 1 E i + 1 , j k − 2 E i , j k + E i − 1 , j k + S 1 ∑ v = 0 k − 1 C v + 1 q − C v q E i + 1 , j k − v − 2 E i , j k − v + E i − 1 , j k − v − S 1 ∑ v = 1 k − 1 C v q − C v − 1 q E i + 1 , j k − v − 2 E i , j k − v + E i − 1 , j k − v ] + [ S 2 E i , j + 1 k + 1 − 2 E i , j k + 1 + E i , j − 1 k + 1 − S 2 E i , j + 1 k − 2 E i , j k + E i , j − 1 k + S 2 ∑ v = 0 k − 1 C v + 1 q − C v q E i , j + 1 k − v − 2 E i , j k − v + E i , j − 1 k − v − S 2 ∑ v = 0 k − 1 C v q − C v − 1 q E i , j + 1 k − v − 2 E i , j k − v + E i , j − 1 k − v S 2 ∑ v = 0 k − 1 C v + 1 q − C v q E i , j + 1 k − v − 2 E i , j k − v + E i , j − 1 k − v − S 2 ∑ v = 0 k − 1 C v q − C v − 1 q E i , j + 1 k − v − 2 E i , j k − v + E i , j − 1 k − v ] . (16)

We assume that the growth factor takes the form of a single Fourier mode as E i , j k = λ k e − 1   σ Δ x i + Δ y j , (17) where σ and Δ x and Δ y are the mode number and step sizes, respectively. Eq. 17 can be the solution of the aforementioned error in Eq. 16. λ k + 1 e − 1   i Δ x σ 1 + j Δ y σ 2 − 2 λ k e − 1   i Δ x σ 1 + j Δ y σ 2 + λ k − 1 e − 1   i Δ x σ 1 + j Δ y σ 2 = [ S 1 λ k + 1 e − 1   i + 1 Δ x σ 1 + j Δ y σ 2 − 2 λ k + 1 e − 1 i Δ x σ 1 + j Δ y σ 2 + λ k + 1 e − 1   i − 1 Δ x σ 1 + j Δ y σ 2 − S 1 λ k e − 1   i + 1 Δ x σ 1 + j Δ y σ 2 − 2 λ k e − 1   i Δ x σ 1 + j Δ y σ 2 + e − 1   i − 1 Δ x σ 1 + j Δ y σ 2 + S 1 ∑ v = 0 k − 1 C v + 1 q − C v q λ k − v e − 1 i + 1 Δ x σ 1 + j Δ y σ 2 − 2 λ k − v e − 1 i Δ x σ 1 + j Δ y σ 2 + λ k − v e − 1 i − 1 Δ x σ 1 + j Δ y σ 2 − S 1 ∑ v = 1 k − 1 C v q − C v − 1 q λ k − v e − 1   i + 1 Δ x σ 1 + j Δ y σ 2 − 2 λ k − v e − 1   i Δ x σ 1 + j Δ y σ 2 + λ k − v e − 1   i − 1 Δ x σ 1 + j Δ y σ 2 ] + [ S 2 λ k + 1 e − 1   i Δ x σ 1 + j + 1 Δ y σ 2 − 2 λ k + 1 e − 1   i Δ x σ 1 + j Δ y σ 2 + λ k + 1 e − 1   i Δ x σ 1 + j − 1 Δ y σ 2 − S 2 λ k e − 1   i Δ x σ 1 + j + 1 Δ y σ 2 − 2 λ k e − 1   i Δ x σ 1 + j Δ y σ 2 + λ k e − 1   i Δ x σ 1 + j − 1 Δ y σ 2 + S 2 ∑ v = 0 k − 1 C v + 1 q − C v q λ k − v e − 1   i Δ x σ 1 + j + 1 Δ y σ 2 − 2 λ k − v e − 1   i Δ x σ 1 + j Δ y σ 2 + λ k − v e − 1   i Δ x σ 1 + j − 1 Δ y σ 2 − S 2 ∑ v = 1 k − 1 C v q − C v − 1 q λ k − v e − 1   i Δ x σ 1 + j + 1 Δ y σ 2 − 2 λ k − v e − 1 i Δ x σ 1 + j Δ y σ 2 + λ k − v e − 1   i Δ x σ 1 + j − 1 Δ y σ 2 ] . (18)

Dividing both sides by e − 1   σ Δ x Δ y i , j and then replacing e − 1   σ Δ x Δ y + e − − 1   σ Δ x Δ y = 2 − 4 sin 2 σ ∆ x ∆ y 2 , we obtain λ k + 1 = 1 1 + μ   λ k 2 + μ − λ k − 1 − μ ∑ v = 0 k − 1 C v + 1 q − C v q λ k − v − ∑ v = 1 k − 1 C v q − C v − 1 q λ k − v , (19) where μ = 4 S 1 sin 2 σ 1 ∆ x 2 + 4 S 2 sin 2 σ 2 ∆ y 2 .

Proposition 1

Suppose λ k + 1 , where k = 0,1 , … , N − 1 , is the solution of Eq. 19; then, we need to prove that λ k + 1 ≤ λ 0 . (20)

Proof. Let us consider k = 0 in Eq. 19 to prove the proposition using the induction method as follows: λ 1 = 1 1 + μ λ 0 2 + μ − λ − 1 . (21)

We get λ − 1 = λ 1 by substituting Eqs 13, 17 into Eq. 21. The provided equation becomes easier after simplification. λ 1 = λ 0 2 + μ 2 + μ . (22)

We obtain the following relation: | λ 1 | ≤ | λ 0 | . (23)

Suppose λ 1 ≤ λ 0 holds true for k = 1 , 2 , … , N − 1 .

Using Eqs 19–, 23 and Lemma 2, we have λ k − 1 | = 1 1 + μ λ k 2 + μ − λ k − 1 | − μ ∑ v = 0 k − 1 C v + 1 q − C v q λ k − v | − ∑ v = 1 k − 1 C v q − C v − 1 q λ k − v | , = 1 1 + μ [ λ 0 2 + μ − λ 0 − μ ∑ v = 0 k − 1 C v + 1 q − C v q − ∑ v = 1 k − 1 C v q − C v − 1 q λ 0 ] , ≤ 1 1 + μ [ 1 + μ − μ ( ∑ v = 0 k − 1 C v + 1 q − C v q − ∑ v = 1 k − 1 C v q − C v − 1 q ) ] λ 0 , ≤ 1 1 + μ 1 + μ − μ C 1 − C k − 1 − C 1 − C k − 2 λ 0 , ≤ 1 1 + μ 1 + μ − μ C k − 2 − C k − 1 λ 0 . (24)

From Lemma 2, the value is 0 < C k − 2 − C k − 1 < 1 , so it is clear that 0 < 1 + μ − μ C k − 2 − C k − 1 1 + μ < 1  and  λ K + 1 ≤ λ 0 . (25)

Here, λ K + 1 ≤ λ 0 and E i , j k + 1 ≤ E i , j 0 , so it is written as E i , j k + 1 2 ≤ E i , j 0 2 . It reveals that the proposed method is unconditionally stable.

Here, following the approach in [41], we assume that w is the closed-form solution, W is the estimated solution, and the function Y W = 0 is the approximated scheme for the proposed equation at the mesh point x i , y j , t k to discover the consistency analysis. The local truncation error at x i , y j , t k was subsequently indicated by Y W = T i k .

Theorem 1

The local truncation error T x , y , t for the suggested scheme is T i , j k = O Δ t 2 + O Δ x 2 + O Δ y 2 . T i , j k + 1 = W i , j k + 1 − 2 W i , j k + W i , j k − 1 − [ S 1 ∑ v = 0 k C v q W i + 1 , j k − v + 1 − 2 W i , j k − v + 1 + W i − 1 , j k − v + 1 − 2 S 1 ∑ v = 0 k − 1 C v q W i + 1 , j k − v − 2 W i , j k − v + W i − 1 , j k − v + S 1 ∑ v = 0 k − 2 C v q W i + 1 , j k − v − 1 − 2 W i , j k − v − 1 + W i − 1 , j k − v − 1 ] − [ S 2 ∑ v = 0 k C v q W i , j + 1 k − v + 1 − 2 W i , j k − v + 1 + W i , j − 1 k − v + 1 − 2 S 2 ∑ v = 0 k − 1 C v q W i , j + 1 k − v − 2 W i , j k − v + W i , j − 1 k − v + S 2 ∑ v = 0 k − 2 C v q W i , j + 1 k − v − 1 − 2 W i , j k − v − 1 + W i , j − 1 k − v − 1 ] . (26)

Using the Taylor series, we obtain T i , j k + 1 = w i , j k + Δ t ∂ u ∂ t | i , j k + Δ t 2 2 ! ∂ 2 u ∂ t 2 | i , j k + Δ t 3 3 ! ∂ 3 u ∂ t 3 | i , j k + Δ t 4 4 ! ∂ 4 u ∂ t 4 | i , j k + … − 2 w i , j k + w i , j k − Δ t ∂ u ∂ t | i , j k + Δ t 2 2 ! ∂ 2 u ∂ t 2 | i , j k − Δ t 3 3 ! ∂ 3 u ∂ t 3 | i , j k + Δ t 4 4 ! ∂ 4 u ∂ t 4 | i , j k − … − [ S 1 ∑ v = 0 k C v q w i , j k − v + 1 + Δ x ∂ u ∂ x | i , j k − v + 1 + Δ x 2 2 ! ∂ 2 u ∂ x 2 | i , j k − v + 1 + Δ x 3 3 ! ∂ 3 u ∂ x 3 | i , j k − v + 1 + x 4 4 ! ∂ 4 u ∂ x 4 | i , j k − v + 1 + … − 2 w i , j k − v + 1 + w i , j k − v + 1 − Δ x ∂ u ∂ x | i , j k − v + 1 + Δ x 2 2 ! ∂ 2 u ∂ x 2 | i , j k − v + 1 − Δ x 3 3 ! ∂ 3 u ∂ x 3 | i , j k − v + 1 + Δ x 4 4 ! ∂ 4 u ∂ x 4 | i , j k − v + 1 − … − 2 S 1 ∑ v = 0 k − 1 C v q   w i , j k − v + Δ x ∂ u ∂ x | i , j k − v + Δ x 2 2 ! ∂ 2 u ∂ x 2 | i , j k − v + Δ x 3 3 ! ∂ 3 u ∂ x 3 | i , j k − v + Δ x 4 4 ! ∂ 4 u ∂ x 4 | i , j k − v + … − 2 w i , j k − v + w i , j k − v − Δ x ∂ u ∂ x | i , j k − v + Δ x 2 2 ! ∂ 2 u ∂ x 2 | i , j k − v − Δ x 3 3 ! ∂ 3 u ∂ x 3 | i , j k − v + Δ x 4 4 ! ∂ 4 u ∂ x 4 | i , j k − v − … + S 1 ∑ v = 0 k − 2 C v q w i , j k − v − 1 + Δ x ∂ u ∂ x | i , j k − v − 1 + Δ x 2 2 ! ∂ 2 u ∂ x 2 | i , j k − v − 1 + Δ x 3 3 ! ∂ 3 u ∂ x 3 | i , j k − v − 1 + Δ x 4 4 ! ∂ 4 u ∂ x 4 | i , j k − v − 1 + … − 2 w i , j k − v − 1 + w i , j k − v − 1 − Δ x ∂ u ∂ x | i , j k − v − 1 + Δ x 2 2 ! ∂ 2 u ∂ x 2 | i , j k − v − 1 − Δ x 3 3 ! ∂ 3 u ∂ x 3 | i , j k − v − 1 + Δ x 4 4 ! ∂ 4 u ∂ x 4 | i , j k − v − 1 − … ] − [ S 2 ∑ v = 0 k C v q w i , j k − v + 1 + Δ y ∂ u ∂ y | i , j k − v + 1 + Δ y 2 2 ! ∂ 2 u ∂ y 2 | i , j k − v + 1 + Δ y 3 3 ! ∂ 3 u ∂ y 3 | i , j k − v + 1 + Δ y 4 4 ! ∂ 4 u ∂ y 4 | i , j k − v + 1 + … − 2 w i , j k − v + 1 + w i , j k − v + 1 − Δ y ∂ u ∂ y | i , j k − v + 1 + Δ y 2 2 ! ∂ 2 u ∂ y 2 | i , j k − v + 1 − Δ y 3 3 ! ∂ 3 u ∂ y 3 | i , j k − v + 1 + Δ y 4 4 ! ∂ 4 u ∂ y 4 | i , j k − v + 1 − … − 2 S 2 ∑ v = 0 k − 1 C v q   w i , j k − v + Δ y ∂ u ∂ y | i , j k − v + Δ y 2 2 ! ∂ 2 u ∂ x 2 | i , j k − v + Δ y 3 3 ! ∂ 3 u ∂ y 3 | i , j k − v + Δ y 4 4 ! ∂ 4 u ∂ y 4 | i , j k − v + … − 2 w i , j k − v + w i , j k − v − Δ y ∂ u ∂ y | i , j k − v + Δ y 2 2 ! ∂ 2 u ∂ y 2 | i , j k − v − Δ y 3 3 ! ∂ 3 u ∂ y 3 | i , j k − v + Δ y 4 4 ! ∂ 4 u ∂ y 4 | i , j k − v − … + S 2 ∑ v = 0 k − 2 C v q { w i , j k − v − 1 + Δ y ∂ u ∂ y | i , j k − v − 1 + Δ y 2 2 ! ∂ 2 u ∂ y 2 | i , j k − v − 1 + Δ y 3 3 ! ∂ 3 u ∂ y 3 | i , j k − v − 1 + Δ y 4 4 ! ∂ 4 u ∂ y 4 | i , j k − v − 1 + … } − 2 w i , j k − v − 1 + w i , j k − v − 1 − Δ y ∂ u ∂ y | i , j k − v − 1 + Δ y 2 2 ! ∂ 2 u ∂ y 2 | i , j k − v − 1 − Δ y 3 3 ! ∂ 3 u ∂ y 3 | i , j k − v − 1 + Δ y 4 4 ! ∂ 4 u ∂ y 4 | i , j k − v − 1 − … ] . (27)