Numerical study for a class of time fractional diffusion equations using operational matrices based on Hosoya polynomial

1.   Introduction

Polynomial approximation serves as a fundamental method across various domains of numerical analysis ^[1,2]. Not only does it provide a robust tool for approximating complex functions, but it also plays a crucial role in numerical integration and solving differential and integral equations. The Lagrange interpolation polynomial at Chebyshev points of the first or second kind has been observed to mitigate the Runge phenomenon ^[3], surpassing interpolants at equally spaced points. Moreover, the accuracy of approximation exhibits rapid enhancement with an increase in the number of interpolation points ^[4,5]. Functions of bounded variation hold significant importance in various branches of mathematical physics, optimization ^[6], free-discontinuity problems ^[7], and hyperbolic systems of conservation laws ^[8]. Additionally, these functions find application in image segmentation and related models ^[9]. However, despite their relevance, the theory of numerical approximations for such functions remains relatively underdeveloped, primarily due to the inherent singularities they exhibit.

A significant body of research has focused on approximating non-smooth functions through decay estimates of series coefficients. Xiang ^[10] explored the decay behavior of coefficients in polynomial expansions of functions with limited regularity, specifically examining Jacobi and Gegenbauer polynomial series. The goal is to derive optimal asymptotic results for the decay of these coefficients, investigating how the decay rate varies for functions with both interior and boundary singularities, across different parameters. Francesco et al. ^[11,12] introduced the constrained mock-Chebyshev least squares (CMCLS) approximation method, which aims to mitigate the Runge phenomenon by interpolating functions on nodes near Chebyshev-Lobatto nodes and using remaining nodes for regression in univariate and bivariate functions. More recently, Wang ^[13,14,15] addressed error localization in Chebyshev spectral methods for functions with singularities. This study begins with a pointwise error analysis for Chebyshev projections of such functions, revealing that the convergence rate away from the singularity is faster than at the singularity itself by a factor of $\frac{1}{x}$. The analysis rigorously explains the observed error localization phenomenon, suggesting that Chebyshev spectral differentiations generally outperform other methods, except near singularities, where the latter exhibit faster convergence.

Liu et al. ^[16] introduced a novel theoretical framework grounded in fractional Sobolev-type spaces, leveraging Riemann-Liouville fractional integrals/derivatives for optimal error estimates of Chebyshev polynomial interpolation for functions with limited regularity. Key components include fractional integration by parts and generalized Gegenbauer functions of fractional degree (GGF-Fs). This framework facilitates the estimation of the optimal decay rate of Chebyshev expansion coefficients for functions with singularities, leading to enhanced error estimates for spectral expansions and related operations. In a separate study, Wang ^[15] focused on deriving error bounds for Legendre approximations of differentiable functions using Legendre coefficients by Hamzehnejad ^[17]. Additionally, Xie ^[14] recently obtained bounds for Chebyshev polynomials with endpoint singularities. Zhang and Boyd ^[18] derived estimates for weak endpoint singularities, while Zhang ^[19,20] obtained bounds for logarithmic endpoint singularities. In ^[21], the focus lied on a specialized filtered approximation technique that generates interpolation polynomials at Chebyshev zeros using de la Vallée Poussin filters. The aim is to approximate locally continuous functions equipped with weighted uniform norms. Ensuring that the associated Lebesgue constants remain uniformly bounded is crucial for this endeavor.

The methodologies discussed above primarily concentrate on Chebyshev interpolation ^{[19,21,22,23]}, yielding results with exact Chebyshev series coefficients. However, computing exact series coefficients poses a general challenge and proves impractical for numerical algorithms, diminishing their utility in practical applications. Furthermore, these approaches usually involve Jacobi, Gegenbauer, and Legendre polynomials on fixed intervals where the respective series' basis functions are defined. Such limitations highlight the need for more versatile and efficient approximation methods in numerical analysis. Addressing this gap necessitates the utilization of efficient approximation techniques. Chebyshev polynomials, renowned for their versatility and effectiveness across diverse fields such as digital signal processing ^[24], spectral graph neural networks ^[25,26,27], image processing ^[22], and graph signal filtering ^[28,29] present a promising avenue for approximating functions of bounded variation. Many physical systems modeled using partial differential equations (PDEs) involve boundary layers or discontinuities, where functions of bounded variation frequently occur. Extending Chebyshev approximations to these functions allows for more accurate error analysis and truncation in numerical simulations of such systems.

