High-order finite difference method for the two-dimensional variable-order fractional cable equation in complex systems and neuronal dynamics

1.   Introduction

The core concepts of an integral and derivative of integer order form the foundation of classical calculus. Many physical phenomena, including displacement, velocity, acceleration, area, and volume, require accurate description through classical calculus. Although classical calculus is very flexible, it is not very good at describing systems with memory effects, where the system's current state depends on its history ^[9]. For instance, the traditional diffusion model describes transport processes well with a mean square displacement that changes with time $\langle X^2 \rangle \sim K t$, where K is the diffusion constant. In more complicated settings, like porous surfaces or biological systems, the mean square displacement changes in a way that is not linear with time $\langle X^2 \rangle \sim K t^\gamma$, where $\gamma$ is the diffusion exponent. The change from standard operations results in "anomalous behavior, " characterized by diffusion occurring at rates that are either slower or faster than those predicted by the conventional diffusion model ^[7]. Therefore, it is clear that classical diffusion theory is not an appropriate way to explain these kinds of processes.

Classical calculus can be expanded by fractional calculus, which lets fundamental and differential operators have orders that are not integers ^[16]. Due to the non-local property of fractional differential operators, they are very useful for modeling complicated systems that traditional calculus cannot handle well. Because of this, fractional differential equations can be used to model different kinds of real-world problems in many areas, such as medicine, finance, science, and engineering ^{[10,11,12,26,32]}. These tools are very important for figuring out and explaining how complicated systems work when their dynamics are off or when their features depend on memory. Fractional-order derivatives have been developed in different forms and are available in the literature. The Riemann–Liouville, Caputo, Grünwald–Letnikov, Riesz, and Hadamard derivatives are among the most frequently used ^[3,13,20,25]. These operators, collectively referred to as classical fractional derivatives of the power-law type, form the foundation of traditional fractional calculus. Specifically, new fractional operators in terms of Mittag–Leffler and exponential kernel types have been defined which offer better adaptability and amenability. Further details on fractional derivatives are given in ^[15,18,21].

Fractional calculus is a very rapidly developing area of applied mathematics due to its unique capabilities in modeling various real phenomena via equations containing fractional derivatives. The proper choice of a fractional derivative for a model remains an open question despite many applications. Matlob and Jamali ^[17] asserted that no definitive guidelines exist for choosing the appropriate fractional derivative for mathematical modeling. In 2020, Tarasov and Tarasova ^[29] established a correspondence between the properties of the kernel of the fractional operator and the relevant physical phenomenon. Although immensely useful, analytical solutions for FPDEs are difficult to obtain, as already commented by ^[6]. Thus, the focus has shifted to numerical methods for solving FPDEs. The finite difference method, finite element method, finite volume method, spectral method, collocation method, reproducing kernel method, mesh-free method, domain decomposition method, and so on, are some of the standard techniques available for the solution of FPDEs ^[22,24].

V-O fractional calculus provides an alternative to the inappropriate use of fixed-order FPDEs for the formulation of phenomena that involve the variable memory effect, a function of time or space. Samko and Ross introduced this in 1993 through a generalization of Riemann–Liouville and Marchaud derivatives; it was later further developed by Lorenzo, Hartley, and Coimbra, among others, to offer diverse definitions of systems with differing memory. The variable-order (V-O) Caputo fractional derivative has an advantage over other forms of fractional derivatives in that it can better simulate dynamic systems with memory effects that change over time and space ^[30]. In contrast to the Riemann-Liouville derivative, the Caputo derivative permits physically accurate initial conditions, rendering it more appropriate for practical applications like neural dynamics and diffusion phenomena. The V-O Caputo derivative offers enhanced flexibility relative to constant-order fractional derivatives, allowing it to characterize intricate events where memory and heredity effects vary across time or geographic domains. This adaptability is especially vital in biological and physical systems, where phenomena like viscoelasticity, anomalous diffusion, and wave propagation display non-uniform behavior. As a result, models in these domains are more accurate and interpretable when the V-O Caputo fractional derivative is used. The V-O FPDEs are accurate; however, they are challenging to solve analytically and thus require numerical methods for solution. Consequently, this study addresses the numerical solution of the non-homogeneous initial-boundary value problem for the 2-D V-O FCE in the following form:

being subject to conditions

where the domain $\Omega = [0, L_0]\times [0, L_1]$, $\partial \Omega$ is the boundary, and $\mu$ is a positive constant. The V-O Caputo fractional derivative ($_0^CD_t^{\beta(x, y, t)})$ of $u(x, y, t)$ is defined as

The FCE is a biophysical model describing the voltage across neuronal cell membranes. It derives from the anomalous electrodiffusion of ions at spiny neuronal dendrites that necessarily escapes the classical cable equation. Originally introduced by Henry and Langlands in 2008 for anomalous diffusion processes in nerve cells, the FCE has found its way into applications as varied as neuronal dynamics, control theory, heterogeneity of neuronal tissue, and viscoelastic materials. The fact that it can model so many complex phenomena from one field to another demonstrates why solving the FCE is such a critical importance when it comes to getting science and engineering to move forward in those areas of research and application.

The progress in the numerical solution of FCEs has been very significant over the last few years because analytical solutions are generally very difficult to obtain. Bhrawy and Zaky ^[4] developed a shifted Jacobi polynomial-based spectral collocation method for fractional nonlinear cable equations in one and two dimensions. The Jacobi operational matrix for Caputo V-O fractional derivatives approximates global spatial and temporal discretization. The numerical results show that the method is accurate and converges, reducing error over finite difference and homotopy-based methods. However, operational matrices are highly computational, especially for large-scale problems or extremely complex geometries, limiting their use. The implicit radial basis function (RBF) meshless approach by Mohebbi and Saffarian ^[19] solves the 2-D V-O FC equation. These methods use thin plate splines for spatial discretization and second-order approaches for temporal derivatives. This numerical simulation method was accurate and efficient in regular and irregular domains. The study provides second-order temporal precision but does not address higher-order approaches or large-scale computing complexity, which may restrict its ability to scale for complex neural or biological systems. Adel et al. ^[2] established the non-standard wighted average finite difference method (NSWAFDM) to solve the 2-D V-O FCE, which is essential for neural dynamics and subdiffusion. According to John von Neumann, an arbitrary weight factor and various discretization schemes make the system stable. As $N$ and step size are refined, numerical results demonstrate the method's accuracy, reducing errors and confirming stability and convergence. The accuracy and computing efficiency are impressive, but it is limited to simple geometries and parameterizations, making it difficult to apply to complex systems. More research must focus on optimization and higher-dimensional problem adaptation. Similarly, Zheng et al.'s finite difference/Legendre spectral approach ^[31] is essential for understanding complex electrophysiology and neural dynamics phenomena, addressing the 2-D V-O time FCE. Finite difference using the midpoint quadrature rule, as well as temporal and Legendre spectral spatial discretization converted the V-O equation into a multi-term fractional system. The stability and convergence study determined that the approach converges with an order of $O(\tau^2 + \sigma^2 + N^{-s})$, and numerical results match predictions from theory. The method is versatile and successful, but high polynomial degrees and fine spatial-temporal resolutions need a lot of computer resources, which may restrict its scalability for large-scale simulation. Zou et al. ^[32] suggested a meshless method for 2-D nonlinear multi-term time FCEs (NM-TTFCEs) that converts governing equations into boundary value problems. The dual reciprocity method (DRM) and improved singular boundary method (ISBM) calculate specific and homogeneous solutions, while finite difference methods discretize time. The method's numerical results are accurate and efficient, especially in irregular geometries and complex computing domains. However, the study does not use high-order finite difference approaches or V-O fractional equations, which may limit its applicability to high-precision neural models. Habibirad et al. ^[8] investigated the direct meshless local Petrov–Galerkin (DMLPG) method for solving the 2-D time FCE. With GMLS spatial discretization and finite difference techniques for fractional derivatives, the system achieved unconditional stability and a convergence rate of $O(\tau^{2 - \max\{\alpha, \beta\}})$. In multiple 2-D geometrical domains, the algorithm proved its versatility and precision. DMLPG simplifies integration, but its complex shape function approximations could raise CPU consumption in larger systems or irregular node distributions, limiting scalability for more complex models. Sweilam et al. ^[27] proposed a higher-order nonstandard implicit compact finite difference technique for solving a 2-D nonlinear V-O FCE using fractional derivatives in the Atangana-Baleanu context. Jon von Neumann stability analysis confirmed the method's stability. Two numerical examples showed better accuracy, fewer errors, and more robust convergence than literature methods. Despite these benefits, the method requires extensive stability computations and is limited to basic geometries and parameterizations, making it difficult to apply to complex domains or higher-dimensional systems.

