Euler Wavelet Method as a Numerical Approach for the Solution of Nonlinear Systems of Fractional Differential Equations

Abstract

In this paper, a numerical approach for solving systems of nonlinear fractional differential equations (FDEs) is presented Using the Euler wavelets technique and associated operational matrices for fractional integration, we try to solve those systems of FDEs. The method’s major objective is to transform the nonlinear FDE into a nonlinear system of algebraic equations that is straightforward to solve with matrix techniques. The Euler wavelets are constructed using Euler polynomials, which have fewer terms than most other polynomials used to construct other types of wavelets, therefore, using Euler wavelets for the numerical approach provides sparse operational matrices. Thanks to the sparsity of those operational matrices, the proposed numerical approach requires less computation and takes less time to evaluate. The approach described here is also applicable to systems of fractional differential equations with variable orders. To illustrate the strength and performance of the method, four numerical examples are provided.

1. Introduction

Fractional calculus has drawn growing attention for decades since it is essential to many fields of research and engineering. Due to the fact that fractional order differential operators are nonlocal operators, compared to integer order differential equations, fractional differential equations have the benefit of being able to characterize some natural physics processes and dynamic system processes better. The use of FDEs allows for the elegant modeling of numerous physics, applied mathematics, and engineering systems [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. However, for most of those cases, analytical solutions to those FDE systems are very hard to obtain. Therefore, there has been a considerable effort to obtain numerical solutions to ordinary and partial FDEs. As a result, a variety of numerical methods have been developed by many researchers recently. Those methods include the spline collocation method [16], B-spline method [17], block-pulse function method [18], Laplace transform method [19], variational iteration method [20], finite difference method [21], Adomian decomposition method [22], Mellin transform method [23], homotopy perturbation method [24], Laplace optimized decomposition method [25], predictor–corrector approach [26], fractional complex transform method [27], computational matrix method [28], spectral collocation methods [29], Legendre spectral-collocation method [30], block-by-block method [31], and finite element method [32].

Most of those methods are applied to equations having only one or a few fractional differential terms in them. Therefore, they become impractical and cumbersome for the multiterm or variable-order FDEs. What is more, the finite elements method approximates the differential equation whereas the finite differences and spectral methods approximate the solution of the differential equation. Spectral base functions are infinitely differentiable, though they are nonzero for all of the regions of interest. On the other hand, the base functions used in the finite differences and finite elements methods have compact support though they are not necessarily continuous. The wavelet methods have the upper hand compared to those methods such that they have both compact support and they are smooth functions.

The wavelet theory has received a lot of attention in recent years. There are a lot of wavelet methods proposed for numerical approximations such as Haar, Legendre, Chebyshev, Euler, sine-cosine, and Gegenbauer wavelets. With the use of the orthogonal basis of those wavelets, it is possible to reduce the current problem to a set of either linear or nonlinear algebraic equations. Large systems of algebraic equations may require more computational load. The simplest wavelet method is the Haar wavelet and thanks to their simplicity they are commonly preferred. The drawback for the Haar wavelets is their lower accuracy. Sine-cosine wavelets produce better results for periodic functions. Therefore, their applications are limited.

In this paper, we focus on the Euler wavelet approximation to explore the advantages of the method. Since the Euler wavelets are produced using Euler polynomials, and since Euler polynomials generally have fewer terms than other polynomials used for the other wavelet methods, they produce sparser operational matrices for the numerical approach. Sparse operational matrices are essential for fast evaluation and low computational requirements in the numerical approach, especially for the large collocation points used to increase accuracy.

The theory of elasticity, dynamics, fluid mechanics, oscillation, and quantum dynamics are only a few examples of the topics where the system of fractional differential equations is frequently utilized to model the system behavior. In this paper, we aim to numerically solve the system of FDEs in the form of:

$$D_{*}^{{\alpha }_i}{u}_i\left(x\right)=f_i\left(x,u_1,u_2,\dots ,u_n\right), u_{{}_i}^{\left(r\right)}\left(0\right)=c_i, 1\le i\le n, 0\le r\le \lceil {\alpha }_i\rceil$$

(1)

where $m_i-1<{\alpha }_i\le m_i, m_i\in {\mathbb{Z}}^{+}, D_{*}^{{\alpha }_i}$ is the Caputo fractional differential operator.

The paper is organized as follows. The brief definitions and basic properties of fractional calculus are presented in Section 2. In Section 3, Euler wavelets are introduced and the operational matrices for the numerical integration for the solution of (one) are constructed. Several numerical examples are presented in Section 4 to demonstrate the validity of the method. The paper concludes by stating the important results of the method in Section 5.

2. Fundamentals of the Fractional Calculus

Definition 1.

The Riemann–Liouville fractional integral operator of order $\alpha$ is given as

$$I^{\alpha }f\left(t\right)=\left\{\begin{array}{l}\frac{1}{\Gamma \left(\alpha \right)} {\displaystyle \underset{0}{\overset{t}{\int }} \frac{f\left(\tau \right)}{{\left(t-\tau \right)}^{1-\alpha }}d\tau }, \alpha >0, t>0\\ f\left(t\right) , \alpha =0\end{array}\right\}$$

(2)

For $\alpha \ge 0$, $\beta \ge 0$, $a\ge 0$ and $\eta \ge -1$ we have the following properties of the Riemann–Liouville fractional integral

$$\left(\mathrm{i}\right) I^{\alpha }{I}^{\beta }=I^{\beta }{I}^{\alpha }$$

(3)

$$\left(\mathrm{ii}\right) I^{\alpha }\left(I^{\beta }f\left(t\right)\right)=I^{\beta }\left(I^{\alpha }f\left(t\right)\right)=I^{\alpha +\beta }f\left(t\right)$$

(4)

$$(\mathrm{iii}) I^{\alpha }{\left(t-a\right)}^{\eta }=\frac{\Gamma \left(\eta +1\right)}{\Gamma \left(\alpha +\eta +1\right)}{\left(t-a\right)}^{\alpha +\eta }$$

(5)

The Riemann–Liouville fractional derivative is defined by

$$D^{\alpha }f\left(t\right)={\left(\frac{d}{dt}\right)}^n\left(I^{n-\alpha }f\left(t\right)\right), 0\le n-1<\alpha \le n$$

(6)

where n is an integer and t > 0. When representing real-world processes, the derivative of the Riemann–Liouville operator has a few drawbacks. As a result, in this study, we employ Caputo’s modified fractional differential operator $D^{\alpha }$, which is described in the formulation that follows.

Definition 2.

The following expression is the fractional derivative operator defined by Caputo:

$$D^{\alpha }f\left(t\right)=\left\{\begin{array}{l}\frac{d^{n}f\left(t\right)}{d{t}^n} , \alpha =n\in R\\ \frac{1}{\Gamma \left(n-\alpha \right)} {\displaystyle \underset{0}{\overset{t}{\int }} \frac{f^{(n)}\left(t\right)}{{\left(t-\tau \right)}^{1-n+\alpha }}}d\tau , 0\le n-1<\alpha \le n\end{array}\right\}$$

(7)

The following two common equations can be used to express the relationship between the Riemann–Liouville operator and the Caputo operator:

$$D^{\alpha }{I}^{\alpha }f\left(t\right)=f\left(t\right)$$

(8)

and

$$I^{\alpha }{D}^{\alpha }f\left(t\right)=f\left(t\right)-{\displaystyle \sum _{k=0}^{n-1}{f}^{(k)}}\left(0^{+}\right)\frac{t^k}{k!}$$

(9)

For more information on fractional differentiation and integration, the reader is directed to [5].

3. Euler Wavelets and Derivation of Operational Matrices for Euler Wavelets

3.1. Euler Wavelets

Localized wavelike functions known as “wavelets” are used in wavelet analysis. A mother wavelet and its enlarged and translated variations make up a family of wavelets. We may create the following family of continuous wavelets as [33] by continually varying the translation parameter b and the dilation parameter a.

$${\psi }_{a,b}\left(t\right)={\left|a\right|}^{-1/2}\psi \left(\frac{t-b}{a}\right), a, b\in R, a\ne 0$$

(10)

When the parameters for translation and dilation are chosen to have discrete values $a=a_0^{-k}$, $b=n{b}_0{a}_0^{-k}$, $a_0>1, b_0>0$, we obtain the discrete wavelets as

$${\psi }_{kn}\left(t\right)={\left|a_0\right|}^{k/2}\psi \left(a_0^{k}t-n{b}_0\right)$$

(11)

Euler wavelets ${\psi }_{nm}=\psi \left(k,\tilde{n},m,t\right)$ have four arguments: $\tilde{n}=n-1$, $n=1,2,3,\cdots ,2^{k-1}$, m is the order for Euler polynomials, t is the normalized time, and k can take any positive integer value. Euler wavelets for $t\in \left[0,1\right)$ results

$${\psi }_{nm}\left(t\right)=\left\{\begin{array}{l}{2}^{\frac{k-1}{2}}{\tilde{E}}_m\left(2^{k-1}t-n+1\right), \frac{n-1}{2^{k-1}}\le t<\frac{n}{2^{k-1}}\\ 0 , otherwise\end{array}\right\}$$

(12)

and

$${\tilde{E}}_m\left(t\right)=\left\{\begin{array}{l} 1 , m=0\\ \frac{1}{\sqrt{\frac{2{\left(-1\right)}^{m-1}{\left(m!\right)}^2}{\left(2m\right)!}{E}_{2m+1}\left(0\right)}} , m>0\end{array}\right\}$$

(13)

where m = 0, 1, 2, …, M − 1 and $n=1,2,3,\dots ,2^{k-1}$ and $E_m\left(t\right)$ are the Euler polynomials defined by [34]

$${\displaystyle \sum _{k=0}^m\left(\begin{array}{c}m\\ k\end{array}\right)E_k\left(t\right)+E_m\left(t\right)}=2t^m$$

(14)

The Euler polynomials start with the polynomials given in (15)

$$E_0\left(t\right)=1, E_1\left(t\right)=t-\frac{1}{2}, E_2\left(t\right)=t^2-t, E_3\left(t\right)=t^3-\frac{2}{3}{t}^2+\frac{1}{4},\dots$$

(15)

3.2. Function Approximation

Euler wavelets can be used to represent a function $f\left(t\right)$, $t\in \left[0,1\right)$ as in the following:

$$f\left(t\right)={\displaystyle \sum _{n=1}^{2^{k-1}}{\displaystyle \sum _{m=0}^{M-1}{c}_{nm}}}{\Psi }_{nm}\left(t\right)=C^T\psi \left(t\right)$$

(16)

where T represents transposition and

$$C={\left[c_{10}, c_{11}, \dots c_{1(M-1)}, c_{20}, c_{21}, \dots c_{2(M-1)}\dots c_{2^{k-1}0}, c_{2^{k-1}1}, \dots c_{2^{k-1}(M-1)}\right]}^T$$

(17)

$$\Psi ={\left[{\Psi }_{10}, {\Psi }_{11}, \dots {\Psi }_{1(M-1)}, {\Psi }_{20}, {\Psi }_{21}, \dots {\Psi }_{2(M-1)}\dots {\Psi }_{2^{k-1}0}, {\Psi }_{2^{k-1}1}, \dots {\Psi }_{2^{k-1}(M-1)}\right]}^T$$

(18)

Here, we define the total number of discrete collocation points $t_i=\frac{i-0.5}{m^{\prime }}$ $i=1,2,3,\dots ,m^{\prime }$ as $m^{\prime }=2^{k-1}M$. Using those points, the Euler wavelet matrix ${\varphi }_{m^{\prime }x{m}^{\prime }}$ becomes

$${\varphi }_{m^{\prime }x{m}^{\prime }}=\left[\Psi \left(t_1\right) \Psi \left(t_2\right) \Psi \left(t_3\right) \cdots \Psi \left(t_{m^{\prime }}\right)\right]$$

(19)

The Euler wavelet matrix for k = 2, M = 3, and $\alpha =0.5$ yields

$${\varphi }_{m^{\prime }x{m}^{\prime }}=\left[\begin{array}{l} 1.4142 1.4142 1.4142 0 0 0\\ -0.9428 0 0.9428 0 0 0\\ -0.4811 -0.8660 -0.4811 0 0 0\\ 0 0 0 1.4142 1.4142 1.4142\\ 0 0 0 -0.9428 0 0.9428\\ 0 0 0 -0.4811 -0.8660 -0.4811\end{array}\right]$$

(20)

3.3. Euler Wavelet Operational Matrix of Fractional Integration

Block Pulse Functions

The block pulse functions (BPFs) of set $m^{\prime }$ are given by

$$b_i\left(t\right)=\left\{\begin{array}{l} 1 , \frac{i-1}{m^{\prime }}\le t<\frac{i}{m^{\prime }}\\ 0 , otherwise\end{array}\right\}$$

(21)

where $i=1,2,3,\dots ,m^{\prime }$. The functions $b_i\left(t\right)$ are both orthogonal and disjoint. For $t\in \left[0,1\right)$

$$b_i\left(t\right)b_j\left(t\right)=\left\{\begin{array}{l} 0 , i\ne j\\ b_i\left(t\right) , i=j\end{array}\right\}$$

(22)

$${\displaystyle \underset{0}{\overset{1}{\int }}{b}_i\left(\tau \right)b_j\left(\tau \right)d\tau }=\left\{\begin{array}{l} 0 , i\ne j\\ 1/m^{\prime } , i=j\end{array}\right\}$$

(23)

We can use the $m^{\prime }$ set of BPFs to represent a function $f\left(t\right)$, $t\in \left[0,1\right)$ provided that $f\left(t\right)$ is square integrable in the represented interval such that

$$f\left(t\right)={\displaystyle \sum _{i=1}^{m^{\prime }}{f}_i} b_i\left(t\right)=f^T{B}_{m^{\prime }}\left(t\right)$$

(24)

where $f={\left[f_1, f_2, \dots ,f_{m^{\prime }}\right]}^T$, $B_{m^{\prime }}\left(t\right)={\left[b_1\left(t\right), b_2\left(t\right), \dots ,b_{m^{\prime }}\left(t\right)\right]}^T$ and $f_i$ are

$$f_i=\frac{1}{m^{\prime }}{\displaystyle \underset{\left(i-1\right)/m^{\prime }}{\overset{i/m^{\prime }}{\int }}f\left(t\right) b_i\left(t\right)dt}$$

(25)

Similarly, the Euler wavelet matrix can also be represented with an $m^{\prime }$ set of BPFs as

$$\psi \left(t\right)={\varphi }_{m^{\prime }x{m}^{\prime }} B_{m^{\prime }}\left(t\right)$$

(26)

Now, we use the definition of the block pulse operational matrix for fractional integration $F^{\alpha }$ given in [35]

$$I^{\alpha }{B}_{m^{\prime }}\left(t\right)\approx F^{\alpha } B_{m^{\prime }}\left(t\right)$$

(27)

where

$$F^{\alpha }=\frac{1}{m^{\alpha }} \frac{1}{\Gamma \left(\alpha +2\right)}\left[\begin{array}{l}1 {\xi }_1 {\xi }_2 {\xi }_3\cdots {\xi }_{m^{\prime }-1}\\ 0 1 {\xi }_1 {\xi }_2\cdots {\xi }_{m^{\prime }-2}\\ 0 0 1 {\xi }_1\cdots {\xi }_{m^{\prime }-3}\\ ⋮ ⋮ \ddots \ddots ⋮ ⋮\\ 0 0 \cdots 0 1 {\xi }_1\\ 0 0 \cdots 0 0 1\end{array}\right]$$

(28)

with ${\xi }_k={\left(k+1\right)}^{\alpha +1}-2k^{\alpha +1}+{\left(k-1\right)}^{\alpha +1}$. Using (26)–(27) we can write

$$I^{\alpha }\psi \left(t\right)\approx I^{\alpha }{\varphi }_{m^{\prime }x{m}^{\prime }} B_{m^{\prime }}\left(t\right)={\varphi }_{m^{\prime }x{m}^{\prime }} I^{\alpha }{B}_{m^{\prime }}\left(t\right)\approx {\varphi }_{m^{\prime }x{m}^{\prime }}{F}^{\alpha }{B}_{m^{\prime }}\left(t\right)\approx {\varphi }_{m^{\prime }x{m}^{\prime }}{F}^{\alpha }{\varphi }_{m^{\prime }x{m}^{\prime }}^{-1}\psi \left(t\right)$$

(29)

$$I^{\alpha }\psi \left(t\right)\approx P_{m^{\prime }\times m^{\prime }}^{\alpha } \psi \left(t\right)$$

(30)

$$P_{m^{\prime }x{m}^{\prime }}^{\alpha }\approx {\varphi }_{m^{\prime }x{m}^{\prime }}{F}^{\alpha }{\varphi }_{m^{\prime }x{m}^{\prime }}^{-1}$$

(31)

where matrix $P_{m^{\prime }\times m^{\prime }}^{\alpha }$ is called the Euler wavelet operational matrix of fractional integration. As an example, the operational matrix for $\alpha =0.5$, k = 2, M = 3 yields:

$$P_{m^{\prime }\times m^{\prime }}^{\alpha }=\left[\begin{array}{l} 0.4616 0.3150 -0.1631 0.5012 -0.1509 0.1404\\ 0.0878 0.2243 0.4203 0.0717 -0.0449 0.0626\\ -0.1305 -0.1591 0.2354 -0.2110 0.0615 -0.0545\\ 0 0 0 0.4616 0.3150 -0.1631\\ 0 0 0 0.0878 0.2243 0.4203\\ 0 0 0 -0.1305 -0.1591 0.2354\end{array}\right]$$

(32)

Let us point out that, considering the two vectors given as $P={\left[p_1,p_2,\dots p_{m^{\prime }}\right]}^T$, $R={\left[r_1,r_2,\dots r_{m^{\prime }}\right]}^T$ we have $P*R={\left[p_1{r}_1,p_2{r}_2,\dots p_{m^{\prime }}{r}_{m^{\prime }}\right]}^T$ and $P^n={\left[p_1^n,p_2^n,\dots p_{m^{\prime }}^n\right]}^T$ for the multiplication and the nth power of those vectors thanks to the properties of the BPFs.

The convergence analysis of the Euler wavelets can be found in [33].

4. Numerical Examples

Here, four nonlinear systems of FDEs are solved using the proposed method to point out the accuracy and efficiency of the method.

4.1. Example 1

First, consider the system of FDE defined as:

$$\begin{array}{l}{D}^{\alpha }u\left(t\right)=\frac{3}{4}{v}^2\left(t\right) , 0<\alpha \le 1\\ D^{\beta }v\left(t\right)=u\left(t\right)v\left(t\right)-\frac{v^4\left(t\right)}{8}+2, 0<\beta \le 1\end{array}$$

(33)

With the initial values $u\left(0\right)=0$, $v\left(0\right)=0$ and the exact solution for $\alpha =1$ and $\beta =1$ is given as $u_{ex}\left(t\right)=t^3$, $v_{ex}\left(t\right)=2t$, $t\in \left[0,1\right]$.

Now, applying the proposed method to the fractional derivatives, we have

$$\begin{array}{l}{D}^{\alpha }u\left(t\right)\approx R_{m^{\prime }}^T \psi \left(t\right)\\ D^{\beta }v\left(t\right)\approx S_{m^{\prime }}^T \psi \left(t\right)\end{array}$$

(34)

where $R_{m^{\prime }}^T =\left[r_1,r_2,\dots ,r_{m^{\prime }}\right]$ and $S_{m^{\prime }}^T =\left[s_1,s_2,\dots ,s_{m^{\prime }}\right]$ are the unknown coefficients. Using (9), (26), (30), (34), and the initial conditions we obtain

$$\begin{array}{l}u\left(t\right)=I^{\alpha }{D}^{\alpha }u\left(t\right)+u\left(0\right)\approx R_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\alpha }\psi \left(t\right)\approx \underset{H_{m^{\prime }}^T}{\underbrace{R_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\alpha }{\varphi }_{m^{\prime }\times m^{\prime }}}}{B}_{m^{\prime }}\left(t\right)\\ v\left(t\right)=I^{\beta }{D}^{\beta }v\left(t\right)+v\left(0\right)\approx S_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\beta }\psi \left(t\right)\approx \underset{K_{m^{\prime }}^T}{\underbrace{S_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\beta }{\varphi }_{m^{\prime }\times m^{\prime }}}}{B}_{m^{\prime }}\left(t\right)\end{array}$$

