An approach to solve fractional optimal control problems via fractional-order Boubaker wavelets

Abstract

This work is devoted to numerical solution to solve fractional optimal control problems (FOCPs). The fractional-order Boubaker wavelets (FOBWs) are introduced. By applying hypergeometric functions, an exact formula for Riemann–Liouville fractional integral operator for FOBW is obtained. By using this formula, the FOCP is reduced into a corresponding optimization problem which can be solved by applying the Lagrange multiplier method. Also, the convergence for the present method is provided. Seven numerical examples are included to demonstrate the applicability of the proposed technique. In Example 6, it will be shown that we can achieve the exact solutions which were not obtained previously in the literature. In addition, in Example 7, a practical example in a cancer model is provided.

Keywords

Boubaker wavelet hypergeometric function fractional optimal control Riemann–Lioville integration numerical method

1. Introduction

These days several problems in engineering and science are modeled in fractional calculus. Fractional calculus plays a critical role in many areas, such as viscoelasticy, the area of signal processing, and electrochemistry (Baleanu et al., 2016; Cattani et al., 2017). Therefore, finding the solution of fractional calculus problems has received considerable interest of many researchers.

Optimal control problems are used to minimize (or maximize) a cost function subjected to certain criteria. In fractional optimal control problems (FOCPs), the cost function or the criteria contain fractional derivative (Tricaud and Chen, 2010). Some numerical methods for solving FOCPs are given in the literature, for example, eigen functions (Agrawal, 2008), Legendre polynomials (Lotfi et al., 2013), Bernoulli polynomials (Behroozifar and Habibi, 2017; Keshavarz et al., 2016), hybrid of block-pulse functions with Bernoulli and Chelyshkov polynomials (Mashayekhi and Razzaghi, 2018; Mohammadi et al., 2018), finite difference method (El-Kady, 2003), modified hat functions (Nemati et al., 2019), and direct transcription method (Salati et al., 2019).

In this paper, Boubaker wavelets are used to find the solution of FOCPs. Boubaker polynomials were presented in 2007. These functions were used for solving the Boltzmann diffusion equation Boubaker (2007, 2011). Love’s equation was solved by Kumar via the Boubaker polynomials in Kumar (2010). These polynomials were applied in Rabiei et al. (2017) for solving FOCPs. These polynomials were used in Rabiei and Ordokhani (2019) for delay system.

Wavelets have many applications in science and engineering. Some of the wavelets in the literature are CAS, Legendre, Bessel, Chebyshev, Haar, and Boubaker (Dehestani et al., 2020; Heydari et al., 2016; Li and Zhao, 2010; Rafiei et al., 2018; Rabiei and Ordokhani, 2020; Saeedi et al., 2011). For solving fractional calculus problems by using wavelet method, the following formula was applied

I^{β} Ψ (t) ≅ P^{β} Ψ (t)

(1)

In this formula, I^β is the Riemann–Liouville fractional integral operator (RLFIO) of order β, and P^β is the operational matrix of the wavelet Ψ (t). For the mentioned wavelets, P^β was calculated approximately. In 2016, the authors in Mashayekhi and Razzaghi (2016) calculated I^βΨ (t) exactly for the first time by using the hybrid of block-pulse functions and Bernoulli polynomials. This exact formula was later obtained by using Taylor wavelets Toan et al. (2019); Vichitkunakorn et al. (2020). For the problems for which the exact solutions contain fractional-order power terms, fractional-order Legendre, generalized Laguerre functions, and Bernoulli wavelets were introduced in Kazem et al. (2013); Bhrawy et al. (2014); Rahimkhani et al. (2016) to obtain a more accurate solution. However, none of these functions and wavelets obtained RLFIO exactly.

In this paper, fractional-order Boubaker wavelets are used to solve FOCPs. Boubaker polynomials are applied for constructing the Boubaker wavelets. Next, we include the advantages of Boubaker polynomials. Shifted Legendre polynomials have the most satisfactory numerical findings among orthogonal functions Razzaghi and Elnagar (1994). In comparison with the shifted Legendre polynomials, Boubaker polynomials have fewer terms. For example, the Boubaker polynomial of degree 6 has 4 terms, while the shifted Legendre polynomial of degree 6 has 7 terms, and this difference will be greater by increasing the degrees of these polynomials. Hence, for function approximation, less CPU time is consumed by applying Boubaker polynomials. Moreover, the operational matrices of Boubaker polynomials are more sparse than that of the shifted Legendre polynomials which decreases computational cost. In this paper, an exact formula in terms of the hypergeometric functions is constructed for the RLFIO of fractional-order Boubaker wavelets. It is noted that the hypergeometric functions used in this paper can be calculated effectively by Mathematica. We obtain very accurate numerical solutions using the exact formula which will be shown in several examples. We will show in Examples 1–5 that our method has more accurate results than existing methods. Also in Example 6, it can be seen that when the exact solutions are of the form of fractional power terms, then the exact solutions can be obtained. These exact solutions were not obtained in the literature before. In addition, in Example 7, a practical example in cancer model is provided. In Section 2, we recall some preliminaries from fractional calculus as well as some properties of hypergeometric functions. In Section 3, the fractional-order Boubaker wavelets are constructed. In Section 4, the exact values of the RLFIO of fractional-order Boubaker wavelets in terms of hypergeometric functions are obtained. In Section 5, our numerical method for solving FOCPs by using fractional-order Boubaker wavelets is provided. The convergence of the method is included in Section 6. Finally, in Section 7, several examples are given to show the accuracy and the applicability of the proposed method.

2. Preliminaries and notations

Definition 1

For μ > − 1, the RLFIO of a function f and order ν ≥ 0 is Podlubny (1999)

