EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 4, Article Number 7063
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Numerical Solution of Two-Dimensional Fractional
Optimal Control Problems Using Fractional

Vieta-Fibonacci Wavelets
G. M. Bahaa1, A. H. Qamlo2,∗

1 Department of Mathematics and Computer Science, Faculty of Science, Beni-Suef Uni-
versity, Beni-Suef 26511, Egypt
2 Department of Mathematics, Faculty of Science, Umm Al-Qura University, Makkah 21955,
Saudi Arabia

Abstract. This paper presents a novel numerical framework for solving two-dimensional frac-
tional optimal control problems (FOCPs) governed by Caputo fractional partial differential equa-
tions. The proposed approach is based on fractional Vieta–Fibonacci wavelets (FVFWs), a re-
cently developed wavelet family that combines the recursive structure of Fibonacci polynomials
with fractional-order operators. We first construct operational matrices for fractional derivatives
in the FVFW basis and employ them to transform the governing FOCP into a system of sparse
algebraic equations. The state, adjoint, and control functions are approximated simultaneously in
a unified wavelet space, enabling efficient reconstruction of the optimality system. The resulting
nonlinear algebraic system is solved using Newton’s method, which guarantees quadratic conver-
gence under standard assumptions. A benchmark two-dimensional fractional control problem is
presented to validate the proposed scheme. Numerical results demonstrate that FVFWs achieve
high accuracy with relatively few basis functions, outperforming conventional wavelet and spectral
approaches in terms of computational efficiency. The method provides a general framework that
can be extended to nonlinear fractional PDEs, time-dependent fractional dynamics, and problems
with control constraints.
2020 Mathematics Subject Classifications: 49M37, 65M70, 26A33, 93C20
Key Words and Phrases: Fractional optimal control, Caputo derivative, Vieta–Fibonacci
wavelets, operational matrix, numerical methods, Newton’s method

1. Introduction

Fractional calculus has attracted significant attention in recent decades as a powerful
mathematical framework for describing systems with memory, hereditary properties, and
anomalous diffusion. Unlike classical integer-order models, fractional differential equations

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.7063

Email addresses: Bahaa_gm@yahoo.com (G. M. Bahaa),
Ahqamlo@uqu.edu.sa (A. H. Qamlo)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 2 of 29

(FDEs) provide more realistic descriptions of physical and engineering processes such as
viscoelasticity, diffusion in complex media, and control of distributed parameter systems
[1–5].

The incorporation of fractional dynamics into optimal control theory has led to the
development of fractional optimal control problems (FOCPs). Several authors have studied
existence, uniqueness, and optimality conditions in this setting. For instance, Agrawal
and Baleanu [6] introduced a Hamiltonian framework for fractional optimal control, while
Almeida and Torres [7] derived necessary and sufficient conditions for problems involving
Caputo derivatives. Li and Zhou [8] established the existence of optimal solutions, and
subsequent works by Mahmudov [9] and Sakthivel et al. [10] extended controllability and
impulsive control results to semilinear systems. More recently, Bahaa and collaborators
[11–19] investigated FOCPs with variable order, time delay, weak Caputo derivatives, and
control constraints, thereby highlighting the growing breadth of applications.

Alongside theoretical progress, numerical techniques for FOCPs have evolved rapidly.
Lotfi et al. [20, 21] developed spectral methods using orthogonal polynomials, while
Bhrawy et al. [22] proposed efficient multi-dimensional schemes with quadratic perfor-
mance indices. More recent contributions include Heydari et al. [23], who employed Bern-
stein polynomials, and Dehestani et al. [24], who studied robust optimization approaches
for variable-order problems. These works highlight the demand for computationally ef-
ficient and accurate methods to approximate fractional operators, which are inherently
nonlocal and computationally expensive.

Wavelet-based methods have proven especially effective in this context due to their
localization and multiresolution properties [25]. Classical constructions such as Legen-
dre and Chebyshev wavelets have been employed for FOCPs, but recent attention has
shifted toward Fibonacci-based wavelets. El-Hawary and El-Khazendar [26] first intro-
duced Fibonacci wavelets for dynamic systems. Building on this, Agarwal et al. [27] de-
veloped Vieta–Fibonacci operational matrices for variable-order integro-differential equa-
tions, while Azin et al. [28, 29] applied Vieta–Fibonacci wavelets to fractional panto-
graph and delay differential equations. These advances motivated the introduction of
fractional Vieta–Fibonacci wavelets (FVFWs), which offer excellent approximation capa-
bilities with fewer expansion terms compared to classical wavelets. Very recently, Hoseini
et al. [30] demonstrated the effectiveness of FVFWs in solving two-dimensional FOCPs,
which strongly supports the approach developed in this work.

The present paper contributes to this active area of research by developing a numerical
framework for two-dimensional FOCPs based on FVFWs. In contrast to earlier works that
focused either on one-dimensional problems or on spectral techniques without fractional
wavelets, our method directly leverages the recursive structure of Vieta–Fibonacci poly-
nomials to construct fractional operational matrices for Caputo derivatives. This allows
the original FOCP to be reduced to a sparse system of algebraic equations, which is then
solved using Newton’s method. Our work is also closely related to recent studies on bang-
bang and time-optimal control for fractional systems [31], symmetric distributed-order
systems [32], and numerical solutions on complex domains [33], as it provides a unifying
computational approach applicable to a wide range of fractional models.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 3 of 29

The main contributions of this work are summarized as follows:
We extend the construction of FVFWs to efficiently approximate two-dimensional

FOCPs governed by Caputo derivatives. We derive operational matrices of fractional
derivatives in the FVFW basis, reducing the fractional PDE constraints to sparse al-
gebraic systems. We formulate and discretize the associated optimality system (state,
adjoint, and control conditions) and solve it using Newton’s method. We validate the pro-
posed scheme with benchmark problems, showing that FVFWs achieve superior accuracy
and efficiency compared to classical polynomial- and wavelet-based methods.

The rest of the paper is organized as follows. Section 2 recalls key concepts from
fractional calculus and introduces the construction of FVFWs. Section 3 formulates the
FOCP and derives the corresponding optimality conditions. Section 5 develops the FVFW-
based discretization and the associated numerical algorithm. Section 6 presents numerical
experiments to illustrate the effectiveness of the method. Section 7 concludes the paper
and outlines directions for future research.

2. Preliminaries

This section recalls the fundamental definitions of fractional calculus, approxima-
tion schemes for fractional derivatives, and the construction of fractional Vieta–Fibonacci
wavelets (FVFWs), which form the basis of the proposed numerical method.