(35)

where $H_{m^{\prime }}^T =\left[h_1,h_2,\dots ,h_{m^{\prime }}\right]$ $K_{m^{\prime }}^T =\left[k_1,k_2,\dots ,k_{m^{\prime }}\right]$ are also the vectors of size $1\times m^{\prime }$

$$\begin{array}{l}{v}^4\left(t\right)\approx {\left(K_{m^{\prime }}^T \right)}^4{B}_{m^{\prime }}\left(t\right)\\ u\left(t\right)v\left(t\right)\approx \left(H_{m^{\prime }}^T *K_{m^{\prime }}^T \right)B_{m^{\prime }}\left(t\right)\end{array}$$

(36)

Finally, using (34)–(36) in (33) we obtain the system of algebraic equations with $2m^{\prime }$ unknowns in the form of $R_{m^{\prime }}^T$ and $S_{m^{\prime }}^T$ coefficients:

$$\begin{array}{l}{R}_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=\frac{3}{4}{\left(K_{m^{\prime }}^T \right)}^2,\\ S_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=H_{m^{\prime }}^T *K_{m^{\prime }}^T-\frac{1}{8}{\left(K_{m^{\prime }}^T \right)}^4+{\left[2,2,\dots ,2\right]}_{1\times m^{\prime }}\end{array}$$