The literature discusses different ways to solve the FCE, but many numerical methods for V-O FCEs often face problems with stability, convergence, and efficiency because of the complicated changes in fractional orders over space and time. This paper presents a solution to these problems by introducing the HFDS that improves accuracy and efficiency while also being stable and converging. The suggested method uses a C-N higher-order finite difference approach, which is known for its precision, dependability, and effectiveness in dealing with complicated problems. The numerical scheme's efficiency is validated by a comprehensive stability and convergence analysis. Numerical simulations, presented in both tabular and graphical forms, further demonstrate the efficiency and accuracy of the suggested method.

The subsequent article's contents are organized in the same way that the C-N high-order finite difference scheme (HFDS) develops in Section 2. Similarly, Sections 3 and 4 describe the theoretical analysis (stability and convergence). Section 5 considers, solves, and discusses three test problems using the HFDS. The final section is the article's conclusion.

2.   The C-N HFDS

This section discretizes the V-O time-fractional and spatial derivatives in Eq (1.1) to develop a CN-HFDS. The following notations are utilized for this purpose:

where $\zeta$ and $\chi$ stand for time and space steps, respectively, and $L = L_0 = L_1$.

Let $\delta_x^2 \rho = \rho_{m+1, n}^k-2\rho_{m, n}^k +\rho_{m-1, n}^k$, and then, using Taylor series

From (2.1) and (2.2),

The use of higher-order correction terms in spatial discretization improves the method's accuracy. Incorporating terms up to $O(h^4)$ reduces truncation error, improving numerical stability and accelerating convergence. The high-order terms improve the approximation of second-order derivatives, guaranteeing that the numerical approach captures the dynamics of variable-order fractional derivatives with more precision.

The numerical approximation of the Caputo V-O time-fractional derivative of order $(0 \leq \beta (x, y, t) \leq 1)$ at the time level $t_{k+1/2}$ is obtained by employing the following discretization scheme ^[14]:

where

and $|\delta^{k+1/2}| \leq C\zeta.$

The C-N discretized form of (1.1) is

By using (2.3)–(2.5), the V-O FCE becomes

After simplification, the following HFDS is obtained:

where

Suppose

The fully discretized scheme can be written as:

where

and

The arrangement of the matrices given above indicates that "A" is strictly diagonally dominating, because $|m_6|\geq 2|m_7|+2|m_8|$. Since the matrix "A" is non-singular, the discrete scheme defined in Eq (2.6) has a unique solution.

Figures 1 and 2 show nine points on the grid and computational molecule for HFDS (2.6), where

Additionally, in Figure 3, the steps for implementing the V-O time-fractional numerical method for the HFDS are outlined clearly.

Lemma 1. ^[14]. Let $\mathcal{H}_s^{i, j, k}$, $0 \leq i \leq M - 1$, $0 \leq j \leq M - 1$, be the coefficients defined in Eq (2.6), then the following hold:

$({\rm{i}})$ $\mathcal{H}_s^{i, j, k} \geq \mathcal{H}_{s+1}^{i, j, k}, \quad 1 \leq s \leq k = 1, 2, \dots, N - 1.$

$({\rm{ii}})$ $\sum_{s = 1}^{k-1} \left(\mathcal{H}_{k-s}^{i, j, k} - \mathcal{H}_{k-s+1}^{i, j, k} \right) = \mathcal{H}_1^{i, j, k} - \mathcal{H}_k^{i, j, k}.$