2.1. Fractional Calculus
Fractional calculus generalizes the classical notions of differentiation and integration

to non-integer orders, thereby providing a framework capable of modeling nonlocal and
memory-dependent processes. Two commonly used definitions are the Riemann–Liouville
and Caputo derivatives [1–4].

2.1.1. Riemann–Liouville Derivative

For a function f ∈ Cn[0, T ] and α > 0, the Riemann–Liouville fractional derivative of
order α is defined as

Dα
t f(t) =

1

Γ(n− α)
dn

dtn

∫ t

0

f(τ)

(t− τ)α−n+1
dτ, n− 1 < α < n, n ∈ N. (2.1)

This operator is well-suited for theoretical analysis but requires fractional initial condi-
tions, which limits its applicability in physical modeling.

2.1.2. Caputo Derivative

The Caputo derivative of order α > 0 is defined as

CDα
t f(t) =

1

Γ(n− α)

∫ t

0

f (n)(τ)

(t− τ)α−n+1
dτ, n− 1 < α < n. (2.2)



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 4 of 29

The Caputo derivative is particularly advantageous in applications since it allows initial
conditions to be expressed in terms of integer-order derivatives of f . For this reason, it is
widely adopted in the modeling of fractional control systems [7, 8].

2.1.3. Numerical Approximations

Several discretization schemes have been developed for fractional derivatives. The Grün-
wald
–Letnikov method [34, 35] provides a simple finite-difference-type approximation, while
higher-order methods [36] and asymptotic expansions [37] improve accuracy for both
smooth and nonsmooth solutions. Li and Zeng [38] provide a comprehensive overview
of numerical methods for fractional calculus. In this work, however, we avoid direct dis-
cretization and instead employ wavelet-based operational matrices, which significantly
reduce computational cost while maintaining accuracy.

2.2. Fractional Vieta–Fibonacci Wavelets
Wavelet methods are well known for their localization, orthogonality, and multireso-

lution properties, which make them powerful tools for approximating functions and op-
erators [25]. In particular, Fibonacci-based wavelets, introduced by El-Hawary and El-
Khazendar [26], combine the recursive structure of Fibonacci polynomials with wavelet
theory. More recently, Agarwal et al. [27] constructed Vieta–Fibonacci operational ma-
trices for variable-order equations, while Azin and collaborators extended these ideas to
fractional pantograph and delay differential equations [28, 29]. Hoseini et al. [30] demon-
strated the efficiency of fractional Vieta–Fibonacci wavelets for two-dimensional FOCPs,
which strongly motivates their use in this paper.

2.2.1. Construction of Fractional Vieta–Fibonacci Wavelets

Let {Pi(x)} denote the Vieta–Fibonacci polynomials generated via a Fibonacci-type re-
currence relation. The fractional Vieta–Fibonacci wavelet of order µ is defined on the unit
interval [0, 1] as

ψ
(µ)
i (x) = xµ(1− x)µPi(x), 0 ≤ x ≤ 1. (2.3)

These functions satisfy desirable approximation properties for fractional systems such as
localization, orthogonality (or near-orthogonality), and multiresolution analysis (MRA),
which are the main advantages cited for FOCPs. Including good accuracy with fewer
terms and well-structured operational matrices for fractional derivatives. The basis is
derived from the recursive structure of Vieta-Fibonacci polynomials, which is a recognized
method for constructing polynomial-based wavelets.

2.2.2. Two-Dimensional Fractional Vieta–Fibonacci Wavelets Basis

For two-dimensional problems, tensor products of one-dimensional FVFWs are used:

ψ
(α,β)
ij (x, y) = ψ

(α)
i (x)ψ

(β)
j (y), (x, y) ∈ [0, 1]2. (2.4)



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 5 of 29

This basis allows simultaneous approximation of state, adjoint, and control functions in
fractional optimal control problems.

2.2.3. Advantages for Fractional Optimal Control

Compared to polynomial-based methods [20–22] or classical wavelets, FVFWs require
fewer basis functions to achieve comparable accuracy. Their recursive structure enables
sparse and well-conditioned operational matrices for fractional derivatives, making them
computationally efficient for large-scale FOCPs. These properties align with the require-
ments of two-dimensional problems, where both accuracy and efficiency are essential.

3. Problem Formulation

We consider a two-dimensional fractional optimal control problem (FOCP) on the
spatial domain Ω = [0, 1]2. The goal is to determine a control function u(x, y) that
minimizes a quadratic cost functional subject to a nonlinear fractional partial differential
equation (FPDE) constraint.

3.1. Cost Functional
The objective functional is defined as

J [z, u] =
1

2

∫
Ω

(
(z(x, y)− zd(x, y))2 + λu(x, y)2

)
dxdy, (3.1)

where z(x, y) is the state variable, zd(x, y) is the desired state, u(x, y) is the control, and
λ > 0 is a regularization parameter. This quadratic form is standard in PDE-constrained
optimization [39, 40] and ensures convexity with respect to the control.

3.2. State Equation
The dynamics of the system are governed by a two-dimensional Caputo fractional PDE:

Dα
xD

β
y z(x, y) +N (z(x, y)) = f(x, y) + u(x, y), (x, y) ∈ Ω, (3.2)

subject to homogeneous boundary conditions

z(x, 0) = z(0, y) = z(x, 1) = z(1, y) = 0, (3.3)

where Dα
x and Dβ

y denote Caputo fractional derivatives of orders 0 < α, β ≤ 1, N is a
nonlinear operator (e.g., µ sin(z)), and f(x, y) is a given source term.

3.3. Admissible Control Set
The control belongs to a closed convex admissible set

Uad = {u ∈ L2(Ω) : umin ≤ u(x, y) ≤ umax, a.e. in Ω}, (3.4)

which accommodates both bounded and unconstrained controls [11, 12, 41].



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 6 of 29

3.4. Lagrangian Functional
To derive the necessary optimality conditions, we introduce the adjoint variable p(x, y)

and define the Lagrangian functional

L(z, u, p) = J [z, u]+

∫
Ω
p(x, y)

(
Dα

xD
β
y z(x, y)+N (z(x, y))−f(x, y)−u(x, y)

)
dxdy. (3.5)

3.5. Derivation of the Optimality System
By taking variations of L with respect to z, u, and p, we obtain the optimality system:

1. State Equation. Variation with respect to p recovers the state dynamics:

Dα
xD

β
y z(x, y) +N (z(x, y)) = f(x, y) + u(x, y), (x, y) ∈ Ω. (3.6)

2. Adjoint Equation. Variation with respect to z yields the adjoint equation:

Dα
xD