Truncated Chebyshev expansions have proven capable of yielding minimax polynomial approximations for analytic functions ^[30]. Our objective is to employ these polynomials not only for approximating functions of bounded variation but also for conducting a comprehensive convergence analysis of Chebyshev polynomial approximation techniques. At the core of our convergence analysis lies the estimation of Chebyshev coefficients' decay. We leverage two recently established decay estimates: Majidian's decay estimate for Chebyshev series coefficients of functions defined on the interval $[-1, 1]$, subject to specific regularity conditions ^[31]; and a sharper decay estimate demonstrated by Xiang ^[10], with a more relaxed smoothness assumption on the function. In our pursuit of convergence results, we take an initial step by extending these decay bounds to encompass Chebyshev series coefficients of functions defined on the interval $[a, b]$.

The main contributions are three folds:

1. Generalization of Chebyshev polynomial approximation: The article extends the traditional Chebyshev polynomial approximation to a broader domain beyond the fixed intervals where basis functions are typically defined. This generalization allows for more flexible application of Chebyshev approximations in various settings.

2. Optimal error estimates: Two new optimal error estimates for Chebyshev polynomial approximations are presented, specifically tailored for functions of bounded variation. These error estimates are derived using approximated Chebyshev series coefficients rather than exact ones, addressing a significant gap in the existing literature.

3. Practical computation with approximated coefficients: By focusing on approximated Chebyshev series coefficients, the article offers a more practical approach for numerical algorithms, overcoming the challenges associated with computing exact series coefficients. The theoretical findings are supported by numerical experiments, providing empirical evidence of the efficacy and accuracy of the proposed error estimates.

These contributions collectively advance the understanding and application of Chebyshev polynomial approximations, particularly for functions with bounded variation, and offer practical solutions for numerical analysis.

While preparing this manuscript we found some very interesting recent works in the domain of approximation theory and Chebyshev polynomials. One of them ^[32] introduced unified Chebyshev polynomials (UCPs) and established their foundational properties, including analytic forms, moments, and inversion formulas. UCPs are shown to be expressible through three consecutive Chebyshev polynomials of the second kind. The authors derive new derivative expressions and connection formulas between different UCP classes, linking them with orthogonal and non-orthogonal polynomials. The second one ^[33] proposed two numerical schemes for solving the time-fractional heat equation (TFHE) using collocation and tau spectral methods. The authors introduce a new basis: Unified Chebyshev polynomials (UCPs) of the first and second kinds, deriving novel theoretical results for these polynomials.

The structure of the article is as follows: Section 2 provides the necessary preliminaries, including the Chebyshev series expansion of a function, the Gauss-Chebyshev quadrature rule, and several lemmas that are utilized to develop the main results. In Section 3, we derive decay bounds for the Chebyshev coefficients for functions of bounded variation and functions with limited smoothness. Section 4 presents $L^1$-error estimates for the Chebyshev approximation of $f$, leveraging the two decay estimates established in Section 3. Section 5 numerically demonstrates that the improved decay estimates of the Chebyshev coefficients and the $L^1$-error estimates of the truncated Chebyshev series approximation obtained in Section 4 are sharper than previously known results. Finally, we conclude the paper in Section 6 by outlining some promising future research directions.

2.   Preliminaries

2.1.   Chebyshev series expansion

The Chebyshev polynomial of the first kind, denoted as $T_j(t)$ for a given integer $j \geq 0$, is defined as:

where $\theta = \cos^{-1}(t)$ and $t \in [-1, 1]$. Notably, $T_j(t)$ is a polynomial of degree $j$ in the variable $t$. These polynomials exhibit orthogonality with respect to the weight function $\omega(t) = \frac{1}{\sqrt{(1-t^2)}},$ within the interval $[-1, 1]$. Specifically, they satisfy the following orthogonality relations:

The Chebyshev series expansion of a function $f:[-1, 1]\rightarrow \mathbb{R}$ is expressed as follows:

and

The norm $\|T_j\|_\omega$ is computed as:

Hence, the Chebyshev coefficients $c_j$ can be obtained using the integral form:

Given the difficulty in accurately evaluating the integral (2.4) for general functions, we resort to employing the Gauss-Chebyshev quadrature rule to approximate $c_j$, the $j$th coefficient of the series.

2.1.1.   Gauss-Chebyshev quadrature formula

Quadrature methods are renowned for numerically computing definite integrals of the type presented in (2.4). The Gauss-Chebyshev quadrature formula, a variant of Gaussian quadrature employing the weight function $\omega$ and $n$ Chebyshev points, provides an explicit formula for numerical integration (see ^[3,34,35]):

where $t_1, t_2, \ldots, t_n$ represent the $n$ roots of a Chebyshev polynomial $T_n(t)$ of degree $n$, also known as Chebyshev points, defined as:

Leveraging the quadrature formula (2.5), we can readily approximate Chebyshev series coefficients (2.4) for any function using the formula provided by Rivlin [35, p. 148] :

Here, $c_{k, n}$ denotes the approximated Chebyshev coefficients using $n$ quadrature points.

2.2.   Notations

We denote the Chebyshev series expansion of a function $f\in L^2_\omega[a, b]$ by $\mathsf{C}_\infty[f](x)$, defined as:

where $\mathsf{G}:[-1, 1]\rightarrow [a, b]$ is a bijection map given by:

The Chebyshev coefficients $c_j$ are calculated as:

Utilizing the change of variable $t = \cos \theta$, we express $c_j$ as:

The $d^{\rm th}$ partial sum, $C_d[f](x)$, approximates the function $f$ at a point $x\in [a, b]$, given by:

In our results, we use:

to denote the corresponding approximation of $f$ using $n$ quadrature points.

Additionally, we represent the Chebyshev series expansion of $f$, approximated with $n$ quadrature points, as:

The following lemmas are used in deriving the required error estimates.

Lemma 2.1. For a given positive integer $n$, we have

for any integer $k$ such that $0\le k < 2n$.

Proof. Since the identity is trivially satisfied for $k = n$, we assume that $k\ne n$.

Using (2.2) in the quadrature formula (2.7), we get

First, consider the case when $k = 0$. In this case, we can write

Using the fact that $T_0\equiv 1$, we see that

The first possibility is:

Thus, we can write

which is the required identity for $k = 0.$

We now assume that $k\ne 0$ (and recall that we have already assumed that $k\ne n$). We can write

Let us first evaluate the second term on the right-hand side. Since $0 < k < 2n$ with $k\ne n$, we see by taking $j = k$ that

and

and hence we see that

Therefore, the above expression can be written as

Let us now consider two cases, namely, $0 < k < n$ and $n < k < 2n$. We skip the proof of $0 < k < n$ and consider only the case when $n < k < 2n$ (note that we have already proved for $k = n$ separately).

(1) For $j = 0, \ldots, k-1$, we write $j = k-\alpha$ for $\alpha = 1, 2, \ldots, k$. Then, for some $p\in {\mathbb{{Z}}}^+$,

Since $\alpha$ ranges from 1 to $k$, we cannot have $p = 0$, for otherwise $\alpha = 2k$ which is not possible. Also, we see that $j+k\ne 2pn$, for any $p = 2, 3, \ldots$, for then $\alpha$ becomes negative. However, for $p = 1$, we have $\alpha = 2k-2n$. Thus,

Thus, we see that one term in the first summation within the brace of (2.12) is nonzero depending on the given value of $k$ between $n$ and $2n$. Note that, for this to happen, we need $n\ge 2$ (because only then $\alpha$ will have a meaningful range). Also note that in order for the present case of $n < k < 2n$ to happen, we need $n\ge 2.$ An interesting point is that this is one of the assumptions in Theorem 3.1.

