On spectral numerical method for variable-order partial differential equations

1.   Introduction

Fractional calculus has got great popularity among the researchers of different areas of science and engineering. It has been reported that dynamical systems typically undergo two stages of development, one is from integer-order dynamical systems to fractional-order systems and techniques in the domains of solid mechanics ^[1], physics ^[2], finance ^[3], population growth ^[4], physiology ^[5], electro-mechanical ^[6], statistical mechanics ^[7], systems of fractional differential equations(FDEs) of different kinds appear frequently. Many engineering, physical and biological problems can be modeled using the applications of fractional calculus for more better results than classical calculus. The results of FDEs and integral equations are more accurate and precise as compared to conventional differential equations, see examples ^{[8,9,10,11,12]}. On the basis of uniqueness and existence, plenty of research articles have been published (we refer ^{[13,14,15,16,17]}).

Here we state that the area where fractional derivative of constant real or complex order has been explored very well. Further huge amount of research work in the said area have been published. As in fractional calculus, the order of differentiation and integration is arbitrary or it can also be a function. Hence the generalization of Reimann-Liouville fractional order to variable-order was founded by Ross and Samko in 1993. Thus Ross and Samko gave the spearheading and initiative work on variable-order operator of differentiation and integration (see ^[18]). A while later, several contributions have been reported including applications of variable-order operators in mathematical models of various phenomenons (we refer ^[19,20,21]). Consequent studies explored that variable-order fractional calculus can be exceptionally valuable in areas like viscous flows, mechanics and modeling many phenomenons. In this regards various articles have been written to investigate the aforesaid area from different aspects including existence theory, numerical and analytical results (see ^{[22,23,24,25]}). Due to these facts said derivatives have the ability to describe real problems in more comprehensive ways. To the best of our knowledge the numerical results based on operational matrices for variable order FDEs have very rarely investigated.

The theory of numerical approximation for solving differential equations of variable order is getting to be increasingly imperative. Since spectral methods based on polynomial basis have been considered in previous time significantly for the computation of numerical results of various problems. The aforesaid methods are stable and powerful as compared to other numerical methods like finite element and residual power series methods, radial base function methods, etc. Also the mentioned methods applied to any linear problems are stable and having faster convergence rate (see some detail in ^[26]). Here we remark that spectral method based on polynomials basis including Legendre, Jacobi and Bernstein polynomials, etc have been reported stable apply to any linear problems (see detail in ^[27]). Further the said methods have been considered very well due to the following points of interest like

● simple to utilize and having clear and evident form;

● offer a suitable framework to approximate the solution of a problem;

● avoiding complex calculations due to simpler weight functions;

● spectral methods use the idea of global representations to find greater order approximations;

● can be easily applied to any differential equations of variable and constant order.

● as compared to finite difference methods, spectral methods are global techniques.

Here it is interesting to say that there is significant difference between difference and spectral methods. The spectral methods utilize basis functions which are nonzero over the entire domain, while difference methods use basis functions that are nonzero only on very small sub domains. So far we know the aforesaid methods have been applied for PDEs with constant real or complex order. But never use in computation of numerical solutions of variable-order PDEs.

Motivated and inspired from the mentioned work above, the objective of this manuscript is devoted to form a scheme for the computation of numerical solutions to a class of variable-order PDEs. To achieve this goal, we utilize operational matrices of integration and differentiation based on shifted Legendre polynomials. We extend operational matrices of differentiation and integration using shifted Legendre polynomials from simple real or complex order to any variable -order. Hence an algorithm is established for a class of multi-terms variable- order PDEs to compute numerical solutions. These operational matrices are based on shifted Legendre polynomials of two variables. Based on these matrices, we reduce the proposed variable order PDE to some algebraic equation. Keeping in mind, that the proposed method omits discretization and collocation which save much more time and memory. We investigate the class of linear multi terms variable order PDE studied earlier in ^[28] under variable order derivative as

with initial conditions

where $\mathrm{c}_1, \mathrm{c}_2, \mathrm{c}_3, \mathrm{c}_4, \mathrm{c}_5$ and $\mathrm{c}_6$ are constants

Also

If we select $\alpha(t) = 2$ and $\beta{(t)} = 1$ in (1.1), it becomes integer order PDEs. Further, the problem (1.1) represents various classes of PDEs, by selecting difference values for coefficient and taking $\alpha(t) = 2$ and $\beta{(t)} = 1$ as:

● If $({\mathrm{c}_2)^{2}}-4({\mathrm{c}_1})({\mathrm{c}_2}) > 0$, then problem (1.1) becomes a class of hyperbolic PDEs.

● If $({\mathrm{c}_2)^{2}}-4({\mathrm{c}_1})({\mathrm{c}_2}) = 0$, then problem (1.1) becomes a class of parabolic PDEs.