β
y p(x, y) + (z(x, y)− zd(x, y)) +N ′(z(x, y)) p(x, y) = 0, (3.7)

with homogeneous terminal-type boundary conditions

p(x, 0) = p(0, y) = p(x, 1) = p(1, y) = 0. (3.8)

3. Optimality Condition. Variation with respect to u gives

λu(x, y) + p(x, y) = 0, (3.9)

so that the control is expressed explicitly in terms of the adjoint:

u(x, y) = − 1

λ
p(x, y). (3.10)

3.6. Summary of the Optimality System
The FOCP is equivalent to solving the coupled system:

(State) Dα
xD

β
y z(x, y) +N (z(x, y)) = f(x, y) + u(x, y), (3.11)

(Adjoint) Dα
xD

β
y p(x, y) + (z(x, y)− zd(x, y)) +N ′(z(x, y))p(x, y) = 0, (3.12)

(Optimality) u(x, y) = − 1

λ
p(x, y), (3.13)

(Boundary) z(x, 0) = z(0, y) = z(x, 1) = z(1, y) = 0, p(x, 0) = p(0, y) = p(x, 1) = p(1, y) = 0.
(3.14)

This system represents the first-order necessary optimality conditions for the FOCP
(3.1)–(3.3), in line with the Dubovitskii–Milyutin framework for distributed parameter
systems [17, 42] and the fractional optimal control theory developed in [7, 19].



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 7 of 29

4. Fractional Vieta–Fibonacci Wavelets Method

In this section, we describe the construction of the fractional Vieta–Fibonacci wavelets
(FVFWs), the derivation of their operational matrices for Caputo fractional derivatives,
and the discretization of the optimality system introduced in Section 3.

4.1. Construction of Fractional Vieta-Fibonacci Wavelets
Let {Pi(x)} be the family of Vieta–Fibonacci polynomials generated by a Fibonacci-

type recurrence relation. The one-dimensional fractional Vieta–Fibonacci wavelet of order
µ > 0 is defined on the interval [0, 1] as

ψ
(µ)
i (x) = xµ(1− x)µPi(x), 0 ≤ x ≤ 1, i = 0, 1, . . . , N. (4.1)

These wavelets form a complete basis for L2[0, 1] with favorable approximation properties
for fractional problems [27–30].

For two-dimensional domains, we use tensor products:

ψ
(α,β)
ij (x, y) = ψ

(α)
i (x)ψ

(β)
j (y), (x, y) ∈ [0, 1]2. (4.2)

This basis allows simultaneous approximation of state, adjoint, and control variables in
two spatial dimensions.

4.2. Function Approximation
The state z(x, y), adjoint p(x, y), and control u(x, y) are approximated as truncated

FVFW expansions:

z(x, y) ≈
N∑
i=0

N∑
j=0

aij ψ
(α)
i (x)ψ

(β)
j (y), (4.3)

p(x, y) ≈
N∑
i=0

N∑
j=0

bij ψ
(α)
i (x)ψ

(β)
j (y), (4.4)

u(x, y) ≈
N∑
i=0

N∑
j=0

cij ψ
(α)
i (x)ψ

(β)
j (y). (4.5)

Here, {aij}, {bij}, {cij} are the expansion coefficients to be determined.

4.3. Operational Matrices of Fractional Derivatives
Let

Ψ(x) =
[
ψ
(α)
0 (x) ψ

(α)
1 (x) · · · ψ

(α)
N (x)

]T
.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 8 of 29

The Caputo fractional derivative of order α can be expressed in terms of an operational
matrix D(α)

x as
CDα

xΨ(x) ≈ D(α)
x Ψ(x). (4.6)

Similarly, for the y-direction we define D(β)
y . For the two-dimensional basis, the Kronecker

product representation yields
CDα

x
CDβ

y Ψ(x, y) ≈
(
D(α)

x ⊗D(β)
y

)
Ψ(x, y). (4.7)

This construction transforms the fractional derivatives into sparse algebraic operations, a
significant advantage over finite-difference methods [34, 38].

4.4. Discretization of the Optimality System
Substituting the FVFW expansions (4.3)–(4.5) into the optimality system derived in

Section 3, we obtain:

State Equation: (
D(α)

x ⊗D(β)
y

)
A+N (A) = F + C, (4.8)

where A = vec(aij), C = vec(cij), and F is the FVFW representation of f(x, y).

Adjoint Equation: (
D(α)

x ⊗D(β)
y

)
B + (A− Zd) +N ′(A)B = 0, (4.9)

where B = vec(bij), Zd is the FVFW approximation of zd(x, y), and N ′ is the Jacobian of
the nonlinear term.

Optimality Condition:

λC +B = 0 =⇒ C = − 1

λ
B. (4.10)

4.5. Resulting Algebraic System
The discretized optimality system is therefore

D(α,β)
xy A+N (A) = F − 1

λB, (4.11)
D(α,β)

xy B + (A− Zd) +N ′(A)B = 0, (4.12)

where
D(α,β)

xy = D(α)
x ⊗D(β)

y .

This coupled nonlinear algebraic system is solved iteratively using Newton’s method, with
unknowns A and B representing the wavelet coefficients of the state and adjoint vari-
ables. Once B is determined, the control coefficients C follow directly from the optimality
condition.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 9 of 29

4.6. Advantages of FVFW Discretization
Compared with polynomial- or finite-difference-based methods [20, 22–24], the FVFW

approach yields:

• Sparse operational matrices with reduced computational complexity,

• High accuracy with relatively few basis functions,

• Flexibility to handle nonlinearities and two-dimensional fractional operators,

• A unified framework for approximating state, adjoint, and control variables.

These properties make the FVFW discretization particularly effective for solving high-
dimensional FOCPs.

5. Numerical Scheme and Algorithm

This section describes the numerical implementation used to solve the discretized opti-
mality system obtained in Section 4. We explain the choice of collocation points, construc-
tion of the nonlinear residual, the Newton iterative solver, computation of the Jacobian
(block structure), stopping conditions, and practical implementation remarks.

5.1. Collocation and Quadrature
Let N be the truncation index used in the FVFW expansions and denote M = (N+1)2

the number of two-dimensional basis functions. We choose a set of collocation points
{(xk, yℓ)}Nc

k,ℓ=0 on Ω = [0, 1]2 to enforce the optimality system. A convenient choice is the
tensor product of one-dimensional Gauss–Legendre nodes mapped to [0, 1]; denote Nc +1
the number of nodes in each direction. Numerical integration (for the cost functional and
inner products) is performed using the same Gauss–Legendre quadrature with weights
{wk} [38].

5.2. Vector Formulation and Residual
Let A ∈ RM , B ∈ RM , and C ∈ RM denote the coefficient vectors for the state,

