A novel numerical approach for solving fractional order differential equations using hybrid functions

1.   Introduction

Fractional order differential equations (FODEs) are the extension of classical integer order differential equations (IODEs) using the concept of fractional calculus (FC) ^[1]. Compared to integer calculus (FC), FC is non-local and long-time memory, which enables FODEs be an alternative powerful mathematical tool for modeling the nature of many physical systems and real processes. For example, the relaxation modulus in viscoelastic material exhibits power-law behavior, which indicates that in the relaxation process the material has history memory. In this case, the relationship between stress and strain can be modeled using a FODE as $\sigma(t) = E\tau^\alpha\frac{d^{\alpha}\varepsilon(t)}{dt^{\alpha}}$ ^[2]. Also, the Warburg impedance has a fractional-power-law dependence on angular frequency, which is modeled as $Z(s) = As^{-1/2}$ ^[3].

In the last decades, linear or nonlinear FODEs have emerged in many scientific and engineering fields such as neural networks ^[4,5], chaos ^[6], diffusion process ^[7], fluid flows in porous media ^[8] and optimal control ^[9]. Now, more and more applications of FODEs can be found in practice, there is a great need of finding the solution of FODEs. However, it is not an easy thing to obtain the exact or analytic solution for most FODEs, pursuing numerical solution of FODEs is an important thing. In the past, the finite difference method ^[10], predictor-corrector method ^[11], Adomian decomposition method ^[12], variational iteration method ^[13] and Homotopy analysis method ^[14] had been proposed to solve FODES.

Though those methods are efficient and can provide good approximation to the real solution, however, they possess high computational complexities, which hampers its real application. So, developing efficient and effective methods for solving FODEs becomes an urgent task. Recently, the operational matrices had been widely adopted by researchers to solve different kinds of FODEs ^[15]. The core idea of operational matrices based methods is to transform the FODEs to a set of algebraic equations. As a result, the problem is dramatically simplified. In the literature, the operational matrices of fractional operators of block pulse functions (BPFs) ^[16], Bernstein polynomials ^[17], Bernoulli polynomials ^[18], Taylor polynomials ^[19] had been developed and adopted by several scholars to solve FODEs numerically. Among the above different polynomials or functions, the BPFs' fractional integral operational matrix is upper triangular matrix, in which the $k-$th row can be obtained by shifting the $(k-1)-$th row to the right. This an appealing property which enable one to solve the set of algebraic equations more easily and drastically reduce the computation burden involved in the computation of matrix inverse. However, BPFs are not enough smooth, more BPFs are required if high approximation accuracy is desired, and in turn, the dimension of operational matrix increases.

In recent years, the hybrid functions, which mainly combines BPFs or Haar function with other polynomials, had been widely adopted by researchers to solve different kind of FODEs. For example, in ^[20], the hybrid of BPFs and Bernoulli polynomials is adopted to find the numerical solution of nonlinear fractional integro-differential equations. In ^[21], the hybrid of Legendre polynomial and Haar function, called Legendre wavelet is used to solve distributed order FODEs. In ^[22], Chebyshev wavelets is adopted to solve nonlinear variable order FODEs. One of the advantage of hybrid functions is that they are piecewise polynomial instead of piecewise constant in each interval as BPFs. Thus, hybrid functions are smoother than BPFs, which leads to a more accurate approximation to a function than BPFs if an equal number of basis functions are used. Therefore, one can obtain more accurate solution for a FODEs using hybrid functions than BPFs.

Observing the aforementioned facts, this paper presents a new numerical method for solving FODEs based on hybrid of BPFs and Taylor polynomials (HBT). The fractional integrals operational matrix of HBT is derived through projecting the HBT functions onto a set of BPFs. Then, the FODEs to be solved is transformed into a set of algebraic equations using the derived operational matrix. Through solving the algebraic equations, one can obtain the numerical solution of FODEs.

The organization of this paper is arranged as follows. Some basics about fractional calculus are briefly reviewed in Section 2. In Sections 3, the basic formulation of HBT functions is given, and the calculation of fractional integral operational matrix of HBT is given in Section 4. The error analysis of the proposed method is given in Section 5. In Section 6, numerical experiments are presented to verify the effectiveness of the proposed method. In the last, the conclusive remarks are presented in Section 7.

2.   Definition of fractional calculus