On the other hand, for any nonnegative integer $s$,

Since $\alpha$ ranges from 1 to $k$ and $n < k < 2n$, we see that the above condition does not hold for any nonnegative integer $s$ and therefore

Thus we see that

Thus, we have

(2) For $j = k+1, k+2, \ldots$, let us write $j = k+\alpha,$ for $\alpha = 1, 2, \ldots$. For some $p\in {\mathbb{{Z}}}^+$,

This is not possible for $p = 0, 1$ because then $\alpha$ becomes negative. However, this is possible for $p = 2, 3, \ldots$, for which we have

On the other hand, for any $q\in {\mathbb{{Z}}}^+$,

This is possible for $q = 1, 2, \ldots,$ for which we have

Note that we have to see if both $j+k = 2pn$ and $|j-k| = 2qn$ hold for some $p, q\in {\mathbb{{Z}}}^+$. If this is so, then we must have the corresponding $\alpha$ values be equal. That is, for some $p$ and $q$,

This shows that both these cases happen if and only if $k$ is a multiple of $n$. But our present case is that $n < k < 2n$ and hence both these cases cannot happen simultaneously.

From the above discussions, we see that

for any $p = 2, 3\ldots$ and $q = 1, 2, \ldots.$. Thus, we have for $j = k+1, k+2, \ldots,$

Using this, we can write

Substituting (2.13) and (2.14) in (2.12), we get

This completes the proof.

Lemma 2.2. For $0\le d < 2n$, we have

for any $f\in L^1[a, b].$

Proof. For any $t\in [a, b]$, we have

By Lemma 2.1, we have

Note that each term of the right-hand-side series can be rewritten as

Substituting this expression in the above inequality, we get

Therefore,

Using the preliminaries and the lemmas presented in Section 2, we establish the decay estimates for the Chebyshev coefficients.

3.   Decay bounds for the Chebyshev coefficients

In this section, we extend the decay bounds established in prior works by Majidian ^[31] and Xiang ^[10]. This generalization is pivotal for numerous applications. In practical scenarios, the function to be approximated may not always reside solely within the domain $[-1, 1]$. Moreover, in various applications, local schemes ^[36] or piecewise approximation ^[37,38] are preferred over global ones. In such cases, decay estimates on a general domain become imperative.

Theorem 3.1. For some integer $k\ge 0$, let $f$, $f'$, $\ldots,$ $f^{(k-1)}$ be absolutely continuous on the interval $[a, b]$. If $V_k: = \|f^{(k)}\|_T < \infty,$ where

then for $j\ge k+1$ and for some $s\ge 0$, we have

where $c_k$, $k = 0, 1, \ldots$, are the Chebyshev coefficients of $f$.

Proof. For a given nonnegative integer $r$, we define:

with the understanding that $c_j^{(0)} = c_j$. In dealing with non-smooth functions, we must utilize the weak derivative (distributional derivative) on the right-hand side of the above expression, if it exists. Employing integration by parts in (3.3), we can express $c_j^{(r)}$ as:

for $j = 1, 2, \ldots$. In order to prove the required estimate, we prove the following general inequality :

for $m = 0, ...., k$ and $j \geq m+1$ using induction on $m$. Then $m = k$ gives our required result. First let us claim that (3.5) holds for $m = 0$. From (3.4) and (3.3), we have

This is precisely the inequality (3.5) for $m = 0.$ Let us assume that the inequality (3.5) holds for $m = 2s$ for some $s\ge 1$. Then for $m = 2s+1$ (odd), we have

Using the assumption that the inequality (3.5) holds for $m = 2s$, we can write

By simplifying the right-hand side, we get

which is precisely the required inequality (3.5) for $m = 2s+1.$ Finally, asssume that the inequality (3.5) holds for $m = 2s+1$. Then for $m = 2s+2$ (even), we have