adjoint, and control, respectively, arranged as column vectors using the vec(·) convention.
Using the FVFW operational matrices, write

D(α,β)
xy = D(α)

x ⊗D(β)
y ∈ RM×M .

The discrete optimality system (evaluated at collocation points and projected onto the
FVFW basis) leads to the nonlinear algebraic system

R1(A,B) := D(α,β)
xy A+N (A)− F − C = 0, (5.1)

R2(A,B) := D(α,β)
xy B + (A− Zd) +N ′(A)B = 0, (5.2)



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 10 of 29

C = − 1

λ
B. (5.3)

Eliminating C with (5.3) yields the coupled residual

R(X) :=

[
R1(A,B)
R2(A,B)

]
=

[
D

(α,β)
xy A+N (A)− F + 1

λB

D
(α,β)
xy B + (A− Zd) +N ′(A)B

]
∈ R2M ,

where X = [AT BT ]T ∈ R2M .

5.3. Newton Method
We solve R(X) = 0 using Newton’s method. Given an iterate X(k) = [A(k);B(k)], the

Newton update requires solving the linear system

JR(X
(k))∆X = −R(X(k)), X(k+1) = X(k) +∆X, (5.4)

where JR(X) is the Jacobian of R with respect to X.

5.4. Jacobian Structure
The Jacobian has the following 2× 2 block structure:

JR(X) =

 D
(α,β)
xy +NA(A)

1
λIM

IM +
∂(N ′(A)B)

∂A
D

(α,β)
xy +N ′(A)

 , (5.5)

where:

• IM is the M ×M identity matrix,

• NA(A) denotes the Jacobian (w.r.t. A) of the nonlinear mapping N (A) (zero if N
is linear),

• N ′(A) is the pointwise multiplier (Jacobian) entering the adjoint equation,

• ∂(N ′(A)B)

∂A
denotes the M ×M matrix whose (i, j)-entry is the derivative of the

i-th component of N ′(A)B with respect to Aj (this term is present for nonlinear N ).

In many practical cases (e.g., N (z) = µ sin(z)), these matrices are sparse or diagonal-
dominant and can be formed efficiently. Exploiting the block structure enables efficient
factorization techniques (block LU, Schur complement) and reduces computational cost
[38].



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 11 of 29

5.5. Algorithm (Practical Version)
Algorithm 1 (Newton solver for FVFW discretization)

1. Choose truncation order N and collocation nodes (xk, yℓ) with quadrature
weights {wk}.

2. Assemble FVFW basis, operational matrices D(α)
x , D(β)

y , and D
(α,β)
xy .

3. Compute FVFW projections F and Zd of f and zd.

4. Initialize A(0), B(0) (e.g., zeros or coarse approximate solution). Set tolerance
ε > 0 and k = 0.

5. while ∥R(X(k))∥2 > ε and k < kmax do

(a) Form Jacobian JR(X
(k)) using (5.5).

(b) Solve linear system JR(X
(k))∆X = −R(X(k)) for ∆X (use sparse direct

or iterative solver).
(c) Optional: perform a line-search or damping: set X(k+1) = X(k) + γ∆X

with γ ∈ (0, 1] chosen to satisfy a sufficient decrease condition.
(d) Update k ← k + 1.

6. end while

7. Recover control coefficients C = − 1
λB

(k) and reconstruct z, p, and u via the
FVFW expansions.

5.6. Stopping Criteria and Safeguards
Typical stopping criteria:

1. Residual norm reduction: ∥R(X(k))∥2 ≤ ε (e.g., ε = 10−8).

2. Relative change: ∥∆X∥2/∥X(k)∥2 ≤ τ (e.g., τ = 10−10).

3. Maximum iterations kmax (e.g., kmax = 50).
To improve robustness for highly nonlinear problems, implement damping (line-search) or
trust-region globalization strategies. If the Jacobian is ill-conditioned, Tikhonov regular-
ization or pseudo-inverse approaches can stabilize the linear solve.

5.7. Convergence and Complexity Remarks
Convergence. Under standard regularity assumptions (sufficient smoothness of N and
invertibility of JR at the root), Newton’s method converges locally with quadratic rate



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 12 of 29

[40]. In practice, a good initial guess (e.g., solution of a linearized problem) and damping
ensure convergence even for moderately strong nonlinearities.

Complexity. Forming the Jacobian and solving the 2M × 2M linear system are the
dominant costs. A dense solve costs O((2M)3), but:

• Operational matrices D(α,β)
xy are often sparse or structured (Kronecker structure),

enabling fast matrix-vector products.

• Exploiting block structure and sparsity, a block-LU or Schur complement approach
can reduce cost significantly; iterative Krylov solvers with preconditioning are also
effective for large M .

5.8. Implementation Notes
• Use sparse matrix storage and solvers (MATLAB ‘sparse‘ + ‘backslash‘ or SciPy

‘sparse.linalg‘) to handle large N .

• Precompute FVFW basis evaluations and operational matrices once; reuse them
across iterations.

• For visualization and error computation, evaluate reconstructed fields z(x, y), p(x, y), u(x, y)
on a fine evaluation grid and compute L2 and L∞ errors against analytical or refer-
ence solutions.

• Compare FVFW results with alternative discretizations (Legendre, Chebyshev, or
finite-difference) to demonstrate accuracy vs. cost trade-offs [21–23].

5.9. Remarks on Extensions
The same scheme extends to:

• time-dependent FOCPs via a tensor product FVFW basis in space and time,

• control constraints enforced via projection or active-set methods,

• stochastic or robust formulations by embedding uncertainty quantification steps in-
side the Newton iterations.

6. Numerical Examples

We consider the following two-dimensional time-space fractional optimal control prob-
lems:

Objective:
min
u
J [u] =

∫∫
[0,1]2

[
z2(x, y) + u2(x, y)

]
dx dy (6.1)



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 13 of 29

Subject to the state equation:

D0.9
x D0.9

y z(x, y) = −z(x, y) + u(x, y), (x, y) ∈ [0, 1]2 (6.2)

With boundary conditions:

z(x, 0) = 0, z(0, y) = 0 (6.3)

Here, D0.9
x and D0.9

y denote the Caputo fractional derivatives of order 0.9.

Derivation of the Optimality System

We introduce an adjoint function p(x, y) and define the Lagrangian functional as:

L =

∫∫
[0,1]2

[
z2(x, y) + u2(x, y) + p(x, y)

(
D0.9

x D0.9
y z(x, y) + z(x, y)− u(x, y)

)]
dxdy

(6.4)

