A numerical approach for 2D time-fractional diffusion damped wave model

1.   Introduction

In recent years, the fractional partial differential equations (PDEs) have attracted a growing attention in modeling in various applications in science and technology. The fractional derivative which is simultaneously possesses memory and non-local property can describe different non-linear phenomena more accurately and efficiently in comparison with the integer order derivative. This makes the fractional PDE a powerful tool for modeling complex dynamical systems ^[43]. Nowadays, fractional calculus is being used to report the non-differentiable problems for various types of fractional PDEs. In this regard, Singh et al. ^[45] find the non-differentiable solutions of local fractional wave equation in fractal strings by applying Local Fractional Homotopy Perturbation Laplace Transform Method (LFHPLTM). The uniqueness and convergence were also proved by the suggested method. Yang et al. ^[46] developed a new Boussinesq-type Model in a fractal domain which was based on the formulation of the local fractional derivative. To achieve exact traveling wave solution which can be obtained from the non-differentiable graph, first they convert the local fractional Boussinesq equation into nonlinear local fractional PDE. Fellah et al. ^[44] reviewed a new concept of the fractional PDE for the sound waves propagation in rigid porous materials. They experienced such a type of expression in which fractional power of frequency is related to the time derivative of fractional order that will describe the behavior of the sound waves. The results show the accuracy and capability. Kumar et al. ^[47] proposed a numerical algorithm to solve the vehicular fractal traffic flow problem. They solved the problem with the help of the local fractional Homotopy Perturbation Sumudu Transform method and the results were computationally sound for similar kinds of fractional PDE occurring in natural sciences. Besides this, Hajipour et al. ^[48] suggested a Variable-Order Fractional (VOF) model for the accurate approximate solution of reaction-diffusion equation in which they utilized a weighted shifted Grunwald derivative to solve the fractional part of the time derivative. This discretization technique was superior when compared with other methods in the literature. The variable order fractional calculus is an efficient tool to predict the compression deformation of amorphous glassy polymers. Meng et al. ^[49] proposed a three-regions- fitting-method and found that the VOF model is very efficient with higher accuracy.

In fractional calculus, the mathematical models are more accurately illustrate the characteristics of the real-world phenomena than the classical calculus, for example, Jajarmi and Baleanu ^[50] formulated a perfect model for describing the pathological behavior of HIV infection in fractional calculus that can never be investigated asymptomatic behaviors during the modeling with the integer-order derivatives. Baleanu et al. ^[51] examined their proposed model regarding the poor nutrition in the life cycle of humans with Mittag-Leffler nonsingular kernel in Caputo's fractional sense. The comparative numerical analysis reveals that the model based on the fractional derivative with Mittag-Leffler kernel has a different asymptotic behavior to the classic integer order derivatives. Shamasneh et al.^[52] proposed a Local Fractional Entropy (LFE)-based model for kidney images enhancement where the quality of MRI images are unpredictable. They tested the model on poor-quality kidney image and found the excellent results when LFE techniques are applied. Thus, the new aspects of fractional calculus provide us with more stretchy mathematical models that help us to formulate more realistic judgments about real-world problems in a better way.

In many aspects, fractional damped wave equations (FDWEs) have taken priority over the ordinary diffusion wave equations with damping due to the numerous implementations in practical natural occurrences. In the literature, much importance has been received to the solution of FDWEs because of the mathematical concept in various branches of applied sciences ^{[1,2,3,4,5,36,39,40,41]}. In several areas of study, it is essential to find the exact or approximate solution of many types of fractional differential equations (FDEs) that will help researchers understand the mechanism and difficulty of the challenge. In the previous era, many investigators, scholars and researchers attempted to formulate the analytical solution to the challenges of various distinct forms of FDEs that were difficult assignments for them. For example, in the solution of time fractional diffusion damped wave equation (TFDWE), Schneider and Wyss ^[6] utilized Greens' functions and their charateristics in aspects of Fox functions to solve the diffusion wave equation, Gorenflo et al. ^[7] used a similar technique to solve the fractional diffusion wave equation employing Laplace transformation and found a scale-invariant solution in forms of weight functions. Many other analytical solutions of FDEs can be found in ^{[8,9,10,11,12]} utilizing the Laplace transformation, the Fourier transformation, the Green function and Homotopy perturbation techniques, etc.