Using the assumption that the inequality (3.5) holds for $m = 2s+1$, we can write

By simplifying the right-hand side, we get

which is precisely the required inequality (3.5) for $m = 2s+2.$ The proof now follows by induction.

Remark 3.1. Note that if $f^{(k)}$ is absolutely continuous, then $V_k$ is precisely the total variation of $f^{(k)}$ and hence in this case, the assumption that $V_k$ is finite implies $f^{(k)}$ is of BV on $[a, b]$. If $f^{(k)}$ involves a jump discontinuity, then one has to necessarily use the distribution derivative of $f^{(k)}$ in computing $V_k$.

The following lemma is the generalization of a result in ^[10].

Lemma 3.1. Let $f$ be a function defined on an interval $[a, b]$ such that for some integer $k\geq 1$, $c_j^{(k)}$ is well-defined and $f^{(k)}$ is of bounded variation on $[a, b]$. Then we have

where $j = p, p+1, \ldots$ and $p = 1, 2, \ldots, k$.

Theorem 3.2. Let $f$ be a function defined on $[a, b]$ such that for some nonnegative integer $k$, $f^{(k)}$ is of bounded variation with $V_k = Var(f^{(k)}) < \infty$. Then we have

for $j = k+1, k+2, \ldots$.

Proof. Since $f^{(k)}$ is of bounded variation, we can write (see Lang ^[39])

where $g_1$ and $g_2$ are monotonically increasing functions on $[a, b]$. Define $u(\theta): = g_i(\mathsf{G}(\cos\theta))$ which is monotonically decreasing for $\theta \in [0, \pi]$. Further, we have

where $v(\theta) = -u(\theta)$. By the second mean value theorem of integral calculus (Apostol [40, Theorem 7.37]), there exists $x_0 \in [0, \pi]$ such that

By the definition of $v$, we have

Substituting in (3.10) and then integrating yields

Using (3.3) and (3.9), we get $|c_j^{(k)}| \leq \dfrac{2V_k}{j\pi},$ for $j = 1, 2, \ldots.$ Consider (3.6) with $p = k$ which gives

Taking the modulus on both sides and using the above inequality, we get

for $j = k+1, k+2, \ldots,$ which leads to the desired result. Note that, in this case, $j = k$ is not defined when $i$ is $0$.

For the well-known decay estimate of the Chebyshev coefficients of a real analytic function, we refer to Rivlin ^[35] (also see Xiang {et al.} ^[41]).

4.   Error estimate for Chebyshev approximation

In this section, we derive $L^1$-error estimates for the Chebyshev approximation of $f$, utilizing two decay estimates provided in Theorems 3.1 and 3.2. Specifically, we establish the error estimate for the truncated Chebyshev series approximation, relying on the decay estimate (3.2) of the Chebyshev coefficients, as presented in the following theorem.

Theorem 4.1. Assume the hypotheses of Theorem 3.1. Then for any given integers $n$ and $d$ such that $n-1\ge k\ge 1$ and $k\le d \le 2n-k-1$, we have

where

(1) if $d = n-l$, for some $l = 1, 2, \ldots, n-k$, then we have

(2) if $d = n+l$, for some $l = 0, 1, \ldots, n-k-1$, then we have

for some integer $s\ge 0$, where

Proof. We have

For estimating the second term on the right-hand side of the above inequality, we use the well-known result (see Fox and Parker ^[42])

for $0\le d < 2n$. Using this property, with an obvious rearrangement of the terms in the series, we can obtain

where

By adding some appropriate positive terms, we can see that (also see Xiang et al. ^[41])

We restrict the integer $d$ to $k\le d \le 2n-k-1$ so that the decay estimate in Theorem 3.1 can be used. Now using the telescopic property of the resulting series (see also Majidian ^[31]), we can arrive at the required estimate.

Remark 4.1. From the above theorem, we see that for a fixed $n$ (as in the hypotheses), the upper bound $\mathbb{T}_{d, n}$ decreases for $d < n$ and increases for $d\ge n$. Further, we see that $\mathbb{T}_{n-l-1, n} = \frac{2}{3}\mathbb{T}_{n+l, n}$, however computationally $\mathsf{C}_{n-l-1, n}[f]$ is more efficient than $\mathsf{C}_{n+l, n}[f].$