3.   Stability

This section presents the stability analysis for the HFDS (2.6). Let $\rho_{m, n}^{k}$ and $V_{m, n}^{k}$ be the approximate and exact solutions for (1.1), respectively, and $\Lambda_{m, n}^{k} = V_{m, n}^{k}-\rho_{m, n}^{k}$. Then from (2.6), we obtain

The function $\Lambda^{k}(x, y)$ is defined as

The Fourier series expansion of error function $\Lambda_{m, n}^{k}$ ^[5] is

where

The $l^{2}$-norm function is used to define the error function $\Lambda_{m, n}^{k}$,

and the $l^{2}$-norm's Parseval equality relationship is

Suppose

where $\varphi_{1} = \frac{2\pi l_{0}}{L_0}$, $\varphi_{2} = \frac{2\pi l_{1}}{L_0}$, and $I = \sqrt{-1}$.

When (3.7) is substituted into (3.1) and after simplification, we obtain

Using Euler's formula, we have

Substituting (3.9) in (3.8), we have

After simplifying (3.10) for $\xi^{k+1}$, we got the simplified form:

where $p_0 = \cos(\varphi_{1} \chi) +\cos(\varphi_{2}\chi)$ and $p_1 = \cos(\varphi_{1}\chi)\cos(\varphi_{2}\chi)$.

Lemma 2. If $\xi^{k}$ $(0\leq k \leq N-1)$ satisfies (3.11), and $Q_2Q_1+Q_0 > 0$, then

Proof. The proof is completed using mathematical induction. Initially, select $k = 0$, resulting in

Taking absolute function and after simplifying (3.12), we have

where $Q_0 = 18 h^2 \mu _0-36 m_1 \left(p_0-2\right)+72 m_2 \left(p_0-p_1-1\right), Q_1 = 5 p_0+p_1+25$.

where $Q_2 = \frac{1}{2} 2^{\beta_{(i, j, k+\frac{1}{2})}} h^2 \sigma _1$. Since, $h, \sigma_1 \in (0, 1)$ and $\mathcal{H}_k^{i, j, 0} > 0$, therefore

Let

and we need to prove for $s = k+1$,

Using (3.13) and Lemma 1:

where $\aleph_0 = 25+ 5p_0 + p_1$.

After simplifing (3.15), we have

If $Q_2Q_1+Q_0 > 0$, then

Hence

□

Theorem 1. The HFDS (2.6) is stable when $Q_2Q_1+Q_0 > 0$.

Proof. Utilizing Lemma 2 and Parseval's equality, we obtain

Therefore,

which shows the HFDS (2.6) is conditionally stable. □

4.   Convergence

This section examines the convergence analysis of the proposed HFDS in relation to solving the 2-D FCE. The analysis quantifies the error bounds and demonstrates that the scheme attains a convergence order of $O(\zeta + \chi^4)$, thereby establishing its reliability and efficiency for numerical simulations of complex systems.

First, let the truncation error at $\rho(x_{m}, y_{m}, t_{k+1})$ be represented by $\Re^{k+\frac{1}{2}}$, and then

where $f_1$ is a constant.

But since $\Lambda_{m, n}^{k} = V_{m, n}^{k}-\rho_{m, n}^{k}$, then from (2.6):

subject to the conditions:

Similarly, the function for truncation error is defined as

where

and $\varphi_{1} = \frac{2\pi l_{1}}{L}$, $\varphi_{2} = \frac{2\pi l_{2}}{L}$.

Substituting (3.7) and (4.4) into (4.2), we get

Simplifying (4.5) for $\xi ^{k+1}$, we obtain

Lemma 3. Let $\xi^{k}$ $(0\leq k \leq N-1)$ satisfy (4.6), and $Q_2Q_1+Q_0 > 0$, and then

Proof. Since, from (4.3) and (3.3), we have

then from (4.1),

The proof is completed through mathematical induction, starting with $k = 0$ in (4.6):

Suppose

and we need to prove for $s = k+1$. Then from (4.6), we obtain