The majority of the FDEs cannot be determined analytically so there is a need for approximations to find the approximate numerical solution by various iterative techniques such as the finite element approach ^{[13,14,15,16,17]}, finite difference method ^[18], spectral approach ^[19] and collocation techniques ^{[20,21,22,23,24]}. Such type of iterative technique reduces the computational complexities of the algorithm. As a result, the following sparse system of linear equations is achieved: $\textbf{Au = b}$ where $\textbf{A}$ is a non-singular matrix that is known, $\textbf{b}$, and $\textbf{u}$ are a constant vector and the unknown column vector. At this stage, the rotated iterative techniques become more significant in order to solve the discretized linear system compared to the conventional (standard) iterative methods. The basic principle of the rotated iterative scheme is the use of nearly the half of interior nodes of the solution domain and the involvement of these interior mesh points in the iterative process to achieve rapid convergence. The values on the corresponding grid points are explicitly evaluated using the finite difference formula until convergence is attained. In this way, the rotated point iterative scheme reduces total arithmetic operations used and consequently there is will be a reduction in computational complexities, execution of CPU timings (per iteration) and number of iterations of the algorithm. As a result of this concept, numerous researchers have proposed a lot of explicit group iterative techniques for solving the various types of partial differential equations (PDEs). Martins et al. ^[25], Othman et al. ^[26], Yousif and Evans ^[27], Ali and Kew ^[28], Othman and Abdullah ^[29], Kew and Ali ^[30], Evans and Haghighi ^[31] and Abu Arqub et al. ^[32,33] contributed their work with most efficient results.

Owing to their good performances to solve integer order derivatives, the efforts are now being created to implement new methodologies to solve fractional PDEs in a better way ^[42]. Preliminary, Balasim and Ali ^[34,35] published recently their work on the 2D time-fractional diffusion, advection and cable differential equations in which promising outcomes were observed for the fractional parameter $\alpha \in (0, 1)$ when the Caputo time-fractional were combined with skewed approximations for the space fractional derivatives. In this article, we have extended our research for the numerical solution of diffusion damped wave equation of fractional order $\alpha$, $(1 < \alpha < 2)$ by applying the rotated finite difference iterative method and detected the proficient results.

Various studies ^[28,30] have proved the rotated five-point finite difference methods provide better and excellent results than the standard five-point difference ones for the numerical solution of multi-dimensional PDEs. More precise to the rotated Crank-Nicolson method, seldom authors utilized this iterative method for the numerical solution of different types of 2D-PDEs for fast convergence and found capable results. Balasim and Ali ^[34,35] have shown promising results for the rotated five-point Crank-Nicolson over the standard five-point Crank-Nicolson for solution of 2D subdiffusion PDE of fractional order. Ali and Kew ^[28] and Kew and Ali ^[30] applied rotated Crank-Nicolson schemes on two and three-dimensional hyperbolic telegraph PDEs respectively and results were much better and efficient than the standard Crank-Nicolson S(C-N) method. Applying the rotated implicit difference approximation in combination with Caputo's fractional derivative on 2D-FDEs is very simple and straightforward in the case of fractional PDEs. But applying the rotated Crank-Nicolson difference R(C-N) approximation together with Caputo's fractional derivative on 2D fractional superdiffusion type PDEs (like fractional wave diffusion, fractional damped wave-diffusion and hyperbolic telegraph equation etc.) is a very difficult task for the researchers. Initially, Balasim and Ali ^[34,35] utilized R(C-N) technique for the numerical solution of 2D fractional subdiffusion PDEs, where the Caputo's fractional derivative is utilized to solve the fractional case. The implementation of rotated Crank-Nicolson method together with Caputo's fractional derivative on 2D fractional superdiffusion PDEs is still an open problem for the investigator in this field of research.

The organization of the article is as follows: We shall derive the $(3-\alpha)$ approximation from the Caputo fractional derivative at the point $(u_{i}, v_{j}, t_{k+1/2})$ and an expression for the standard finite difference S(FPCN) in section 2. The derivation of finite difference expression R(FPCN) will be given in section 3. The stability and convergence of the R(FPCN) technique will be discussed respectively in sections 4 and 5. The numerical experiments will be performed in section 6, followed by computational complexity. In section 7, we present some final remarks about our derived results in the conclusion.

2.   S(FPCN) iterative technique for the TFWDE

Consider the second-order 2D TFDWE as:

where, $1 < \alpha < 2,$ with initial conditions

and boundary conditions