(37)

Solving (37) for the $R_{m^{\prime }}^T$ and $S_{m^{\prime }}^T$ coefficients also provides the numerical solution for $u\left(t\right)$ and $v\left(t\right)$, as indicated in (35).

Table 1. Absolute errors $E_u$ and $E_v$ corresponding to $u\left(t\right)$ and $v\left(t\right)$ for several $m^{\prime }$ values of example 1.

t	${\mathit{m}}^{\prime }=12$		${\mathit{m}}^{\prime }=24$		${\mathit{m}}^{\prime }=48$		${\mathit{m}}^{\prime }=96$
	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$
0	1.09 × 10−3	2.51 × 10−6	1.36 × 10−4	7.85 × 10−8	1.70 × 10−5	2.45 × 10−9	2.12 × 10−6	7.67 × 10−11
0.1	1.89 × 10−4	2.75 × 10−6	9.92 × 10−5	6.36 × 10−7	2.42 × 10−5	1.48 × 10−7	5.23 × 10−6	3.64 × 10−8
0.2	7.94 × 10−4	2.03 × 10−5	1.93 × 10−4	4.74 × 10−6	4.19 × 10−5	1.16 × 10−6	1.05 × 10−5	2.90 × 10−7
0.3	9.47 × 10−4	6.46 × 10−5	2.42 × 10−4	1.57 × 10−5	6.69 × 10−5	3.91 × 10−6	1.66 × 10−5	9.75 × 10−7
0.4	1.56 × 10−3	1.51 × 10−4	3.38 × 10−4	3.70 × 10−5	8.52 × 10−5	9.21 × 10−6	2.21 × 10−5	2.30 × 10−6
0.5	2.87 × 10−3	2.89 × 10−4	5.81 × 10−4	7.14 × 10−5	1.28 × 10−4	1.78 × 10−5	2.99 × 10−5	4.45 × 10−6
0.6	2.04 × 10−3	4.88 × 10−4	5.60 × 10−4	1.21 × 10−4	1.39 × 10−4	3.03 × 10−5	3.40 × 10−5	7.56 × 10−6
0.7	2.76 × 10−3	7.52 × 10−4	6.84 × 10−4	1.87 × 10−4	1.64 × 10−4	4.67 × 10−5	4.12 × 10−5	1.17 × 10−5
0.8	3.12 × 10−3	1.08 × 10−3	7.83 × 10−4	2.68 × 10−4	2.02 × 10−4	6.69 × 10−5	5.04 × 10−5	1.67 × 10−5
0.9	4.05 × 10−3	1.44 × 10−3	9.57 × 10−4	3.60 × 10−4	2.40 × 10−4	8.99 × 10−5	6.07 × 10−5	2.25 × 10−5