I^{ν} f (t) = {\begin{array}{l} \frac{1}{Γ (ν)} \int_{0}^{t} \frac{f (s)}{{(t - s)}^{1 - ν}} d s & ν > 0, t > 0 \\ f (t) & ν = 0 \end{array}

(2)

This operator has the following properties

1. I^ν (υ₁f₁ (t) + υ₂f₂ (t)) = υ₁I^νf₁ (t) + υ₂I^νf₂ (t), for constants υ₁ and υ₂,

2. $I^{ν} t^{β} = \frac{Γ (β + 1)}{Γ (β + ν + 1)} t^{ν + β}, β > - 1$ .

Definition 2

For $n - 1 < ν \leq n, n \in ℕ$ , the Caputo derivative of order ν of the function f is Podlubny (1999)

D^{ν} f (t) = {\begin{matrix} \frac{1}{Γ (n - ν)} \int_{0}^{t} \frac{f^{(n)} (s)}{{(t - s)}^{n - 1 - ν}} d s & i f ν \neq n \\ f^{(n)} (t) & i f ν = n \end{matrix}

(3)

This definition satisfies the following properties

1. D^νc = 0,

2. $I^{ν} D^{ν} f (t) = f (t) - \sum_{i = 0}^{n - 1} f^{(i)} (0) \frac{t^{i}}{i!}$

3. $D^{ν} t^{β} = {\begin{array}{l} 0 & ν \in N_{0}, β < ν \\ \frac{Γ (β + 1)}{Γ (β + 1 - ν)} t^{β - ν} & o t h e r w i s e \end{array}$

Definition 3

The hypergeometric function ${{}_{2}F}_{1} (a, b, c; z)$ is given by Andrews et al. (1999)

{{}_{2}F}_{1} (a, b, c; z) = \frac{Γ (c)}{Γ (b) Γ (c - b)} \int_{0}^{1} t^{b - 1} {(1 - t)}^{c - b - 1} {(1 - z t)}^{- a} d t

Let c be replaced by b + 1, so

{{}_{2}F}_{1} (a, b, b + 1; z) = b \int_{0}^{1} t^{b - 1} {(1 - z t)}^{- a} d t

Letting

t = s / c

, then

{{}_{2}F}_{1} (a, b, b + 1; z) = \frac{b}{c^{b}} \int_{0}^{c} s^{b - 1} {(1 - \frac{z s}{c})}^{- a} d s

The relationship between the operator I^β and the function ₂F₁ is derived as follows. The unit step function μ_c (t) is defined by

μ_{c} (t) = {\begin{array}{l} 1 & t \geq c \\ 0 & t < c \end{array}

If t < c, then

I^{β} (t^{α} μ_{c} (t)) = 0

. Assume that t ≥ c. Then we have

\begin{array}{l} I^{β} (t^{α} μ_{c} (t)) = \frac{1}{Γ (β)} \int_{0}^{t} \frac{s^{α} μ_{c} (s)}{{(t - s)}^{1 - β}} d s \\ = \frac{Γ (α + 1) t^{α + β}}{Γ (α + β + 1)} - \frac{t^{β - 1} c^{α + 1}}{Γ (β) (α + 1)} {{}_{2}F}_{1} (1 - β, α + 1, α + 2; \frac{c}{t}) \end{array}

(4)

3. The Boubaker wavelets

Here, we give the Boubaker wavelets as well as construct the fractional-order Boubaker wavelets. Boubaker wavelets b_nm(t) = b (k, n, m, t) are defined over [0, 1) by

b_{n, m} (t) = {\begin{array}{l} \frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} B_{m} (2^{k + 1} t - 4 n + 2) & \frac{n - 1}{2^{k - 1}} \leq t < \frac{n}{2^{k - 1}} \\ 0 & o t h e r w i s e \end{array}

(5)

with

w_{m} (t) = {\begin{array}{l} 4 & m = 0 \\ \sum_{r = 0}^{⌊ m / 2 ⌋} \frac{2^{2 m - 4 r + 2} α_{m, r}^{2}}{2 m - 4 r + 1} + 2 \sum_{k = 0}^{⌊ m / 2 ⌋ - 1} \sum_{s = k + 1}^{⌊ m / 2 ⌋} \frac{2^{2 m - 2 (k + s) + 2} α_{m, k} α_{m, s}}{2 m - 2 (k + s) + 1} & m \geq 1 \end{array}

(6)

where for

α_{i, r} = {(- 1)}^{r} (\begin{matrix} i - r \\ r \end{matrix}) \frac{i - 4 r}{i - r}

The definition of Boubaker polynomials over [ − 2, 2] is

B_{0} (t) = 1, B_{i} (t) = \sum_{r = 0}^{⌊ i / 2 ⌋} α_{i, r} t^{i - 2 r}, i \geq 1

where m = 0, 1, 2, …, M; k = 1, 2, 3, …; n = 1, 2, …, 2^k−1. Here M is the degree of the Boubaker polynomials and w_m is the normality coefficient.

3.1. Fractional-order Boubaker wavelets

To obtain the fractional-order Boubaker wavelets, t is replaced by t^α in Boubaker wavelets, where the real number α is an additional parameter. In particular, the formal definition of fractional-order Boubaker wavelets is as follow

ψ_{n, m}^{α} (t) = {\begin{matrix} \frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} B_{m} (2^{k + 1} t^{α} - 4 n + 2) & {(\frac{n - 1}{2^{k - 1}})}^{1 / α} \leq t < {(\frac{n}{2^{k - 1}})}^{1 / α} \\ 0 & otherwise \end{matrix}

3.2. Function approximation

A function f (t) with t ∈ [0, 1) can be expanded by

f (t) = \sum_{n = 1}^{\infty} \sum_{m = 0}^{\infty} c_{n, m} ψ_{n, m}^{α} (t)

This series can be truncated as

f (t) ≃ \sum_{n = 1}^{2^{k - 1}} \sum_{m = 0}^{M} c_{n, m} ψ_{n, m}^{α} (t) = C^{T} Ψ^{α} (t)

where

C^{T} = [c_{1,0}, \dots, c_{2^{k - 1}, 0}, c_{1,1}, \dots, c_{2^{k - 1}, 1}, \dots, c_{1, M}, \dots, c_{2^{k - 1}, M}]

and

Ψ^{α} {(t)}^{T} = [ψ_{1,0}^{α} (t), \dots, ψ_{2^{k - 1}, 0}^{α} (t), ψ_{1,1}^{α} (t), \dots, ψ_{2^{k - 1}, 1}^{α} (t), \dots, ψ_{1, M}^{α} (t), \dots, ψ_{2^{k - 1}, M}^{α} (t)]

(7)

4. The Riemann–Liouville fractional integral operator for fractional-order Boubaker wavelets

Here, a formula is provided for computing the RLFIO of the fractional-order Boubaker wavelets. We will establish I^βΨ^α (t) where Ψ^α (t) is given in equation (7) and we write

I^{β} Ψ^{α} (t) = {\bar{Ψ}}^{α} (t, β)

(8)

where

{\bar{Ψ}}^{α} (t, β) = {[I^{β} ψ_{1,0}^{α} (t), \dots, I^{β} ψ_{2^{k - 1}, 0}^{α} (t), \dots, I^{β} ψ_{1, M}^{α} (t), \dots, I^{β} ψ_{2^{k - 1}, M}^{α} (t)]}^{T}

To obtain the elements of vector

{\bar{Ψ}}^{α} (t, β)

, we have

\begin{array}{l} ψ_{n, m}^{α} (t) = \frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} μ_{{(\frac{n - 1}{2^{k - 1}})}^{1 / α}} (t) B_{m} (2^{k + 1} t^{α} - 4 n + 2) \\ - \frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} μ_{{(\frac{n}{2^{k - 1}})}^{1 / α}} (t) B_{m} (2^{k + 1} t^{α} - 4 n + 2) \end{array}

(9)

By integrating both sides of equation (9), we obtain

\begin{array}{l} I^{β} {ψ_{n, m}^{α} (t)} = I^{β} (\frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} μ_{{(\frac{n - 1}{2^{k - 1}})}^{1 / α}} (t) B_{m} (2^{k + 1} t^{α} - 4 n + 2)) \\ - I^{β} (\frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} μ_{{(\frac{n}{2^{k - 1}})}^{1 / α}} (t) B_{m} (2^{k + 1} t^{α} - 4 n + 2)) \end{array}

(10)

Now for m ≥ 1, by considering the following relations

B_{m} (2^{k + 1} t^{α} - 4 n + 2) = \sum_{r = 0}^{⌊ m / 2 ⌋} {(- 1)}^{r} (\begin{matrix} m - r \\ r \end{matrix}) \frac{m - 4 r}{m - r} {(2^{k + 1} t^{α} - 4 n + 2)}^{m - 2 r}

and

{(2^{k + 1} t^{α} - 4 n + 2)}^{m - 2 r} = \sum_{j = 0}^{m - 2 r} (\begin{matrix} m - 2 r \\ j \end{matrix}) {(2 - 4 n)}^{m - 2 r - j} 2^{(k + 1) j} t^{α j}

and using equation (10) we get

\begin{array}{l} F (t) = \frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} \sum_{r = 0}^{⌊ m / 2 ⌋} \sum_{j = 0}^{m - 2 r} {(- 1)}^{r} (\begin{matrix} m - r \\ r \end{matrix}) \frac{m - 4 r}{m - r} (\begin{matrix} m - 2 r \\ j \end{matrix}) {(2 - 4 n)}^{m - 2 r - j} 2^{(k + 1) j} \\ [t^{j α + β} \frac{Γ (j α + 1)}{Γ (j α + β + 1)} - \frac{t^{β - 1}}{Γ (β) (j α + 1)} {(\frac{n - 1}{2^{k - 1}})}^{(j α + 1) / α} {{}_{2}F}_{1} (1 - β, j α + 1, j α + 2; t {(\frac{n - 1}{2^{k - 1}})}^{1 / α})] \end{array}