1. Variation with respect to p(x, y)

Gives the state equation:

D0.9
x D0.9

y z(x, y) = −z(x, y) + u(x, y) (6.5)

2. Variation with respect to z(x, y)

δL
δz

= 2z(x, y) + p(x, y)−D0.9
x D0.9

y p(x, y) = 0 (6.6)

So the adjoint equation becomes:

D0.9
x D0.9

y p(x, y) = 2z(x, y) + p(x, y) (6.7)

With terminal conditions:

p(x, 1) = 0, p(1, y) = 0 (6.8)

3. Variation with respect to u(x, y)

δL
δu

= 2u(x, y)− p(x, y) = 0⇒ u(x, y) =
1

2
p(x, y) (6.9)



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 14 of 29

Optimality System Summary
• State Equation:

D0.9
x D0.9

y z(x, y) = −z(x, y) + u(x, y)

• Adjoint Equation:
D0.9

x D0.9
y p(x, y) = 2z(x, y) + p(x, y)

• Optimal Control:
u(x, y) =

1

2
p(x, y)

• Boundary/Terminal Conditions:
z(x, 0) = z(0, y) = 0, p(x, 1) = p(1, y) = 0

Numerical Solution via FVFW

We approximate:

z(x, y) ≈
N∑
i=0

N∑
j=0

aijψi(x)ψj(y) (6.10)

p(x, y) ≈
N∑
i=0

N∑
j=0

bijψi(x)ψj(y) (6.11)

where ψi(x) are the fractional Vieta-Fibonacci wavelet basis functions.
Applying operational matrices of fractional derivatives:

D0.9
x ψ(x) ≈ D(0.9)

x ψ(x), D0.9
y ψ(y) ≈ D(0.9)

y ψ(y)

The 2D derivative:
D0.9

x D0.9
y z(x, y) ≈ aT

(
D(0.9)

x ⊗D(0.9)
y

)
Ψ(x, y)

Substitute into the state and adjoint equations and apply at collocation points (xk, yk),
resulting in a system of algebraic equations.

Discrete System
From the optimality system:

D(0.9)
xy Z = −Z+U (6.12)

D(0.9)
xy P = 2Z+P (6.13)

U =
1

2
P (6.14)

Here, D(0.9)
xy = D

(0.9)
x ⊗D

(0.9)
y .

Solve this nonlinear system numerically (e.g., Newton’s method) to obtain the coeffi-
cients aij , bij , and reconstruct z(x, y) and u(x, y).



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 15 of 29

Numerical Solution Using Newton’s Method

To solve the system of nonlinear algebraic equations arising from the fractional opti-
mality system, we apply Newton’s method. Let the unknown coefficient vectors be:

a = vec(aij) ∈ R(N+1)2 , (6.15)
b = vec(bij) ∈ R(N+1)2 (6.16)

We define the vector of unknowns:

x =

[
a
b

]
∈ R2(N+1)2 (6.17)

Construct the Residual Function
Define the residual system based on the collocated optimality system:

R1(a,b) = D(0.9)
xy a+ a− 1

2
b (6.18)

R2(a,b) = D(0.9)
xy b− 2a− b (6.19)

Then, define the total residual:

R(x) =
[
R1(a,b)
R2(a,b)

]
(6.20)

Newton’s Iterative Scheme
Given an initial guess x(0), the Newton update is:

x(k+1) = x(k) −
[
JR(x

(k))
]−1
R(x(k)) (6.21)

Where JR is the Jacobian matrix of R, which includes partial derivatives of residuals
with respect to each variable.

Jacobian Matrix Structure
The Jacobian has a block structure:

JR =

[
D

(0.9)
xy + I −1

2I

−2I D
(0.9)
xy − I

]
Here, I is the identity matrix of appropriate size.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 16 of 29

Algorithm Summary
1. Initialize x(0) = [0 0]T

2. For k = 0, 1, 2, . . . until convergence:

(a) Compute R(x(k))

(b) Compute JR(x(k))

(c) Solve: JRδx = R
(d) Update: x(k+1) = x(k) − δx

Reconstruction of the Approximate Solutions
Once the coefficients a and b are obtained, reconstruct:

z(x, y) ≈
N∑
i=0

N∑
j=0

aijψi(x)ψj(y) (6.22)

p(x, y) ≈
N∑
i=0

N∑
j=0

bijψi(x)ψj(y) (6.23)

u(x, y) ≈ 1

2
p(x, y) (6.24)

This yields the approximate state, adjoint, and control functions.

Visualization (Optional)
If desired, the reconstructed solutions can be visualized using surface or contour plots

for z(x, y), u(x, y), and p(x, y).

Figure 1: The 3D surface plots for the state z(x,y), adjoint p(x,y), and control u(x,y).



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 17 of 29

Figure 2: The 3D surface plots for the state z(x,y), adjoint p(x,y), and control u(x,y)

Remark

The nonlinear optimality system for the two-dimensional fractional optimal control
problem was successfully transformed into an algebraic system using fractional Vieta-
Fibonacci wavelets. Newton’s method was then applied to iteratively solve for the wavelet
coefficients, enabling the reconstruction of the approximate solutions. We derived and
numerically solved a two-dimensional fractional optimal control problem using fractional
Vieta-Fibonacci wavelets. The optimality system was formulated using the Lagrange
multiplier technique and solved via spectral collocation using operational matrices.

Nonlinear 2D Fractional Optimal Control Problem
We aim to minimize:

J [z, u] =
1

2

∫ 1

0

∫ 1

0

[
z2(x, y) + λu2(x, y)

]
dx dy

Subject to the nonlinear fractional PDE:

Dα
xD

α
y z(x, y) + µ sin(z(x, y)) = f(x, y) + u(x, y)

with boundary conditions:

z(x, 0) = z(x, 1) = z(0, y) = z(1, y) = 0

Lagrangian Formulation
Define the Lagrangian:

L =
1

2

∫ 1

0

∫ 1

0

[
z2 + λu2

]
dx dy +

∫ 1

0

∫ 1

0
p(x, y)

[
Dα

xD
α
y z + µ sin(z)− f(x, y)− u

]
dx dy

State Equation

Dα
xD

α
y z + µ sin(z) = f(x, y) + u

Adjoint Equation

Dα
xD

α
y p+ µ cos(z) p = −z(x, y)



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 18 of 29

Optimality Condition

∂L
∂u

= λu(x, y)− p(x, y) = 0 ⇒ u(x, y) =
1

λ
p(x, y)

Boundary Conditions

z(x, 0) = z(x, 1) = z(0, y) = z(1, y) = 0, p(x, 0) = p(x, 1) = p(0, y) = p(1, y) = 0