summarizes the absolute errors obtained from the proposed method for $u\left(t\right)$ and $v\left(t\right)$ for several $m^{\prime }$ values ($\alpha =1$, $\beta =1$). $E_u$ and $E_v$ represent the absolute errors in $u\left(t\right)$ and $v\left(t\right)$, respectively. We present the solution graphs $u\left(t\right)$ and $v\left(t\right)$ for integer orders with the exact solution in Figure 1 and the fractional orders $\alpha =0.7, 0.8, 0.9$, $\beta =0.7, 0.8, 0.9$ with the integer orders $\alpha =1$, $\beta =1$ in Figure 2. As Table 1 and Figure 1 suggest, the proposed method follows the exact solution closely for the integer-orders of $\alpha =1$, $\beta =1$. Also, as can be seen from the table, the absolute errors decrease with the larger $m^{\prime }$ values, as expected. The absolute errors are approximately on the order of 10⁻⁴ for $m^{\prime }=24$, on the order of 10⁻⁵ for $m^{\prime }=48$, and on the order of 10⁻⁵–10⁻⁶ for $m^{\prime }=96$. In Figure 2, the fractional orders $\alpha ,\beta$ are chosen to be equal, though they can be chosen arbitrarily on the interval $\left(0,1\right]$. Several different $\alpha ,\beta$ values are chosen in the following examples. As Figure 2 demonstrates, the solution approaches to the exact solution when $\alpha ,\beta$ approach one. As a result, Table 1 and Figure 1 and Figure 2 suggest that the proposed method is a valid approximation to the system of FDE in question.