where $\Omega = \{(u, v, t)/ \; \; 0 < u, v < L, \; \; 0\leq t \leq T \}$ and f is the general source term. In Eq (2.1), the term $\frac{\partial z(u, v, t)}{\partial t}$ is called the damped factor which usually leads initial directionality propagation of wave in a string. The motion of wave become completely random if the damped factor become zero and wave motion behavior like fractional wave-diffusion equation. We utilize the following expression of Caputo fractional derivative to solve the fractional part of Eq (2.1)

Define $t_{k} = k\tau, k = 1, 2, \cdots, l$, $u_{i} = ih_{u}, \; i = 1, 2, \cdots, m$ and $v_{j} = jh_{v}, \; j = 1, 2, \cdots, n$ where, $\tau = \frac{T}{l}$, $h_{u} = \frac{L}{m}$ and $h_{v} = \frac{L}{n}$. The knots in the interval are defined by $(iu, jv, k\tau) = (u_{i}, v_{j}, t_{k}),$ where

Let $Z_{i, j}^{k}$, and $z_{i, j}^{k}$ be exact and the approximate solution of Eq (2.1), take $f(u_{i}, v_{j}, t_{k})$ as $f_{i, j}^{k}$. We employ S(C-N) technique which is consistently based on space and time, and the Caputo fractional derivative of order $\alpha$, which correctly handles initial value problems for the solution of Eq (2.1).

The second order differential operator is discretized as follows,

Initially, to derive the required expression for Caputo time-fractional derivative $\frac{\partial^{\alpha} z(u, v, t)}{\partial t^{\alpha}}$, $(1 < \alpha < 2)$ at knot $(u_{i}, v_{j}, t_{k+1})$ for Eq (2.1), consider the Eq (2.2) for $p = 2$, we have

Hence,

where, $b_{s} = (s+1)^{2-\alpha}-s^{2-\alpha}, s = 0, 1, 2, \cdots, k$ and $\xi = (t_{k+1}-\zeta$).

In particular, for $s = 0$, we have

Next, we find the truncation error $R^{k+1}_{\tau}$ of above iterative scheme with the help of Taylor's series expansion, since

The truncation error $R^{k+1}_{\tau}$ evaluated as,

Taking the absolute value over the truncation error $R^{k+1}_{\tau}$, we have

where, $C$ is a constant and $M$ is the upper bound of $| \frac{\partial^{3} z(u_{i}, v_{j}, t_{k}) }{\partial t^{3}}|$. Therefore,

Therefore, the discrete form of time-fractional derivative $\frac{\partial^{\alpha} z(u, v, t)}{\partial t^{\alpha}}$ at the knot $(u_{i}, v_{j}, t_{k+1})$ is

where, $|R_{1}^{k+1}| \leq C\tau$.

Similarly, the discrete form of time-fractional derivative $\frac{\partial^{\alpha} z(u, v, t)}{\partial t^{\alpha}}$ at $(u_{i}, v_{j}, t_{k+1/2})$ is

where, $|R_{2}^{k+1/2}| \leq C\tau^{3-\alpha}$.

By utilizing Crank-Nicolson difference technique for the space derivatives, Caputo time-fractional derivative from Eq (2.3) and central difference approximation for first order differential operator leads Eq (2.1) to the following expression,

for $i = 1, 2, \cdots, m-1$, $j = 1, 2, \cdots, n-1$ and $k = 0, 1, 2, \cdots, l$, where

Rewrite the above equation as,

for $i = 1, 2, \cdots, m-1$, $j = 1, 2, \cdots, n-1$ and $k = 0, 1, 2, \cdots, l.$

Finally, the $S(FPCN)$ can be obtained as follows:

for $i = 1, 2, \cdots, m-1$, $j = 1, 2, \cdots, n-1$ and $k = 0, 1, 2, \cdots, l$, with initial conditions: $z_{i, j}^{0} = g(u_{i}, v_{j}),$ and boundary conditions:

$0\leq u, v\leq L, \; \; \; 0 \leq t \leq T.$

The technique defined in Eq (2.5) is a point-dependent technique on the computational domain that utilizes all the points for the knot convergence as illustrated in Figure 1. All knots of the solution rectangular domain are considered as iterative knots in the iterative procedure at each temporal level until the knot converged solutions are obtained for some predetermined convergence criterion as well as for some tolerance factor. After achieving grid converged solutions, the predicted knot solutions are used as a beginning prediction for the next time level.

3.   R(FPCN) iterative scheme for TFDWE

