Foreword to the special issue "Contemporary PDEs between theory and modeling"

1.   Introduction

Orthogonal polynomials have significant roles in various mathematical and applied sciences. This is due to the polynomials' unique properties and characteristics. For example, these polynomials are very useful in approximation theory, numerical integration, and the solution of differential equations; see, for example, ^[1,2]. The Jacobi polynomials (JPs) are fundamental classical orthogonal polynomials used extensively in numerical analysis and approximation theory. These polynomials contain six essential classes. Four classes are the celebrated four kinds of Chebyshev polynomials (CPs). The first and second kinds of CPs were widely employed in many applications. The authors of ^[3,4] employed CPs of the first kind to handle some BVPs that arise in some applications. The authors in ^[5] numerically solved some singular Emden-Fowler equations using some combinations of the second-kind CPs. The authors of ^[6] utilized the CPs of the second kind to solve the population balance equation. Higher-order Emden-Fowler equations were solved via the CPs of the third-kind in ^[7]. There are many contributions regarding general JPs; see, for example, ^[8,9,10]. These polynomials involve two parameters. The derivative formulas of these polynomials involve a terminating hypergeometric function of the type $\ _3F_{2}(1)$. Generally, we cannot sum this hypergeometric function in a closed form, but we can do so for some specific choices of the two parameters of the JPs. The authors of ^[11] obtained reduced derivative expressions of certain Jacobi polynomials that generalize CPs of the third kind. This paper will utilize these polynomials to obtain new solutions for the fractional cable problem.

Fractional differential equations (FDEs) are powerful tools in many branches of the practical sciences. They shed light on many phenomena that standard DEs fail to address. The capacity to model genetic and memory processes is crucial. These equations can describe models related to biological and physiological processes, such as tumor growth and neuron action. Some applications for FDEs can be found in ^[12,13]. In most cases, the analytical solutions for the FDEs are unavailable, so it is vital to employ numerical analysis. Many approaches exist in the literature to handle different DEs. For example, in ^[14], the splines method was utilized for certain time-fractional diffusion equations. The Adomian decomposition method was applied in ^[15]. A collocation algorithm was applied in ^[16] for a type of KdV equation. In ^[17], the inverse Laplace transform methods were utilized to treat some FDEs. A neural network numerical method was applied in ^[18] for generalized Caputo FDEs. Based on an extreme learning machine, the authors in ^[19] followed another approach for FDEs. Matrix methods were used in ^[20,21,22] to treat some FDEs. A Haar wavelets method was applied in ^[23] for certain pantograph FDEs. Another wavelet approach was used in ^[24] to treat a system of nonlinear FDEs.

The cable equation is one of the most fundamental equations for modeling neural dynamics. Under simplifying assumptions, the Nernst-Planck electro-diffusion equation for ion transport in neurons has also been shown to be identical to the cable equation ^[25]. It has been explained that if ions are experiencing anomalous subdiffusion, then it is best to employ models that account for anomalous diffusion rather than those that assume normal diffusion since the latter is likely to provide deceptive diffusion coefficient values ^[26]. A fractional form of the Nernst-Planck equation was constructed by Langlands et al. ^[27] to represent the anomalous subdiffusion of the ion. Due to the importance of the different types of these equations, many authors were interested in dealing with them numerically. The authors of ^[28] proposed a numerical approximation to TFCE. Atta in ^[29] proposed two spectral methods for treating linear and non-linear TFCE. The authors of ^[30] proposed some approximate solutions for the two-dimensional TFCE. Another approach based on the finite element method was presented in ^[31] to deal with the distributed-order TFCE in two dimensions. A specific interpolation method was used in ^[32] for solving two-dimensional FCDE. Another method was used in^[33] for the two-dimensional fractional cable equations. In ^[34], a finite difference scheme was utilized to solve specific distributed-order cable equations. The author in ^[35] numerically treated the variable-order fractional cable equation.