● If $({\mathrm{c}_2)^{2}}-4({\mathrm{c}_1})({\mathrm{c}_2})\leq0$, then the problem (1.1) becomes a class of elliptic PDEs.

● If $\mathrm{c}_1 = 1, \ {\mathrm{c}_2} = 0, \ \mathrm{c}_3 = -c^2, \ {\mathrm{c}_4} = {\mathrm{c}_5} = {\mathrm{c}_6} = 0, \ g(t, x) = x$, the problem (1.1) becomes the famous wave PDE with source term in space variable been studied in ^[35] by using HWCs method.

In same way fixing the coefficients and orders, the famous Poisson PDE and Laplace PDE become special cases of the problem (1.1). The mentioned PDEs have wide range applications in the field of mechanics and electro-magnetics, solitons and turbulent flow theory. We utilize our adopted procedure to present the numerical results graphically. Also comparison with HWCs method is given. Further, absolute errors are recorded accordingly. For some recent work in numerical analysis of various problems, we refer ^[31,33] and the references therein.

This article is organized as: In Section 2, basics results are recalled from literature of fractional calculus. Section 3 is devoted to the development of operational matrices for variable order integration and differentiation. Section 4 is devoted to the establishment of the proposed method for variable-order PDEs and also supporting examples along with the graphs are provided. Last section is devoted to brief conclusion.

2.   Basic results

In this section, we recall some basic results which can be found in ^[18,29,30].

Definition 2.1. Let $\alpha > 0$ be continues and bounded function, then the variable- order Riemann–Liouville integral of a function $\mathrm{h}\in L[0, 1]$ is defined by

provided that integral exists on right side.

Definition 2.2. The Caputo derivative of a function $\mathrm{h}\in C[0, 1]$ with variable- order $n-1 < \alpha(t)\leq n$ can be defined as

From (2.1), the following results hold

Definition 2.3. ^[32] The Legendre polynomials on the interval $[-1, 1]$ are recalled as below

where ${Ł}_{0}(x) = 1$ and ${Ł}_{x}(x) = x.$

By the use of transformation $x = \frac{2t+1}{L}-1,$ we obtain the recursive relation for shifted Legendre polynomials denoted by ${Ł}_{L, r}(t)$ on $t\in [0, L]$ as

Further the shifted Legendre polynomial of degree $n$ in $[0, 1]$ is expressed as

Further, the orthogonality condition for shifted Legendre polynomials over [0, 1] is given by

Further in two variables, the said polynomials are expressed as

The orthogonality is computed as

Let $\mathbf{B} = C([0, 1]\times[0, 1])$ be the Banach space of square integrable functions. Then any function $U\in \mathbf{B}$ can be approximated in terms of two variables shifted Legendre polynomials as

where

On using $\mathbf{C}_{r, q} = \mathbf{C}_n$, then in analytical form of the said polynomial (2.6) can be written in matrix form as

where ${\mathbb{X}_{M}}$ is the coefficient vector of dimension $1\times M^2$ and ${\phi}_{M^2}(x, t)$ function vector of dimension $M^2\times 1$.

2.1.   Convergence and error bounds

Here we can prove the convergence of the proposed method in same line as done in ^[36].

Theorem 2.4. Suppose that $U (t, x)$ and $U_n (t, x)$ are the exact and the approximated solution respectively computed by the proposed algorithm, then

where $\Omega_1, \Omega_2, \Omega_3$ are terminated constants given by

Proof. Let $U_n(x, t)$ is the approximate solution computed in terms of interpolated polynomial for $U(x, t)$ at point $(x_i, t_i)$ by means of aforesaid method such that $x_i = \frac{i}{n}$ and $t_i = \frac{i}{n}$. Then using (2.8) and after computation, briefly we can write the error bounds as

Therefore

Now from (2.9), the error bound for the maximum absolute value

Since the space $\mathbf{B}$ bounds and the absolute error bound coincides, Therefore when $n\to \infty$, the right side of (2.9) goes to zero. Hence one has

This show the convergence of the procedure.

3.   Operational matrices

This section is devoted to the construction of operational matrices of variable- order differentiation and integration based on shifted Legendre polynomials. For the required purpose we follow the same rules as followed in ^[34].

Theorem 3.1. If ${\phi}_{M^2}(x, t)$ be the function vector, thenvariable- order integration with respect to $t$ of ${\phi}_{M^2}(x, t)$ yields

where

where $r = Mi+j+1$, $v = Ma+b+1,$ ${\zeta}_{v, r, z} = {\sigma}_{i, j, b, z}$ for $i, j, a, b = 0, 1, 2, \cdots, m$

and

Theorem 3.2. The variable-order differentiation of ${\phi}_{M^2}(x, t)$ with respect to $t$ is given as

where ${R}_{M^2 \times M^2}^{{\alpha(t)}, t}$ is the operationalmatrix of differentiation and is given as