Hence, by using equation (4), for m ≥ 1, we have

I^{β} ψ_{n, m}^{α} (t) = {\begin{array}{l} 0 & t \in [0, {(\frac{n - 1}{2^{k - 1}})}^{1 / α}) \\ F (t) & t \in [{(\frac{n - 1}{2^{k - 1}})}^{1 / α}, {(\frac{n}{2^{k - 1}})}^{1 / α}) \\ F (t) - G (t) & t \in [{(\frac{n}{2^{k - 1}})}^{1 / α}, 1) \end{array}

where

\begin{array}{l} G (t) = \frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} \sum_{r = 0}^{⌊ m / 2 ⌋} \sum_{j = 0}^{m - 2 r} {(- 1)}^{r} (\begin{matrix} m - r \\ r \end{matrix}) \frac{m - 4 r}{m - r} (\begin{matrix} m - 2 r \\ j \end{matrix}) {(2 - 4 n)}^{m - 2 r - j} 2^{(k + 1) j} \\ [t^{j α + β} \frac{Γ (j α + 1)}{Γ (j α + β + 1)} - \frac{t^{β - 1}}{Γ (β) (j α + 1)} {(\frac{n}{2^{k - 1}})}^{(j α + 1) / α} {{}_{2}F}_{1} (1 - β, j α + 1, j α + 2; t {(\frac{n}{2^{k - 1}})}^{1 / α})] \end{array}

and

\begin{array}{l} I^{β} {ψ_{n, m}^{α} (t)} = \frac{2^{(k + 1) / 2}}{\sqrt{w_{m}}} \sum_{r = 0}^{⌊ m / 2 ⌋} \sum_{j = 0}^{m - 2 r} {(- 1)}^{r} (\begin{matrix} m - r \\ r \end{matrix}) \frac{m - 4 r}{m - r} (\begin{matrix} m - 2 r \\ j \end{matrix}) {(2 - 4 n)}^{m - 2 r - j} 2^{(k + 1) j} \\ (I^{β} (μ_{{(\frac{n - 1}{2^{k - 1}})}^{1 / α}} (t) t^{α j}) - I^{β} (μ_{{(\frac{n}{2^{k - 1}})}^{1 / α}} (t) t^{α j})) \end{array}

Because the formula of Boubaker polynomials is valid for m ≥ 1, the same process for m = 0, will be considered separately. Since

ψ_{n, 0}^{α} (t) = 2^{(k + 1) / 2} μ_{{(\frac{n - 1}{2^{k - 1}})}^{1 / α}} (t) - 2^{(k + 1) / 2} μ_{{(\frac{n}{2^{k - 1}})}^{1 / α}} (t)

We have

I^{β} {ψ_{n, 0}^{α} (t)} = 2^{(k + 1) / 2} (I^{β} (μ_{{(\frac{n - 1}{2^{k - 1}})}^{1 / α}} (t)) - I^{β} (μ_{{(\frac{n}{2^{k - 1}})}^{1 / α}} (t)))

Therefore, for m = 0, we obtain