The R(FPCN) can be developed by utilization of R(C-N) difference technique, the Caputo time fractional derivative of order $\alpha$ from Eq (2.3) and central difference approximation for first order derivative. This will leads Eq (2.1) to the following expression,

for $i = 1, 2, \cdots, m-1$, $j = 1, 2, \cdots, n-1$ and $k = 0, 1, 2, \cdots, l.$

Finally, we get the R(FPCN) as follows,

for $i = 1, 2, \cdots, m-1$, $j = 1, 2, \cdots, n-1$ and $k = 0, 1, 2, \cdots, l,$ with initial conditions: $z_{i, j}^{0} = g(u_{i}, v_{j}),$ and boundary conditions:

$0\leq u, v\leq L, \; \; \; 0 \leq t \leq T.$

The main purpose of the rotated grid scheme is to evaluate the values of the half of nodes of the discretized domain to get quick convergence. On the other hand, the remaining grid points can be determined by the standard grid scheme.

The technique in Eq (3.1) is a half grid dependent technique on the solution rectangular domain that employs half of the knots for the grid convergence. In the half grid dependent technique, we divide knots of entire discretized domain into sets of two knots, one type of knot is recognized as iterative knots $\bullet$ while the other types of knots are recognized as direct knots $\circ$ as shown in the Figure 2. This procedure divides in two parts, in first part the approximation in Eq (3.1) iterates only on iterative knots $\bullet$ of the suggested domain until a certain convergence is attained. In the second part, when convergence is attained, the values on direct points $\circ$ calculated by approximation defined in Eq (2.5). By sweeping half of the knots, we decrease approximately half of the algorithm's computational difficulty, resulting in lower CPU timings for each iteration.

Lemma 1. The coefficients $b_{s}$, $s = 0, 1, 2, \cdots$ shown in Eq (2.3) fulfill the following characteristics:

(1) $b_{0} = 1$, $b_{s} > 0$, $\forall\; s = 0, 1, 2, \cdots$.

(2) $b_{s} > b_{s+1}$, $\forall\; s = 0, 1, 2, \cdots$.

4.   Stability analysis of the R(FPCN) iterative technique for TFDWE

The stability of the suggested finite difference technique using the Fourier method will be discussed in this part. Consider $\xi_{i, j}^{k}$ be the approximated solution of Eq (3.1) and the error is defined as:

for $i = 1, 2, \cdots, m-1$, $j = 1, 2, \cdots, n-1$ and $k = 0, 1, 2, \cdots, l$, and

This error satisfies Eq (3.1) and we have,

for all $1\leq i\leq m-1$, $1\leq j \leq n-1$, $0 \leq k \leq l$, with the initial and boundary conditions are

For $k = 0, 1, 2, \cdots, l$, define the grid function as,

The $\xi_{(u, v)}^{k}$ can be expressed in the Fourier series as follows:

where,

By the $l^{2}$ error norm definition and Parseval's equality, we have

Let the error Eq (4.1) has the solution of the following form

where, $\sigma_{1} = 2\pi l_{1}/L$ and $\sigma_{2} = 2\pi l_{2}/L.$

Using the above expression in Eq (4.1) we obtain,

where

Proposition 1. The R(FPCN) define in Eq (3.1) is unconditionally stable.

Proof. First, we utilize the mathematical induction to satisfy the following result.

If $\lambda^{k} (k = 1, 2, \cdots, l)$ satisfy Eq (4.2), then $|\lambda^{k}| \leq |\lambda^{0}|.$ For $k = 0$, from Eq (4.2), we have

Since, $\mu \geq 0$ and $0 < b_{1} < 1$, therefore, we have

Now assume that $|\lambda^{m}| \leq |\lambda^{0}|$ for $m = 1, 2, \cdots, k$, we have

i.e.,

Now utilizing the preceding result and Parseval's equality, we obtain the following inequality,

that is,

for all $k = 1, 2, \cdots, l.$

Therefore, the R(FPCN) in Eq (3.1) is unconditionally stable.

5.   Convergence of R(FPCN) iterative technique for TFWDE

To investigate the convergence of the R(FPCN) iterative technique defined in Eq (3.1). Suppose that $z(u_{i}, v_{j}, t_{k+1/2})$ is the exact solution illustrated by the Taylor series. Then the truncation error at point $z(u_{i}, v_{j}, t_{k+1/2})$ is denoted by $R_{i, j}^{k+1/2}$ and is defined as,