Unlike integer calculus, the definition of FC is not unique. There are several definitions for fractional integration and derivatives. Among these definitions, the widely used one are the Riemann-Liouville (R-L) definition and the Caputo definition.

In R-L definition, the fractional integral of function $f(t)$ is given as ^[23]

where $\alpha > 0$ is the order, 0 and $t$ are the low and upper limits of fractional integrals, $\Gamma(\cdot)$ denotes the Gamma function. The fractional derivative of a function $f(t)$ in the sense of Caputo definition is

R-L fractional integral and Caputo fractional derivative has the following relationship

and

Since the derivation is the inverse operation of integration, the derivation of a function can easily be achieved from the result of integration.

3.   Properties of HBT functions

3.1.   BPFs and HBT functions

The Taylor polynomials defined on the interval $[0, T]$ is ^[24]

A set of BPFs $\psi_i(t), i = 1, 2, \cdots N$ is given as follows ^[25]

The BPFs $\psi_i(t)$ are disjoint and orthogonal, that is to say

A function $f(t)$ in $C^{n+1}[0, T)$ can be expanded into an $N$-term BPF series as

where $\Psi(t) = [\psi_1(t), \psi_2(t), \ldots, \psi_N(t)]^{\mathsf{T}},$ the constant coefficients $c_i$ in Eq (3.2) are given by

On the interval $[0, T)$, the hybrid of BPFs and Taylor polynomials are defined as ^[24]

where $i$ and $j$ are the order of BPFs and the degree of Taylor polynomials respectively.

3.2.   Function approximation

Let $H = C^{n+1}[0, T), \; \{h_{10}(t), h_{11}(t), \cdots, h_{N(M-1)}(t)\}\subset H,$ be a set of HBT functions, $X = span\{h_{10}(t), \cdots, h_{1(M-1)}(t), h_{20}(t), \cdots, h_{2(M-1)}(t), \cdots, h_{N0}(t), \cdots, h_{N(M-1)}(t)\},$ and $f(t)$ be a function in $H.$ Since $X$ is a finite dimensional subspace of $H$, $f(t)$ has the unique best approximation in $X$, such as $f^*(t)\in X,$ that is

Since $f^*(t)\in X,$ there exist the unique coefficients $c_{10}, c_{11}, \cdots, c_{N(M-1)}$ such that

where ${\mathsf{T}}$ denotes transpose of vector or matrix, $C$ and $H(t)$ are $N\times M$ column vectors

and

To evaluate $C$, firstly let

Using Eq (3.8) one gets

where

Therefore,

or

where

and

is a matrix of order $NM\times NM$ and is given by

For HBT functions $D$ has the following form

where

and $T(t) = [T_0(t), T_1(t), \ldots, T_{M-1}(t)]^{\mathsf{T}}.$ Hence, $C^{\mathsf{T}}$ in Eq (3.8) is given by

4.   Fractional integral operational matrix of HBT

Here, the fractional integral operational matrix of HBT is derived. To simplify the derivation process, the HBT functions are first projected onto a set of BPFs, and the projection matrix is calculated. The HBT operational matrix is obtained by using matrix operations of the projection matrix and BPF operational matrix.

The R-L fractional integration operator of hybrid function vector $H(t)$ can be written as

where matrix $P^\alpha$ is the fractional integral operational matrix of HBT. To simplify the derivation process, the hybrid functions are projected onto a set of BPFs $\psi_i(t), i = 1, 2, \cdots, m$, and $m = N\times M$. Specifically, the hybrid function vector defined in Eq (3.10) can be expressed as

where $\Psi(t) = [\psi_1(t), \psi_2(t), \cdots, \psi_m(t)]^\mathrm{T}$ and $\Phi$ is the projection matrix which transform the hybrid functions onto BPFs.

In general, for arbitrary $M$ and $N$, one has

The expressions of $a_{ij}$, $b_{ij}$ and $d_{ij}$ are given in Appendix A. As an example, let $N = 2, M = 3$ and $T = 1$, one has

the function vector (4.4) can be expressed as

In this case,

where

Hence, we have

The projection matrix is derived in detail in Appendix A and omitted here. It is easy to see from the above example that the obtained projection matrix is in block-diagonal form. This elegant property is useful for simplifying the computation of the operational matrices.