By using (4.11), (4.8), and Lemma 1, we have

After simplifying (4.12), we have

If $Q_2Q_1+Q_0 > 0$, then

□

Theorem 2. The HFDS (2.6) is conditionally convergent with convergence order $O (\zeta + h^4)$.

Proof. Using Lemma 3 and (3.6), we get

Therefore,

Thus, the proof is complete.□

5.   Numerical experiements

This section includes numerical examples of the2-D V-O FCE that are solved and examined using the HFDS. Numerical simulations are performed to assess the feasibility, accuracy, and efficiency of the proposed approach. Computational results regarding CPU time ($time$), iteration count ($Iter.$), average absolute error ($AAE$), and maximum absolute error ($MAE$) are presented in both tabular and graphical formats for comparative analysis. Moreover, the numerical examples were conducted on a personal computer with 4 GB of RAM and a Core i3 processor at 2.40 GHz, utilizing Mathematica (Version 14.1) with varying mesh sizes and time intervals. Furthermore, numerical simulations were conducted with maximal tolerance and error $(l_\infty)$. The convergence order of the HFDS is determined using the subsequent formulas ^[1].

Furthermore, the absolute maximum error is approximated using the formula

Test Problem 1: Let the source function in (1.1) be

with analytical solution

Test Problem 2: Suppose the source function for (1.1) is

with analytical solution

The numerical findings in Tables 1 and 2 validate the precision and robustness of the HFDS technique for addressing 2-D FCEs. The AAE and MAE for Test Problem 1 decrease continuously as $\chi^{-1}/\zeta^{-1}$ increases, as shown in Table 1. For example, with

the AAE decreases from $1.5896E-2$ at $\chi^{-1}/\zeta^{-1} = 5$ to $6.106 E-5$ at $\chi^{-1}/\zeta^{-1} = 30$, indicating convergence with smaller grid sizes. Similarly, the AAE of $\beta(x, y, t) = \frac{15+\sin(xt)}{16}$ decreases from $1.8464E-3$ to $7.4085E-5$ as the grid refinement improves, as shown in Table 3. These patterns emphasize the method's ability to effectively manage V-O derivatives while achieving error reduction, which is consistent with theoretical estimations. The ability to comply with von Neumann stability criteria and uniform performance across discretizations confirms the HFDS as a reliable and accurate method for addressing the V-O FCE.

Similarly, Tables 4–6 illustrate the computational advantages of the HFDS over Salama's approach ^[23] for the different grid sizes. In Table 4, the HFDS used $\chi = \frac{1}{18}$ and $\beta(x, y, t) = \frac{1 - (xyt)^3 + \sin(xyt)^2}{10}$ to minimize the maximum error from $2.5741E-3$ to $1.2631 E-5$ which shows the proposed scheme's accuracy in solving for the V-O FCE be demonstrating a gain of more than 200 times. While comparing the computational iterations, Salama required 150 iterations, but the HFDS required 45 which indicates the computing efficiency of the proposed scheme. Similarly, as illustrated in Table 5, the HFDS significantly reduced the error to $1.7315E-5$ for Test Problem 1 with $\beta(x, y, t) = 5t^{-3}$, surpassing Salama's ^[23] $1.5878 E-3$. Additionally, the number of iterations decreased from 435 to 51. This demonstrates, HFDS's capability to address the V-O FCE effectively at a lower computing cost. In Table 6, the HFDS had a $2.3955E-5$ error for

and $\chi = \frac{1}{18}$, but Salama's method had a $3.2975 E-4$ error for Test Problem 2. The accuracy of the HFDS in the test problem is demonstrated by this approximately 100-times improvement. These findings imply that the HFDS is computationally efficient and reduces errors. Also, increased accuracy and scalability are achieved through the optimization of sparse computational matrices, which demonstrates that problems with variable memory effects could be solved by the HFDS when accuracy and processing cost are important considerations.

Figure 4 illustrates the comparative accuracy of the proposed HFDS against Salama's method ^[23] for solving the 2-D V-O FCE. The results highlight a significant reduction in the MAE achieved by the HFDS, demonstrating its superior numerical precision. The figure further reinforces that as the grid is refined, the proposed method consistently outperforms the existing approach, reducing errors by several orders of magnitude.