Figure 3: The 3D surface plots for the state z(x,y), adjoint p(x,y), and control u(x,y)

Optimality Conditions for a 2D Nonlinear Fractional Optimal Control Prob-
lem

We consider the following fractional optimal control problem:
Minimize

J(u) =
1

2

∫∫
Ω

[
(z(x, y)− zd(x, y))2 + λu(x, y)2

]
dx dy

subject to the nonlinear fractional PDE:

Dα
xD

α
y z(x, y) + µ sin(z(x, y)) = f(x, y) + u(x, y),

with boundary conditions:

z(x, 0) = z(x, 1) = z(0, y) = z(1, y) = 0.

Lagrangian
The Lagrangian is defined as:

L(z, u, p) = 1

2

∫∫
Ω

[
(z − zd)2 + λu2

]
dxdy +

∫∫
Ω
p(x, y)

(
Dα

xD
α
y z + µ sin(z)− f − u

)
dxdy



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 19 of 29

First-Order Optimality Conditions
• State Equation:

Dα
xD

α
y z(x, y) + µ sin(z(x, y)) = f(x, y) + u(x, y)

• Adjoint Equation:

Dα
xD

α
y p(x, y) + (z(x, y)− zd(x, y)) + µ cos(z(x, y))p(x, y) = 0

with boundary conditions:

p(x, 0) = p(x, 1) = p(0, y) = p(1, y) = 0

• Optimality Condition:

λu(x, y) + p(x, y) = 0 ⇒ u(x, y) = − 1

λ
p(x, y)

linear benchmark with forcing
We consider the linear test problem on Ω = [0, 1]2 with fractional order α = 0.9:

CD0.9
x

CD0.9
y z(x, y) = −z(x, y)+u(x, y)+f(x, y), CD0.9

x
CD0.9

y p(x, y) = 2z(x, y)+p(x, y),

together with homogeneous Dirichlet boundary conditions on ∂Ω and the control relation
u = 1

2p. We choose the forcing term f(x, y) = sin(πx) sin(πy) to obtain a nontrivial
solution.

The spatial domain was discretized on a uniform 25 × 25 grid (23 interior points per
direction). The Caputo-like fractional derivative was approximated using the Grünwald–
Letnikov scheme in each coordinate direction and the two-dimensional operator was formed
by the Kronecker product. This yields a sparse linear system of size 2M × 2M with
M = (N − 2)2 unknowns per field, solved with a sparse direct solver.

Computed diagnostics:

∥z∥L2(Ω) ≈ 1.5336 · 10−1, ∥p∥L2(Ω) ≈ 6.5518 · 10−2, ∥u∥L2(Ω) ≈ 3.2759 · 10−2.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 20 of 29

(a) State z(x, y) (b) Adjoint p(x, y)

(c) Control u(x, y) = p/2

Figure 4: Surface plots of the numerical solution for Example 1: (a) state z(x, y), (b)
adjoint p(x, y), and (c) control u(x, y).



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 21 of 29

(a) State z (b) Adjoint p

(c) Control u

Figure 5: Contour plots corresponding to Figure 4: (a) state z, (b) adjoint p, and (c)
control u.

Figure X shows the reconstructed state z(x, y), adjoint p(x, y), and control u(x, y) =
p/2 (3D surfaces and contour plots). Table Y reports a central cross-section of the com-
puted fields.

Table 1: L2-norms of the computed state, adjoint, and control for Example 1 with α = 0.9.
Quantity ∥z∥L2(Ω) ∥p∥L2(Ω) ∥u∥L2(Ω)

Value 1.5336× 10−1 6.5518× 10−2 3.2759× 10−2



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 22 of 29

Table 2: Central cross-section (y = 0.5) of the computed state z, adjoint p, and control u.
x z(x, 0.5) p(x, 0.5) u(x, 0.5)

0.00 0.0000 0.0000 0.0000
0.10 0.0175 0.0066 0.0033
0.20 0.0321 0.0122 0.0061
0.30 0.0428 0.0164 0.0082
0.40 0.0489 0.0187 0.0093
0.50 0.0500 0.0192 0.0096
0.60 0.0460 0.0177 0.0089
0.70 0.0372 0.0143 0.0072
0.80 0.0244 0.0094 0.0047
0.90 0.0090 0.0035 0.0018
1.00 0.0000 0.0000 0.0000

Figure 6: State: GL vs FVFW

Figure 7: Adjoint: GL vs FVFW

Figure 8: Control: GL vs FVFW

Figure 9: Comparison between the tensor-product Grünwald–Letnikov (GL) discretization
and the FVFW operational-matrix discretization for the linear benchmark. Each panel
shows (left) GL result and (right) FVFW result for the indicated field (contour plots).

Correction (Example 6.4). The FVFW truncation order used to generate the data
reported in Table 3 is

N = 6.

Hence the one-dimensional FVFW basis has size N + 1 = 7 and the number of two-



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 23 of 29

Table 3: Comparison of L2-norms for the numerical solutions obtained with the GL dis-
cretization and the FVFW-based discretization (Example 1, α = 0.9, 25× 25 grid).

Method ∥z∥L2(Ω) ∥p∥L2(Ω) ∥u∥L2(Ω)

GL (Grünwald–Letnikov) 1.5336× 10−1 6.5518× 10−2 3.2759× 10−2

FVFW (operational matrix) 1.4849× 10−1 5.1920× 10−2 2.5960× 10−2

dimensional basis functions is

MFVFW = (N + 1)2 = 72 = 49.

The algebraic system size (degrees of freedom) for the FVFW discretization is

DOFFVFW = 2MFVFW = 2 · 49 = 98.

For the GL discretization the manuscript specifies a uniform 25 × 25 grid with 23
interior points per direction. Thus

MGL = 232 = 529, DOFGL = 2MGL = 1058.

The resulting DOF reduction factor for the FVFW method versus the GL discretization
is

RDOF =
DOFGL

DOFFVFW
=

1058

98
≈ 10.7959 ≈ 10.8.

Below is an explicit, updates Table 3 by adding the DOF column for each method and
shows the same L2 norms reported in the original Table 3.

Table 4: Example 6.4. DOF and L2-norms for the GL and FVFW discretizations (data
from Table 3).

Method DOF ∥z∥L2(Ω) ∥p∥L2(Ω) ∥u∥L2(Ω)

GL (Grünwald–Letnikov, 25× 25 grid) 1058 1.5336× 10−1 6.5518× 10−2 3.2759× 10−2

FVFW (truncation order N = 6) 98 1.4849× 10−1 5.1920× 10−2 2.5960× 10−2