$r = Mi+j+1$, $v = Ma+b+1$, $\mathbf{\Phi}_{v, r, z} = {\mathbf{\Pi}}_{i, j, a, b, z},$ for $i, j, b, z = 0, 1, 2, \cdots, n$

and

Theorem 3.3. The fractional order derivative of ${\phi}_{M^2}(x, t)$ w.r.t $x$ and $t$ is given by

and

where $r = Mi+j+1$, $v = Ma+b+1$, $\mathbf{\Omega}_{v, r, z} = \emptyset_{i, j, a, b},$ for $, a, b, i, j = 0, 1, 2, \cdots, k$

and

Remark 1. Operational matrices based on shifted Legendre polynomials play important roles in numerical analysis.

4.   General algorithm

On the bases of operational matrices developed previously, our considered problem (1.1) is reduced into simple Sylvester type equation. Then the algebraic equation is solved for the unknown coefficient matrix to compute the required numerical result for $U(x, t)$. Consider the approximation as

By using 3.1), one has

Now to find the values of $l_1$ and $l_2$, we use initial conditions ${U}(0, x) = \theta(x)$ and $U_{t}(0, x) = \phi(x)$ in (4.2). Hence one has $l_1 = \theta(x)$ and $l_2 = \phi(x).$ By inserting these values in (1.1), one has

We can write as

where

and

Now, from (4.3), we calculate the approximation of remaining terms of (1.1) as

So our considered class of variable -order PDEs (1.1) by using (4.3) and (4.4) reduces to the following form

By further simplification, (4.6) yields

Further simplification of (4.7) yields

Using $B = {G_{M^2}{A}-{Q}_{M^2}}$ and (4.8) in (4.7), one has

Hence (4.9) is of the form $X+XA+B = 0.$ Which can be solved using Matlab by Gauss elimination method to get $X$ for the required numerical solution.

5.   Numerical examples

The following examples are provided to demonstrate the method.

Example 1.

At $\alpha(t) = 2$, ${\beta(t)} = 1$, the exact solution is given by

Now we are approximating the solution through the suggested method and also represented graphically in Figure 1 respectively. Now in Table 1, we give absolute errors at $\alpha(t) = 1+\exp(-t), \ \beta = 1$

Example 2.

For ${\alpha(t)} = 2$ and ${\beta(t)} = 1$, the integer order solution is given by $U(t, x) = \sin(x+t)$ and $g(t, x) = 8\cos(x+t)-9\sin(x+t).$ Now we approximate the solution through the suggested method and presents graphically in Figure 3 respectively. Now in Table 2, we give absolute errors at various fractional order.

Example 3. Let choose the values of constants in $(1.1)$ such that $\mathrm{c}_1 = 1$, $\mathrm{c}_2 = 0$, $\mathrm{c}_3 = -c^2$, $\mathrm{c}_4 = \mathrm{c}_5 = \mathrm{c}_6 = 0$ and $g(t, x) = x$, we get the following wave equation with source term or in-homogeneous differential equation ^[35]

with

and

For $\alpha = \beta = 2$ $(3)$ gives exact solution of the following form

Now in Table 3, we give absolute errors at various fractional order. Further in Table 4, we compared the maximum absolute error at various fractional order of our solution and compared with HWCs method ^[35]. Let the maximum absolute error be $\|U-\bar{U}\|_{\infty}$, then Also we compared the exact solution with the given approximate solution in Figures 5 and $(6)$ respectively. We have compared our results with the results of HWCs method for the Example 3 by taking scale level $10$ and for HWCs method taking collocation points $32$. We see that the adopted spectral method is much more better than the HWCs method. Here if we increase the scale level further, the accuracy can be further enhanced.

6.   Conclusion and discussion

We have presented an algorithm for the numerical solution of variable- order PDEs by using the concept of operational matrices. The concerned matrices have been created by using shifted Legendre polynomials under variable -order Caputo derivative and Riemann-Liouville integration. The aforesaid matrices have the ability to reduce the proposed problem to some algebraic equation. With the help of Matlab, we have evaluated the mentioned equation to compute numerical results. Three test problems have been investigated by using different kinds of variable orders. We see that the proposed method provides more better results. The concerned differential and integral operators have greatest degree of freedom as compared to classical fractional order. We have also used different spaces points and computed the maximum absolute error for various scale levels. By increasing the scale level, the efficiency of the method can be enhanced and the accuracy be improved further. Further, we have compared our results in Example 3 with the results of HWCs. We see that spectral method adopted in this work provides much more better results than the aforementioned numerical method. Further the procedure can be extended to other nonlinear problems also. Further variable -order PDEs can also be used as power full tools to investigate various real world problems like diffusion, blood flow, etc. For the mentioned problems such type numerical methods are more powerful and helpful.

Acknowledgments

Authors Kamal Shah and Thabet Abdeljawad would like to thank Prince Sultan University for paying the APC and support through the research lab TAS.

Conflict of interest

There is no competing interest regarding this work.