Spectral methods have grown to the point where they are indispensable for improving numerical solutions of differential equations (DEs) in all fields. Compared with more conventional numerical approaches, these methods stand out due to their many advantages. The pinpoint accuracy of their solutions is one of their most vital points. It is also crucial to note that solutions are expressed as combinations of specific special functions or polynomials. There are three main methods of classifying spectral techniques. The Galerkin method has some conditions in choosing the trial and test functions; see, for example, ^{[11,36,37,38]}. One benefit of the tau technique over the Galerkin method is its greater flexibility when choosing the basis functions; see, for instance, ^[39,40,41]. Any differential equation can be treated numerically using the collocation approach, so many contributions exist in this area; see, for example, ^{[42,43,44,45,46]}.

This article proposes a tau approach to solving the TFCE using the GCPs3 as basis functions. The key formulas to design the desired method are the expressions of the polynomials' integer and fractional derivatives. We will explicitly transform the TFCE with its underlying conditions into a matrix system, yielding an algebraic system of equations that we can numerically treat. We can summarize our goals as follows:

● Developing some inner product formulas involving the derivatives of the GCPs3.

● Obtaining the fractional derivatives of the GCPs3.

● Applying the tau method to the TFCE to convert the problem into a matrix system that can be handled.

● Obtaining an error bound for the proposed numerical method.

● Testing the algorithm numerically by displaying examples supported by comparisons with other techniques.

The paper is organized as follows. Section 2 displays some fundamentals and necessary formulas. Section 3 focuses on implementing the tau algorithm designed for the TFCE's numerical treatment. Section 4 studies the error bound of the method used. Section 5 gives some illustrative examples accompanied by comparisons with other methods. We end the article with some conclusions in Section 6.

2.   Fundamental and some used definitions and formulas

This section presents some useful fundamentals in the sequel. The Caputo fractional derivative is given. We provide a brief account of JPs. Some characteristics of the GCPs3 are also taken into account.

2.1.   The Caputo sense of fractional derivative

Definition 1. ^[47] In the sense of Caputo, the fractional-order derivative is given by

The following identities hold:

where ${\mathbb{N}} = \{1, 2, \ldots\}$ and ${\mathbb{N}}_0 = \{0, 1, 2, \ldots\}$, and $\lceil \nu \rceil$ is the ceiling function.

2.2.   An account on the shifted normalized JPs

The shifted normalized JPs on $[0, 1]$ are a sequence of polynomials $\left\{\phi_k^{(\nu, \theta)}(\zeta)\right\}_{k\ge 0}$ that can be defined as

where $P_k^{(\nu, \theta)}(2\, \zeta-1)$ are the classical shifted JPs on $[0, 1]$.

$\left\{ \phi_k^{(\nu, \theta)}(\zeta) \right\}_{k \ge 0}$ is an orthogonal set on $[0, 1]$ in the sense that

where

and $\delta_{n, m}$ is the well-known Kronecker delta.

$\phi_k^{(\nu, \theta)}(\zeta)$ can be expressed as

where

Moreover, the inversion formula of $\phi_k^{(\nu, \theta)}(\zeta)$ is

2.3.   The shifted GCPs3

The shifted GCPs3 that were studied in ^[11] are specifically those of the normalized shifted JPs $\phi_k^{(\nu, \theta)}(\zeta)$ with $\theta = \nu+1$. Many formulas concerned with these polynomials were given in ^[11]. Here, we provide some formulas that will be used in this paper.

The orthogonality relation of $\phi_k^{(\nu, \nu+1)}(\zeta)$ is given by

where $w(\zeta) = (1-\zeta)^{\nu}\, \zeta^{\nu+1}$, and

$\phi_k^{(\nu, \nu+1)}(\zeta)$ may be expressed as

where

In addition, $\zeta^r$ can be represented as

In the following, an expression for the derivatives of the shifted GCPs3 is given.

Theorem 1. ^[11] For $s, q\in \mathbb{Z^+}$, with $s\ge q$, we have

where $G_{i, s, q}$ are given by

3.   Tau approach for the TFCE

This section analyzes an approach based on applying the Tau method to solve the TFCE. The integer and fractional derivatives of the basis functions of the GCPs3 will be used. Now, consider the TFCE ^[48]

governed by the conditions

where $K > 0$ and $\hat{\nu}$ are constants, $g(\zeta), \, h_{1}(t), \, h_{2}(t)$ are known continuous functions, and the source term is $f(\zeta, t)$.

Now, one may set

where $\Omega = (0, 1)\times(0, 1).$