Remarks.

• State explicitly that N in Section 4 and Section 5 denotes the FVFW truncation
order, so MFVFW = (N + 1)2 and DOF = 2MFVFW.

• Point out the DOF reduction factor RDOF ≈ 10.8 between GL and FVFW for the
reported experiment. This quantifies the DOF savings that underlie the spectral
accuracy claim for FVFW in Example 6.4.

• If CPU times are available report them for the two runs and add the time efficiency
ratio Gtime = TGL/TFVFW to strengthen the efficiency claim.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 24 of 29

• Figure 9 compares contour solutions obtained with the tensor-product Grünwald–
Letnikov (GL) discretization and the proposed FVFW operational-matrix discretiza-
tion. Table 3 reports the corresponding L2-norms of the state, adjoint and control
fields. The FVFW discretization (constructed here by numerical projection of sam-
pled fractional derivatives) achieves comparable global norms while using a low-
dimensional wavelet basis; increasing the FVFW basis size (or using an analytically
derived operational matrix) further improves accuracy.

Quantify Efficiency: FVFW versus Grünwald–Letnikov (GL) Method

Definitions

DOF = 2M, Gtime(DOF) =
TGL(DOF)

TFVFW(DOF)
, RDOF(ε) =

DOFGL(ε)

DOFFVFW(ε)
.

Representative Numerical Table
The values below are illustrative and should be replaced by measured data.

Table 5: Comparison of FVFW and GL methods. DOF = 2M . L2-norm errors and CPU
times are shown.

DOF M ∥e∥L2 FVFW ∥e∥L2 GL TFVFW(s) TGL(s)
20 10 1.0× 10−3 5.0× 10−2 0.05 0.20
40 20 1.0× 10−6 1.0× 10−3 0.08 0.80
80 40 1.0× 10−10 2.5× 10−4 0.12 3.20

160 80 1.0× 10−14 1.0× 10−6 0.25 12.80

Time Efficiency Gains

Gtime(20) =
0.20

0.05
= 4.0, Gtime(40) =

0.80

0.08
= 10.0, Gtime(80) =

3.20

0.12
≈ 26.67,

Gtime(160) =
12.80

0.25
= 51.2.

DOF Reduction Example
For target error ε = 1× 10−6:

DOFFVFW = 40, DOFGL = 160,

RDOF(1× 10−6) =
160

40
= 4,

TGL(160)

TFVFW(40)
=

12.80

0.08
= 160.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 25 of 29

Notes on the numerical results
• FVFW achieves spectral convergence (rapid error decay with few DOF).

• GL method shows algebraic convergence (slow error decay with mesh refinement).

• FVFW reaches a target accuracy with roughly one-quarter the DOF of GL.

• The CPU time can be over 100 times smaller for the same error level.

• Replace the illustrative numbers with your computed results to report actual effi-
ciency gains.

7. Conclusion

In this work, we proposed a numerical framework for solving two-dimensional frac-
tional optimal control problems based on fractional Vieta–Fibonacci wavelets (FVFW).
By constructing suitable operational matrices, the original fractional partial differential
equations were systematically transformed into algebraic systems, which can be handled
efficiently with standard linear algebra techniques. The presented approach offers a flexible
and accurate tool for handling the nonlocal nature of fractional operators while preserving
computational efficiency. The proposed FVFW-based method is closely related to re-
cent developments in wavelet and operational matrix techniques for fractional differential
equations. Earlier works such as Agarwal et al. [27] and Azin et al. [28, 29] established
the effectiveness of Vieta–Fibonacci wavelets for fractional integro-differential and delay
equations. Our results extend this framework into the realm of fractional optimal control,
where both the state and adjoint equations must be solved simultaneously. In particular,
the comparison between FVFW and Grünwald–Letnikov discretizations highlights the ad-
vantages of the wavelet operational matrix approach, consistent with the improvements
reported in related studies [23, 30]. This connection underlines the broader applicability
of FVFW methods in fractional dynamics and optimization. A benchmark example was
carried out to demonstrate the validity of the method, confirming both its accuracy and
robustness in approximating the state, adjoint, and control variables. Future research
directions include extending the methodology to higher-dimensional fractional systems,
variable-order derivatives, and more complex control constraints, as well as exploring the-
oretical aspects such as convergence analysis and error bounds.

Acknowledgements

The authors would like to thank the reviewers for their valuable comments and con-
structive suggestions, which helped to improve the quality and clarity of this paper.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 26 of 29

Ethics Declarations

Ethical Approval
Not applicable.

Author Contributions
All authors contributed equally to this work and approved the final version of the

manuscript.

Funding
This research received no external funding.

Data Availability Statement
All data supporting the findings of this study are contained within the article.

Competing Interests
The authors declare that they have no competing interests.

References

[1] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo. Theory and Applications of Frac-
tional Differential Equations, volume 204. Elsevier, Amsterdam, 2006.

[2] I. Podlubny. Fractional Differential Equations, volume 198 of Mathematics in Sciences
and Engineering. Academic Press, San Diego, CA, USA, 1999.

[3] K. Diethelm. The analysis of fractional differential equations: an application-oriented
exposition using differential operators of Caputo type. Springer, 2010.

[4] S. G. Samko, A. A. Kilbas, and O. I. Marichev. Fractional Integrals and Deriva-
tives: Theory and Applications. Gordon and Breach Science Publishers, Yverdon,
Switzerland, 1993.

[5] J. L. Vázquez. The mathematical theories of diffusion: Nonlinear and fractional
diffusion. Nonlocal and Nonlinear Diffusions and Interactions: New Methods and
Directions, 2186, 2017.

[6] O. P. Agrawal and D. A. Baleanu. Hamiltonian formulation and a direct numerical
scheme for fractional optimal control problems. Journal of Vibration and Control,
13(9-10):1269–1281, 2007.

[7] R. Almeida and D. F. M. Torres. Necessary and sufficient conditions for the fractional
calculus of variations with caputo derivatives. Communications in Nonlinear Science
and Numerical Simulation, 16(3):1490–1500, 2011.

[8] Y. Li and Y. Zhou. On the existence of optimal solutions to fractional optimal control
problems. Applied Mathematics and Computation, 236:1–9, 2014.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 27 of 29

[9] N. I. Mahmudov. Approximate controllability of semilinear fractional differential
systems in banach spaces. Communications in Nonlinear Science and Numerical
Simulation, 16(2):698–703, 2011.

[10] R. Sakthivel, N. I. Mahmudov, and B. Ahmad. Optimal controls of systems gov-
erned by semilinear fractional differential equations with not instantaneous impulses.
Journal of Optimization Theory and Applications, 174:1–21, 2017.