for all $i = 1, 2, \cdots, m-1$, $j = 1, 2, \cdots, n-1$ and $k = 0, 1, 2, \cdots, l.$

From Eq (2.1), we have

There is a constant $C_{1}$ for all $i = 1, 2, \cdots, m-1$, $j = 1, 2, \cdots, n-1$ and $k = 0, 1, 2, \cdots, l$ such that,

where,

Define the function,

From Eq (5.1), we have

In order to get the error equation subtract Eq (5.3) from Eq (3.1), we have

for all

This equation verifies the boundary conditions

and initial conditions

For $k = 0, 1, 2, \cdots, l$, define the grid function

and

for all $1\leq i\leq m-1$, $1\leq j \leq n-1$, $0 \leq k \leq l.$

The functions $\phi^{k}(u, v)$ and $R^{k}(u, v)$ can be defined in the Fourier series,

where,

By the $l^{2}$ error norm and Parseval's equality, we have

where,

Suppose that the solutions have the following forms,

where, $\sigma_{1} = 2\pi l_{1}/L$ and $\sigma_{2} = 2\pi l_{2}/L.$

Using the above equation in Eq (5.4) we have,

where, $\mu$ is defined in section (4).

Proposition 2. Assume that $\eta^{k} (k = 1, 2, \cdots, l)$ is the solution of Eq (5.9), then there is a positive constant $C_{2}$, so that $|\eta^{k}|\leq C_{2} |\zeta^{1/2}|, k = 1, 2, \cdots, l$.

Proof. From Eq (5.5), noticing that $\phi^{0} = 0$, we have, $\eta^{0}(l_{1}, l_{2}) = 0$.

For convergence the series, there is a positive constant $C_{2} > 0$ of the right side of Eq (5.8) such that,

where, $C_{2} = \underset{1 \leq k \leq l}{ \max}\{{C}^{k}\}$. We will verify this by mathematical induction.

For $k = 0$, from Eq (5.9) we have,

Since, $\mu \geq 0$ and $m_{o} < 2+m_{1}+\mu$, from Eq (5.10), we get thus,

Assume that $|\eta^{m}|\leq C_{2} |\zeta^{1/2}|, m = 1, 2, \cdots, k$. We prove that it is true for $n = k+1$. Again from Eq (5.9) we have,

From Lemma 1(2), we have $1 = b_{0} > b_{1} > b_{2} > \cdots > b_{k-1} > b_{k}$ and also we have $m_{o} = (\tau)^{\alpha} \Gamma(3-\alpha)$, $m_{1} = (\tau)^{\alpha-1} \Gamma(3-\alpha)$. In both cases, as $\tau \rightarrow 0$ (i.e. increase the number of time levels) implies that $b_{k-1}, \; b_{k}, \; m_{o}, \; m_{1} \rightarrow 0$ thus, we have

This completes the proof.

Theorem 1. The R(FPCN) defined by Eq (3.1) is $L_{2}$ convergent of order $\tau^{3-\alpha}+(h_{u})^{2}+(h_{v})^{2}$.

Proof. Apply the Proposition 2 and using Eqs (5.7) and (5.8), we have

There is a positive constant $C_{1}$ as in inequality Eq (5.2) such that

Therefore, $\|\phi^{k}\|_{2} \leq C(\tau^{3-\alpha}+(h_{u})^{2}+(h_{v})^{2})$ where $C = C_{1}C_{2}.$

6.   Numerical examples and discussions

To demonstrate the effectiveness and simplicity of the R(FPCN) technique in solving the 2D TFDWE, the numerical examples were performed on a PC with a Core 2 Duo 2.8 GHz processor and 2GB of RAM running with Window XP SP3 operating system in Mathematica 11. To get accurate results, the numerical examples were shown at three time levels for various mesh sizes 4, 8, 12 and 16. For the space discretization, we assume that $h = h_{u} = h_{v}$ in both $u$ and $v$ directions and different time steps 1/4, 1/8, 1/12 and 1/16 are considered for the time discretization $(0 < t < 1)$ in the suggested domain. The Gauss-Seidel technique was used throughout our computations for convenience with relaxation factor $\omega_{e}$ set to 1. The $l_{\infty}$ norm was utilized for the convergence with the tolerance factor $\varepsilon = 10^{-5}$.