The graphical representation of the error convergence is depicted in Figure 5 to further validate the theoretical convergence order $O(\zeta + \chi^4)$. The plot illustrates the reduction in the maximum absolute error (MAE) as the grid size ($\chi^{-1}$) and time-step size ($\zeta^{-1}$) are optimized. The convergence behavior is consistent across multiple lines that correspond to different fractional orders $\beta(x, y, t)$. The robustness and efficacy of the HFDS approach for solving V-O FCEs are confirmed by these results.

Figures 6–9 present 3D error graphs illustrating smoother distributions and reduced errors as the grid is refined, demonstrating the consistency and stability of the approach. As can be observed in Figure 4(a), where the numerical solutions correspond to the analytical solutions in Figure 6, where the maximum errors for Test Problem 1 decrease. Likewise, for Test Problem 2, errors are decreased with finer grids (Figure 6(b)). Furthermore, Figures 10 and 11 show heatmaps and contour plots of the HFDS's ability to handle the V-O FCE. The localized error peaks observed in coarser grids significantly decrease with finer discretizations, as illustrated in Figure 10(c) for Test Problem 1 and Figure 11(c) for Test Problem 2. The method's precision and adaptability in complex V-O FCEs are demonstrated by consistent patterns in analytical and numerical solutions (Figures 9 and 11) and the decrease in errors, which concludes the precision and computational effectiveness of the HFDS in handling the V-O FCE.

The spatial order of convergence for three fractional orders is illustrated in Table 3 and Figure 12, which illustrates the convergence of the HFDS for various $\beta$ values. The spatial order of the HFDS converges with a spatial order approaching 4 at early temporal orders for different tested $\beta$ values, as demonstrated in Figure 12, which shows that the theoretical convergence order of the approach is $O(\zeta + \chi^4)$. Moreover, the convergence trends are consistent across a variety of definitions of $\beta$, including $\beta = \frac{1 - (xyt)^3 + \sin(xyt)^2}{10}$, $\beta = \frac{5 + (xy)^5}{50}$, and $\beta = \frac{1 - (xyt) + \sin(xyt)}{10}$. Regardless of the complexity of fractional orders or the refinement of the mesh size, the HFDS's efficacy and versatility in managing the V-O FCE guarantee consistent performance.

From the above discussion, we conclude that the HFDS offers V-O FCE solutions that are precise, reliable, and effective, as well as flexibility for complex systems like neural dynamics. The results facilitate its application in more intricate higher-dimensional problems and complex geometries.

6.   Conclusions

In this work, a novel HFDS for solving the 2D V-O FCE—a mathematical model essential for modelling the dynamics of neurons and other biological systems with memory-dependent characteristics is presented. The proposed HFDS is developed using the C-N approach, which guarantees improved computational efficiency and accuracy. A thorough theoretical examination confirms that the approach is conditionally stable and convergent, attaining a convergence order of $O(\zeta^{\beta} + \chi^4)$, along with the existence of a unique solution. Numerical simulations confirm the efficacy and resilience of the proposed strategy, showcasing substantial enhancements compared to traditional methods. The comparative results in Tables 4–6 demonstrate significant decreases in maximum absolute errors and computational iterations, affirming the enhanced accuracy and efficiency of the HFDS. The results indicate that the suggested method is effective for addressing complex 2D FCEs, especially in the modeling of neural processes. Future research could aim to improve computing power by using parallel computing and adaptive mesh refinement.

Author contributions

Muhammad Asim Khan: Writing-original draft, Software, Formal analysis; Majid Khan Majahar Ali: Supervision, Validation; Saratha Sathasivam: Writing-review & editing, Supervision. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was supported by the Ministry of Higher Education Malaysia (MOHE) through the Fundamental Research Grant Scheme (FRGS), FRGS/1/2022/STG06/USM/02/11, and Universiti Sains Malaysia.

Conflict of interest

The authors declare that they have no conflict of interest.