The following theorem states the error estimate for the truncated Chebyshev series approximation based on the decay estimate (3.8).

Theorem 4.2. Assume the hypotheses of Theorem 3.2 for an integer $k\ge 1$. Let, for given integers $d$ and $n$ such that $n\geq k$ and $k\leq d \leq 2n-k-1$, $\mathsf{C}_{d, n}[f]$ be the truncated Chebyshev series of $f$ with approximated coefficients. Then we have

where

(1) if $d = n-l$, for some $l = 1, 2, \ldots, n-k$, we have

(2) if $d = n+l$, for some $l = 0, 1, 2, \ldots, n-k-1$, we have

Proof. Recall the $L^1$-error estimate for the truncated Chebyshev series expansion of $f$ with approximated coefficients given by (4.4) and (4.5):

Case 1. Now let us take the first case in the above estimate and apply Theorem 3.2 for $d = n-l, l = 1, 2, \ldots, n-k$, to get

Hence using the telescopic property of the above series, we have

Case 2. Similarly, for $d = n+l, l = 0, 1, \ldots, n-k-1$, we apply Theorem 3.2 to get

By substituting (4.10) and (4.11) in (4.9), we get

Hence, we have the required results.

5.   Numerical comparison

In this section, we numerically illustrate that the improved decay estimate of the Chebyshev coefficients and the $L^1$-error estimate of the truncated Chebyshev series approximation, obtained in Section 4, are sharper than the earlier ones obtained in ^[31].

Example 5.1. Let us consider the following example:

The function $g$ is absolutely continuous and

which is not continuous. Therefore, we have to take $k = 1$. Let us check the other hypothesis of Theorems 3.1 and 3.2. Using the weak derivative of $g'$, we can compute $V_1$ in Theorems 4.1 and 4.2 as

and the bounded variation of $g'$ is approximately equal to 2.7778, which is taken as the value of $V_1$ in Theorems 3.2 and 4.2.

The decay estimates (bounds given in (3.2) and (3.8)) of the Chebyshev series coefficients $c_j$ of $g$, for $j = 2, 3, \ldots$, as a function $j$ given in Theorems 4.1 and 3.2 are depicted in Figure 1(a). The error estimates for the truncated Chebyshev series that we obtained in Theorems 4.1 and 4.2 are demonstrated in Figure 1(b). It can be seen that the improved estimates we derived by using the results of Xiang ^[10] are sharper than the earlier ones we obtained using the results of Majidian ^[31].

6.   Conclusions

In conclusion, this study extends the applicability of Chebyshev series to functions defined beyond the traditional $[-1, 1]$ domain, broadening the scope of Chebyshev approximations for a variety of real-world applications. By introducing generalized decay bounds and truncation error results for Chebyshev approximations, we provide a more efficient and accessible framework for approximating functions, particularly in situations where exact computation of Chebyshev coefficients is not feasible. The results are highly relevant for fields such as spectral graph neural networks (GNNs), where Chebyshev approximations are commonly used to analyze graph signals and compute spectral filters efficiently. Additionally, these findings can enhance approximation techniques in image processing, particularly in tasks like edge detection and image compression, where rapid, accurate approximations are crucial. Overall, the theoretical advancements presented here offer a promising path to improve computational efficiency and approximation accuracy across a wide range of applications, particularly for functions with bounded variation.

Author contributions

S. Akansha: Conceptualization, methodology, writing-original draft preparation, visualization; Aditya Subramaniam: Investigation, validation, writing and proofreading, reviewing and editing. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

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

Acknowledgments

Dr. Akansha expresses her sincere gratitude to her Ph.D. supervisor, Prof. S. Baskar, for his invaluable insights and constructive feedback on the research problem addressed in this article. His thoughtful guidance greatly enhanced the clarity and overall quality of this work.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.