From Eq (4.1) and Eq (4.2), one has

According to Ref.^[25],

where

is the BPFs operational matrix of fractional integration with $\xi_k = (k+1)^{\alpha+1}-2k^{\alpha+1}+(k-1)^{\alpha+1}.$ Substituting Eq (3.6) into Eq (4.6), one can get

Therefore,

5.   Error analysis

Denote the HBT approximation to function $f(t)$ at the level $l$ as $f^*(t) = \sum_{i = 1}^N\sum_{j = 0}^{M-1}c_{ij}h_{ij}(t),$ where $l$ is determined by $N$ and $M \; (l = l_{N, M}).$ We replace $f(t)$ with $f^*(t)$ in Eq (2.1) and call the resulting integral HBT approximation of the $\alpha-$ order R-L fractional integral of $f(t).$

then the absolute error between $I^\alpha f(t)$ and $I^\alpha f^*(t)$ is

Theorem 1. Let $f(t)\in C^M[0, 1)$, $f^*(t) = \sum_{i = 1}^N\sum_{j = 0}^{M-1}c_{ij}h_{ij}(t)$ is the best approximation of $f(t)$ in $X$, then

1). $|f(t)-f^*(t)|\leq\frac{T^M}{M!N^M}|f^{(M)}(t)|, \quad t\in[\frac{i-1}{N}T, \frac{i}{N}T), \quad i = 1, 2, \ldots, N.$

2). $\varepsilon_l\leq\frac{KT^{\alpha+M}}{M!N^M\Gamma(\alpha+1)}, \; \mid f^{(M)}(t)\mid\leq K, \quad K$ is a finite positive value.

Proof. $1).$ Using Taylor formula, $f(t)$ can be approximated in the $i$th interval $[\frac{i-1}{N}T, \frac{i}{N}T)$ with

the truncation error is

where $\xi_i$ lies between $\frac{i-1}{N}$ and $t$. Since the best approximation is unique ^[26], for all $t\in[\frac{i-1}{N}T, \frac{i}{N}T),$ the absolute error between $f(t)$ and $f^*(t)$ is given as

$2).$ The absolute error between $I^\alpha f(t)$ and $I^\alpha f^*(t)$ is

according to the assumption $\mid f^{(M)}(t)\mid\leq K,$ the absolute error between $I^\alpha f(t)$ and $I^\alpha f^*(t)$ can be estimated as

where $K$ is a finite positive value.

In order to verify the maximum absolute error arised by HBT approximation is smaller than the upper bound in theory given in Theorem 1, the function $f(t) = t$ is selected as an example. The analytic expression of $\alpha$ order fractional integral of $f(t) = t$ is given as

Using Eq (4.9) with $\alpha = 0.5, T = 1$ and $N = 2, M = 3,$ the HBT estimation of $(_0I^\alpha_t)(t)$ is obtained as

Using BPFs, $(_0I^\alpha_t)(t)$ is approximated as ^[27]

Table 1 presents the absolute errors obtained by the approximation of BPFs and HBT functions. The piecewise-polynomial property of HBT functions lead to a more accurate approximation of the fractional integral even with the equal number of basis functions as BPFs. Therefore, HBT approximation of R-L fractional integral is effective.

6.   Results

To show the effectiveness of our presented method, the fractional integral operational matrix of HBT is used to solve several FODEs. The solutions obtained by our proposed method are compared with their exact solutions, those by block-pulse operational matrix (BPOM) method and fractional Taylor operational matrix (FTOM) method ^[19].

6.1.   Example 1

A fractional Riccati equation as ^[28,29]

subject to the initial condition $x(0) = 0$ is considered. Let

together with the initial condition, one has

Since $H(t) = \Phi\Psi(t),$ from Eq (6.3) one has

Let

and using Eqs (3.3) and (3.4) one can obtain

Substituting Eqs(4.2), (6.2), (6.4) and (6.6) into Eq (6.1), one has

Equation (6.7) is a set of algebraic equations. In this paper, the Matlab function $\verb|fsolve|$ is used to solve Eq (6.7). The numerical solution, for $N = 4, M = 6$ is shown in Figure 1. When $\alpha = 1$, the exact solution of this equation is