Consequently, any function $\xi^{\mathcal{L}}(\zeta, t)\in \mathcal{P}^{\mathcal{L}}(\Omega)$ can be written in the form

where $\mathit{\boldsymbol{\psi}}(\zeta) = [\phi_0^{(\nu, \nu+1)}(\zeta), \phi_1^{(\nu, \nu+1)}(\zeta), \ldots, \phi_{\mathcal{L}}^{(\nu, \nu+1)}(\zeta)]$, $\left(\mathit{\boldsymbol{\psi}}(t)\right)^{T} = [\phi_0^{(\nu, \nu+1)}(t), \phi_1^{(\nu, \nu+1)}(t), \ldots, \phi_{\mathcal{L}}^{(\nu, \nu+1)}(t)]^{T},$ and $\mathit{\boldsymbol{C}} = (c_{ij})_{0\leq i, \, j\leq \mathcal{L}}$ is a matrix of unknowns with order $(\mathcal{L}+1)^2$.

The residual $\mathbf{Res}(\zeta, t)$ of Eq (3.1) has the following form:

The Tau technique, when used, leads to

with $\hat{w}(\zeta, t) = w(\zeta)\, w(t).$

Assume that

Therefore, Eq (3.7) can be rewritten as

or in the following matrix form:

The governing conditions in (3.2) and (3.3) give

Now, a suitable approach may be used to solve the system of algebraic equations of order $(\mathcal{L}+1)^{2}$, which includes Eqs (3.15) and (3.16).

In the following theorem, we will give explicitly the entries of the matrices $\mathit{\boldsymbol{\mathcal{A}}}, \, \mathit{\boldsymbol{\mathcal{B}}}, \, \mathit{\boldsymbol{\mathcal{H}}}, \, \mathit{\boldsymbol{\mathcal{K}}}$, and $\mathit{\boldsymbol{\mathcal{Q}}}$.

Theorem 2. The elements $a_{i, r}, \, b_{j, s}, \, {h}_{i, r}, \, k_{j, s}$, and $q_{j, s}$ are given by

where

and $\bar{B}_{r, j}$ is given in (2.11).

Proof. The first part of Theorem 2 is clear from the orthogonality relation (2.9). To show the second and third parts of Theorem 2, we will give a closed form for the integral

Based on Eq (2.13), we can write

and $G_{\ell, i, q}$ are as given in (2.14).

After applying the the orthogonality relation (2.9), relation (3.23) reduces to the following form:

and thus

and accordingly, we have

Now, if we set $q = 2$ in (3.26), then the elements of $h_{i, r}$ can be obtained as in (3.19). Furthermore, setting $q = 1$ in (3.26) yields $b_{j, s}$ as in (3.20).

To find $k_{j, s}$, we use Eq (2.10) to get

which leads to the following result:

Following similar procedures to those given to obtain elements of the matrix $\mathit{\boldsymbol{\mathcal{K}}}$, the entries of the matrix $Q$ can be obtained.

3.1.   Comments on computational complexity of the resulting system

This part is confined to describing the structure of the matrices $\mathit{\boldsymbol{\mathcal{A}}}, \, \mathit{\boldsymbol{\mathcal{B}}}, \, \mathit{\boldsymbol{\mathcal{H}}},$ and $\mathit{\boldsymbol{\mathcal{K}}}$ that appear in system (3.15). In addition, we comment on the resulting system and its numerical solution.

Here, we give the structure of the matrices $\mathit{\boldsymbol{\mathcal{A}}}, \, \mathit{\boldsymbol{\mathcal{B}}}, \, \mathit{\boldsymbol{\mathcal{H}}}, \mathit{\boldsymbol{\mathcal{K}}}$, and $\mathit{\boldsymbol{\mathcal{Q}}}$ taking the following forms for $\rho_1 = 0.2, \, \rho_2 = 0.8$, $\nu = 2,$ and $\mathcal{L} = 5$. The other choices for $\rho_1, \rho_2,$ and $\nu$ lead to the same structure of these matrices.