I^{β} ψ_{n, 0}^{α} (t) = {\begin{array}{l} 0 & t \in [0, {(\frac{n - 1}{2^{k - 1}})}^{1 / α}) \\ F_{0} (t) & t \in [{(\frac{n - 1}{2^{k - 1}})}^{1 / α}, {(\frac{n}{2^{k - 1}})}^{1 / α}) \\ F_{0} (t) - G_{0} (t) & t \in [{(\frac{n}{2^{k - 1}})}^{1 / α}, 1) \end{array}

where

F_{0} (t) = 2^{(k + 1) / 2} [t^{β} \frac{1}{Γ (β + 1)} - \frac{t^{β - 1}}{Γ (β)} {(\frac{n - 1}{2^{k - 1}})}^{1 / α} {{}_{2}F}_{1} (1 - β, 1,2; t {(\frac{n - 1}{2^{k - 1}})}^{1 / α})]

and

G_{0} (t) = 2^{(k + 1) / 2} [t^{β} \frac{1}{Γ (β + 1)} - \frac{t^{β - 1}}{Γ (β)} {(\frac{n}{2^{k - 1}})}^{1 / α} {{}_{2}F}_{1} (1 - β, 1,2; t {(\frac{n - 1}{2^{k - 1}})}^{1 / α})]

5. Numerical method

Consider the following FOCP Nemati et al. (2019):

Minimize (maximize)

J (t, x (t), u (t)) = \int_{0}^{1} ϕ (t, x (t), u (t)) d t

(11)

subjected to

D^{β_{0}} x (t) = g (t, x (t), D^{β_{1}} x (t), \dots, D^{β_{r}} x (t), u (t))

(12)

with the initial conditions

x^{(k)} (0) = d_{k}, k = 0, \dots, ⌈ β_{0} ⌉ - 1

(13)

where

x (t) = {[x_{1} (t), x_{2} (t), \dots, x_{l} (t)]}^{T}

is the vector of state functions, ϕ and g are given continuous functions, β₀ ≥ β₁ ≥…≥ β_r ≥ 0, while ⌈β₀⌉ indicates the ceiling function.

Without loss of generality, equations (12) and (13) can be written as

D^{β_{0}} x_{s} (t) = P_{s}^{1} (t) D^{β_{1}} x (t) + \dots + P_{s}^{r} (t) D^{β_{r}} x (t) + Q_{s} (t) u (t) + F_{s} (t), s = 1,2, \dots, l

(14)

and

x_{s}^{(k)} (0) = d_{s, k}, k = 0,1,2, \dots, ⌈ β_{0} ⌉ - 1

(15)

where

P_{s}^{i} (t), i = 1,2, \dots, r

, is a vector of continuous functions, Q_s (t), F_s (t) are continuous functions and d_sk, k = 0, …, ⌈β₀⌉ − 1, are constants.

In this paper, we want to find u (t) and the corresponding vector x (t) satisfying equations (14) and (15) while minimizing (maximizing) J in equation (11). For this purpose, assume that Q_j (t) ≠ 0, then from equation (14), u (t) can be written as

u (t) = \frac{1}{Q_{j} (t)} (D^{β_{0}} x_{j} (t) - P_{j}^{1} (t) D^{β_{1}} x (t) - \dots - P_{j}^{r} (t) D^{β_{r}} x (t) - F_{j} (t))

(16)

By substituting equation (16) to equations (11) and (14), the following optimization problem is obtained

m i n J (x, u) = \int_{0}^{1} ϕ (t, x, \frac{1}{Q_{j}} (D^{β_{0}} x_{j} - P_{j}^{1} D^{β_{1}} x - \dots - P_{j}^{r} D^{β_{r}} x - F_{j})) d t

subject to the following constraints for s = 1, …, l

\begin{array}{l} D^{β_{0}} x_{s} (t) = A_{s}^{1} (t) D^{β_{1}} x (t) + \dots + A_{s}^{r} (t) D^{β_{r}} x (t) + G_{s} (t) \\ x_{s}^{(k)} (0) = d_{s, k} \end{array}

(17)

where

A_{s}^{i} (t), i = 1,2, \dots, r

, and G_s (t) are the new but known vector and functions due to eliminating u (t) from the previous equations

To solve this optimization problem for s = 1, …, l, $D^{β_{0}} x_{s} (t)$ is expanded by the fractional-order Boubaker wavelets as

D^{β_{0}} x_{s} (t) ≃ \sum_{n = 1}^{2^{k - 1}} \sum_{m = 0}^{M} c_{n, m}^{s} ψ_{n, m}^{α} (t) = C_{s}^{T} Ψ^{α} (t)

(18)

Applying

I^{β_{0}}

to equation (18) results

x_{s} (t) = I^{β_{0}} D^{β_{0}} x_{s} (t) + \sum_{k = 0}^{⌈ β_{0} ⌉ - 1} \frac{t^{k}}{k!} x_{s}^{(k)} (0) ≃ C_{s}^{T} {\bar{Ψ}}^{α} (t, β_{0}) + \sum_{k = 0}^{⌈ β_{0} ⌉ - 1} \frac{t^{k}}{k!} d_{s, k}

(19)

where

d_{s, k} = x_{s}^{(k)} (0)

. From the above equation, we get

D^{β_{i}} x_{s} (t) ≃ C_{s}^{T} {\bar{Ψ}}^{α} (t, β_{0} - β_{i}) + \sum_{k = 0}^{⌈ β_{0} ⌉ - 1} \frac{D^{β_{i}} (t^{k})}{k!} d_{s, k}, for i = 1,2, \dots, r

(20)

By using equations (18)–(20) in equations (11) and (17), we obtain an optimization problem. By solving this problem, the values for C_s can be determined.

6. Convergence analysis

Theorem 1

Assume that D^(kα)f ∈ C(0, 1], k = 0, 1, …, M + 1. Then ∀t ∈ (0, 1]

f (t) = \sum_{k = 0}^{M} \frac{t^{k α}}{Γ (k α + 1)} D^{(k α)} f (0^{+}) + \frac{t^{M α + α}}{Γ (M α + α + 1)} D^{(M α + α)} f (ξ)

with 0 < ξ ≤ t. Moreover

| f (t) - \sum_{k = 0}^{M} \frac{t^{k α}}{Γ (k α + 1)} D^{(k α)} f (0^{+}) | \leq χ_{α} \frac{t^{M α + α}}{Γ (M α + α + 1)}

(21)

where χ_α ≥ sup_ξ∈(0,1]|D^(Mα+α)f(ξ)|.

Proof

See Odibat and Shawagfeh (2007).

Theorem 2

Suppose that f ∈ L²[0, 1), D^(kα)f ∈ C(0, 1], k = 0, 1, …, M + 1 and C^TΨ^α (t) is the best approximation of f out of $S = s p a n {ψ_{n, m}^{α} (t) | n = 1, \dots, 2^{k - 1}, m = 0, \dots, M}$ , then

{‖ f - C^{T} Ψ^{α} (t) ‖}_{L^{2} [0,1)} \leq \frac{1}{Γ (M α + α + 1) \sqrt{2 (M + 1) α + 1}} s u p_{ξ \in (0,1]} | D^{(M α + α)} f (ξ) |

Proof

The interval [0, 1) is devided to the 2^k−1 smaller intervals $I_{k, n} = [{(\frac{n - 1}{2^{k - 1}})}^{1 / α}, {(\frac{n}{2^{k - 1}})}^{1 / α}), n = 1,2, \dots, 2^{k - 1}$ . On each of the interval I_k,n, the wavelet Ψ^α(t) coincides with a polynomial in t^α of degree at most M. Using the generalized fractional-order Taylor’s expansion in equation (21), on the interval I_k,n for every n, we have

Hence

\begin{array}{l} {‖ f - C^{T} Ψ^{α} (t) ‖}_{L^{2} [0,1)}^{2} \leq \sum_{n = 1}^{2^{k - 1}} {‖ f - C^{T} Ψ^{α} (t) ‖}_{L^{2} I_{k, n}}^{2} \\ \leq {(s u p_{ξ \in (0,1]} | D^{(M α + α)} f (ξ) | \frac{1}{Γ (M α + α + 1)})}^{2} . \sum_{n = 1}^{2^{k - 1}} {‖ t^{M α + α} ‖}_{L^{2} I_{k, n}}^{2} \\ = {(s u p_{ξ \in (0,1]} | D^{(M α + α)} f (ξ) | \frac{1}{Γ (M α + α + 1)})}^{2} {‖ t^{M α + α} ‖}_{L^{2} [0,1)}^{2} \\ = {(\frac{1}{Γ (M α + α + 1) \sqrt{2 (M + 1) α + 1}} s u p_{ξ \in (0,1]} | D^{(M α + α)} f (ξ) |)}^{2} \end{array}

Taking the square root of the above inequality completed the proof.

Now regarding equation (18), we know that c_n,m are chosen in such a way that C^TΨ^α (t) be the best approximation of $D^{β_{0}} x (t)$ , so by considering theorem 1, it can be concluded that by increasing M the $D^{β_{0}} x (t) - C^{T} Ψ^{α} (t)$ , tends to zero.

Regarding the equation (19), since the RLFIO is achieved using the new idea of beta function and without any approximation, the error occurs while I^β is taken on the best approximation of function not on the function itself. In other words, in the second term of following relation

I^{β_{0}} D^{β_{0}} x (t) + x (0) ≃ I^{β_{0}} C^{T} Ψ^{α} (k) + x (0) = C^{T} {\bar{Ψ}}^{α} (t, β_{0}) + x (0)

The following theorem shows that this error tends to zero as well, if the number of basis functions becomes greater.

Theorem 3

Suppose that f ∈ L² [0, 1) and C^TΨ^α (t) is the best approximation of f. Then, we have

\lim_{M \to \infty} {‖ I^{β} f - I^{β} C^{T} Ψ^{α} ‖}_{L^{2} (0,1)} = 0

Proof

In Theorem 2, it is shown that C^TΨ^α(t) converges to f uniformly because ‖f − C^TΨ^α‖₂ = 0 as M tends to infinity. By using the fact that limit and integration are commutative for uniformly convergence function over finite intervals Rudin (1994), we can write

\begin{array}{l} \lim_{M \to \infty} | I^{β} f - I^{β} C^{T} Ψ^{α} | = \frac{1}{Γ (β)} \lim_{M \to \infty} \int_{0}^{t} | {(t - s)}^{β - 1} f (s) - C^{T} Ψ^{α} (s) | d s \\ = \frac{1}{Γ (β)} \int_{0}^{t} \lim_{M \to \infty} | {(t - s)}^{β - 1} (f (s) - C^{T} Ψ^{α} (s)) | d s \\ = 0 \end{array}

\begin{array}{l} \int_{0}^{1} \lim_{M \to \infty} {| I^{β} f - I^{β} C^{T} Ψ^{α} |}^{2} = \lim_{M \to \infty} \int_{0}^{1} {| I^{β} f - I^{β} C^{T} Ψ^{α} |}^{2} \\ = \lim_{M \to \infty} {‖ I^{β} f - I^{β} C^{T} Ψ^{α} ‖}_{2}^{2} \\ = 0 \end{array}

The convergence of equation (20) is similar to equation (19) and so we skip it.

6.1. Illustrative examples

We give six examples and compare the numerical solutions obtained by using our numerical methods with those from several other existing methods. The computation of our numerical solutions is performed by using Mathematica 12.1.1.

Example 1

Minimize Barikbin and Keshavarz (2020)

J (x, u) = \int_{0}^{1} (0.625 x^{2} (t) + 0.5 u^{2} (t) + 0.5 x (t) u (t)) d t

subject to

D^{β} x (t) = 0.5 x (t) + u (t), t \in [0,1], 0 < β \leq 1

and

x (0) = 1

For β = 1, the exact optimal control is given by El-Kady (2003)

u (t) = \frac{- (\tanh (1 - t) + 0.5) \cosh (1 - t)}{\cosh (1)}

For this u (t), the exact cost function J is given by J = 0.380 797 077.

This example was solved for β = 1 in El-Kady (2003) by using Chebyshev finite difference technique with M₁ = 7, in Rabiei et al. (2017) by applying Boubaker polynomials with M₂ = 5, and in Barikbin and Keshavarz (2020) by using Bernoulli wavelets with k = 1 and M₃ = 4. These results are given in Table 1 together with the value of J obtained by using our method with β = k = 1 and M = 3. In this table, M₁, M₂, and M₃ are the degrees of the Chebyshev, Boubaker, and Bernoulli polynomials.

In Figure 1, the numerical and the exact solutions of u (t) and x (t) are given for β = 0.5, 0.8, 0.9, 1 with k = 2, M = 3. The figure shows that when β tends to 1, the numerical solutions approach the exact solutions.

Table 1.

Comparison of J for different methods for Example 1.

Method	J
Chebyshev finite difference
M₁ = 7, β = 1	0.380 797
Boubaker polynomials
M₂ = 5, β = 1	0.380 797
Bernoulli wavelet
k = 1, M₃ = 4, β = 1	0.380 797
Present method
k = 1, M = 3, β = 1	0.380 797 07

Figure 1.

Graphs of the numerical and exact values for x (t) and u (t) with β = 0.5, 0.8, 0.9 and 1 in Example 1.

Example 2

Minimize Behroozifar and Habibi (2017)

J = \int_{0}^{1} (x^{2} (t) + u^{2} (t) + e^{2 t} + \frac{1}{4} e^{- 2 t + 2 t^{3 / 2}} - \frac{3}{8} e^{- 2 t + t^{3 / 2}} \sqrt{π} + \frac{9}{64} e^{- 2 t} π + t^{3} + e^{- t + t^{3 / 2}} u (t) - \frac{3}{4} e^{- t} \sqrt{π} u (t) - 2 t^{3 / 2} x (t)) d t

(22)

with

D^{1.5} x (t) = 2 e^{t} u (t) + e^{x (t)}

and

x (0) = x^{'} (0) = 0

The exact state, control, and cost functions are

x (t) = t^{3 / 2}

u (t) = \frac{1}{2} e^{- t} (- e^{t^{3 / 2}} + \frac{3 \sqrt{π}}{4})

, and J = 3.194 528 049 46. In Table 2, the numerical solutions obtained by using the present method by selecting k = 1, M = 3, α = 0.5 are compared with those obtained by using the Bernoulli polynomial method Behroozifar and Habibi (2017) by selecting M₃ = 6 and 8. This problem is also solved for k = 1, M = 3, α = 0.5, and we report the absolute errors of the obtained x (t) and u (t) in Table 3.

Table 2.

Comparison of J for Example 2.

Method	Error x (t)	Error u (t)	J
Bernoulli polynomials
M₃ = 6	Multiples of 10⁻⁴	Multiples of 10⁻⁴	−
M₃ = 8	Multiples of 10⁻⁴	Multiples of 10⁻⁵	3.194 53
Present method
k = 1, M = 3, α = 0.5	Multiples of 10⁻¹³	Multiples of 10⁻¹²	3.194 528 049 5

Table 3.

AEs for k = 1, M = 3 and α = 0.5, for Example 2.

t	Errors of x(t)	Errors of u(t)
0.1	1.48 × 10⁻¹³	2.01 × 10⁻¹²
0.2	2.49 × 10⁻¹³	6.84 × 10⁻¹³
0.3	2.59 × 10⁻¹³	5.84 × 10⁻¹³
0.4	1.76 × 10⁻¹³	1.21 × 10⁻¹²
0.5	3.75 × 10⁻¹⁴	1.21 × 10⁻¹²
0.6	9.15 × 10⁻¹⁴	7.06 × 10⁻¹³
0.7	1.22 × 10⁻¹³	1.54 × 10⁻¹³
0.8	5.07 × 10⁻¹⁴	1.21 × 10⁻¹²
0.9	1.53 × 10⁻¹⁵	2.32 × 10⁻¹²

Example 3

Minimize the following nonlinear problem Dehestani et al. (2020)

J = \int_{0}^{1} ({(x (t) - t^{2})}^{2} + {(u (t) - t e^{- t} + \frac{1}{2} e^{t^{2} - t})}^{2}) d t

with

D^{β} x (t) = e^{x (t)} + 2 e^{t} u (t), 0 < β \leq 1

and

x (0) = 0

For β = 1, the exact x (t) and u (t) are given by x (t) = t² and

u (t) = t e^{- t} - \frac{1}{2} e^{t^{2} - t}

. For these values, we get the cost function as J = 0.

In Tables 4 and 5, the AEs of the control and state functions are given with those obtained by Bessel wavelets method reported in Dehestani et al. (2020) for k = 1, M₄ = 5 with α = 1, and α = 0.5, respectively. In both tables, M₄ gives the degree of the Bessel polynomials.

In Figure 2, the graphs of the exact and calculated state and control are plotted for k = 2 and M = 3 and several values of β. Similarly to Example 1, the numerical solutions tend to the exact one as β approaches 1.

Table 4.

AEs of x (t) and u (t) for α = 1, k = 1 and M = M₄ = 5 in Example 3.

t	Present method		Bessel wavelet
	$\bar{x (t)}$	u (t)	$\bar{x (t)}$	u (t)
0	0	1.640 × 10⁻¹²	5.145 × 10⁻²⁷	2.754 × 10⁻⁴
0.1	6.848 × 10⁻¹⁵	6.601 × 10⁻¹³	5.935 × 10⁻⁵	3.336 × 10⁻⁵
0.2	9.889 × 10⁻¹⁴	7.022 × 10⁻¹⁴	4.962 × 10⁻⁵	4.223 × 10⁻⁵
0.3	4.178 × 10⁻¹⁴	4.661 × 10⁻¹³	2.364 × 10⁻⁵	9.709 × 10⁻⁵
0.4	8.604 × 10⁻¹⁴	3.179 × 10⁻¹³	7.165 × 10⁻⁶	5.704 × 10⁻⁵
0.5	1.256 × 10⁻¹³	1.557 × 10⁻¹³	6.183 × 10⁻⁶	1.375 × 10⁻⁵
0.6	2.503 × 10⁻¹⁴	4.154 × 10⁻¹³	1.442 × 10⁻⁵	3.822 × 10⁻⁵
0.7	1.144 × 10⁻¹³	1.975 × 10⁻¹³	2.072 × 10⁻⁵	3.723 × 10⁻⁶
0.8	1.175 × 10⁻¹³	2.845 × 10⁻¹³	1.648 × 10⁻⁵	5.719 × 10⁻⁵
0.9	5.628 × 10⁻¹⁴	3.401 × 10⁻¹³	3.060 × 10⁻⁶	2.902 × 10⁻⁵
1	1.987 × 10⁻¹⁴	9.996 × 10⁻¹³	8.043 × 10⁻⁷	1.098 × 10⁻⁴

Table 5.

AEs of x (t) and u (t) for α = 0.5, k = 1 and M = M₄ = 5 in Example 3.

t	Present method		Bessel wavelet
	$\bar{x (t)}$	u (t)	$\bar{x (t)}$	u (t)
0	0	1.683 × 10⁻¹³	1.769 × 10⁻²⁵	9.352 × 10⁻²
0.1	1.580 × 10⁻¹⁵	2.303 × 10⁻¹⁴	3.875 × 10⁻²	4.421 × 10⁻²
0.2	3.205 × 10⁻¹⁵	1.415 × 10⁻¹⁴	4.401 × 10⁻²	9.830 × 10⁻³
0.3	5.329 × 10⁻¹⁵	5.828 × 10⁻¹⁶	3.914 × 10⁻²	1.192 × 10⁻²
0.4	4.274 × 10⁻¹⁵	9.547 × 10⁻¹⁵	2.813 × 10⁻²	1.080 × 10⁻²
0.5	1.110 × 10⁻¹⁵	1.178 × 10⁻¹⁴	1.463 × 10⁻²	1.065 × 10⁻²
0.6	2.553 × 10⁻¹⁵	8.160 × 10⁻¹⁵	8.052 × 10⁻⁴	1.246 × 10⁻²
0.7	4.996 × 10⁻¹⁵	6.175 × 10⁻¹⁶	1.070 × 10⁻²	1.547 × 10⁻²
0.8	3.996 × 10⁻¹⁵	8.909 × 10⁻¹⁵	1.759 × 10⁻²	1.806 × 10⁻²
0.9	1.998 × 10⁻¹⁵	1.801 × 10⁻¹⁴	1.769 × 10⁻²	1.834 × 10⁻²
1	1.509 × 10⁻¹⁴	2.453 × 10⁻¹⁴	8.939 × 10⁻²	1.457 × 10⁻²

Figure 2.

Graphs of the numerical and exact values for x (t) and u (t) with β = 0.5, 0.8, 0.9, 1 for Example 3.

Example 4

Consider the following time varying problem Sahu and Saha Ray (2018)

m i n J (x, u) = \frac{1}{2} \int_{0}^{1} (x^{2} (t) + u^{2} (t)) d t

subject to

D^{β} x (t) = t x (t) + u (t), x (0) = 1

In this example, we compare our method with several wavelets that were compared to each other in Sahu and Saha Ray (2018). These wavelets are Legendre, Chebyshev, Laguerre and CAS wavelets. In that paper, these wavelets were compared by using k = 1 and M₅ = 4 for Legendre, M₁ = 4 for Chebyshev, M₇ = 4 for Laguerre and M₆ = 4 for CAS wavelets. These wavelets are compared with our method with k = 1, 2 and different M in Table 6. In this table, M₇ is the degree of the Laguerre polynomials.

Table 6.

Comparison of J between our method and the methods in Sahu and Saha Ray (2018) for Example 4.

Method	J
Legendre wavelets
k = 1, M₅ = 4	0.484 268
Chebyshev wavelets
k = 1, M₁ = 4	0.484 232
Laguerre wavelets
k = 1, M₇ = 4	0.483 476
CAS wavelets
k = 1, M₆ = 4	0.500 511
Present method
k = 1, M = 2	0.484 422 222 671 866 8
k = 1, M = 4	0.484 267 697 101 286 7
k = 1, M = 7	0.484 267 696 228 727 2
k = 1, M = 8	0.484 267 696 228 727 1
k = 1, M = 9	0.484 267 696 228 727 1
k = 2, M = 4	0.484 267 696 228 727 4
k = 2, M = 5	0.484 267 696 228 727 2

Example 5

Consider minimizing the following nonlinear problem Nemati et al. (2019)

J = \int_{0}^{1} (e^{t} {(x (t) - t^{4} + t - 1)}^{2} + (1 + t^{2}) {(u (t) + t^{4} - t + 1 - \frac{8000 t^{21 / 10}}{77 Γ (1 / 10)})}^{2}) d t

(23)

with

D^{1.9} x (t) = x (t) + u (t)

and

x (0) = 1, x^{'} (0) = - 1

The exact state and control variables are given by x(t) = t⁴ − t + 1 and

u (t) = - t^{4} + t - 1 + (8000 t^{21 / 10} / 77 Γ (0.1))

. With these values, the cost function is J = 0. This example was solved in Lotfi et al. (2013) by using the Legendre orthonormal functions with M₅ = 4 and M₅ = 8, in Keshavarz et al. (2016) via Bernoulli polynomials with M₂ = 4 and M₂ = 8, and in Nemati et al. (2019) by modified hat functions with M₈ = 4 and M₈ = 8. These methods are compared to our method with k = 1, 2 and different M in Table 7. In this table, M₈ is the number of subintervals in the hat function method.

Table 7.

Comparison of J for Example 5.

Method	J
Legendre basis
M₅ = 4	5.42 × 10⁻⁷
M₅ = 8	8.22 × 10⁻¹⁰
Bernoulli basis
M₂ = 4	5.16 × 10⁻⁷
M₂ = 8	4.23 × 10⁻¹¹
modified hat basis
M₈ = 4	9.64 × 10⁻⁷
M₈ = 8	1.00 × 10⁻⁸
Present method
k = 1, M = 4	6.05 × 10⁻⁷
k = 1, M = 6	1.64 × 10⁻¹¹
k = 1, M = 8	3.71 × 10⁻¹⁴
k = 2, M = 4	2.08 × 10⁻⁸
k = 2, M = 6	8.12 × 10⁻¹³
k = 2, M = 8	1.05 × 10⁻¹⁵

Example 6

We consider the two dimensional FOCP given as Lotfi and Yousefi (2014)

\min J = \int_{0}^{1} ({(x_{1}^{2} (t) - 1 - t^{3 / 2})}^{2} + {(x_{2} (t) - t^{5 / 2})}^{2} + {(u (t) - \frac{3 \sqrt{π}}{4} t + t^{5 / 2})}^{2}) d t

(24)

with

D^{0.5} x_{1} (t) = x_{2} (t) + u (t)

(25)

D^{0.5} x_{2} (t) = x_{1} (t) + \frac{15 \sqrt{π}}{16} t^{2} - t^{3 / 2} - 1

(26)

and

x_{1} (0) = 1 x_{2} (0) = 0

The exact state and control functions are

(x_{1} (t), x_{2} (t)) = (1 + t^{3 / 2}, t^{5 / 2})

and

u (t) = ((3 \sqrt{π}) / 4) t - t^{5 / 2}

. For these values, the cost function is J = 0.

This problem was solved by Epsilon-Ritz method in Lotfi and Yousefi (2014), and the best reported solution is 8.002 7 × 10⁻⁶ by using 10 basis functions. Here by choosing M = 3, k = 1 and $α = 1 / 2$ , we consider

D^{1} x_{2} (t) = C^{T} Ψ^{α} (t)

Then by using equation (8), we have

x_{2} (t) = C^{T} {\bar{Ψ}}^{α} (t, 1)

(27)

By substituting equation (27) into equation (26), we get

x_{1} (t) = C^{T} {\bar{Ψ}}^{α} (t, \frac{1}{2}) - \frac{15 \sqrt{π}}{16} t^{2} + t^{3 / 2} + 1

(28)

Then by substituting the above equation into equation (25), u (t) is obtained as

u (t) = C^{T} {\bar{Ψ}}^{α} (t, 0) - C^{T} {\bar{Ψ}}^{α} (t, 1) + \frac{3 \sqrt{π}}{4} t - \frac{5}{2} t^{3 / 2}

(29)

By substituting the equations (27)–(29) into equation (24) and solving the obtained optimization problem, C^T is calculated as

C^{T} = {- \frac{\sqrt{\frac{313805}{273}}}{16}, \frac{11 \sqrt{\frac{313805}{273}}}{64}, \frac{\sqrt{\frac{941415}{91}}}{32}, \frac{\sqrt{\frac{313805}{273}}}{64}}

(30)

By using equation (30), we get

(x_{1} (t), x_{2} (t)) = (1 + t^{3 / 2}, t^{5 / 2})

and

u (t) = ((3 \sqrt{π}) / 4) t - t^{5 / 2}

which are the exact state and control functions. These values give the exact cost function J = 0. These exact values were not obtained previously in the literature.

Example 7

We consider the following cancer model Hassani et al. (2021)

\min J (u) = \int_{0}^{1} (20 T (t) - N (t) + u^{2} (t)) d t

(31)

subject to the nonlinear fractional differential system

D^{β} T (t) = 1.5 T (t) (1 - T (t)) - 0.5 T (t) I (t) - T (t) N (t) + 1.5 T (t) F (t) - 0.08 D_{1} (t) T (t)

(32)

D^{β} I (t) = 0.33 + 0.01 \frac{T (t) I (t)}{0.3 + T (t)} - 0.5 T (t) I (t) - 0.2 I (t) - 2 \times 10^{- 11} D_{1} (t) I (t)

(33)

D^{β} N (t) = N (t) (1 - N (t)) - T (t) N (t) - 0.008 D_{1} (t) N (t)

(34)

D^{β} F (t) = 0.75 F (t) (1 - 1.5 F (t)) - 0.1 T (t) F (t) - 0.008 D_{1} (t) F (t)

(35)

D^{β} D_{1} (t) = u (t) - 0.1 D_{1} (t)

(36)

where T (t) stands for tumor cells, I (t) immune cells, N (t) normal cells, and F (t) fat cells around tumor cells at time t. D₁ (t) is the chemotherapeutic drug, and u (t) is the dose of this drug that is injected (for more details about the model, see Hassani et al. (2021)). The initial conditions are

T (0) = 1, I (0) = 0.001, N (0) = 4, F (0) = 0.25, D_{1} (0) = 0.5

(37)

This problem is solved for β = 0.95 that is used in Hassani et al. (2021), α = 1, k = 1, and M = 4. First, D^βT(t), D^βI(t), D^βN(t), D^βF(t), and D^βD₁(t) are expanded by fractional-order Boubaker wavelets as follows

\begin{array}{l} D^{β} T (t) = C_{1}^{T} Ψ^{1} (t) \\ D^{β} I (t) = C_{2}^{T} Ψ^{1} (t) \\ D^{β} N (t) = C_{3}^{T} Ψ^{1} (t) \\ D^{β} F (t) = C_{4}^{T} Ψ^{1} (t) \\ D^{β} D_{1} (t) = C_{5}^{T} Ψ^{1} (t) \end{array}

(38)

where C_i for i = 1, …, 5 are unknown coefficient vectors and Ψ¹ is introduced in equation (7). Next, by applying the operator I^β in equation (38), we get

T (t) = C_{1}^{T} {\bar{Ψ}}^{1} (t, β) + 1

(39)

I (t) = C_{2}^{T} {\bar{Ψ}}^{1} (t, β) + 0.001

(40)

N (t) = C_{3}^{T} {\bar{Ψ}}^{1} (t, β) + 4

(41)

F (t) = C_{4}^{T} {\bar{Ψ}}^{1} (t, β) + 0.25

(42)

D_{1} (t) = C_{5}^{T} {\bar{Ψ}}^{1} (t, β) + 0.5

(43)

Then, by substituting D^βD₁(t) from equation (43) in equation (36), u (t) is achieved as

u (t) = C_{5}^{T} Ψ^{1} (t) + 0.1 C_{5}^{T} {\bar{Ψ}}^{1} (t, β) + 0.5

(44)

Finally, after substituting equations (38)–(43) in the performance index J in equation (31), and differential equations in equations (32)–(36), the problem is reduced to a parameter optimization problem. The approximate functions obtained for the states and control variables are shown in Figure 3. From this figure, it is shown that over time, the number of tumor cells (T (t)) decreases, the number of normal cells starts to grow, and the immune and fat cells population increase, while the drug concentration and its dose decline after destroying most of the tumor cells.

Figure 3.

Numerical values of T (t), I (t), N (t), F (t), D₁ (t), and u (t) obtained by the present method for k = 1 and M = 4 for Example 7.

7. Conclusion

This article is devoted to an effective scheme for solving FOCPs by using fractional-order Boubaker wavelets. The convergence analysis of the method is provided. Several numerical examples were given to indicate the validity of the method.

Some advantages of the proposed method are as follows:

1. We could obtain the exact value of the RLFIO for fractional-order Boubaker wavelets.

2. The proposed method could obtain the exact solutions of some problems whose exact solutions are fractional-order power terms. These exact solutions were not obtained previously in the literature.

3. Our numerical method gives more accurate solutions than those shown in the literature.

Footnotes

Acknowledgments

The authors wish to express their sincere thanks to the anonymous referees and the editor for valuable suggestions that improved the final version of the manuscript.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Mohsen Razzaghi

References

Agrawal

(2008) Fractional optimal control of a distributed system using eigenfunctions. ASME J Comput Nonlinear Dyn. 3(2): 021204.

Andrews

Askey

Roy

(1999) Special Functions. Cambridge University Press, Cambridge.

Baleanu

Diethelm

Scalas

, et al. (2016) Fractional Calculus: Models and Numerical Methods. Singapore: World Scientific.

Barikbin

Keshavarz

(2020) Solving fractional optimal control problems by new Bernoulli wavelets operational matrices. Optim Control Appl Meth 41: 1188–1210.

Behroozifar

Habibi

(2017) A numerical approach for solving a class of fractional optimal control problems via operational matrix Bernoulli polynomials. J Vib Control 24 (12): 2494–2511

Bhrawy

Alhamed

Baleanu

(2014) New spectral techniques for systems of fractional differential equations using fractional-order generalized Laguerre orthogonal functions.Fract Calc Appl Anal 17: 1138–1157.

Boubaker

(2007) On modified Boubaker polynomials: some differential and analytical properties of the new polynomial issued from an attempt for solving bi-varied heat equation.Trends in Applied Sciences Research 2: 540–544.

Boubaker

(2011) Boubaker polynomials expansion scheme (BPES) solution to Boltzmann diffusion equation in the case of strongly anisotropic neutral particles forward backward scattering. Annals of Nuclear Energy 38: 1715–1717.

Cattani

Guariglia

Wang

, et al. (2017) On the critical strip of the Riemann zeta fractional derivative. Fundam Inf 151(1–4): 459–472.

10.

Dehestani

Ordokhani

Razzaghi

(2020) Fractional-order Bessel wavelet functions for solving variable order fractional optimal control problems with estimation error. International Journal of Systems Science 51(6): 1032–1052.

11.

El-Kady

(2003) A Chebyshev finite difference method for solving a class of optimal control problems. Intern J Computer Math 7: 883–89.

12.

Hassani

Machado

Mehrabi

(2021) An optimization technique for solving a class of nonlinear fractional optimal control problems: application in cancer treatment. Mathematical Modelling 93: 868–884.

13.

Heydari

Hooshmandasl

Maalek Ghaini

, et al. (2016) Wavelets method for solving fractional optimal control problems. Applied Mathematics and Computation 286: 139–154.

14.

Kazem

Abbasbandy

Kumar

(2013) Fractional-order Legendre functions for solving fractional-order differential equations. Appl Math Model 37(7): 5498–5510.

15.

Keshavarz

Ordokhani

Razzaghi

(2016) A numerical solution for fractional optimal control problems via Bernoulli polynomials. J Vib Control. 22(18): 3889–3903.

16.

Kumar

(2010) An analytical solution to applied mathematics-related Love’s equation using the Boubaker polynomials expansion scheme. Journal of the Franklin Institute. 347(9): 1755–1761.

17.

Zhao

(2010) Haar wavelet operational matrix of fractional order integration and its applications in solving the fractional order differential equations. Appl Math Comput 216: 2276–2285.

18.

Lotfi

Yousefi

Dehghan

(2013) Numerical solution of a class of fractional optimal control problems via the Legendre orthonormal basis combined with the operational matrix and the Gauss quadrature rule. J Comput Appl Math 250: 143–160.

19.

Lotfi

Yousefi

(2014) Epsilon − Ritz method for solving a class of fractional constrained optimization problems. J Optim Theory Appl 163: 884–899.

20.

Mashayekhi

Razzaghi

(2016) Numerical solution of the fractional BagleyTorvik equation by using hybrid functions approximation. Mathematical Methods in the Applied Sciences 39(3): 353–365.

21.

Mashayekhi

Razzaghi

(2018) An approximate method for solving fractional optimal control problems by hybrid functions. J Vib Control. 24(9): 1621–1631

22.

Mohammadi

Moradi

Baleanu

, et al. (2018) A hybrid functions numerical scheme for fractional optimal control problems: application to nonanalytic dynamic systems, Journal of Vibration and Control 24(21): 5030–5043.

23.

Nemati

Lima

Torres

DFM

(2019) A numerical approach for solving fractional optimal control problems using modified hat functions. Commun Nonlinear Sci Numer Simulat 78: 104849.

24.

Odibat

Shawagfeh

(2007) Generalized Taylors formula. Appl Math Comput 186: 286–293.

25.

Podlubny

(1999) Fractional Diferential Equations. New York: Academic Press.

26.

Rabiei

Ordokhani

Babolian

(2017) The Boubaker polynomials and their application to solve fractional optimal control problems. Nonlinear Dyn 88(2): 1013–1026.

27.

Rabiei

Ordokhani

(2019) Solving fractional pantograph delay differential equations via fractional-order Boubaker polynomials. Engineering with Computers 35: 1431–1441.

28.

Rabiei

Ordokhani

(2020) A new operational matrix based on Boubaker wavelet for solving optimal control problems of arbitrary order. Transactions of the Institute of Measurement and Control 42(10): 1858–1870.

29.

Rafiei

Kafash

Karbassi

(2018) A new approach based on using Chebyshev wavelets for solving various optimal control problems. Computational and Applied Mathematics 37(1): 144–157.

30.

Rahimkhani

Ordokhani

Babolian

(2016) Fractional-order Bernoulli wavelets and their applications. Appl Math Model 40: 8087–8107.

31.

Razzaghi

Elnagar

(1994) Linear quadratic optimal control problems via shifted Legendre state parametrization. Int J Syst Sci 25: 393–9.

32.

Rudin

(1994) Principles of Mathematical Analysis. 3rd edition. McGraw-Hill, New York. 151. Theorem 7.16.

33.

Saeedi

Moghadam

Mollahasani

, et al. (2011) A CAS wavelet method for solving nonlinear Fredholm integro-differential equations of fractional order. Communications in Nonlinear Science and Numerical Simulation 16(3): 1154–1163.

34.

Sahu

Saha Ray

(2018) Comparison on wavelets techniques for solving fractional optimal control problems. Journal of Vibration and Control 24(6): 1185–1201.

35.

Salati

Shamsi

Torres

DFM

(2019) Direct transcription methods based on fractional integral approximation formulas for solving nonlinear fractional optimal control problems. Communications in Nonlinear Science and Numerical Simulation 67: 334–350.

36.

Toan

Razzaghi

(2019) Taylor wavelet method for fractional delay differential equations. Eng Comput 37(1): 231–240.

37.

Tricaud

Chen

(2010) An approximate method for numerically solving fractional order optimal control problems of general form. Comput. Math Appl 59(5): 1644–1655.

38.

Vichitkunakorn

Razzaghi

(2020) A numerical method for fractional pantograph differential equations based on Taylor wavelets. T I Meas Control 42: 1334–1344.