Table 2 shows comparisons of the approximate solutions at different values of $t$ obtained by the present method with $N = 32, M = 6$, BPOM method with $N = 32$, FTOM method with $M = 6$ and the exact solution for $\alpha = 1$. In spite of the large amount of calculation, the solutions obtained by the present method are the same as exact solutions.

On the other hand, the numerical solutions obtained by the present method with $N = 4, M = 6$ for $\alpha = 0.7, 0.8, 0.9, 1$ are shown in Figure 1. It is demonstrate that the approximate solution obtained by the present method is in good agreement with the exact solution when $\alpha = 1$. Moreover, as indicated in ^[28], the solution continuously depends on the time-fractional derivative.

6.2.   Example 2

A fractional differential equation as

and

subject to $x(0) = x'(0) = 0$ is considered ^[30]. The exact solution of this equation is $x(t) = \frac{1}{3}t^3.$

Let

and take the initial states into consideration, one has

Since $H(t) = \Phi \Psi(t),$ from Eq (6.9) we have

Let

then

Expanding $f(t)$ onto HBT functions, one has

where $f^\mathrm{T}$ is a known constant vector. Substituting Eq (6.9)–(6.11) and Eq (6.14) into Eq (6.8), together with $H(t) = \Phi \Psi(t)$, then

Eq (6.16) is a set of algebraic equations, which is solved via Matlab function $\verb|fsolve|$. Figure 2 shows the numerical solution obtained by the present method with $N = 4, M = 6$ when $a = 1, b = 1, c = 1, e = 1, \alpha_1 = 0.333, \alpha_2 = 1.234$. It can be seen that the numerical solution is almost the same as the exact solution.

The comparison of absolute errors in $x(t)$ obtained by the BPOM method with $N = 4$, FTOM method with $M = 6$ and present method with $N = 4, M = 6$ at given points for different values of $N$ and $M$ are shown in Table 3. The computational results illustrate that, compared with BPOM method and FTOM method, the present method can get more accurate approximate solutions although it has a large amount of calculation. Moreover, for the present method, the error is smaller and smaller when $N$ and $M$ increase. Therefore for higher accuracy of the approximation, lager $N$ and $M$ are recommended.

7.   Conclusions

In this paper, a new numerical method is presented to solve FODEs by using the hybrid functions consisting of BPFs and Taylor polynomials. The fractional integral operational matrix of the hybrid functions is derived. The FODEs to be solved is transformed into a set of algebraic equations via the operational matrix. Finally, the numerical solutions of FODEs are get by solving the algebraic equations. Illustrative examples verify that the proposed algorithm can get more accurate numerical solutions of FODEs than BPOM method and FTOM method although it has a large amount of calculation.

Acknowledgments

This work is partially supported by Public Welfare Project of Zhejiang Province of China (LGN20C050002).

Conflict of interest

The authors declare no conflict of interest.

A.   The projection matrix of HBT functions

In this section, the derivation process of the projection matrix from HBT functions to BPFs is given. Note that Eq (3.10) is the HBT function vector. First, elements $h_{1j}(t)$ of $H(t)$ are considered and can be expanded into BPFs as follows:

where $j = 0, 1, \ldots, (M-1)$ and $m = N\times M.$ The coefficient $a_{(j+1)k}$ can be computed according to Eq (3.6)

Since $h_{1j}(t)$ is nonzero when $0\leq t < \frac{T}{N},$ and $\psi_k(t)$ is defined on the interval $[\frac{k-1}{NM}T, \frac{k}{NM}T), h_{1j}(t)\psi_k(t)$ is nonzero only when $1\leq k\leq M.$ Observing this fact, Eq (A.2) can be written as

After some manipulations, it is easy to obtain

Define

Then

Thus, when $0\leq t < \frac{1}{N}$

where

similarly, when $\frac{T}{N}\leq t < \frac{2}{N}T, \quad h_{2j}(t), j = 0, 1, \ldots, (M-1)$ can be also expanded into BPFs, and

where

Let $x = t-\frac{T}{N}$, one has

Rewriting $b_{(j+1)k}$ into vector form,

Then

similarly, when $\frac{N-1}{N}T\leq t < T,$

where $j = 0, 1, \ldots, (m-1)$ and $m = N\times M.$

And finally, take the same process as above, the matrix $\Phi_N$ can be obtained, where $\Phi_N = [d_{(j+1)k}]_{M\times NM}$

Then

Therefore