Remark 1. It is clear that the matrices $\mathit{\boldsymbol{\mathcal{A}}}, \, \mathit{\boldsymbol{\mathcal{B}}},$ and $\, \mathit{\boldsymbol{\mathcal{H}}}$ have special structures, however, the matrices $\mathit{\boldsymbol{\mathcal{K}}}$ and $\mathit{\boldsymbol{\mathcal{Q}}}$ are full matrices. This is because these matrices resulted from the fractional derivatives in (3.1).

Remark 2. It should be stressed that the numerical treatment is not effective alone when the resulting linear systems of algebraic equations are large, dense, and ill-conditioned, making the system potentially large and computationally intensive. In such case, we have to use suitable numerical solvers to treat these systems; one can refer to ^{[49,50,51,52]}.

4.   Error bound

This section examines the error analysis of the numerical solution $\xi^{\mathcal{L}}(\zeta, t)$ when compared to the exact solution $\xi(\zeta, t)$ of Eq (3.1), divided into the following two cases:

$1)$ Error analysis in the $L^{\infty}$ norm.

$2)$ Error analysis in the $L^{2}$ norm.

4.1.   Error analysis in $L^{\infty}$- norm

Assume $\mathcal{P}^{\mathcal{L}}(\Omega)$ as defined in (3.4). Then, for each $\hat{\xi}^{\mathcal{L}}(\zeta, t)\in \mathcal{P}^{\mathcal{L}}(\Omega)$ there exists a unique best approximation $\xi^{\mathcal{L}}(\zeta, t)\in \mathcal{P}^{\mathcal{L}}(\Omega)$, which is given by

The final inequality remains valid when $\hat{\xi}^{\mathcal{L}}(\zeta, t)$ represents the polynomials that interpolates $\xi(\zeta, t)$ at the points $(\zeta_i, t_j)$, where $\zeta_i$ are the zeros of $\phi_i^{(\nu, \nu+1)}(\zeta)$ and $t_j$ are the zeros of $\phi_j^{(\nu, \nu+1)}(t)$. Subsequently, by performing the same procedures outlined in ^[53,54], we obtain

where $\eta, \hat{\eta}, \mu, \hat{\mu}\in (0, 1)$.

Now, we have

Since $\xi(\zeta, t)$ is a smooth function on $\Omega$, then we can assume the existence of the positive constants $L_1, L_2,$ and $L_3$, such that

Using the one-to-one mapping $\zeta = \frac{1}{2}\, (z+1)$ between $[-1, 1]$ and $[0, 1]$, we can minimize the factor $\left\|\prod_{i = 0}^{\mathcal{L}}(\zeta-\zeta_{i})\right\|_{\infty}$. More precisely, we have

where $V_\mathcal{L}^{\nu} = \frac{2^{-\mathcal{L}}\, \Gamma (\nu +1) \, \Gamma (2 \, \mathcal{L}+2\, \nu +2)}{\Gamma (\mathcal{L}+\nu +1) \, \Gamma (\mathcal{L}+2\, \nu +2)}$ is the leading coefficient of $\phi_{\mathcal{L}+1}^{(\nu, \nu+1)}(z)$, and $z_i$ are the zeros of $\phi_{\mathcal{L}+1}^{(\nu, \nu+1)}(z)$. Moreover, $\left|\prod_{j = 0}^{\mathcal{L}}(t-t_{j})\right|_{\infty}$ can be minimized with the aid of the mapping: $t = \frac{1}{2}(\bar{t}+1)$.

where $\hat{V}_\mathcal{L}^{\nu} = \frac{2^{-\mathcal{L}}\, \Gamma (\nu +1) \, \Gamma (2 \, \mathcal{L}+2\, \nu +2)}{\Gamma (\mathcal{L}+\nu +1) \, \Gamma (\mathcal{L}+2\, \nu +2)}$ is the leading coefficient of $\phi_{\mathcal{L}+1}^{(\nu, \nu+1)}(\bar{t})$ and $\bar{t}_j$ are the roots of $\phi_{\mathcal{L}+1}^{(\nu, \nu+1)}(\bar{t}).$

Now, we have

Therefore, inequality (4.4) together with Eqs (4.5) and (4.6) leads to

This gives an estimation of the absolute error.

4.2.   Error analysis in $L^{2}$- norm