4.2. Example 2

For the second example, consider the system of FDEs given as:

$$\begin{array}{l}{D}^{\alpha }u\left(t\right)=-1002u\left(t\right)+1000v^2\left(t\right), 0<\alpha \le 1\\ D^{\beta }v\left(t\right)=u\left(t\right)-v\left(t\right)-v^2\left(t\right) , 0<\beta \le 1\end{array}$$

(38)

where the initial values are given as $u\left(0\right)=1$, $v\left(0\right)=1$ and the exact solution for $\alpha =1$, $\beta =1$ becomes $u_{ex}\left(t\right)=e^{-2t}$, $v_{ex}\left(t\right)=e^{-t}$, $t\in \left[0,1\right]$.

As in Example 1, applying the proposed method to the fractional derivatives using (9), (26), (30), (34) and initial conditions for this system of FDEs yield

$$\begin{array}{l}u\left(t\right)=I^{\alpha }{D}^{\alpha }u\left(t\right)+u\left(0\right)\approx R_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\alpha }\psi \left(t\right)+1\approx \underset{H_{m^{\prime }}^T }{\underbrace{R_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\alpha }{\varphi }_{m^{\prime }\times m^{\prime }}}}{B}_{m^{\prime }}\left(t\right)+1\\ v\left(t\right)=I^{\beta }{D}^{\beta }v\left(t\right)+v\left(0\right)\approx S_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\beta }\psi \left(t\right)+1\approx \underset{K_{m^{\prime }}^T}{\underbrace{S_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\beta }{\varphi }_{m^{\prime }\times m^{\prime }}}}{B}_{m^{\prime }}\left(t\right)+1\end{array}$$

(39)

and

$$v^2\left(t\right)\approx {\left(K_{m^{\prime }}^T +1\right)}^2{B}_{m^{\prime }}\left(t\right)\approx \left[{\left(K_{m^{\prime }}^T\right)}^2 +2K_{m^{\prime }}^T +1\right]B_{m^{\prime }}\left(t\right)$$

(40)

Finally, substituting (34), (39), and (40) into (38) we obtain the system of algebraic equations with $2m^{\prime }$ unknowns in the form of $R_{m^{\prime }}^T$ and $S_{m^{\prime }}^T$coefficients:

$$\begin{array}{l}{R}_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=-1002\left(H_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\right)+1000\left[{\left(K_{m^{\prime }}^T\right)}^2 +2K_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\right],\\ S_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=\left(H_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\right)-\left(K_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\right)-\left[{\left(K_{m^{\prime }}^T\right)}^2 +2K_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\right]\end{array}$$

(41)