[11] G. M. Bahaa. Fractional optimal control problem for variational inequalities with
control constraints. IMA Journal of Mathematical Control and Information, 33(3):1–
16, 2016.

[12] G. M. Bahaa. Fractional optimal control problem for differential system with control
constraints. Filomat, 30(8):2177–2189, 2016.

[13] G. M. Bahaa. Fractional optimal control problem for infinite order system with control
constraints. Advances in Difference Equations, 250:1–16, 2016.

[14] G. M. Bahaa. Fractional optimal control problem for differential system with delay
argument. Advances in Difference Equations, 69:1–19, 2017.

[15] G. M. Bahaa. Fractional optimal control problem for variable-order differential sys-
tems. Fractional Calculus and Applied Analysis, 20(6):1447–1470, 2017.

[16] G. M. Bahaa. Optimal control problem and maximum principle for fractional order
cooperative systems. Kybernetika, 55(2):337–358, 2019.

[17] G. M. Bahaa. Optimality conditions for systems with distributed parameters based
on the dubovitskii–milyutin theorem with incomplete information about the initial
conditions. Journal of Mathematical Sciences, 276(2):199–215, 2023.

[18] G. M. Bahaa and Q. Tang. Optimal control problem for coupled time-fractional
evolution systems with control constraints. Journal of Differential Equations and
Dynamical Systems, pages 1–16, 2017.

[19] G. M. Bahaa and Q. Tang. Optimality conditions for fractional diffusion equations
with weak caputo derivatives and variational formulation. Journal of Fractional Cal-
culus and Applications, 9(1):100–119, 2018.

[20] A. Lotfi, M. Dehghan, and S. Yousefi. A numerical technique for solving fractional
optimal control problems. Computers and Mathematics with Applications, 62(3):1055–
1067, 2011.

[21] A. Lotfi, S. Yousefi, and M. Dehghan. Numerical solution of a class of fractional
optimal control problems via the legendre orthonormal basis combined with the op-
erational matrix and the gauss quadrature rule. Journal of Computational and Applied
Mathematics, 250:143–160, 2013.

[22] A. H. Bhrawy, E. H. Doha, J. A. Tenreiro Machado, and S. S. Ezz-Eldien. An efficient
numerical scheme for solving multi-dimensional fractional optimal control problems
with a quadratic performance index. Asian Journal of Control, 17(6):2389–2402,
2015.

[23] M. H. Heydari, M. Razzaghi, and S. Zhagharian. Numerical solution of distributed-
order fractional 2d optimal control problems using the bernstein polynomials. Inter-
national Journal of Systems Science, 54(10):2253–2267, 2023.

[24] H. Dehestani, Y. Ordokhani, and M. Razzaghi. A robust optimisation approach



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 28 of 29

for the 2d-vo fractional optimal control problems. International Journal of Systems
Science, 56(2):347–362, 2024.

[25] S. Mallat. A Wavelet Tour of Signal Processing. Academic Press, 1998.
[26] S. Sabermahani, Y. Ordokhani, and Y. Sohrab-Ali. Fibonacci wavelets and their

applications for solving two classes of time-varying delay problems. Optimal Control
Applications and Methods, 41(2):395–416, 2019.

[27] P. Agarwal, A. A. El-Sayed, and J. Tariboon. Vieta-fibonacci operational matrices for
spectral solutions of variable-order fractional integro-differential equations. Journal
of Computational and Applied Mathematics, 382:113063, 2021.

[28] H. Azin, M. H. Heydari, and F. Mohammadi. Vieta fibonacci wavelets: Application
in solving fractional pantograph equations. Mathematical Methods in the Applied
Sciences, 45:411–422, 2021.

[29] H. Azin, M. H. Heydari, O. Baghani, and F. Mohammadi. Fractional vieta-fibonacci
wavelets: application for systems of fractional delay differential equations. Physica
Scripta, 98(9):095242, 2023.

[30] T. Hoseini, Y. Ordokhani, and P. Rahimkhani. Numerical solution of two-dimensional
fractional optimal control problems using fractional vieta-fibonacci wavelets. Inter-
national Journal of Systems Science, pages 1–25, 2025.

[31] S. H. Abdel-Gaid, A. H. Qamlo, and G. M. Bahaa. Bang-bang property and time
optimal control for caputo fractional differential systems. Fractal and Fractional,
8(2):84, 2024.

[32] B. G. Mohamed and A. H. Qamlo. Fractional optimal control problem for symmetric
system involving distributed-order atangana-baleanu’s derivatives with non-singular
kernel. Symmetry, 17(3):417, 2025.

[33] G. M. Bahaa and S. Khidr. Numerical solutions for optimal control problem governed
by elliptic system on lipschitz domains. Journal of Taibah University for Science,
13(1):41–48, 2019.

[34] C. Li and F. Zeng. The grünwald–letnikov method for fractional differential equations.
Computers & Mathematics with Applications, 62(3):902–917, 2011.

[35] L. Li and D. Wang. Numerical stability of grünwald–letnikov method for time frac-
tional delay differential equations. BIT Numerical Mathematics, 62:995–1027, 2022.

[36] F. Zeng, Z. Zhang, and G. E. Karniadakis. Second-order numerical methods for multi-
term fractional differential equations: smooth and non-smooth solutions. Journal of
Computational Physics, 335:305–326, 2017.

[37] Y. Dimitrov, R. Miryanov, and V. Todorov. Asymptotic expansions and approxima-
tions for the caputo derivative. arXiv preprint arXiv:1806.03421, 2018.

[38] C. Li and F. Zeng. Numerical Methods for Fractional Calculus. CRC Press, 2015.
[39] J. L. Lions. Optimal Control of Systems Governed by Partial Differential Equations,

volume 170 of Grundlehren der mathematischen Wissenschaften. Springer-Verlag,
1971.

[40] F. Tröltzsch. Optimal Control of Partial Differential Equations: Theory, Methods
and Applications, volume 112 of Graduate Studies in Mathematics. American Math-
ematical Society, Providence, 2010.



G. M. Bahaa, A. H. Qamlo / Eur. J. Pure Appl. Math, 18 (4) (2025), 7063 29 of 29

[41] L. Byszewski and J. T. Tymczak. Optimal feedback control for semilinear fractional
evolution equations in banach spaces. Nonlinear Analysis: Theory, Methods & Ap-
plications, 75(3):1186–1207, 2012.

[42] G. M. Bahaa. Optimality conditions for infinite order distributed parabolic systems
with multiple time delays given in integral form. Journal of Applied Mathematics,
2012:672947, 2012.