Theorem 3. Given that $\frac{\partial^{i+j}\, \xi(\zeta, t)}{\partial\, \zeta^i\, \partial\, t^j}\in \mathbf{C}(\Omega)$, $i, j = 0, 1, 2, ..., \mathcal{L}+1$, and let $\xi^{\mathcal{L}}(\zeta, t)$ be the proposed numerical solution belonging to $\Delta^{\mathcal{L}}$, and

where $\Omega = (0, 1)\times(0, 1)$. Then, the following estimation holds:

where $\mathfrak{\hat{a}}\lesssim \mathfrak{\bar{a}}$ means that there exist a generic constant $\mathfrak{n}$ such that $\mathfrak{\hat{a}} \leq \mathfrak{n}\, \mathfrak{\bar{a}}$.

Proof. Assume that

is the Taylor expansion of $\xi(\zeta, t)$ about the point $(0, 0)$, and

Since $\xi^{\mathcal{L}}(\zeta, t)$ is the best approximate solution of $\xi(\zeta, t)$, then according to the concept of the best approximation, we obtain

According to the inequality ^[55]

where $\mathrm{z}\geq 1,$ $\mathrm{z}+\mathrm{a} > 1$, $\mathrm{z}+\mathrm{b} > 1,$ and $\mathrm{a}, \mathrm{b}$ are any constants, and

We can rewrite Eq (4.13) as

Consequently, we get the following estimation

This completes the proof of this theorem.

Theorem 4. Suppose that $\xi(\zeta, t)$, $\xi^{\mathcal{L}}(\zeta, t),$ and $\frac{\partial^{i+j}\, \xi(\zeta, t)}{\partial\, \zeta^i\, \partial\, t^j}$ satisfy the condition of Theorem 3 and

Then, the following estimation holds:

Proof. Assume that $\frac{\partial^m\, v^{\mathcal{L}}(\zeta, t)}{\partial\, \zeta^m}$ is the Taylor expansion of $\frac{\partial^m\, \xi(\zeta, t)}{\partial\, \zeta^m}$ about the point $(0, 0)$. Then, the residual between $\frac{\partial^m\, \xi(\zeta, t)}{\partial\, \zeta^m}$ and $\frac{\partial^m\, v^{\mathcal{L}}(\zeta, t)}{\partial\, \zeta^m}$ can be written as

Since $\frac{\partial^m\, \xi^{\mathcal{L}}(\zeta, t)}{\partial\, \zeta^m}$ is the best approximation of $\frac{\partial^m\, \xi(\zeta, t)}{\partial\, \zeta^m}$, by the concept of the best approximation, we get

The desired result can be obtained by repeating similar procedures as in Theorem 3.

Theorem 5. Suppose that $\xi(\zeta, t)$, $\xi^{\mathcal{L}}(\zeta, t),$ and $\frac{\partial^{i+j}\, \xi(\zeta, t)}{\partial\, \zeta^i\, \partial\, t^j}$ satisfy the condition of Theorem 3 and

Then, the following estimation holds:

Proof. Similar to the proof of Theorem 4.

Theorem 6. Suppose that ${D_{t}^{\alpha}}\, \xi(\zeta, t)\in \mathbf{C}(\Omega)$, $\alpha\in (0, 1),$ satisfy the conditions of Theorem 3. Then, the following estimation holds:

Proof. According to Eq (4.12) and the properties of the Caputo operator in (2.3), one gets

Now, taking ${\|.\|}_{2}$ for both sides and follwoing similar steps as in Theorem 3 and 4, we get the desired result.

Corollary 1. The following estimation holds:

Proof. Special case from Theorem 6 after replacing $\alpha$ with $1-\alpha$.

Corollary 2. The following estimation holds:

Proof. The proof of Corollary 2 can be easily obtained from the application of Corollary 1 along with Theorem 4.

Theorem 7. The norm ${\|\mathbf{Res}(\zeta, t)\|}_{2}$ will be sufficiently small for sufficiently large values of $\mathcal{L}.$

Proof. Equations (3.6) and (3.1) enable us to write $\mathbf{Res}(\zeta, t)$ as

Taking the $L^2$-norm and using Theorem 5 along with Corollaries 1 and 2, we get

Lastly, it is clear from the final equation that for large enough values of $\mathcal{L}$, ${\|\mathbf{Res}(\zeta, t)\|}_{2}$ will be small enough. We have finished proving the theorem.