$$\begin{array}{l}{R}_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=-1002H_{m^{\prime }}^T +1000{\left(K_{m^{\prime }}^T \right)}^2+2000K_{m^{\prime }}^T -{\left[2,2,\dots ,2\right]}_{1\times m^{\prime }},\\ S_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=H_{m^{\prime }}^T -{\left(K_{m^{\prime }}^T \right)}^2-3K_{m^{\prime }}^T -{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\end{array}$$

(42)

Solving (42) for the $R_{m^{\prime }}^T$ and $S_{m^{\prime }}^T$ coefficients also provides the numerical solution for $u\left(t\right)$ and $v\left(t\right)$, as indicated in (35).

Similar to the first example results,

Table 2. Absolute errors $E_u$ and $E_v$ corresponding to $u\left(t\right)$ and $v\left(t\right)$ for several $m^{\prime }$ values of Example 2.

t	${\mathit{m}}^{\prime }=12$		${\mathit{m}}^{\prime }=24$		${\mathit{m}}^{\prime }=48$		${\mathit{m}}^{\prime }=96$
	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$
0	6.71 × 10−4	7.06 × 10−4	3.35 × 10−4	1.96 × 10−4	1.17 × 10−4	5.16 × 10−5	3.70 × 10−5	1.32 × 10−5
0.1	1.46 × 10−3	7.59 × 10−4	3.27 × 10−4	1.82 × 10−4	8.40 × 10−5	4.55 × 10−5	2.08 × 10−5	1.15 × 10−5
0.2	9.32 × 10−4	6.03 × 10−4	2.45 × 10−4	1.52 × 10−4	6.32 × 10−5	3.88 × 10−5	1.60 × 10−5	9.68 × 10−6
0.3	8.05 × 10−4	5.27 × 10−4	1.97 × 10−4	1.31 × 10−4	4.68 × 10−5	3.20 × 10−5	1.17 × 10−5	8.01 × 10−6
0.4	5.01 × 10−4	4.10 × 10−4	1.44 × 10−4	1.08 × 10−4	3.70 × 10−5	2.70 × 10−5	8.81 × 10−6	6.65 × 10−6
0.5	8.53 × 10−5	2.51 × 10−4	6.64 × 10−5	7.47 × 10−5	1.89 × 10−5	2.03 × 10−5	5.62 × 10−6	5.28 × 10−6
0.6	3.51 × 10−4	3.00 × 10−4	7.60 × 10−5	7.03 × 10−5	1.86 × 10−5	1.76 × 10−5	4.96 × 10−6	4.49 × 10−6
0.7	2.08 × 10−4	2.21 × 10−4	5.39 × 10−5	5.58 × 10−5	1.47 × 10−5	1.45 × 10−5	3.66 × 10−6	3.62 × 10−6
0.8	1.70 × 10−4	1.89 × 10−4	4.34 × 10−5	4.70 × 10−5	9.75 × 10−6	1.13 × 10−5	2.46 × 10−6	2.82 × 10−6
0.9	9.89 × 10−5	1.30 × 10−4	2.94 × 10−5	3.61 × 10−5	7.67 × 10−6	8.98 × 10−6	1.74 × 10−6	2.19 × 10−6

summarizes the absolute errors obtained from the proposed method for $u\left(t\right)$ and $v\left(t\right)$ for several $m^{\prime }$ values ($\alpha =1$, $\beta =1$). $E_u$ and $E_v$ represent the absolute errors in $u\left(t\right)$ and $v\left(t\right)$, respectively. As can be seen from the table, the absolute errors decrease with the larger $m^{\prime }$ values, as expected. The exact and proposed method solutions for $\alpha =1$, $\beta =1$ are also plotted in Figure 3. The absolute errors are roughly on the order of 10⁻⁴–10⁻⁵ for $m^{\prime }=24$, on the order of 10⁻⁵ for $m^{\prime }=48$, and on the order of 10⁻⁶ for $m^{\prime }=96$. Table 2 demonstrates that the proposed method follows the exact solution closely for the integer-orders of $\alpha =1$ and $\beta =1$. The solution graphs $u\left(t\right)$ and $v\left(t\right)$ for several fractional orders $\alpha ,\beta$ with integer orders $\alpha =1$, $\beta =1$ are given in Figure 4. Here, the fractional orders are chosen to be not equal. Numerical simulation is done for the fractional order pairs $\alpha =0.4, \beta =0.5$, $\alpha =0.6, \beta =0.7$, $\alpha =0.8, \beta =0.9$ and integer orders. As Figure 4 demonstrates, the solution approaches to the exact solution when $\alpha ,\beta$ approach one, which thus validates the fractional solutions.

4.3. Example 3

Here, we consider the following system of FDEs

$$\begin{array}{l}{D}^{\alpha }u\left(t\right)=\frac{u\left(t\right)}{2} , 0<\alpha \le 1\\ D^{\beta }v\left(t\right)=u^2\left(t\right)+v\left(t\right), 0<\beta \le 1\end{array}$$

(43)

with $u\left(0\right)=1$, $v\left(0\right)=0$. The exact solution for the integer orders $\alpha =1$, $\beta =1$ is given as $u_{ex}\left(t\right)=e^{t/2}$, $v_{ex}\left(t\right)=t{e}^t$, $t\in \left[0,1\right]$.

Applying the proposed method to this system of FDEs by following the same steps of the previous two examples results in the following system of algebraic equations (44), which is then solved for the unknown coefficients $R_{m^{\prime }}^T$ and $S_{m^{\prime }}^T$.

$$\begin{array}{l}{R}_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=\frac{1}{2}\left(H_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\right),\\ S_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}={\left(H_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\right)}^2-\left(K_{m^{\prime }}^T \right)\\ ={\left(H_{m^{\prime }}^T \right)}^2+2H_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}-K_{m^{\prime }}^T\end{array}$$

(44)

The results obtained for $u\left(t\right)$ and $v\left(t\right)$ for several $m^{\prime }$ values ($\alpha =1$, $\beta =1$) are summarized in

Table 3. Output values $u\left(t\right)$ and $v\left(t\right)$ of the proposed method for $\alpha =1$, $\beta =1$ and several $m^{\prime }$ values with the exact solution of example 3.