The computational difficulty of any iterative numerical technique is entirely determined by the number of arithmetic operations performed per iteration. The higher the computational difficulty of iterative technique, the more operation performs by the algorithm denoting the slow convergence rate. To assess the computational complexity of both techniques suppose that the suggested domain is discretized with grid size $n$, and that the internal grid knots will be $\lambda^{2}$ where $\lambda = n-1$. Table 1 describes the computational complexity of the $S(FPCN)$ and the $R(FPCN)$ iterative methods before and after convergence, whereas Table 2 summarizes the total number of arithmetic used. The results of conducted numerical experiments of both the iterative methods are shown in Table 3 where the promising results of $R(FPCN)$ indicates that this method is more efficient than $S(FPCN)$ in the forms of execution time (in seconds), number of iterations (Ite.) and the total number of mathematical operations at various values of $h$, $\tau$ and $\alpha$. Moreover, $R(FPCN)$ technique has 40–45$\%$ less computational difficulty than $S(FPCN)$, whereas the calculating time and the number of iterations were observed almost 30–35$\%$ and 20–25$\%$ less, with the same degree of accuracy.

Example 1. Consider the following TFDWE ^[37],

where,

The initial conditions are, $z(u, v, 0) = g(u, v) = 0$, $z_{t}(u, v, 0) = 0$, and boundary boundary conditions

The exact analytical solution is $z(u, v, t) = t^{2} (u+v)(2-u-v).$

The diagram in Figure 3 represents the graphs of total number of operations, number of iterations and execution of timings (in sec.) of S(FPCN) and R(FPCN) iterative methods against the mesh sizes at different values of $\alpha$. It can easily be seen from these graphs, the quantity of elapsed time, number of iterations and total number of operations of R(FPCN) iterative scheme are considerable less as compare to S(FPCN) iterative scheme.

Table 4 shows the temporal convergence order of the numerical data obtained from the proposed rotated Crank-Nicolson scheme for different time nodes $\tau$ and $\alpha$, and for fixed space node $h = \frac{\pi}{60}$, when end time is taken $T_{end} = 1$. It is observed that R(FPCN) iterative scheme generates $(3-\alpha)$ temporal convergence order. Table 5 demonstrates spatial convergence order for different space nodes $h$ and $\alpha$, and for fixed time node $\tau = 0.01$, when space end point is chosen $L_{end} = 1$. This shows that R(FPCN) iterative scheme attains second order spatial accuracy. Therefore convergence rate of both tables verify the theoretical derivation of convergence order of our proposed R(FPCN) iterative scheme as derived in Theorem 1.

We calculate the error $e(\tau, h)$ in $L^{2}$ norm as $e(\tau, h) = \| Z_{exact}-z_{approx.} \|_{2}$ and temporal and spatial convergence order of the proposed method by the following formulas ^[38]:

where, $e_{1}$ and $e_{2}$ are the errors correspond to the time step with size $\tau_{1}$ and $\tau_{2}$, and

where, $e_{1}^{\prime}$ and $e_{2}^{\prime}$ are the errors correspond to the space step with size $h_{1}$ and $h_{2}$.

Figure 4 represents graphs of temporal and spatial convergence rates of our proposed scheme at different values of $\alpha$. It can be observed from both graphs as the time or space steps (i.e. $h$ or $\tau$) decreases the temporal and spatial parameters converges to second order accuracy in time and space respectively, at different values of $\alpha$. Figure 5 illustrates the graph of absolute errors at different values of $\alpha$. It can easily be observed that the absolute error is significantly small at $\alpha = 1.75$ compared with $\alpha = 1.25$.

7.   Conclusions

In this paper, we formulate two numerical iterative schemes S(FPCN) and R(FPCN) by utilizing the $(3-\alpha)$-Caputo approximation for the time fractional derivative, and standard and rotated five points Crank-Nicolson approximations for the space derivatives in solving the diffusion damped wave equation of fractional order. The developed schemes are tested with a numerical example and results are compared. The numerical results revealed that the computational cost of R(FPCN) method is reduced as compared to S(FPCN) method and it also reduces the total number of operations, iteration numbers and CPU-timings which assist the algorithm in attaining rapid convergence. The derived numerical data sports our theoretical derivations. There are a lot of open problems that need to be solved by grouping strategies based on rotated finite difference approximations for several types of FDEs especially for three-dimensional case.

Acknowledgements

This research project is supported by Thailand Science Research and Innovation (TSRI) Basic Research Fund: Fiscal year 2023 under project number FRB660073/0164.

Conflict of interest

The authors declare that they have no conflicts of interest.