5.   Illustrative examples

The method discussed in Section 3 is used for solving a few illustrative examples to demonstrate the viability and effectiveness of the suggested generalized third-kind Chebyshev tau method (G3KCTM).

Test Problem 1. ^[56,57] Consider the equation

governed by

whose exact solution is $\xi(\zeta, t) = t^2\, \sin(\pi\, \zeta)$.

Table 1 presents a comparison of the $L^{\infty}$ error between our method and the methods in ^[56,57] when $\rho_1 = \rho_2 = 0.5$, and $\nu = 3$ with $\mathcal{L} = 13$ and $\mathcal{L} = 14$. At $\rho_1 = \rho_2 = 0.3$, $\nu = 1$, and various values of $\mathcal{L}$, the absolute errors (AEs) are shown in Figure 1. The $L^{\infty}$ error at $\mathcal{L} = 14, \, \nu = 1$ compared to the method in ^[56] at various $\rho_1$ and $\rho_2$ values is also shown in Table 2. Also, at various $\rho_1$ and $\rho_2$ values, Table 3 compares our method to the method in ^[57] in terms of the $L^{\infty}$ error at $\mathcal{L} = 14, \, \nu = 2$. Table 4 shows the maximum absolute errors (MAEs) for various $\nu$ values for $0 < t < 1$ and $\mathcal{L} = 14$ with $\rho_1 = \rho_2 = 0.4.$

Test Problem 2. ^[57] Consider the equation

governed by

whose exact solution is $\xi(\zeta, t) = \left(t^3+1\right) \left(\zeta^2-\zeta\right)$.

At $\mathcal{L} = 3, \nu = 3$ and with different $\rho_1$ and $\rho_2$, we compare our method to the method in ^[57] in terms of $L^{2}$ and $L^{\infty}$ errors in Table 5. This demonstrates the accuracy of our approach. With $\mathcal{L} = 3$ and $\rho_1 = \rho_2 = 0.4$, the AEs are displayed in Figure 2.

Test Problem 3. Consider the equation

governed by

whose exact solution is $\xi(\zeta, t) = t^3 \, \cos (\pi\, \zeta)$.

Figure 3 displays the AEs at $\rho_1 = 0.3, \, \rho_2 = 0.9$, $\nu = 1,$ and different values of $\mathcal{L}$. Table 6 displays the MAEs when $0 < t < 1$ at $\mathcal{L} = 16$ and $\rho_1 = 0.1, \, \rho_2 = 0.7,$ and different values of $\nu$. Figure 4 shows the AEs at $\rho_1 = \rho_2 = 0.5$, $\mathcal{L} = 16$, and different $\nu$.

Test Problem 4. Consider the equation

governed by

where $f(\zeta, t)$ is chosen such that the exact solution of this problem is $\xi(\zeta, t) = t^{5/2}\, \zeta^3$.

Table 7 displays the AEs at different values of $t$ when $\rho_1 = \rho_2 = 0.5, \, \nu = 3,$ and $\mathcal{L} = 8$. Also, Table 8 displays the AEs at different values of $t$ when $\rho_1 = 0.3, \, \rho_2 = 0.9, \, \nu = 2,$ and $\mathcal{L} = 8$. Finally, Figure 5 shows the AEs (left), and approximate solution (right) at $\rho_1 = 0.3, \, \rho_2 = 0.9, \, \nu = 2$, when $\mathcal{L} = 8$.

6.   Concluding remarks

This article analyzed an algorithm based on the tau method utilizing the GCPs3 to solve the TFCE that arises in neuronal dynamics. The philosophy of the derivation of the tau method was based on enforcing the residual of the equation to vanish. The algorithm converted the equation with its underlying conditions into a solvable matrix system. The entries of the matrices are explicitly computed using the explicit formulas of the integer and the fractional derivatives of the GCPs3. The error analysis is discussed in detail. Additionally, some examples and comparisons were displayed. The presented examples showed the applicability and high accuracy of the proposed algorithm. We believe the proposed approach in this paper may be applied to other significant problems in the applied sciences.

Use of AI tools declaration

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

Conflict of interest

The authors declare there is no conflict of interest.