t	Exact Solution		${\mathit{m}}^{\prime }=24$		${\mathit{m}}^{\prime }=48$		${\mathit{m}}^{\prime }=96$		${\mathit{m}}^{\prime }=192$
	$\mathit{u}(\mathit{t})$	$\mathit{v}(\mathit{t})$	$\mathit{u}(\mathit{t})$	$\mathit{v}(\mathit{t})$	$\mathit{u}(\mathit{t})$	$\mathit{v}(\mathit{t})$	$\mathit{u}(\mathit{t})$	$\mathit{v}(\mathit{t})$	$\mathit{u}(\mathit{t})$	$\mathit{v}(\mathit{t})$
0	1.0000	0.0000	1.0001	0.0005	1.0000	0.0001	1.0000	0.0000	1.0000	0.0000
0.1	1.0513	0.1105	1.0513	0.1111	1.0513	0.1107	1.0513	0.1106	1.0513	0.1105
0.2	1.1052	0.2443	1.1052	0.2450	1.1052	0.2444	1.1052	0.2443	1.1052	0.2443
0.3	1.1618	0.4050	1.1619	0.4058	1.1619	0.4052	1.1618	0.4050	1.1618	0.4050
0.4	1.2214	0.5967	1.2215	0.5977	1.2214	0.5970	1.2214	0.5968	1.2214	0.5967
0.5	1.2840	0.8244	1.2841	0.8257	1.2840	0.8247	1.2840	0.8244	1.2840	0.8244
0.6	1.3499	1.0933	1.3499	1.0947	1.3499	1.0936	1.3499	1.0934	1.3499	1.0933
0.7	1.4191	1.4096	1.4192	1.4114	1.4191	1.4101	1.4191	1.4097	1.4191	1.4097
0.8	1.4918	1.7804	1.4919	1.7825	1.4919	1.7810	1.4918	1.7806	1.4918	1.7805
0.9	1.5683	2.2136	1.5684	2.2161	1.5683	2.2143	1.5683	2.2138	1.5683	2.2137

. As can be seen from the table, the outputs obtained in the proposed method are in agreement with the exact solution, and the larger the $m^{\prime }$ value, the better the numerical approximation becomes. The exact and proposed method solutions for $\alpha =1$, $\beta =1$ are also plotted in Figure 5, whereas the solution graphs $u\left(t\right)$ and $v\left(t\right)$ for the fractional orders $\alpha ,\beta$ with integer orders are given in Figure 6. Again, the fractional orders are chosen to be not equal in the simulation. The numerical simulation is done for the fractional order pairs $\alpha =0.25, \beta =0.35$, $\alpha =0.5, \beta =0.6$, and $\alpha =0.75, \beta =0.85$, and the integer orders. As Figure 5 and Figure 6 demonstrate, the solution approaches to the exact solution when $\alpha ,\beta$ approach one, which thus validates the fractional solutions.

4.4. Example 4

For our last example, we consider the following system of FDEs

$$\begin{array}{l}{D}^{\alpha }u\left(t\right)=u\left(t\right) , 0<\alpha \le 1\\ D^{\beta }v\left(t\right)=2u^2\left(t\right) , 0<\beta \le 1\\ D^{\gamma }w\left(t\right)=3u\left(t\right)v\left(t\right) , 0<\gamma \le 1\end{array}$$

(45)

with the initial conditions $u\left(0\right)=1$, $v\left(0\right)=1$, $w\left(0\right)=1$ and the exact solution for $\alpha =1$, $\beta =1$ and $\gamma =1$ is given as $u_{ex}\left(t\right)=e^t$, $v_{ex}\left(t\right)=e^{2t}$, and $w_{ex}\left(t\right)=e^{3t}$, $t\in \left[0,1\right]$.

Applying the proposed method to this system of FDEs, we have:

$$\begin{array}{l}{D}^{\alpha }u\left(t\right)\approx R_{m^{\prime }}^T \psi \left(t\right)\\ D^{\beta }v\left(t\right)\approx S_{m^{\prime }}^T \psi \left(t\right)\\ D^{\gamma }w\left(t\right)\approx T_{m^{\prime }}^T \psi \left(t\right)\end{array}$$

(46)

where $R_{m^{\prime }}^T =\left[r_1,r_2,\dots ,r_{m^{\prime }}\right]$, $S_{m^{\prime }}^T =\left[s_1,s_2,\dots ,s_{m^{\prime }}\right]$, and $T_{m^{\prime }}^T =\left[t_1,t_2,\dots ,t_{m^{\prime }}\right]$ are the unknown coefficients.

$$\begin{array}{l}u\left(t\right)=I^{\alpha }{D}^{\alpha }u\left(t\right)+u\left(0\right)\approx R_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\alpha }\psi \left(t\right)+1\approx \underset{H_{m^{\prime }}^T}{\underbrace{R_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\alpha }{\varphi }_{m^{\prime }\times m^{\prime }}}}{B}_{m^{\prime }}\left(t\right)+1\\ v\left(t\right)=I^{\beta }{D}^{\beta }v\left(t\right)+v\left(0\right)\approx S_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\beta }\psi \left(t\right)+1\approx \underset{K_{m^{\prime }}^T}{\underbrace{S_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\beta }{\varphi }_{m^{\prime }\times m^{\prime }}}}{B}_{m^{\prime }}\left(t\right)+1\\ w\left(t\right)=I^{\gamma }{D}^{\gamma }w\left(t\right)+w\left(0\right)\approx T_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\gamma }\psi \left(t\right)+1\approx \underset{L_{m^{\prime }}^T}{\underbrace{T_{m^{\prime }}^T P_{m^{\prime }\times m^{\prime }}^{\gamma }{\varphi }_{m^{\prime }\times m^{\prime }}}}{B}_{m^{\prime }}\left(t\right)+1\end{array}$$

(47)

Substituting (46)–(47) in (45) we obtain the system of algebraic equations with $3m^{\prime }$ unknowns in the form of $R_{m^{\prime }}^T$, $S_{m^{\prime }}^T$, and $T_{m^{\prime }}^T$coefficients:

$$\begin{array}{l}{R}_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=H_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }},\\ S_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=2\left[{\left(H_{m^{\prime }}^T\right)}^2+2H_{m^{\prime }}^T +{\left[1,1,\dots ,1\right]}_{1\times m^{\prime }}\right],\\ T_{m^{\prime }}^T {\varphi }_{m^{\prime }\times m^{\prime }}=3\left(H_{m^{\prime }}^T *K_{m^{\prime }}^T \right)+3H_{m^{\prime }}^T +3K_{m^{\prime }}^T +{\left[3,3,\dots ,3\right]}_{1\times m^{\prime }}\end{array}$$

(48)

The absolute errors obtained from the proposed method for $u\left(t\right)$, $v\left(t\right)$, and $w\left(t\right)$ for $m^{\prime }=48$ and $m^{\prime }=96$ ($\alpha =1$, $\beta =1$, and $\gamma =1$) are summarized in

Table 4. Output values $u\left(t\right)$, $v\left(t\right)$, and $w\left(t\right)$ of the proposed method for $\alpha =1$, $\beta =1$, and $\gamma =1$, and several $m^{\prime }$ values with the exact solution of Example 4.

t	${\mathit{m}}^{\prime }=48$			${\mathit{m}}^{\prime }=96$
	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$	${\mathit{E}}_{\mathit{w}}$	${\mathit{E}}_{\mathit{u}}$	${\mathit{E}}_{\mathit{v}}$	${\mathit{E}}_{\mathit{w}}$
0	5.72 × 10−5	2.41 × 10−4	5.70 × 10−4	1.39 × 10−5	5.71 × 10−5	1.32 × 10−4
0.1	6.44 × 10−5	2.78 × 10−4	7.14 × 10−4	1.60 × 10−5	6.82 × 10−5	1.74 × 10−4
0.2	7.48 × 10−5	3.42 × 10−4	9.76 × 10−4	1.87 × 10−5	8.57 × 10−5	2.45 × 10−4
0.3	8.82 × 10−5	4.39 × 10−4	1.40 × 10−3	2.20 × 10−5	1.09 × 10−4	3.50 × 10−4
0.4	1.02 × 10−4	5.40 × 10−4	1.90 × 10−3	2.57 × 10−5	1.37 × 10−4	4.86 × 10−4
0.5	1.24 × 10−4	7.53 × 10−4	3.05 × 10−3	3.04 × 10−5	1.80 × 10−4	7.15 × 10−4
0.6	1.39 × 10−4	8.76 × 10−4	3.78 × 10−3	3.45 × 10−5	2.15 × 10−4	9.23 × 10−4
0.7	1.60 × 10−4	1.08 × 10−3	5.07 × 10−3	4.00 × 10−5	2.70 × 10−4	1.27 × 10−3
0.8	1.86 × 10−4	1.37 × 10−3	7.14 × 10−3	4.64 × 10−5	3.42 × 10−4	1.78 × 10−3
0.9	2.13 × 10−4	1.69 × 10−3	9.56 × 10−3	5.35 × 10−5	4.28 × 10−4	2.44 × 10−3

. $E_u$, $E_v$, and $E_w$ represent the absolute errors in $u\left(t\right)$, $v\left(t\right)$, and $w\left(t\right)$, respectively. The absolute errors are roughly on the order of 10⁻³–10⁻⁴ for $m^{\prime }=48$, and on the order of 10⁻⁴–10⁻⁵ for $m^{\prime }=96$. As can be seen from the table, the absolute errors decrease with the larger $m^{\prime }$ values, as expected. The numerical solutions and exact solutions of $u\left(t\right)$, $v\left(t\right)$, and $w\left(t\right)$ for the integer orders $\alpha =1$, $\beta =1$, and $\gamma =1$ are plotted in Figure 7. Figure 7 shows that the numerical solution is very close to the exact solution even for the relatively small value of $m^{\prime }=12$. Also, the numerical solutions of $u\left(t\right)$, $v\left(t\right)$, and $w\left(t\right)$ for the fractional orders $0.75, 0.85, 0.95$ with the integer orders one are plotted in Figure 8. As is the case with the previous examples, the solution approaches to the exact solution when $\alpha ,\beta ,\gamma$ approach one.

The above four examples demonstrate that the proposed algorithm is effective and accurate to apply to the systems of FDEs.

5. Conclusions

In this paper, we proposed a numerical solution method for the systems of FDEs. The proposed method employs discrete Euler wavelets to obtain the so-called operational matrices for the fractional integration. With the help of the block pulse functions and the small number of terms in the Euler polynomials, the operational matrices become very sparse. This is the essence of the proposed method since sparse operational matrices used for the mapping between fractional terms in the FDEs and the discrete terms in the numerical approach determine the required memory and also determine the speed of the numerical method. The system of FDEs is converted into a system of algebraic equations and the algebraic equation system is solved using the Newton–Raphson method to obtain the unknown coefficients. The accuracy of the proposed method increases for the larger number of collocation points, as expected. This is verified in each of the considered examples shown in Table 1, Table 2, Table 3 and Table 4 and Figure 1, Figure 3, Figure 5 and Figure 7. The maximum errors are generally on the order of 10⁻⁴–10⁻⁶ for the collocation points up to $m^{\prime }=96$. For higher accuracy, the number of the collocation points must be increased. The numerical solutions obtained for fractional orders, which are the main target of this study, are also in agreement with the solutions obtained for integer orders. One can see from Figure 2, Figure 4, Figure 6 and Figure 8 that the method provides precise solutions for the fractional orders owing to the fact that when the fractional orders approach the integer values, the solution also approaches the exact solution obtained for the integer orders.

For future work, the proposed method can also be applied to the systems of variable-order FDEs, systems of fractional partial differential equations, systems of fractional integral equations, and systems of fractional delay differential equations.