Development and Evaluation of Numerical Solvers for Mathematical Models with Fractional Derivatives (Разработка и оценивание численных решателей для математических моделей с дробными производными) тема диссертации и автореферата по ВАК РФ 00.00.00, доктор наук Хенди Ахмед Саид Абделазиз

• Хенди Ахмед Саид Абделазиз
• доктор наук
• 2026

Оглавление диссертации доктор наук Хенди Ахмед Саид Абделазиз

Contents

1 Introduction

2 Structure preserving numerical schemes for some fractional order mathematical models

2.1 A pseudo energy-invariant method for relativistic wave equations with Riesz space-fractional derivatives

2.1.1 Mathematical model

2.1.2 Discrete energy

2.1.3 Computer simulations

2.2 On the solution of hyperbolic two-dimensional fractional systems via discrete variational schemes

of high order of accuracy

2.2.1 Numerical method

2.2.2 Numerical results

2.3 A numerically efficient and conservative model for a Riesz space-fractional Klein-Gordon-Zakharov system

2.3.1 Preliminaries

2.3.2 Numerical model

2.3.3 Auxiliary lemmas

2.3.4 Structural properties

2.3.5 Numerical properties

2.4 Theoretical analysis of an explicit energy-conserving scheme for a fractional Klein-Gordon-Zakharov system

2.4.1 Numerical model

2.4.2 Physical properties

2.4.3 Numerical properties

2.4.4 Simulations

2.5 On the dissipativity of some Caputo time-fractional subdiffusion models in multiple dimensions: Theoretical and numerical investigations

2.5.1 Preliminaries

2.5.2 Single-term systems

2.5.3 Multi-term systems

2.5.4 Numerical simulations

3 Galerkin-Legendre spectral schemes for nonlinear time-space fractional partial differential equations

3.1 Semi-implicit Galerkin-Legendre spectral schemes for nonlinear time-space fractional diffusion-

reaction equations with smooth and nonsmooth solutions

3.1.1 Preliminaries

3.1.2 The Numerical Schemes

3.1.3 Some Auxiliary Lemmas

3.1.4 Theoretical Analysis

3.1.5 Numerical Simulations

3.2 Global consistency analysis of Ll-Galerkin spectral schemes for coupled nonlinear time-space fractional Schrodinger equations

3.2.1 Preliminaries

3.2.2 The numerical algorithms

3.2.3 Uniform L1-spectral scheme implementation for smooth solution

3.2.4 Nonuniform L1-spectral scheme implementation for the nonsmooth solution

3.2.5 Well-posedness of the numerical solution

3.2.6 Convergence analysis

3.2.7 Convergence of the semidiscretizaion in space

3.2.8 Convergence analysis of the full discretization scheme (3.91)

3.2.9 Convergence analysis of the full discretization nonuniform scheme (3.102)

3.2.10 Numerical examples

3.3 Theoretical and numerical aspects for the long time behaviour of nonlinear delay time Caputo fractional reaction-diffusion equations

3.3.1 Energy estimate for delay reaction-diffusion equations

3.3.2 Stability and contractivity of the solution

3.3.3 The numerical scheme

3.3.4 Numerical contractivity rates

3.3.5 Convergence analysis

3.3.6 Numerical simulations

3.4 An energy-preserving computational approach for the semilinear space fractional damped KleinGordon equation with a generalized scalar potential

3.4.1 Preliminaries

3.4.2 Crank-Nicolson implicit Galerkin-Legendre spectral methodology

3.4.3 Parallel implementation

3.4.4 Convergence and stability estimates

3.4.5 Discrete energy properties

3.4.6 Numerical simulations of the scheme

4 Numerical methods for variable order time-fractional reaction-diffusion equations and their applications

4.1 An L1 type difference / Galerkin spectral scheme for variable-order time-fractional nonlinear diffusion-reaction equations with fixed delay

4.1.1 The linearized numerical scheme

4.1.2 Theoretical analysis of the fully discrete scheme

4.1.3 Numerical verification

4.2 On the existence and uniqueness of solutions to a nonlinear variable order time-fractional reaction-diffusion equation with delay

4.2.1 Uniqueness of a solution

4.2.2 Local existence of a solution if s ^ T

4.2.3 Rothe-Galerkin spectral method

4.2.4 Numerical results

4.3 Space-dependent variable-order time-fractional wave equation: Existence and uniqueness of its weak solution

4.3.1 Uniqueness of a solution

4.3.2 A priori estimates

4.3.3 Existence of a solution

4.3.4 Existence of a weak solution for the wave equation with fractional derivative of constant

order

4.4 On the Rothe-Galerkin spectral discretization for a class of nonlinear variable-order fractional

wave equations

4.4.1 Weak formulation

4.4.2 Uniqueness of a solution

4.4.3 Time discretization

4.4.4 Existence of a solution

4.4.5 Numerical implementation

5 Time-dependent source reconstruction in a time-fractional diffusion equation with non-smooth solutions

5.1 A solely time-dependent source reconstruction in a one-term time-fractional diffusion equation with non-smooth solutions

5.1.1 Weak formulation

5.1.2 Existence of a unique solution

5.1.3 Rothe time discretization based on graded meshes

5.1.4 Numerical experiments

5.2 A solely time-dependent source reconstruction in a multiterm time-fractional order diffusion equation with non-smooth solutions

5.2.1 Uniqueness of a solution to the ISP

5.2.2 Reconstruction of the unknown source

5.2.3 Numerical experiments

5.2.4 Experiment 1: M =

5.2.5 Experiment 2: M =

6 Numerical analysis of time-fractional nonlinear subdiffusion equations with delay

6.1 A novel discrete Gronwall inequality in the analysis of difference schemes for time-fractional multi-delayed diffusion equations

6.1.1 Numerical model

6.1.2 Numerical analysis

6.1.3 Simulations

6.2 Convergence and stability estimates in difference setting for time-fractional parabolic equations with functional delay

6.2.1 Numerical model

6.2.2 Efficiency analysis

6.2.3 Numerical simulations

6.3 An easy to implement linearized numerical scheme for fractional reaction-diffusion equations with

a prehistorical nonlinear source function

6.3.1 Problem formulation

6.3.2 Preliminaries

6.3.3 The Numerical Scheme

6.3.4 Technical lemmas

6.3.5 Theoretical analysis of the semi discrete scheme

6.3.6 Theoretical analysis of the fully discrete scheme

6.3.7 Numerical simulations

6.4 Numerical analysis of multi-term time-fractional nonlinear subdiffusion equations with time delay:

What could possibly go wrong?

6.4.1 Applications to some numerical schemes in literature and improvements

6.4.2 First approach

6.4.3 Second approach

6.4.4 Improvement of other approaches

7 Grid schemes for integro differential problems with weakly singular kernels

7.1 A two-grid temporal second-order scheme for the two-dimensional nonlinear Volterra integro-differential equation with weakly singular kernel

7.1.1 Preliminaries

7.1.2 Analysis of the two-grid difference scheme

7.1.3 Numerical experiment

7.2 Optimal error analysis of space-time second-order difference scheme for semi-linear non-local Sobolev-type equations with weakly singular kernel

7.2.1 Time discretization

7.2.2 Truncation error

7.2.3 Stability and convergence analyses

7.2.4 Space disretization

7.2.5 Derivation of the fully-discrete scheme

7.2.6 Stability and convergence analyses

7.2.7 Numerical experiments

7.3 Theta-type convolution quadrature OSC method for nonlocal evolution equations arising in heat conduction with memory

7.3.1 The Q-type convolution quadrature rule

7.3.2 Case of p(t) with a >

7.3.3 Fully discrete scheme

7.3.4 Stability and convergence

7.3.5 Numerical experiments

7.4 Spatial two-grid (STG) compact difference method for nonlinear Volterra integro-differential equation with Abel kernel

7.4.1 Preliminaries

7.4.2 STG compact difference scheme

7.4.3 Analysis of STG compact difference method

7.4.4 Extension to two-dimensional problem

7.4.5 Two-dimensional STG compact difference algorithm

7.4.6 Stability and convergence analysis

7.4.7 Numerical experiments

7.5 BDF2 ADI orthogonal spline collocation method for the fractional integro-differential equations

of parabolic type in three dimensions

7.5.1 Preliminary knowledge

7.5.2 The BDF2 ADI OSC scheme

7.5.3 Theoretical analysis

7.5.4 Numerical experiments

8 Software packages for the solutions of models described by fractional partial differential equations

8.0.1 Software packages description

8.0.2 Numerical implementations of the software package

8.0.3 Structure of the software packages

9 Conclusions and outlook 311 References

The author's publications on the topic of dissertation

The author's certificates of state registration for computer programs on the topic of dissertation

Введение диссертации (часть автореферата) на тему «Development and Evaluation of Numerical Solvers for Mathematical Models with Fractional Derivatives (Разработка и оценивание численных решателей для математических моделей с дробными производными)»

Introduction

Relevance of the Research Topic. The research topic addressed in this dissertation-the numerical analysis and computation of nonlocal problems governed by fractional and integro-differential equations-is of profound and increasing relevance in modern science, engineering, and applied mathematics. At the heart of this relevance is the fact that many real-world phenomena are inherently nonlocal, meaning that their current state is influenced by a continuum of past states or spatial interactions that cannot be effectively captured by classical, local models based on integer-order derivatives.

Fractional calculus, which generalizes classical differentiation and integration to non-integer orders, provides a mathematically rigorous and physically meaningful framework for modeling systems with memory, hereditary properties, and spatial nonlocaliiy. These features are ubiquitous in complex systems: for example, anomalous diffusion in disordered media, viscoelastic and thermoelastic materials, subsurface fluid transport, neural dynamics, financial markets with memory effects, and various processes in quantum mechanics and biological transport.

In recent years, fractional and variable-order differential equations have emerged as powerful tools for describing such systems, as they allow for the inclusion of long-range temporal and spatial correlations. However, the analytical intractability of these models, due to their inherent nonlocality and nonlinearity, presents significant computational challenges. Closed-form solutions are rare or nonexistent, especially for nonlinear and multidimensional cases. As such, the need for robust, accurate, and efficient numerical methods has become critically important-not only for advancing theoretical understanding but also for enabling predictive simulations and engineering applications.

This dissertation responds directly to this need by developing a suite of innovative numerical schemes tailored specifically to fractional and nonlocal models. These include structure-preserving methods that conserve physical invariants like energy and mass, spectral Galerkin methods for high accuracy in spatial domains, and difference schemes for variable-order and delayed systems that better reflect real-life dynamics. Importantly, the work also addresses models with nonsmooth solutions, weakly singular kernels, and multi-scale behavior, which are common in practice but often neglected in numerical literature due to their complexity.

Furthermore, the topic is especially timely given the recent surge in interest in fractional modeling in cutting-edge fields, such as quantum gravity, soft matter physics, biomedical engineering, and control theory. The development of variable-order fractional operators reflects the recognition that system memory and nonlocality may not remain constant over time or space-introducing another layer of realism that this dissertation tackles through innovative numerical treatment and rigorous analysis.

In addition to its strong theoretical grounding, the work has high practical value. By implementing the proposed methods into modular software tools and validating them through a series of carefully designed numerical experiments, the dissertation ensures that its results are not just mathematically elegant, but computationally feasible and ready for application in diverse scientific and industrial contexts.

In sum, the relevance of the research lies in its ability to bridge a critical gap between real-world complexity and computational capability, providing the mathematical foundation, numerical tools, and practical insights necessary to advance the modeling and simulation of systems where nonlocality and memory are not exceptions, but

fundamental characteristics.

The Degree of Development of the Research Topic. The research domain concerning nonlocal problems-especially those governed by fractional and integro-differential equations-has evolved considerably, driven by the growing realization that many natural and engineered systems exhibit memory, hereditary characteristics, and long-range interactions. These features are fundamentally incompatible with the assumptions underlying classical differential equations, making fractional calculus an essential mathematical tool for modeling such behavior.

From a theoretical standpoint, the development of fractional calculus is well-established. Foundational contributions by Riemann, Liouville, Caputo, Riesz, and others have provided rigorous definitions of fractional integrals and derivatives, and their applications have expanded into disciplines as diverse as viscoelasticity, anomalous transport, diffusion in porous media, control theory, and quantum physics. Over time, the theory has matured to accommodate not only constant-order derivatives but also variable-order and distributed-order operators, which more realistically model systems whose dynamics evolve with time, space, or internal state.

Despite this theoretical progress, the field remains in a dynamic phase of development, especially regarding the computational treatment of these equations. Nonlocality introduces a significant departure from classical models: the system's current state depends on an entire history of prior states or on spatial domains far from the point of interest. This complexity poses serious challenges for the development of accurate and efficient numerical methods. Standard tools from classical numerical analysis-such as finite difference or finite element methods-must be adapted or completely re-engineered to deal with the nonlocal memory integrals and singularity of fractional kernels.

Over the last two decades, researchers have proposed a wide variety of numerical schemes, including L1-type finite difference methods, spectral Galerkin methods, and convolution quadrature rules. While effective for certain classes of problems-especially linear and smooth cases-many of these methods struggle with more realistic scenarios involving nonsmooth initial data, nonlinearity, delay, and variable-order derivatives. Moreover, the ability of numerical schemes to preserve intrinsic physical or structural properties of the original models-such as conservation laws, positivity, or dissipativity-remains limited in scope.

In this context, the dissertation makes a timely and original contribution. It addresses important open problems by developing structure-preserving, energy-conserving, and contractive-preserving numerical schemes for both time- and space-fractional partial differential equations. These include extensions to variable-order models and systems with delays-areas that are significantly underdeveloped in the existing literature. Furthermore, the dissertation provides comprehensive theoretical analyses of the proposed methods, establishing results related to convergence, stability, long-time behavior, and error estimates, even under challenging conditions such as non-smooth solutions or weakly singular kernels. In details, we go across the development of each element of our research topic.

The design of numerical methods that preserve the discrete energy of conservative systems governed by partial differential equations (PDEs) has become a central theme in computational physics and numerical analysis [78, 430]. Various approaches have been developed to achieve numerical schemes that maintain the energy characteristics of the original continuous models. These include classical finite difference methods [497], mimetic finite differences [245], finite element techniques [145], and Galerkin-based formulations [108]. This trend gained substantial momentum with the influential contributions of D. Furihata and collaborators around the early 2000s [134], laying the groundwork for structure-preserving algorithms. Subsequent advancements have extended these ideas to preserve other key physical invariants such as mass [185] and momentum [417], with broad applications across many physical systems. In recent years, the field has been further enriched by the incorporation of fractional calculus, which provides a versatile framework for modeling memory effects, hereditary properties, and nonlocal interactions. Fractional derivatives have found widespread application in modeling anomalous transport and diffusion processes [373, 155], long-range interacting particle systems [218, 421], Hamiltonian oscillator chains [422], electron-Induced Charging of Ferroelectrics [300], the theory of differential two-player games [338] and time-Fractional Dual-Phase-Lag Heat Conduction Model [270].

Multidimensional structure-preserving methods have also emerged as a significant research direction [39, 105, 416, 360]. The term "structure preservation" broadly refers to the ability of a numerical scheme to retain important qualitative features of the continuous system, such as positivity [292, 90, 271, 328], boundedness [354, 498], and monotonicity [278, 273], which are essential in fields ranging from population dynamics to thermodynamics.

The study of fractional PDEs has introduced new complexities and opportunities. These equations have become increasingly relevant in diverse fields such as viscoelasticity [204], thermoelasticity [349], quantitative finance [382], quantum mechanics [311], and hydrology [404].

A fractional-differential approach to numerical simulation of electron-induced charging of ferroelectrics, employed a fractional-differential model and an implicit finite-difference scheme to simulate the charging of ferroelectric materials under electron beam irradiation [299]. This approach allowed for improved agreement with experimental data by incorporating memory effects via the fractional derivative, which could be adjusted to capture different dynamic behaviors. Iln [282], the authors proposed a time-fractional modification of the Landau-Ginzburg-Devonshire-Khalatnikov model to describe the dynamics of ferroelectric polarization switching. To solve the resulting time-fractional cubic-quintic partial differential equation numerically, they developed a computational scheme combining an iterative procedure with an implicit finite-difference scheme based on an approximation of the Caputo derivative.

Despite their usefulness, fractional models often lack a well-defined Euler-Lagrange formulation [7], which complicates the design of energy-preserving schemes. Nevertheless, classical models have been generalized using fractional derivatives [381, 314, 424], allowing global and historical effects to be incorporated into local dynamic systems. Some approaches were trying to propose a method for constructing approximate solutions of ordinary fractional differential equations with Riemann-Liouville derivatives in indirect way [269]. The method was based on the two-scale technique, treating the fractional part of the derivative order as a small parameter and introducing two different scales. As a result, the fractional differential equation was reduced to a sequence of integer-order differential equations, all of which were linear, except possibly the first one.

Notably, Riesz space-fractional derivatives have been instrumental in modeling anomalous diffusion [156] and formulating conservation laws and Hamiltonian structures [425]. Also the authors in [37] investigated a fractional Helmholtz equation involving the fractional Laplacian. They constructed fundamental solutions and their factorized representations in terms of H-functions using Fourier and Mellin integral transforms. A multipole expansion for the integral representation of the equation's solution was derived, and a technique for evaluating H-functions from this expansion was proposed. Additionally, they developed a modified multipole method for solving the equation. This motivates the exploration in Sections 2.1 and 2.2 of whether numerical discretizations for nonlinear hyperbolic equations involving Riesz derivatives can be designed to conserve known physical quantities. The Zakharov system, originally formulated to model Langmuir wave propagation in plasmas [499], remains one of the most important coupled systems in nonlinear wave theory. It models the interaction between high-frequency electric fields and low-frequency ion-density fluctuations. Applications extend beyond plasma physics to include shallow water waves [54] and nonlinear optics [92]. The mathematical analysis of Zakharov-type systems has included global existence [432, 496, 407, 495], uniqueness [216], regularity [20], and exact solutions [435, 429]. Well-posedness has also been studied in various functional settings [36, 197, 129, 117].

With the rise of fractional calculus [238, 200], Zakharov-type equations have been extended into the fractional domain, giving rise to fractional versions of the Zakharov-Kuznetsov, Benjamin-Bona-Mahony, and Klein-Gordon-Zakharov systems [372, 32, 362]. Many of these models preserve fractional energy-like invariants, prompting the need for numerical methods that can respect such conservation properties. Sections 2.3 and 2.4 address this need by presenting both implicit and explicit schemes for Riesz space-fractional Klein-Gordon-Zakharov systems.

Despite the success of conservative schemes for integer-order PDEs, long-time behavior and energy dissipation in fractional systems remain relatively unexplored. Anomalous diffusion models, characterized by algebraic decay and long-tail dynamics, are central to understanding such behavior [173, 423]. Dissipative systems arise naturally in engineering and physics [168], requiring schemes that capture complex limit structures.

Historically, early efforts in conservative discretizations include the nonlinear Klein-Gordon equation [400], sine-Gordon equation [119], nonlinear Schr"odinger equation [420], and wave equations in polar coordinates

[333]. This legacy has been carried forward by researchers like Luis V'azquez and colleagues, leading to developments in conservative and dissipative Galerkin and finite difference methods [118, 284, 135, 283].

The Discrete Variational Derivative Method (DVDM) has emerged as a robust approach for constructing such structure-preserving schemes [138, 181, 484, 137, 136]. In parallel, the rapid expansion of fractional calculus [348] has spurred interest in conservative schemes for space-fractional PDEs using Riesz derivatives [34]. Riesz derivatives are particularly well-suited for such models due to their self-adjoint and negative-definite properties.

Recent developments have applied structure-preserving discretizations to space-fractional wave equations [275, 276], sine-Gordon equations [474], and quantum-corrected Zakharov models [471]. Caputo time-fractional derivatives, while less physically grounded than Riesz derivatives [422], are widely used in modeling memory effects in time. Numerical methods for Caputo derivatives include fast convolution algorithms [486], high-order central difference approximations [494], modified integral schemes [71], and predictor-corrector techniques for delay systems [44, 289].

The long-time behavior of Caputo-type systems has recently been analyzed, with studies showing the presence of absorbing sets and asymptotic contractivity [445,447]. The authors in [214] considered a homogeneous Dirich-let initial-boundary value problem for a quasilinear parabolic equation with a time-fractional derivative and elliptic coefficients depending on the gradient of the solution. Conditions on the coefficients ensured the monotonicity and Lipschitz continuity of the elliptic operator on the set of functions with uniformly bounded spatial gradients. For this problem, they constructed and analyzed a linear regularized mesh scheme, and derived a sufficient condition on the regularization parameter to guarantee the local correctness of the scheme. Inspired by these findings, Section 2.5 of this dissertation presents a theoretical and numerical investigation of dissipative mechanisms in nonlinear fractional PDEs with reaction terms.

Before proceeding with the description of Galerkin-Legendre spectral schemes for nonlinear time-space fractional partial differential equations with smooth and nonsmooth solutions, as given in section 3.1, let us give an overview of those advocated in the literature. The Lagrange-Galerkin scheme is studied for degenerate parabolic variational inequality arising in connection with the pricing of American options. This scheme is constructed using a combination of characteristic method for approximating the material derivative and the finite element method for approximating the diffusion part of the equation [96]. The design of an efficient technique to treat numerically the fractional differential operators is not an easy task. The main difficulty is given by the nonlocality of the fractional derivative. This salient feature imposes a serious challenge for the numerical simulation of nonlinear space-time fractional diffusion equations. There are two predominant categories of numerical methods for time stepping, i.e., the interpolation and the fractional linear multistep methods. The former approximates the fractional derivative directly by piecewise polynomials, whereas the latter relies on approximating the fractional derivative in the Laplace domain. For example, the piecewise linear interpolation yields the widely applied L1 method or Diethelm's finite difference method [244, 215]. High-order interpolations can also be applied; see [80]. The fractional linear multistep methods [265] provide another general framework for constructing high-order methods to discretize the fractional integral and derivative operators. These two approaches have their pros and cons: The fractional linear multistep methods inherit the stability properties of linear multistep methods for the initial value problem, which are quite flexible and often much easier to analyze, in a way often strikingly opposed to standard quadrature formulas [265], but it is often restricted to uniform grids. The interpolation methods are of finite difference type and very flexible in construction and implementation and can generalize to nonuniform grids, but often challenging to analyze, especially for nonsmooth data. Generally, these schemes are only first-order accurate when implemented straightforwardly, unless restrictive compatibility conditions are fulfilled. Hence, suitable corrections to the straightforward implementation are needed in order to restore the desired high-order convergence. We refer interested readers to [190] for an updated overview of these numerical methods and relevant convergence rate results.

The convergence analysis of the aforementioned algorithms is generally carried out on the condition that the underlying solution is sufficiently smooth [102]. However, too much regularity might impose unrealistic assumptions on the initial value [402]. When solving fractional differential equations, the singularity requires special

attention to enhance the accuracy of numerical schemes. Recently, many numerical approaches have been developed to deal with the weak singularity, such as using nonuniform/refined grids to capture the initial information and recover an optimal convergence order [241, 403], or using the correction terms to remedy the loss of accuracy and recover high-order schemes [508, 66], or employing non-polynomial (or singular) basis functions to include the correct singularity index [377, 501], or considering fractional coordinate transformations to regularize the problem [503]. Furthermore, fast and parallel algorithms were proposed to overcome the computational cost and speed up the approximation schemes in temporal direction for time fractional diffusion equation [388]. Up to now, most references were concerned with numerical analysis of the numerical scheme for linear fractional problems, e.g, [102,514]. A few researchers considered stability and convergence of the Ll-scheme for nonlinear time-dependent fractional differential equations [361]. However, the results were just held locally in time. In order to overcome the difficulties, fractional Gronwall type inequalities have been introduced [242]. It is worth mentioning that due to these discrete fractional Gronwall-type inequalities, numerical analysis of the L1-type schemes can be obtained without the local assumptions.

Another approach to design an efficient numerical scheme is to discretize the nonlocal differential operators with nonlocal numerical methods. Hence, the global behavior of the solution can be naturally taken into account and the computational cost is not substantially increased when moving from a second-order to a fractional-order diffusion model. Compared to the extensive amount of work put into developing finite difference schemes in the literature, only a little effort has been put into developing global and high order spectral methods. Historically, Galerkin methods of spectral type have been considered efficient schemes to solve partial differential equations. In fact, under suitable conditions, spectral Galerkin approximations usually gain an exponential convergence rate compared to the finite-difference approach. Moreover, one of the limitations of finite differences is that it is hard to raise their computational accuracy in space, especially when the properties of stability and convergence analysis rely on local assumption [541]. To overcome this difficulty, several spectral schemes have been developed. Also, suitable fractional-derivative spaces comparable to fractional-dimensional Sobolev spaces have been proposed [110]. In the way, a variational solution of fractional advection-dispersion equations on bounded domains was proposed in that work. Lin and Xu [244] developed a hybrid scheme for linear time-fractional diffusion problem, treating the time-fractional derivative using L1-scheme on uniform mesh and discretizing the integer-order spatial derivative by a Legendre spectral method. Zhang et al. [514] developed nonuniform L1-type spectral Galerkin approximation for linear time-space fractional equations. However, the theoretical understanding of these formulations remains rather limited in the context of nonlinear time-space fractional differential equations with nonsmooth solutions, and there is a huge demand on developing rigorous theoretical analysis for the nonlinear case.

A fractional generalization of quantum mechanics comes out when the Brownian-like quantum paths substitute with the Levy-like ones in the Feynman path integral. This concept was discovered by Laskin [219]. A space-fractional quantum framework has been introduced into quantum mechanics by extending the Feynman path integral to Lévy one [217]. A time-fractional quantum framework has been introduced into quantum mechanics [309] in which a new version of the space-time fractional Schrödinger equation has been launched. The introduced space-time fractional Schrödinger equation has a new scale parameter, which is a time-fractional generalization of Planck's constant in quantum physics. Recently, many efforts have been devoted to developing effective numerical methods for solving time and/or space fractional differential equations equations [490, 527, 121, 232, 297, 23, 297]. As an efficient tool for numerically solving nonlinear time and/or space fractional differential equations, Galerkin's method of finite element formulation has been introduced in literature. Li et al. [227] proposed a linearized L1-type numerical schemes for the nonlinear time-fractional Schrödinger equations and derived unconditionally optimal error estimates by using a novel discrete fractional Gronwall-type inequality. A nodal discontinuous Galerkin method for solving the nonlinear fractional Schrödinger equation and the strongly coupled nonlinear fractional Schrödinger equations has been proposed in [5]. An efficient Galerkin finite element method for solving the two-dimensional nonlinear time-space fractional Schrödinger equation with time-dependent potential function on an irregular convex domain was analyzed in [116].

Global methods such as spectral methods for fractional partial differential equations attract high attention due to the non locality for this kind of equations. Under some suitable conditions, the spectral Galerkin approximation

usually gains an exponential convergence rate comparing to the finite difference methods. Also, it is hard to raise the computational accuracy in space of a finite difference scheme in which the properties of stability and convergence results of high order schemes rely on some local assumption [541]. To overcome this difficulty, several researches have proposed some spectral schemes [45, 510]. An impressive attempt was made to improve the spectral Galerkin methods for fractional partial differential equations on bounded domains. In [110], fractional derivative function spaces were defined and demonstrated to be equivalent to the fractional Sobolev space, excluding integer multiples of half. Following [110], the authors in [162,244] proposed finite difference/spectral approximations for the time-fractional diffusion equation, two-dimensional Riesz-space fractional nonlinear reaction-diffusion equation and the distributed-order time-space fractional reaction-diffusion equation.

A combination between Galerkin-Legendre spectral method in space and the Chebyshev collocation method in time for solving the fractional Cattaneo equation was introduced in [231]. Recently, Bhrawy et al. [47,48,46] considered the multi-dimensional space and/or time constant-, variable, and distributed-order fractional Schrodinger equations using Jacobi spectral methods. Wang and Li [462] proposed a conservative spectral Galerkin method based on the Crank-Nicolson method for the temporal discretization and the Legendre-Galerkin method for the spatial discretization to solve the coupled nonlinear space-fractional Schrodinger equations. Moreover, a rigorous analysis of the unique solvability and optimal error estimate in the L2-norm were derived. A split-step spectral Galerkin method for the two-dimensional nonlinear space-fractional Schrodinger equation was presented in [463]. An approach with time-splitting spectral method for the coupled Schrodinger-Boussinesq equations with Riesz Fractional Derivative is proposed in [371]. To the best of our knowledge, there is no research content combining difference approximations and spectral methods in space dealing with time-space fractional CNLSEs. This gap in the literature is filled in section 3.2, in addition to the ignorance of most finite difference schemes to the fact of the nonsmooth solutions of fractional differential equations and the high demand of special classes of discrete fractional Gronwall inequalities to study the analysis of their discrete forms, motivates us to consider this case of study in spite of the known considerable challenges such as non-locality in time and space of the problem under consideration.

Differential equations with time delay du = f (t, u(x, t), u(x, t - s)), have been used a lot in the modelling of many observations in different sorts of applied sciences as biology, ecology, control systems medicine, etc. The parameter s > 0 denotes the delay parameter. They are concerned with providing satisfactory meanings of modelling phenomena in which its rate of variation does not rely only on the present affairs but also invokes their history. To give a wide description for delay differential equations and their fields of applications, we refer to the books [38, 58], the review papers [59, 28, 6] and the references therein. To give satisfactory modelling of real-world problems, the spatial operator should be invoked side by side with the temporal one. This yields the delay reaction-diffusion equations of the following form ^ = D1 ^ + f (t, u(x, t), u(x, t - s)), where D1 > 0 is the diffusion coefficient. The following diffusive Mackey-Glass equation (see e.g. [457,470, 228]) is an example of a reaction-diffusion equation with delay: -¡f = D1 ^ - a1u(x, t) + ^a^xf-s). This equation is used to model the hematopoiesis, where the spatial diffusion is used to express the spatial movement of substances effectively from high to low concentration; a1 , a2 and n are known parameters. Control problems related to a class of time fractional reaction-advection-diffusion systems have been discussed in [82, 83, 81]. A priori estimates to solutions of the time-fractional reaction diffusion equation coupled with the Darcy system was discussed in [542]. Another example of reaction-diffusion equations with delay is the diffusive Nicholson's blowflies equation (see e.g. [42]) du = D1 dx! - a1 u(x, t) + a2u(x, t)e-u"(x,i-s), which is a laboratory fly population model; a1, a2 and n are given parameters. Traveling waves describe transition processes from the physical point of view. However, if the reaction goes faster, the propagation speed will become rather large. This behavior is not clarified physically but can be introduced mathematically, cf. [53]. The analogous nonphysical phenomena were appeared in the heat conduction phenomena by Vernotte [52]. Memory effects have been considered in mathematical models to avoid the nonphys-ical behavior. Accordingly, the following non-Fickian delay reaction-diffusion equation has been studied (see e.g. [53, 21]): df = D0 + D2 /¿e-^§(x,s)ds + f (t,u(x,t),u(x,t -s)), where D2 > 0 and 8 > 0. The long time behavior of non-Fickian delay reaction-diffusion equations was studied in [229]. Energy estimates, dissipativity, asymptotic stability, and contractivity of the problems were discussed. A numerical method was presented which

can preserve the stability and contractivity of the underlying systems. Alternatively, fractional-derivative models for non-Fickian transport in a single fracture and its extension were addressed in [236]. For the non-Fickian transport observed in the fracture-matrix system, which is due to commonly observed multiscale physical/chemical heterogeneity within a natural fracture, the mass exchange rate may not remain constant. It is rather a random variable changing with the size and properties of the stagnant zones. This results in a multi-rate mass transfer process and a heavy (i.e. power-law) late-time tail in tracer breakthrough curves. This common sub-diffusion process can be efficiently captured by time-nonlocal transport models, especially the fractional differential equation with Caputo-type temporal fractional derivative with order 0 < a < 1 [521]. The movement of organisms and cells can be governed by occasional long-distance runs in accordance with an approximate Levy walk. Hence, the fractional diffusion can be also tackled for cell movement models with delay [111]. Accordingly, section 3.3 is concerned with the long time behavior of the time-fractional diffusion model.

The intensive work on Quantum gravity enhanced the possibility of the fractal structure of spacetime. In that direction, the causal dynamical triangulations theory which proposed in [17, 18, 16] can be considered as a non-perturbative approach to quantum gravity. It gives evidence about the fractal structure of spacetime such that the number of four dimensions at large scales goes to two dimensions at the Planck scale smoothly. Another theory of quantum gravity that has the same feature is named the asymptotically safe quantum Einstein gravity [220, 315]. It proves the evolving of spacetime from a lower-dimensional scale-invariant structure at high energy or micro scales to that of normal four-dimensional geometry at large scales or low energy. In these two theories, the fractal structure of spacetime and its geometry cannot be explained correctly by a single metric and so introducing a different metric at each energy scale is needed. A similar technique proposed by Modesto [295] used in asymptotically safe quantum gravity to be applied to loop quantum gravity. This technique can also give a similar fractal-like spacetime structure. In the fractal spacetime, a different suggestion of a Lorentz invariant theory needs fixing the dimensions of the spectral and ultraviolet Hausdorff to be two [62]. In the same manner, we can find out a technique to solve the ultraviolet problem by reducing the spacetime dimensions. The common feature of all these different techniques is the dependence of spacetime dimensions on scale and the possibility of having fractal properties at small scales. Fractional calculus can be considered as a useful mathematical tool to deal with fractals and fractal-related phenomena [380, 542, 19, 88, 515]. The realization of the close connection between fractional calculus and fractal geometries is an essential advance in the stimulation of fractional calculus application to physical science [370, 347]. This relation gives the fractional-order differential equations the possibility to be a natural tool for describing transport processes with fractal properties, such as anomalous diffusion, non-Debye relaxation processes, and other fractal phenomena. For a long time, the trajectories of the well-known fractal process, the Brownian motion, have been invoked in Feynman path integral approach to quantum mechanics [125]. The paths of the Brownian motion are known to be fractal curves with Hausdorff dimension two [4]. Accordingly, the employment of fractional calculus and fractional differential equations in Quantum theory is reasoned by the possible existence of fractal structure of quantum spacetime [290]. Anomalous transport systems in fractal media can be described by replacing the integer-order diffusion equation with the fractional order diffusion equation. In the same manner, it is reasonable to consider quantum mechanics and quantum field theory which satisfy fractional generalizations of Schrödinger equation, Klein-Gordon equation, and Dirac equation for quantum theories. The formulation of quantum mechanics through various types of fractional Schrödinger equations and their applications is discussed in [219]. In particular, the fractional Levy path integral approach of fractional quantum mechanics has been considered in [63,433]. It is interesting to see whether it is possible to use the path integral representation of fractional Brownian motion [385] and fractional oscillator processes [107] in the formulation of fractional quantum theory. Fractional Levy path integral approach of fractional quantum mechanics has been discussed in [63, 433]. The possibility of using the path integral representation of fractional Brownian motion [385] and fractional oscillator processes [107] is considered in the fractional quantum theory formulation. According to [317, 318], it has been noticed that the foamy or fractal nature of spacetime is observed at very short distances. So, it is inevitable to incorporate non-locality features in Klein-Gordon equations instead of locality in the fractal spacetime. To investigate how the curved spacetime affect the stability of propagating waves, we study the semilinear space fractional field equation of the Higgs boson type and providing the numerical analyses of the solutions and their discrete energy,

in section 3.4. The motivation of this study yields from the lack of computational structure-preserving approach based on spectral Galerkin scheme to solve Riesz space fractional Higgs boson equations in the de Sitter spacetime.

Fractional theory of calculus of variable-order has recently attracted considerable attention as a powerful mathematical tool for modelling a wide range of discontinuities and nonlinear phenomena [411]. Because many physical systems involve dynamics with memory effects whose behavior changes over time, even while transiting from one fractional-order to another, interest in fractional operators moved increasingly to their variable-order counterparts. The powerful practical implications of these variable-order objects, of course, come at the price of a more complicated mathematical characterization. In contrast to integer-order and fixed fractional-order operators, the variable-order of the fractional operators can be considered as a function of internal or external system variables such as, for example, time, space, state of stress, system energy, temperature, or even a combination of the various variables. Such an extension from fixed-order to variable-order operators also allows the formulation of mathematical models in which the order of the underlying governing equations can be modified using either the system's instantaneous state or its history. As a result, the corresponding model can evolve seamlessly to capture widely dissimilar dynamics without changing the structure of the governing equations which describe the system's response [60]. This emphasizes the evolutionary nature of the variable-order fractional calculus formalism, which indeed can play a critical role in the simulation of nonlinear dynamical models. Recognizing this untapped potential and unique capability of variable-order operators, the scientific community has been intensively investigating applications of variable-order fractional calculus to the modelling of physical and engineering systems [140, 195, 139,25, 141].

It is worth noting that viscoelasticity is certainly the playground of the most widespread applications of variable-order fractional calculus since its appearance [359, 383, 147]. Although some examples of fractional operators with constant fractional-order successfully mimic experimental test findings, these operators may not always be suited for characterizing viscoelastic behavior. Indeed, nonlinear and complex phenomena occur in several materials, such as biological tissues, polymers, and rubbers. These phenomena imply significant variations in the material's mechanical characteristics that cannot be described by fractional viscoelastic models of constant-orders [60]. Similar problems arise to characterize the mechanical behavior of viscoelastic materials under variable environmental conditions, for instance, changes in temperature or viscoelastic aging effects. As a result, variableorder fractional operators have recently been proposed and utilized in a wide range of applications. For more comprehensive literature review on variable-order fractional operators, associated differential equations and latest applications in natural sciences, we refer the reader to Patnaik et al. [330], Samko [378], Sun et al. [411] and Ortigueira et al. [321].

Time-fractional diffusion equations serve as effective tools for modelling and simulation of various processes with hereditary or memory features, such as anomalous diffusion transport, and hence attract considerable interest [290, 267, 41]. Moreover, the theory, application, and implementation of several numerical approaches for the solution of time-fractional differential equations have been proposed in [1, 2, 543, 45, 46]. However, the singularity of the solutions of constant-order time-fractional diffusion equations at the beginning of time seems to be physically irrelevant to the subdiffusive transport of the model. The fundamental reason why this phenomenon occurs lies in the incompatibility between the global nature of the power law decaying tails and the locality of the classical initial condition at the beginning of time [26, 187]. It was pointed out in [450] that the regularity of the solution to the variable-order Riemann-Liouville fractional diffusion equation depends on the behavior of the variable order and its derivatives at time t = 0, in addition to the standard smoothness assumptions. It was demonstrated that the solution to the problem exhibits full regularity like its integer-order analogue if the variable order has an integer limit at t = 0, or it has a singularity at t = 0 like in the case of the constant-order time-fractional diffusion equations if the variable order has a non-integer value at time t = 0. Accordingly, the variable-order fractional operators can eliminate the nonphysical singularity of the solutions to constant-order time-fractional diffusion equations and produce solutions with full regularity. To the best of our knowledge, up to now, there are no closed form solutions to variable-order fractional diffusion equations. The main reason is that, due to the impact of variable fractional order, analytic techniques such as the Laplace transform cannot be utilized to solve variable-order fractional mod-

els. Recently, there are several investigations on numerical schemes to approximate the solutions of variable-order time-fractional diffusion equations in the literature. We mention here the works by Tavares et al. [426], Chen et al. [87], Karniadakis et al. [506, 504, 523], Zheng et al. [534, 529, 186], Pang and Sun [325], Wei et al. [466, 467] and Liu et al. [248]. Nevertheless, rigorous numerical and mathematical analysis of variable-order fractional differential equations remains wide-open.

Most studies in the literature have been absorbed in the linear variable-order fractional diffusion equations without delay [160, 187, 439]. In section 4.1, we develop a linearized explicit finite difference/spectral Galerkin approach for nonlinear variable fractional-order diffusion-reaction equations with a fixed time delay and a drift term, and discuss its stability and the convergence rate. The first time that a variable order fractional differential derivative appeared as a generalization of a constant order fractional differential operator was in an article of Samko and Ross in 1993 [379]. An employment of fractional partial differential equations (FPDEs) with variable orders have been encountered in different physical and dynamical systems [14]. A successful application of variable order differential operators to a wide range of real-world problems has been verified due to the ability of constructing governing equations of evolutionary type [329]. These applications have been efficiently used in fields of science in which investigating memory properties that vary in time and space is needed, as in mechanics, transport processes, control theory and biology [330, 411]. New difficulties from mathematical, numerical and computational viewpoint raised up when dealing with FPDEs. For example, an exhibition of singularity near the initial time is noticed for the first-order time derivative of the solutions to the time-fractional diffusion equations. This weakly singularity gives an evidence that the error estimates in the literature are inappropriate if their proofs were done under full regularity assumptions of the true solutions, which is discussed in e.g. [403]. According to the works [26, 450], the root of occurrence of solutions with non-physical sense for constant order FPDEs near t = 0 comes out from the incompatibility between the non-locality of the governing FPDEs (of the power-law decaying tails) and the locality of the initial conditions, respectively. Elimination of that non-physical singularity can be done as in [450] if we use FPDEs in which the order varies smoothly to an integer value at t = 0. It can be clearly noticed in the literature that variable-order FPDEs arise in many applications, as in [506, 411] in which the variable order indicates the fractal dimension of the porous media. Different numerical approximations were developed and analyzed for FPDEs, see [540, 389, 408, 525, 506, 468, 167]. Well-posedness results for time-dependent variable order problems where the derivative order has not the memory of its history can be found in e.g. [434, 280, 450, 532], whilst for a variable order FPDEs with hidden memory we refer to [533]. All these contributions consider linear sources and do not consider time delay.

Rothe's method (which is also called the method of semi-discretization in time) was initially developed as a discretization method in time for partial differential equations [194]. The work in [194] introduced Rothe's method as an accurate theoretical tool for solving a wide scale of evolution problems. A theoretical and numerical treatment of initial boundary value problems of parabolic type with Volterra operators was constructed in [193], endowed by an integro-differential equation and a natural boundary condition. A constructive proof of the uniqueness and existence of its variational solution under weak assumptions on the data was done by invoking Rothe's method. We also give some recent well-posedness results illustrating the importance of this method. In [440], the authors study the existence and uniqueness of a solution to anisotropic thermoelastic systems. A consideration of a weak solution for a fractional order diffusion equation with Volterra differential operator and with fractional integral condition is shown in [74]. It has been proved that the solution existed and was unique as well as its regularity by designing an appropriate Rothe scheme. Recently, the paper [288] studies analytically and numerically a Rothe-Galerkin finite element method for approximating the solutions of a Boussinesq-type system to model water wave propagation over a time-dependent variable topography. Moreover, in [439], the author studied the well-posedness of a fractional diffusion equation with space-dependent variable order with the aid of Rothe's method. Section 4.2 is devoted to study the existence and uniqueness of variable order time-fractional reaction-diffusion equation with delay by means of Rothe scheme.

The time-fractional diffusion or wave equation is the fractional-order generalization of the classical diffusion or wave equation, respectively. Some real-world applications for different kinds of fractional order wave equations can be found in [175, 24, 266]. For available applications of variable-order fractional diffusion and wave operators

in the general area of scientific and engineering modelling, mathematical foundations and numerical methods, please see [259, 411, 330]. Section 4.3 is concerned with variable-order fractional wave equation with space-dependent variable order. The solution to autonomous (time-independent elliptic part) constant-order fractional wave equations is proved to be existed and unique in [438, 436]. A priori estimates for solutions of boundary value problems for constant-order fractional wave equations are derived in [12]. A fundamental solution to the fractional wave equation of constant-order is determined in [266]. The existence and uniqueness of the solution to generalized fractional differential equations with variable-order operators have been discussed in [363, 479, 511]. Analytical solutions to specific fractional time-dependent variable-order differential equations are obtained in [280]. Interesting studies about the numerical analysis of the variable and distributed order fractional wave equation can be found in [103, 340]. However, to the best of our knowledge, the contribution in section 4.3 deals with the well-posedness of problem (5.46),which is he first in literature by means of Rothe scheme.

One of the mean features of Caputo time fractional subdiffusion equations and their applications is the initial weak singularity [267]. Similar to its subdiffusion analogue, the time-fractional wave equation exhibits an initial weak singularity that does not seem physically relevant. Consequently, assuming that the solutions are smooth is deemed inappropriate in many results proved in literature. The exhibition of two time scales is the main reason behind the inappropriateness of that assumption. More clearly, the Brownian motion of the diffusive process is initially noticed in the heterogeneous material at the first timescale and then goes to anomalous diffusion by time evolving. The vibration of a perfectly elastic membrane can be considered as the motivation of model problem of variable order wave equation. This membrane is tautly stretched along the boundary d Q of the physical domain in the plane Q as the equilibrium position of the membrane, and construct a u axis perpendicular to the x plane [535]. For more explanations about the applications of variable-order fractional differential problems in mechanics, viscoelasticity, transport processes and control theory, we refer to [330]. A survey of the recent relevant literature and findings in primary definitions, models, numerical methods and their applications for variable-order differential equations can be found in [411].

A glimpse on variable-order fractional derivatives was given by Samko and Ross in 1993 [380] initially. Theoretical findings about the well-posedness and regularity of variable-order fractional wave problems are limited. This observation can be due to multiple reasons. The first is related to the impossibility of analysing through giving the analytical formulation in terms of special functions, which happens more easily in their constant-order counterparts. The next reason for not extending the technique in [531] to a class of nonlinear variable order fractional wave equations as will be given later in section 4.4, is the heavy relying on the exponential decaying of the solution operator of the variable-order diffusion equation. This can be helpful to balance the exponential growth of the solutions to the variable-order ordinary differential equations in the spectral decomposition in case of considering diffusion problems. The solution operator of that problems does not decay exponentially, so that technique cannot be invoked here. The false assumptions of sufficient smoothness of the solutions to them is another reason. More recently, a novel technique has been proposed to prove a greatly improved polynomial growth of the solutions to the variable-order wave ordinary differential equations in the spectral decomposition of a linear case of variable-order fractional wave problems with only time-dependent coefficients [535]. This approach is made with respect to the eigenvalues from the previously proven exponential growth. This technique can accordingly prove the well-posedness and regularity estimates of a particular case of variable-order fractional wave problems.

For the existence and uniqueness of the weak forms of variable-order fractional differential problems, we invoke recent studies from the literature. Concerning the space-dependent variable-order, we mention [201, 439]. In [201], the existence of a weak solution to the fractional diffusion equation with initial and boundary conditions has been studied. The governing elliptic operator is supposed to be autonomous. K. Van Bockstal researched a similar problem in [439]. However, the coefficients accompanying the problem have a temporal and spatial variable dependency. Armed with the strongly positive definiteness of the governing kernel and the assumption of the belonging of the initial data to H0(Q), the existence of a unique weak solution has been confirmed in [439]. Concerning the time-dependent variable-order, we state the contributions [532, 535]. In [450], the analogue of (5.46) has been considered for the linear diffusion equation but with a variable-order Riemann-Liouville fractional derivative (with j € C([0, T]) satisfying 0 < jSm ^ j(t) ^ 1 and limt^0(j6(t) - jS(0))ln(t) = 0) and solely space-

dependent coefficients. The authors establish the well-posedness of the problem in multiple space dimensions (using eigenfunction expansion on a smooth domain Q) and show that the regularity of the solution depends on the value of ^ (0). Then, in [532], the authors study the well-posedness of the diffusion problem with current-stated based variable-order operator (4.3). Also, in this situation, the regularity depends on the value of ^(0). In [533], the authors study problem (5.46) for k(x, t) = K > 0 and f = 0. Again using spectral decomposition, the authors prove the well-posedness of the problem if ^ € C1([0, T]) satisfies 0 < ^(t) < 1. A full discretisation of this problem has been studied in [535]. The new aspect of our contribution in section 4.4 is that we formulate the problem on a more general domain and consider a time-dependent diffusion coefficient and nonlinear source term. Consequently, the spectral decomposition approach is not applicable, and for this reason, we tackle the problem with the aid of Rothe's method. Rothe's method was initially developed as a discretization method in time for partial differential equations [365, 194]. The work [194] introduced Rothe's method as an accurate theoretical tool for solving a wide scale of evolution problems. Note that Rothe's methods has been also successfully applied when solving fractional evolution problems of constant order [438, 436].

Unlike to the classical diffusion equation, the subdiffusion model contains a time fractional-order derivative. This is supported by a macroscopic level observation, the diffusion process comes from the random motion of individual particles, and the use of the first-order time derivative in the canonical model rests on a Brownian motion assumption. It is shown experimentally that anomalous diffusion in which the mean-square variance grows slower (subdiffusion) than in a Gaussian process can give a satisfactory fitting to some experimental data [290]. In the direction of mechanically loaded structures, the authors in [209] addressed an inverse problem for determining the strength characteristics of a thin-walled structure subjected to combined mechanical and thermal loads. The problem was approximated using the finite element method with a super-element model naturally corresponding to the structure under study. The subdiffusion model as in [183] can be used to handle a setup in a microwave heating process used in various applications in industry, e.g. in ceramics and in food processing. A supply of external energy is targeted at a controlled level by the microwave generating equipment. However, a spatial and temporal variation for the dielectric constant of the target material is leading to spatially heterogeneous conversion of electromagnetic energy to heat. This can be due to a source term in the subdiffusion model. It is noticed that this spatial variation of absorbing material has no high effect on the thermal diffusivity, this is reasoned to another material at higher concentration.

There is a high demand to mention that the temperature is not so high that temperature dependence of the dielectric constant is important, as in thermal runaway studies [369]. If u(x, t) denotes the concentration of the absorbed energy in this model, then its integral over all volume of the material determines the time dependence absorbed energy. The above inverse problem mentioned here for such a model gives an idea of how total energy content might have an external control. The integral condition that can be accompanied by the subdiffusion model, arises naturally and can be used as supplementary information in the determination of the source term. Such type of condition can model various physical phenomena in the context of many fields [64]. Before starting with discussing the most relevant articles concerning that, we note some interesting studies about the recognition of the space-dependent part of the source for fractional equations, e.g. [203, 183, 492, 374, 487, 396]. Identification of the time-dependent part of the source term in a fractional diffusion equation assuming various boundary conditions and additional measurements has been studied in [469, 192, 394]. The well-posedness of the problem in a one-dimensional setting with nonlocal boundary conditions, F = 0 and a measurement in the form of integral over the domain is studied in [183]. In [203] the problem is addressed assuming the unit square domain and a nonlocal boundary condition. In [492] the authors studied the recognition of the space-dependent source in a 1D case from final data; a numerical scheme based on a local discontinuous Galerkin method is investigated. In [374] the problem is assumed for the fractional diffusion equation with only space-dependent coefficients, uniqueness is proven and a regularization method is applied. In [487] the authors proposed the Landweber regularization technique for the backward problem in a fractional diffusion equation. The Tikhonov regularization method is applied on the problem considering a simple fractional diffusion equation in the one-dimensional case in [446]. In [394] an ISP for a semilinear time-fractional diffusion equation with a solely time-dependent unknown source was

studied. The existence and uniqueness of a weak solution for the ISP were proved. The missing source term was recovered from an integral-type measurement over the domain. A numerical algorithm based on Rothe's method was established and the convergence of approximations towards the exact solution was demonstrated. However, the analysis in this paper gives that the solution is continuously differentiable on the closed time interval. An essential feature of the time Caputo fractional subdiffusion problem is that the solution lacks the smoothness near the initial time although it would be smooth away from t = 0. The related discussions are referred to [56, 267]. Stynes et al [403] showed that under proper regularity and compatibility assumptions, the one-dimensional subdiffusion problem dtSu - puxx + c(x)u = f (x, t) with homogeneous Dirichlet boundary condition has a unique classical solution, and there exists a constant Cw such that

|u''(t)| < Cu(1 + tS-2), 0 < t < T. (1.1)

In [206, Section 6.1], this result is generalized to the case Q C Rd for d € {2,3}. A rigorous numerical analysis of time fractional nonlinear parabolic partial differential equations, with fractional derivative of order a € (0,1) in time is presented in [191]. A complete solution theory is done for the time fractional nonlinear parabolic partial differential equations with a Lipschitz nonlinear source term. Nonetheless, as pointed out by [56, 304, 306], most of existing works on L1 or high order L1-type approximations of the Caputo fractional derivative assumed that the solution is smooth at the initial time t = 0. The corresponding error estimates are always restrictive since often the solution does not have the required regularity. The Rothe's method in [394] was based on a compatibility condition between the initial data, sources and boundary condition, and so the corresponding error estimates are also restrictive. More recently, the authors in [403] analyzed the L1-formula on graded time grids by taking into account the initial singularity in the Caputo's fractional subdiffusion problems. They obtained the error estimate under the regularity assumption (1.1) by using the discrete maximum principle and direct analysis of the local truncation error with the basic assumption Ti-1 ^ Ti, 2 ^ i ^ n, where n is the number of time discretization intervals. This assumption has a reasonableness because graded meshes that concentrate the grid points near t = 0 would be necessary to recompense the lack of smoothness there. With the aid of the discrete fractional Gronwall inequality and the global consistency error analysis. The work in section 5.1 focuses on the global (in time) solvability of the ISP (5.2)-(5.3) under low regularity assumptions and in the proposition of two approximation schemes based on different approximations of the Caputo fractional derivative.

Replacing the standard time derivative with a time fractional derivative produced a time-fractional order diffusion equation. The produced diffusion equation can be handled to explain accurately the superdiffusion and subdiffusion phenomena [291, 120]. Further theoretical investigations are required to incorporate adequate tools for the description of more complicated (and more realistic) random processes, which are described by a set of characteristic exponents, and are therefore of multi-term time-fractional diffusion equation, i.e., a finite linear combination of the fractional derivative with the positive weight coefficients. This kind of model could describe the diffusion phenomenon of a solute in a multi-scale medium. Such processes are believed to provide useful models for a crowd of non-homogeneous and non-stationary processes, for example, it is shown to be efficient models for describing some anomalous diffusion processes in the highly heterogeneous media by [171] in which the authors indicated that diffusion equation with time-fractional derivative was well- performed in describing the long-tailed profile of a particle diffuses in a highly heterogeneous medium.

In case of direct problems of time-fractional diffusion equation (initial value problems and initial boundary value problems), the literature is fruitful with many studies about them. For some theoretical results concerned with uniqueness and existence, we refer to [268, 344, 212]. Numerical implementations for that direct problems are tackled in the literature also by using different numerical methods. Finite difference schemes are invoked in [247, 305], Galerkin spectral methods are utilized in [418, 461] and finite element methods are implemented in [188, 207]. Practically, it can be happened that part of initial conditions, or boundary conditions, or source term may be unknown. Recovering them may need some additional measurements, which leads to fractional diffusion inverse problems. Different numerical techniques are formulated to reconstruct various parameters in time fractional diffusion problems [115, 391, 488]. For time dependent source function reconstruction for time Caputo fractional diffusion equations, which is the main focus of our work here, we mention the following contributions

to the two scales of theoretical analysis (existence and uniqueness) and the numerical identification of the source functions by means of numerical methods. On the scale of the numerical identification, we invoke some recent works [469, 154, 22, 254, 488]. A determination of a time-dependent source term by means of initial and boundary data and an additional measurement data at an inner point for a time-fractional diffusion equation [469]. By using separation of variables and Duhamel's principle, the problem transformed to Volterra integral equation of the first kind, then a combination between boundary element method and a generalized Tikhonov regularization is utilized to solve the resulted integral equation. In the same manner, the identification of the time-dependent source term according to an additional solute concentration distribution measured at a point on the boundary for a time-fractional diffusion-wave equation was done in [154]. A profitable meshless method based on radial basis functions is used for determining a time-dependent factor of an unknown source on some subboundary conditions for a time-fractional diffusion equation from the non-local measurement data in [22]. An inverse source problem of time-fractional diffusion wave equation on spherically symmetric domain was considered in [488]. Landwe-ber's iterative method is used to solve this inverse source problem and derive a priori and a posteriori convergence estimates. On the scale of theoretical analysis, Rothe's semi discretization scheme was firstly formulated for time fractional inverse source problems in [394]. The missing solely time-dependent source is recovered from an additional integral measurement. The existence, uniqueness and regularity of a weak solution is addressed. A numerical algorithm based on Rothe's method is designed, a priori estimates and convergence of iterates towards the exact solution are derived. Although the results are accurate and correctly derived, it was assumed that the solution of the problem is smooth. This naturally contradicts the behavior of the solution near t = 0 as a lack of smoothness occurs due to the singular kernel of time Caputo fractional derivative. For completeness, we also refer to some reference works for the reconstruction of time-dependent sources in classical heat conduction and other settings (from theoretical and numerical side), see [65, 350, 158, 437, 159]. Inspiring by this, we are concerned here with the multiterm time fractional inverse source diffusion problem with smooth or non-smooth solutions, in section 5.2.

The numerical methods for time-fractional differential equations are classified into two categories: indirect and direct methods. The indirect methods are based on the reformulation of the time-fractional differential equations into integro-differential equations. On the other hand, the direct methods are based on direct approximations for time-fractional derivatives [89]. Examples of these methods are the L1-type techniques [320] and the L2-1a schemes [13]. The L1 schemes are based on piecewise polynomial approximations for time-fractional derivatives. It is worth mentioning that most of the reports on L1-type methods study the efficiency of numerical schemes for linear fractional-order problems [539, 189]. A few of those works [226] have discussed the stability and the convergence for nonlinear time-fractional differential equations. However, we must point out that those results were reached locally in time. In order to improve those achievements, Li et al. [227] developed a discrete form of a fractional Gronwall-type inequality which can be used in the numerical analysis of L1-type schemes without local assumptions. However, that discrete fractional Gronwall-type inequality cannot be directly applied to time-fractional differential equations with time delay. Recently, Li et al. [233] established a new discretized fractional Gronwall-type inequality to analyze fractional reaction-diffusion equations with fixed temporal delay using nonlocal assumptions. Motivated by the fact that L1 schemes have only order of approximation O(A2—a), Gao et al. [144] developed a new difference analog of the Caputo fractional derivative L1—2. In turn, Alikhanov [13] considered ath order approximations for Caputo fractional derivatives at the time t = tj+CT with a = 1 — 2, and constructed a new difference analog. This new approach has been called the L2-1a scheme. Motivated by these achievements and some recent applications [509], we present, in section 6.1, a discrete fractional Gronwall inequality which is consistent with the L2-1a. Our Gronwall inequality has the advantage over previous works that it can be used to provide optimal error estimates for multi-delayed fractional problems like (3.71)-(3.74). Helped with such tool, we construct high-order compact linear difference schemes based on the L2-1a formula, and we will discuss in detail the numerical analysis of the proposed scheme. Both the analytical and the numerical results obtained here will show the superiority of the current approach when compared to similar methods available in the literature.

The L1-type schemes have a wide range of applicability in solving differential equations of fractional order

in time [319]. Most reports on L1-type methods refer to the efficiency analysis of numerical schemes for linear fractional-order problems [539, 189]. A few of those works have discussed the stability and the convergence of L1-type schemes for nonlinear time-fractional differential equations. That is the case for articles like [230]. However, the results were controlled by the locality in time in those cases. Recently, Li et a/. [237, 226] developed a fractional Gronwall-type inequality in order to overcome those difficulties. Using that inequality, the numerical analysis of L1-type schemes may be established without local assumptions. It is important to note that the fractional Gronwall-type inequality mentioned above cannot be directly applied to time-fractional differential equations with time delay. Recently, Li et a/. [226] developed a new fractional Gronwall-type inequality for fractional problems, in order to analyze fractional reaction-diffusion equations with fixed temporal delay. This novel form of Gronwall's inequality helps in obtaining an optimal error estimate for fixed time-delay fractional parabolic equations. The present work aims at extending the results in [227] to study time-fractional differential equations with functional delay. Such systems can effectively model physical problems for which the evolution does not only depend on the present state of the system but also on the past history. Those equations provide more realistic models for phenomena that display time-lags or memory-effects. As examples of those phenomena in the sciences, we may cite various problems in automatic control [180], traffic models [97] and population dynamics [343]. Gronwall-type inequalities are crucial in the qualitative analysis of fractional systems in differential and difference settings [537]. In the differential setting, the existence and uniqueness of positive solutions for a class of nonlinear fractional delay differential equations was considered in [239] using a nonlinear alternative of the Leray-Schauder type. The stability and the dissipation of Caputo nonlinear fractional functional differential equations with order 0 < a < 1 were discussed in [445], based on new Gronwall-type inequalities. Moreover, some attractiveness results for fractional functional differential equations were obtained in [77] using a fixed-point theorem. Finally, a robust finite-time stability problem of fractional-order systems with time-varying delay and nonlinear perturbation was investigated in [339] based on a generalized Gronwall inequality. A delay-dependent sufficient condition for robust finite-time stability of such systems was provided there in terms of the Mittag-Leffler function. On the other hand, optimal error estimates of some numerical schemes for multidimensional nonlinear time-fractional Schrödinger equation have been obtained using a discrete form of a fractional Gronwall-type inequality introduced in [227]. The convergence and the stability of the proposed scheme were obtained in terms of a new fractional Gronwall-type inequality. In section 6.2, we extend the role of the discrete form of the fractional Gronwall-type inequality proposed firstly in [227], in order to obtaining optimal error estimates for time fractional parabolic equations with functional delay. Moreover, we study the convergence and stability of the difference method using those results.

Investigating analytically and numerically fractional differential equations had been of wide concern and significant interest for decades. The reason behind that is its availability to model some phenomena more effectively than the integer-order equations. They have been used in many fields of engineering and science as physics, chemistry, biology, viscoelasticity, finance and many different areas of science and technology [410, 272, 3]. Similarly, delay differential equations owe an extendable range of science and engineering applications. The prehistory feature of these equations can be explained in the sense of how the change rate of the time-dependent process of their mathematical modeling is determined through past and present states in the meantime. A vital role in mathematical modeling is played and extensively used in different kinds of applications. Physiological systems [33], population dynamics [249], HIV-contamination modeling [95], chemical reaction modeling [43] are examples of these applications. Including fractional derivatives and delays allow modeling real-world problems more accurately. Fractional delay partial differential equations (FDPDEs) are equations involving fractional derivatives and time delays. Non-local in nature is a feature of fractional derivatives that can represent memory effects whereas the history of an earlier state feature is effectively formulated by time delays. The evolution of a dependent variable of FDPDEs impacts solving FDPDEs. The difficulty of solving FDPDEs comes from the dependence of that evolution at time t on their value at t — s where s is the time delay and all previous solutions due to the character of history dependence of a fractional derivative. In [260], Lu developed a mono iterative scheme for the solution of the finite-difference system derived from a class of nonlinear delay reaction-diffusion systems. Sun and Zhang in [413] used a linear compact difference scheme to solve the scalar delay parabolic equations. Liu and Zhang

in [246] constructed a fully discrete scheme based on combining the Crank-Nicolson method and the Legendre spectral Galerkin method for the two-dimensional nonlinear delay diffusion-reaction equations. Two classes of finite difference methods are constructed in [519] to solve the space fractional semilinear delay reaction-diffusion equations. The fractional-centered finite difference method in space and backward differential formula in time are employed. A second-order accurate implicit scheme based on the Alikhanov formula for the temporal variable and the fractional centered difference formula for spatial discretization is established in [526] to solve a class of timespace fractional diffusion equations with a time drift term and a non-linear delayed source function. A spectrally accurate Petrov-Galerkin spectral scheme for fractional delay differential equations is developed in [505]. This scheme is developed based on a new spectral theory for fractional Sturm-Liouville problems. Recent research has focused on developing alternative approaches to finding the approximate solution to fractional partial differential equations, as well as the appropriateness of the Legendre spectral technique for other nonlinear models [544, 538, 512, 101]. By using the semigroup theory of operators and the monotone iterative technique, the existence and uniqueness of mild solutions for the abstract time-space evolution equation with delay under some quasimonotone conditions are obtained [234].

Despite numerous numerical schemes employed for FDPDEs, convergence and stability analysis for direct numerical methods of FDPDEs is limited especially for nonlinear problems [169, 166]. This is because of the unsuitable application of traditional techniques of discrete energy inequalities and the improper invoking of the traditional versions of discrete Gronwall inequalities. The nonlocality feature of time Caputo fractional derivatives and their singular kernels induced the interpolated numerical approximation for such kinds of derivatives. This makes the use of usual discrete Gronwall inequalities inappropriate. New versions of discrete Gronwall inequalities are recently constructed to fill this gap. These kinds of discrete inequalities are called discrete fractional Gronwall inequalities. They were firstly constructed by Liao and his coauthors in [242] by introducing a new concept of discrete complementary convolution kernels corresponding to Riemann Liouville fractional integral kernels. Inspiring by these recent inequalities, convergence and stability estimates for the proposed scheme introduced here are investigated using discrete energy inequalities.The key to Liao analysis was to establish a new Gronwall type inequality for a positive sequence &n, n = 1, 2, ■ ■ ■ , satisfying 0D&n < + n2&n-1 + gn, where 0D denotes an L1 approximation to the time Caputo fractional derivative operator, n and n2 are both positive constants. In case of delay problems, the analysis is built on establishing a discrete Gronwall inequality for a sequence satisfying the

& > 0 V i > 0, O0 is known and & = 0 if i < 0,

< &j + gj, Vj < n,

0D?< Mj + n2&j-1 + n3&j-n-1 + n4&j-n + , Vj > n.

So, we are extending the use of such novel discrete inequalities to the combined numerical schemes such as the finite difference/Galerkin spectral method constructed in section 6.3. For continuous time-fractional differential equations with delay, the fractional Halanay inequalities [448] are important in the study of stability and dissi-pativity for that kind of problem. Inspired by these inequalities, we can state and prove the convergence and stability estimates for semi-discrete schemes proposed later. The main contribution of that piece of work is to extend these numerical and theoretical results for time delay problems. That extension is not trivial due to the new challenges that appeared in the formulation of a linearized scheme and also in the theoretical analysis (convergence and stability) in the case of semi and fully discrete schemes. The first challenge is solved by invoking the recent discrete inequalities of Halanay type and the second one by recalling a suitable form of discrete fractional Gronwall inequalities which can deal with prehistory functions.

There are two primary implementation forms for numerical methods considered by the experts to approximate time Caputo fractional differential equations: indirect and direct approaches. The indirect ones are built on the transformation of the time-fractional differential equations to integro-differential equations. A very interesting paper is [346] where an iterative technique has been applied for non-Lipschitzian time fractional porous medium equation. The approach based on a transformation of the governing equation to the equivalent nonlin-

ear Volterra equation. The classical convergence theory cannot be applied since the nonlinearity of the equation is non-Lipschitzian and so some specific quadratures used to discretize the corresponding Volterra integral equation. On the other hand, the direct schemes are proposed depending on direct approximations for time-fractional derivatives. The first example of these methods is the L1 technique [319] which based on piecewise polynomial approximations for time Caputo fractional derivatives. Motivated by the fact that L1 schemes have only order of approximation O(A2-a), a € (0,1), if the exact solution of the considered problem have the requisite regularity around the initial point and some lose of accuracy would be reached if the initial singularity is observed. Alikhanov [13] considered a high order approximation for Caputo fractional derivatives at the time t = tj+a with a = 1 - "2, and constructed a novel difference analog. It is worth mentioning that most of the reports on L1 and L2 - 1a methods study the efficiency of numerical schemes for linear fractional-order problems [539, 189]. Some efforts in [227, 242] are exerted to develop a discrete form of a fractional Gronwall-type inequality which can be used in the numerical analysis of these schemes avoiding any local assumptions. However, that discrete fractional Gronwall-type inequality can't be directly applied to time-fractional differential equations with time delay. Recently, Li et al. [233] established a new discretized fractional Gronwall-type inequality to analyze fractional reaction-diffusion equations with fixed temporal delay using nonlocal assumptions. Following these updates and some recent applications [509, 345], we present in section 6.4 a discrete fractional Gronwall inequality which is consistent with the L2 - 1a to deal with the analysis of multi-term time-fractional partial differential equations. The main aim of the proposed discrete Gronwall inequality over previous works that it can be used to provide optimal error estimates for multi-term fractional problems with nonlinear delay. Next, the use of the deduced inequality to improve and correct the analysis of some numerical approaches in literature is the other target of that consideration.

Considering a two-dimensional nonlinear Volterra integro-differential equation with a weakly singular kernel is the scope of section 7.1. Such integro-differential equations with Riemann-Liouville integral operators appear frequently in various mathematical and physical models. Such problems is a commonly used model for studying physical phenomena related to elastic forces. This model is mainly used in the problems of heat conduction, viscoelasticity and population dynamics of materials with memory [130, 164, 294, 357, 491]. In viscoelastic problems, the parameter 1 in this model represents the Newtonian contribution to viscosity, and the integral term represents the viscosity part of the equation.

Semilinear Sobolev-type integro-differential equation (STIDE) is an abstract formulation of partial integro-differential equation (PIDE), which usually appears in various applications of physical processes, for instance, the thermodynamics [86], the propagation of small-amplitude long waves [40], the fluid movement through fissured rocks [31] and so on. The need for practical applications forms the necessity to study such equations. A number of papers, devoted to proposing numerical methods to solve fractional PIDEs, have applied finite difference methods for time discretization. We refer to [476, 419, 356] for more details. Although the methods presented in these papers are effective methods for solving PIDEs, all of them are based on uniform temporal meshes, which ignores the fact that the exact solution may not be smooth enough near t = 0, thus none of them achieves second-order convergence in time. From the ideas of [57, 384], it can be known that nonuniform temporal meshes can provide an efficient approach to obtain reliable numerical approximations of solutions with the singularity at t = 0, which is applicable to (7.45) as well. Indeed, many existing studies have used numerical methods with nonuniform temporal meshes to solve PIDEs. For PIDEs with the R-L fractional integral, the product integration (PI) method is considered widely, see [308] for details. However, to our knowledge, there are still theoretical gaps in error analysis and approximate properties of PI rule employed to discretize R-L fractional integral and equation with practical physical applications has not been considered in any published works. This motivates us to establish a new and rigorous analysis of a finite difference scheme for it. Similar to [308], the main interest of section 7.2 in nonuniform temporal meshes (graded temporal meshes), but the equations, aims and error analysis studied in these papers differ from our work. Furthermore, the primary contributions of the contribution in this section are outlined below. We construct a Crank-Nicolson difference scheme for the approximate solution of problem on a mesh that is graded in time and uniform in space, where a Crank-Nicolson time-stepping method is used to discretize the temporal derivatives and the PI rule (the product trapezoidal integration rule) proposed by [148]

is considered to approximate R-L fractional integral while the standard second-order central difference formula is utilized to approximate Am. Actually, this scheme elaborated just now is a three-level scheme in this paper, which also includes a linearization process and a pre-correction process. For the PI rule presented by [148], we derive some approximate properties. In addition, invoking some ideas from subsection 5.2 in [401] and taking into account the actual situation of the regularity condition, we deduce the approximate error of approximating R-L fractional integral with this PI rule. The error analysis covers both temporal and spatial directional errors, with a focus on the temporal directional error. Combining the regularity condition with the approximate error yields sharp pointwise-in-time error bounds for graded temporal meshes with the mesh grading exponent y > 1. The result indicates that the temporal convergence order is y(1 + a) with 1 < y < y+a and the optimal temporal convergence rate of order 2 can be attained with y > .

The nonlocal evolution equation with usual and tempered convolution kernels can be used to mathematically model the heat transfer theory with memory, population dynamics, and viscoelastic materials; see, e.g., [172] and references therein. An overview of models that refer to materials requiring non-classical memory kernels was introduced in [69] and its references. Tempered kernels were considered in [399] in the context of tempered subdiffusion. The generalized Langevin equation which takes into account the memory effect of the friction, and hence defines a non-Markovian stochastic process was studied in [69]. A generalized Langevin equation was explored for a free particle in the presence of a truncated power-law and tempered Mittag-Leffler memory kernel. There exists some significant work devoted to the numerical analysis for nonlocal evolution equations with tempered kernels (FETs). Chen and Deng [84] proposed the fractional linear multi-step approach with the Fourier transform to analyze the properties of the discretized scheme for FETs, also in [85], they introduced a new second-order convolution quadrature scheme. Zaky [502] proved the existence and uniqueness of solutions for tempered fractional boundary value problems. Then, Li et al. [224] analyzed the well-posedness of the ordinary tempered fractional differential equations. Sultana et al. [406] proposed the quadratic and quadratic-linear schemes for tempered fractional integrodifferential equations. Fernandez and Ustaoglu [124] derived certain significant analytic properties of FETs. Guo et al. [161] proposed efficient fractional multistep methods for FETs and designed a fast algorithm for computing the quadrature weights. Recently, some work has also considered FETs, see [355, 442]. Closer to this research point, we refer to the work of Lubich [261]. We are devoting some efforts to constructing and analyzing a second-order Crank-Nicolson (0 = 1) convolution quadrature product rule for the integral term, for the convolution integrodifferential equations. To the best of our knowledge, there is no published work devoted to considering the 0 € (1, 1) convolution quadrature product rule numerical scheme for the integral term with a = 0 for the convolution integrodifferential equations and a > 0 for the tempered integrodifferential equations. Accordingly, we are planning to fill this gap in the coming context. It is well known that introducing fully discrete schemes for problems containing integral terms with weakly singular kernels is of concern in practical application and theoretical analysis. Armed with that fact and in order to construct a fully discrete scheme for (7.111), we propose a high-order OSC method [51, 50] for the approximation of the spatial operators. Bearing these ideas in mind, the main aspects of section 7.4 can be summarized as follows: The numerical scheme of Crank-Nicolson convolution quadrature product rule in [261] are expanded to the 0-type (0 € (2, 1)) convolution quadrature product rule in case of usual convolution structure a = 0 and tempered convolution structure a > 0. Besides, two different kinds of 0-type convolution quadrature product rules are also proposed for the convolution integro-differential equations and the tempered integro-differential equations, respectively. Detailed error estimates for the proposed 0-type numerical schemes are derived according to both cases of convolution structures (usual and tempered), respectively. Novel theoretical analysis skills are provided to achieve our goal.

The nonlinear Volterra integro-differential equation (VIDE) with Abel kernel appears in many physical models, such as heat conduction problems in materials with memory and problems involving viscoelastic force correlations; see, e.g., [172]. In viscoelastic force problems, for instance, ^ denotes the Newtonian contribution to the viscosity, and the integral term in it represents the viscous part of the equation. For the linear case of problem, many scholars have constructed a variety of numerical methods, such as the Galerkin finite-element scheme [286, 355], the h-p Continuous Petrov Galerkin methods [493], the finite difference schemes [476], the ADI schemes [442], and so on. Furthermore, for the nonlinear case of problem, Lopez [213] designed and analyzed a backward Euler method.

Mustapha [308] studied the generalized extrapolated Crank-Nicolson scheme for semilinear integral-differential equations with weakly singular kernels on graded meshes. Cao et al. [67] utilized a local meshless technique to solve nonlinear integro-differential equations with multi-term kernels. However, the computational complexity of solving nonlinear problems increases dramatically as the mesh is encrypted and the spatial dimensionality increases. To improve the computational efficiency, Xu [477, 478] proposed the two-grid method, which has the advantage of saving computational cost without loss of accuracy. Inspired by this idea, some researchers started to use the time two-grid method to solve nonlinear problems [76, 75] and in particular VIDEs. For example, a time two-grid finite difference method was designed and the stability and convergence of the algorithm were analyzed in [475]. Wang et al. developed the spatial two-grid (STG) finite element methods for model (7.144) with 1 = 1 and either a nonlinear memory [456] or a nonlinear source term [455]. To our best knowledge, the STG finite difference method for that problem has not been considered, especially for the case that 1 = 0. Since the finite difference method obtains discrete solutions on a grid, it is necessary to map the solution on a coarse grid to a fine grid to maintain spatial accuracy when performing the STG finite difference method. Most of the existing research works use piecewise linear interpolation to construct second-order two-grid difference algorithms, see e.g. [98]. To improve the accuracy of the STG finite difference method, Fu et al.[133] proposed and analyzed a higher-order mapping operator between two grids that can be applied to the corresponding higher-order two-grid difference scheme of nonlinear partial differential equations for theoretical analysis. Inspired by [133], we construct the STG compact difference method on a uniform time grid to solve the problem in section 7.5, which incorporates the PI rule to overcome the singularity of the solutions at t = 0. The STG compact difference algorithm can be divided into two steps as follows: Step I, on the coarse grid, we establish a nonlinear compact difference method to obtain the coarse solution UH by solving a small nonlinear system; Step II, on the fine grid, based on the rough solution UH in Step I and the higher-order mapping operator nH, we obtain the accurate numerical solution Uh by solving a large-scale linear system. Furthermore, we rigorously prove that the proposed two-grid algorithm is stable and convergent with accuracy 0(t2 + H8 + h4) under appropriate regularity conditions and the discrete L2 norm, where H and h are the spatial steps on the coarse and fine grids, respectively. Then, the construction and analysis of the algorithm can be extended to the two-dimensional (2D) case. Finally, several numerical experiments are given to verify the effectiveness and efficiency of the algorithm.

While the foundational theory of fractional calculus is mature, the numerical analysis and practical simulation of nonlocal models remain an area of active and ongoing development. This research advances that frontier, addressing critical limitations in existing methods and enabling the accurate and efficient simulation of a broader class of nonlocal phenomena with direct implications across physics, engineering, and the life sciences.

Purpose of the Work. The primary purpose of this thesis is to develop, analyze, and implement robust numerical methods for solving nonlinear time-space fractional partial differential and integro-differential models with weakly singular kernels. These models, which incorporate fractional derivatives and integral terms with weakly singular kernels, naturally arise in the mathematical modeling of various complex phenomena such as anomalous diffusion, quantum transport, viscoelasticity, heat conduction, and biological pattern formation. The nonlocal nature of fractional operators introduces substantial analytical and computational challenges that demand specially tailored numerical approaches capable of capturing memory effects, spatial nonlocality, and long-time dynamics.

The main objectives of this thesis are as follows:

• To ensure the numerical preservation of fundamental physical properties such as energy, mass, and momentum, as well as boundedness and positivity, depending on the nature of the modeled system.

• To construct and validate hybrid spectral-difference numerical schemes that accurately approximate solutions to a wide class of time-space fractional PDEs and integro-differential equations, even in the presence of low regularity or singularities.

• To discuss the existence and uniqueness for different sorts of variable order fractional models by adapting the well known Rothe method.

• To address the existence and uniqueness for fractional order inverse source problems by novel techniques based on nonuniform mesh Rothe techniques.

• To address and overcome challenges associated with time delays, distributed-order derivatives, and singular kernels using adaptive grid techniques, correction strategies, and discrete variational formulations.

• To apply the proposed numerical methods to physically relevant systems including:

- Nonlinear space-time fractional Schrödinger and Klein-Gordon-Zakharov equations,

- Fractional reaction-diffusion equations with delay and memory,

- Semilinear space fractional damped Klein-Gordon equation with a generalized scalar potential,

- Multi-dimensional fractional integro-differential equations with weakly singular kernels.

• To provide rigorous mathematical analysis for the stability, convergence, and efficiency of the proposed schemes using energy-based techniques, fractional Sobolev spaces, and fractional Grönwall inequalities.

• To extend the computational framework to the efficient treatment of large-scale three-dimensional problems, where second-order accurate schemes and alternating direction implicit (ADI) methods are used to handle fractional Volterra-type integro-differential equations with weakly singular memory kernels.

Ultimately, this thesis offers a unified computational framework that blends analytical rigor with practical algorithm design for a broad range of fractional and nonlocal models. The proposed methods are not only mathematically sound but also applicable to real-world problems in physics, engineering, and applied sciences where memory, delay, and spatial nonlocality play a crucial role in system behavior.

Scientific Novelty. The dissertation contains new scientific results in all areas of the specialty.

• Results in the field of mathematical modeling. It offers new theoretical insights and practical tools for the computational modeling of fractional and nonlocal phenomena. Its findings extend the boundaries of existing methods and provide a reliable computational infrastructure for solving a wide spectrum of fractional order models encountered in modern physics, biology, and engineering. Our computational modeling techniques applied to subdiffusion and suberdiffusion models with and without delay effect, time and space fractional order Schrödinger and Klein-Gordon models, spatial fractional order Klein-Gordon models with Higgs potential, spatial fractional order Klein-Gordon-Zakharov models, heat conduction with memory and multi-dimensional fractional integro-differential models with weakly singular kernels. We computationally succeeded in modeling computationally some physical properties of the considered models such as mass and energy in a discrete sense. This gives us an insight about how to study the physical behaviour of the continuous models by showing such behavior computationally. We conducted a computational investigation of several key phenomena, including the emergence of Turing patterns in hyperbolic fractional reaction-diffusion models, the process of nonlinear supratransmission in Riesz space-fractional sine-Gordon equations, and the mechanisms of energy dissipation in a fractional-order Higgs boson model within the de Sitter spacetime framework.

• Results in the field of numerical methods. The scientific novelty of this work lies in its comprehensive and original contributions to the numerical analysis of nonlinear fractional partial differential equations and fractional integro-differential equations. It presents a unified framework for the design and analysis of structure-preserving numerical schemes tailored for systems characterized by memory effects, spatial nonlocality, and singular integral kernels. These schemes are developed with a strong emphasis on maintaining the intrinsic physical and mathematical properties of the underlying models, such as conservation laws, dissipation, boundedness, and long-time stability. One of the key innovations is the construction of hybrid numerical methods that integrate Legendre spectral discretizations in space with high-order temporal schemes, such as

graded L1-type and fractional multistep methods. These techniques enable accurate simulation of nonlinear time-space fractional PDEs with varying regularity, while ensuring the preservation of energy and other integral invariants in both conservative and dissipative contexts. Another major contribution is the incorporation of fractional models with weakly singular kernels and time-delay effects, which are common in realistic applications such as viscoelasticity, anomalous diffusion, and biological transport. The thesis introduces graded temporal meshes and singularity-aware basis functions to tackle the challenges posed by non-smooth data and singular memory terms. Furthermore, it provides a detailed treatment of multi-dimensional problems through the development of second-order accurate schemes, including novel two-grid, alternating direction implicit (ADI) orthogonal spline collocation methods and fast schemes-particularly for fractional Volterra-type integro-differential equations. The work is further distinguished by its rigorous analytical foundation. Stability and convergence analyses are conducted using discrete energy estimates, discrete Sobolev embeddings, and generalized Gronwall-type inequalities adapted for the fractional setting. The thesis also adapted Rothe methods to be applied effectively for the existence and uniqueness of variable order fractional models and time fractional inverse source problems. These results are achieved without assuming high regularity of the exact solutions, thereby enhancing the robustness and applicability of the proposed schemes.

• Results in the field of software development. The compatibility of the proposed schemes with widely used software environments such as MATLAB, Python, Mathematica, and C++ libraries ensures that the developed methods can be independently tested, reproduced, and validated across a range of computational platforms. The accompanying algorithmic implementation guidelines may serve as methodological prototypes for subsequent software development, thereby promoting transparency, reproducibility, and cross-platform integration in computational fractional-order modeling. Furthermore, based on the outcomes of the computational research, eight computer programs were registered as intellectual property objects, underscoring the practical relevance and applicability of the proposed numerical methodologies.

Provisions to be Defended. This thesis puts forward several significant provisions that form the foundation of its contributions to the field of computational mathematics and fractional modeling.

Firstly, the thesis advances the field by extending the structure preserving numerical schemes to fractional-order equations involving Riesz and Caputo fractional derivatives. This adaptation preserves the variational structure in the discrete setting and contributes a framework for deriving structure-preserving methods applicable to a broad class of physical systems with fractional dynamics. This is done side by side with preserving high order of convergence and stability to the constructed scheme, which represents some challenges to preserve both physical and numerical properties.

Secondly, it defends the novel construction of hybrid numerical schemes that integrate Legendre spectral methods for spatial discretization with high-order time integration strategies based on nonuniform grids. This aids to achieve spectral accuracy in space and coping better near the initial singularities of Caputo fractional derivatives, when constructing such hybrid schemes for different sorts of time and space fractional order problems. Also, it has been used successfully to investigate the long time behaviour for time Caputo fractional subdiffusions with delay. These schemes are tailored to preserve key physical properties such as energy, momentum, and mass in conservative systems, as well as to accurately represent dissipative behavior in non-conservative models. The capacity of these methods to maintain structure and achieve long-time stability distinguishes them from traditional approaches. We applied such schemes for Riesz space fractional Klein-Gordan equations with Higgs potential. Thirdly, the existence and uniqueness Theorems for the solutions of some direct and inverse fractional order problems, based on semidiscretization (Rothe) method, are given full consideration. We state and prove such theorem for direct problems such as nonlinear variable order subdiffusion equations with delay, space dependent variable order superdiffusion equations and wave equations with time dependent variable order derivatives. Also, we state and prove such theorems for one term and multiterm time fractional inverse source problems. At which, we construct Rothe methods based on graded meshes to cope better near the initial singularities of time Caputo fractional derivatives, for the first time in literature.

Fourthly, A particularly novel and robust component of the thesis is the rigorous numerical analysis carried out using discrete fractional Grönwall-type inequalities. These inequalities serve as powerful tools for establishing a priori bounds on numerical solutions and are indispensable for proving the stability and convergence of schemes involving nonlocal time derivatives. Their use in this thesis demonstrates a deep integration of analytical techniques with algorithmic development, enabling reliable long-time simulation of fractional systems without resorting to overly restrictive smoothness assumptions. We construct such discrete inequalities and successfully implement to wide range of fractional differential equations with different sorts of delay.

Also, the thesis introduces new computational techniques for multidimensional fractional integro-differential equations with weakly singular kernels such as nonlinear Volterra integro differential equations, Sobolev type equations, nonlocal evolution equations arising in heat conduction with memory and Volterra equations with Abel kernel. The temporal two-grid, orthogonal spline collocation, spatial two grid and BDF2 alternating direction implicit (ADI) methods developed here are demonstrated to achieve high order of convergence and stability with respect to optimal estimates, which is a substantial contribution to handling realistic memory-dependent models. Finally, the work provides a validated methodology for dealing with nonsmooth data and delay terms, which are often present in real-world applications but difficult to treat numerically. The use of adaptive temporal meshes, correction techniques, and singularity-aware basis functions ensures that the proposed methods remain accurate and stable under challenging conditions. In addition, the thesis establishes theoretical guarantees through the derivation of optimal error estimates and stability bounds using advanced techniques. These results strengthen the mathematical rigor of the proposed schemes. The thesis successfully applies the constructed methods to various physically significant problems-such as the fractional Schrödinger, Klein-Gordon, Klein-Gordon-Zakharov systems, non-Fickian fractional reaction-diffusion equations, Volterra and Sobolev integro equations and nonlocal evolution problems, demonstrating their practical relevance and adaptability. The ability to generalize and extend these methods across a spectrum of complex models underlines the robustness and versatility of the computational framework presented.

Methodology and Research Methods. The methodology employed in this thesis is grounded in the meticulous design, development, and rigorous analysis of numerical schemes specifically tailored for time-space fractional partial differential equations (PDEs) and integro-differential equations with singular memory kernels. The core objective of these methods is to preserve the intrinsic physical and mathematical structures of the underlying continuous models, such as conservation laws, dissipation, energy invariants, and variational principles, in their discrete analogs. At the foundation of the proposed computational framework is the development of hybrid numerical schemes. These schemes combine the spectral accuracy of Legendre-based spatial discretization with various high-order temporal approximation strategies, including the L1 method, fractional linear multistep methods, and time-splitting procedures. This hybridization is strategically designed to accommodate the nonlocal and memory-intensive nature of fractional derivatives, enabling the simulation of both smooth and nonsmooth solutions with high fidelity. A central feature of the methodology is the extension of the discrete variational techniques to fractional systems. That techniques offers a systematic approach to deriving numerical schemes that inherently preserve energy or other conserved quantities. In this work, it has been adapted to handle equations with Riesz and Caputo fractional derivatives, ensuring that the numerical approximation respects the variational structure and conservation laws intrinsic to the continuous models. To tackle problems involving weakly singular kernels, delay terms, or low-regularity solutions, the thesis introduces specialized discretization techniques. These include memory-preserving quadrature formulas, time grading techniques that adaptively refine the temporal mesh to handle initial singularities, and orthogonal spline collocation methods particularly effective in multidimensional settings. Such tools are essential in maintaining accuracy without incurring excessive computational costs. A significant methodological advancement of the thesis is the integration of discrete fractional Grönwall-type inequalities into the stability and convergence analysis. These inequalities serve as indispensable tools in deriving uniform-in-time a priori bounds for discrete solutions of nonlocal-in-time equations. Their inclusion allows for rigorous proofs of stability and convergence for the proposed numerical schemes under realistic regularity assumptions. The theoretical framework is further reinforced by analytical tools from fractional Sobolev space

theory. These tools enable the derivation of sharp error estimates in both spatial and temporal norms, accommodating the intrinsic complexity of fractional operators and accounting for singular behaviors that arise naturally in physical applications. The research methodology also includes comprehensive numerical experimentation. The developed methods are validated against benchmark problems with known solutions and are applied to a variety of real-world models, such as fractional Schrödinger equations, nonlinear Klein-Gordon-Zakharov systems, and time-fractional reaction-diffusion equations with memory and delay. These experiments confirm the theoretical findings and demonstrate the robustness, efficiency, and wide applicability of the proposed computational strategies. In summary, the methodology adopted in this thesis offers a cohesive and innovative framework that merges spectral accuracy, structure-preserving discretization, rigorous fractional analysis, and robust numerical implementation for the reliable solution of complex nonlocal and memory-dependent models.

Theoretical and Practical Significance of the Work. The theoretical and practical significance lies in its comprehensive contribution to the numerical analysis and computational modeling of fractional partial differential equations and integro-differential systems with weakly singular kernels. These classes of problems, governed by nonlocal dynamics and memory effects, pose substantial challenges in mathematical modeling and numerical simulation. The thesis addresses these challenges by developing finite difference techniques and spectral methods to incorporate fractional-order derivatives and weakly singular integrals, thereby offering a robust computational framework suitable for modern nonlocal models. These computational frameworks are applied successfully to wide range of mathematical models built on such abstract equations in different disciplines of science. It has been shown that retaining some physical properties of the original continuous models, can be preserved through the constructed numerical schemes. The existence and uniqueness of variable order fractional order models can be theoretically achieved by adapting Rothe methods. Reconstruction of time dependent parameters in fractional order inverse models, has been done through novel adaption of Rothe methods over graded meshes. The advancing of grid schemes for integro-differential models with weakly singular kernels through novel techniques in construction and theory, has been achieved. Also, fast schemes have been developed for fractional order models. Furthermore, the incorporation of discrete fractional Grönwall-type inequalities into the stability and convergence proofs of the proposed schemes, especially for models with time delay, adds a powerful tool to the existing arsenal of numerical analysis for models that incorporate fractional order derivatives and delay effects, allowing researchers to derive sharp a priori bounds even in the absence of high regularity. These discrete inequalities strengthen the theoretical foundations for solving real-world fractional models and highlight the deep interplay between abstract functional analysis and practical numerical computation. From a practical perspective, the significance of this work is equally profound. The developed methods offer efficient, stable, and highly accurate algorithms for simulating a wide range of physical systems characterized by memory and hereditary properties. These include, but are not limited to, anomalous transport phenomena in heterogeneous media, quantum field dynamics with fractal structures, time-fractional diffusion in biological tissues, and wave propagation in viscoelastic materials. The presented algorithms are not only theoretically sound but also computationally viable. By adaptive time-stepping strategies, graded meshes, and parallelizable spline-based collocation methods, the schemes exhibit excellent scalability and robustness in practical applications. These features make them ideal for integration into simulation tools and software packages used in engineering, physics, and biomedical research. Moreover, the numerical algorithms and techniques proposed are well-suited for incorporation into open-source and commercial scientific computing platforms. The modular structure and compatibility of the schemes with modern programming paradigms facilitate their integration into packages such as MATLAB, Python (with NumPy and SciPy), Mathematica, and C++-based simulation libraries. This allows for the broader dissemination and adoption of the developed methods by researchers and practitioners working on real-time modeling and simulation tasks in various disciplines. By combining deep theoretical insights with implementable numerical methods, this work serves as a bridge between abstract mathematical modeling and the practical demands of scientific computing. The results obtained are poised to inform the design of future algorithms, stimulate further investigations into structure-preserving discretizations for nonlocal systems, and enhance our understanding of fractional order models in applied science and technology.

The Degree of Reliability of the Research Results. The reliability of the results presented is ensured through a rigorous combination of theoretical analysis, numerical validation, and methodological consistency. Each numerical scheme proposed in the work is supported by a solid mathematical foundation, including detailed proofs of stability, convergence, and error estimates. These analytical guarantees are derived using advanced theoretical tools such as discrete fractional Gronwall inequalities, and variational formulations tailored for nonlocal models. In addition to theoretical rigor, the reliability of the results is reinforced by extensive computational experiments conducted across a wide variety of test cases. These include problems with known exact solutions, which allow for the direct assessment of convergence behavior and numerical accuracy, as well as realistic physical models where the preservation of energy, mass, or dissipation serves as a qualitative benchmark. The alignment of numerical outcomes with expected physical properties provides further evidence of the correctness and robustness of the proposed methods. All computational experiments were carried out using reproducible numerical implementations, and results were verified for consistency under grid refinement and parameter variation. The stability and long-time behavior of the methods were carefully analyzed, particularly in the context of nonlinear and memory-dependent problems, which are known to present significant numerical challenges. Moreover, the modular nature of the implemented algorithms makes them amenable to reuse and extension, thus facilitating verification by independent researchers. The compatibility of the schemes with widely used software packages such as MATLAB, Python, Mathematica, and C++ libraries ensures that the methods can be independently tested, reproduced, and validated in diverse computing environments. Algorithmic implementation guidelines developed may also serve as templates for future software development, encouraging a culture of transparency, collaboration, and cross-platform integration in computational fractional order modeling.

Personal Contribution of the Author. The author has made significant original and primary contributions in both the theoretical and computational aspects of the research results. The entire theoretical analysis-encompassing the derivation of discrete energy identities, a priori estimates, convergence and stability estimates was primarily carried out by the author. The author also carried out a comprehensive review of the literature, identifying existing gaps in the field and formulating the research objectives to address them effectively. The author contributed substantially to designing and implementing most of the numerical schemes proposed in the thesis, including hybrid finite difference and spectral methods. Most of the computational experiments presented in the thesis were designed, executed, and interpreted by the author. These include validation with known analytical solutions, investigation of long-time behavior, and simulation of realistic physical scenarios governed by fractional PDEs. We here specify the contribution of the author with some details.

• In [541-549] and [573-580], the conceptual framework for the extension of structure preserving schemes to mathematical models with fractional order derivatives was formulated and developed by the author. This includes the rigorous adaptation of the preserving energy methods to fractional models and the integration of discrete energy inequalities for the theoretical analysis of stability and convergence. The theorems of preserving physical properties such as energy preserving and numerical properties such as consistency, stability and convergence, were conducted by the author. The author also takes a part in the describing of the numerical results and how they coincide with the theoretical achievements. The partial contribution of coauthor Macias-Diaz, to the joint publications in that direction consisted of handling the concept of energy functions for the continuous fractional models, construction of the discrete energy function accordingly, techniques of involving average operators to guarantee the numerical properties (convergence and stability) in the meantime and numerical implementation of the constructed schemes. The coauthor Martinez in [538] contributes partially to the solvability of constructed numerical scheme. The coauthor Abbaszadeh in [576] contributes to the numerical implementation of the constructed scheme to show the long time behaviour for the fractional order model under consideration. The coauthor Flores in [577] contributes partially to the numerical implementation and the barrier functions in the proofs of convergence and stability estimates. The coauhor Markov in [607] contributes to the high-performance implementation of the bounded numerical solver for a fractional FitzHugh-Nagumo equation. The coauthor Zaky in [550,576, 578,581, 615] contributes partially to the theory of the continuous models under study, stability analysis and numerical implementation.

• In [550-553], [569-572], [582-587], [589-593], [566] , [597] and [600], the conceptual framework of the creation of hybrid numerical schemes that can combines the advantages of Galerkin spectral and finite difference methods for time and space fractional order models, was formulated and developed by the author. The author takes a part in the novel construction of such hybrid numerical schemes that integrate Legendre spectral methods for spatial discretization with high-order time integration strategies based on uniform and nonuniform grids. The theorems of convergence and stability estimates were conducted by the author. The author also takes a part in the computer implementation of the constructed schemes and describing of the numerical results and how they coincide with the theoretical results. The coauthor Zaky in [550-553], [569571], [581-587], [589-593], [562] , [597] and [600] was involved partially in the construction of Galerkin spectral and hybrid schemes, algorithmic basement and the numerical implementation. The coauthors Ab-baszadeh and Dehghan in [592, 597] were contribution to the construction of the free element stabilization technique and rational RBF partition of unity approach mechanism and their numerical implementations. Taha and Suragan in [553] contributed partially to the theory behind the mathematical model under study and to the explanation of the numerical findings. Van Bockstal in [552] contributed partially to the theoretical findings of the longtime behavior for the mathematical model under consideration. Manimaran and Shangerganesh in [590, 591] contributed to the construction and analysis of BDF2 schemes and its basement for the mathematical models with nonlocal diffusion.

• In [558-561], the conceptual framework of the extension of hybrid numerical schemes to variable order fractional problems and using Rothe methods to discuss the existence and uniqueness of such problems, was formulated and developed by the author. The author takes a part in the construction of the semidiscretization schemes to be adaptable for variable order problems and contributed to stating and proving the theorems related. The author contributes to the computer implementation of the constructed schemes and their explanations. Van Bockstal contributed to the adaption of Rothe methods for variable order mathematical models and the existence theorems [554-557]. Zaky contributed partially to the numerical implementation of the full discretisation models [554-557].

• In [558], [559], the conceptual framework of construction of Rothe methods based on graded meshes for time Caputo fractional inverse source problems was formulated by the author. The author contributed to the proofs of the related existence and uniqueness theorems and the explanation of the numerical results. The coauthor Van Bockstal contributed partially to the theory of the fractional inverse source models, the basement of Rothe methods and the numerical implementations in [558, 559].

• In [560-563], the conceptual framework of construction of convenient discrete fractional Gronwall inequalities for the numerical analysis of time delayed fractional partial differential equations, was formulated by the author. The author contributed to the proofs of the theoretical results and the implementation of numerical results. More exclusively, the author contributed to the convergence and stability theorems and their proofs for some numerical methods constructed to time delayed fractional order models [562, 600, 601]. The partial contribution of other authors V. G. Pimenov, Omran and Zaky to the joint publications consisted of constructing the numerical techniques handling different sorts of delays, creation of the algorithmic implementation to Galerkin spectral methods and theoretical and numerical implementations.

• In [562-568], [602-605] and and [594], the author contributed also to the conceptual framework of the construction of different kinds of grid schemes for weakly singular kernals integro problems. The author contributed to the proofs of the theoretical findings (consistency, convergence and stability) of the constructed schemes and the explanation of numerical results. The coauthors Chen, Qiu, Wang and Cao were involved partially in the creation of the grid schemes and their numerical implementation [564-568], [602], [603]. The coauthors Talaei and Zaky in [604, 594] were contributing partially to the theory of the problem under study and its numerical implementations.

In summary, the work reflects the author's deep engagement with all phases of the research process-from theoretical formulation and numerical design to computational implementation and analytical verification-underscoring a

high degree of scholarly independence and originality.

Approval of the Results. The approval of the results presented in this thesis is carried out through a rigorous and multi-faceted approach that includes analytical verification, benchmark comparisons, and systematic computational experimentation. Theoretical models and numerical methods are validated against known exact solutions where available, ensuring that the implemented schemes produce results consistent with established mathematical behavior. The main results of the work were presented at list of major conferences and seminars: Seminar of the department of computational mathematics and computer science, Institute of natural sciences and mathematics, Ural Federal university, Russia; Seminar of N.N. Krasovskii Institute of Mathematics and Mechanics of the Ural Branch of the Russian Academy of Sciences under the guidance of Professors V.N. Ushakov and A.M. Tarasyev; Seminar of School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China; Seminar of department of mathematical analysis, Ghent university, Ghent, Belgium; International Conference on Computational and Mathematical Methods in Science and Engineering (CMMSE 20182019) in Spain; 10th International Conference "Inverse Problems: Modeling and Simulation"held on May 22 - 28,2022, Malta; International Workshop on Computing Technologies and Applied Mathematics, Blagoveshchensk, 2023; Fractional calculus seminars, SISSA, International School of Advanced Studies, Italy, 2024.

Publications. A total of 75 peer-reviewed publications [541-615] and 8 registered software certifications [616-623] authored and co-authored by Hendy resulted directly from the research presented in this dissertation. These contributions are not only reflective of the scope and depth of the theoretical investigations undertaken but also demonstrate the applicability and robustness of the developed numerical techniques in solving complex real-world problems. The publications collectively validate the methodological innovations proposed in the thesis. They include rigorous numerical analysis supported by convergence proofs and benchmark simulations, aligning closely with the theoretical framework laid out in the core chapters of this research.

Structure and Scope of the Work. This dissertation is structured into seven meticulously developed chapters, each designed to build upon the preceding work and collectively provide a coherent and in-depth exploration of structure-preserving numerical methods for fractional and memory-dependent partial differential equations. The scope integrates rigorous mathematical theory, algorithmic innovation, and detailed computational implementation.

Chapter 1 lays the theoretical and historical groundwork by introducing the mathematical basis of fractional calculus and the significance of structure-preserving methods. It outlines the scientific motivation, literature survey, and sets the objectives of the research.

Chapter 2 focuses on the construction of energy-conserving finite difference schemes for fractional generalizations of the nonlinear Klein-Gordon-Zakharov system. It provides a careful derivation of implicit and explicit discretizations, and theoretical analysis of their energy-conserving properties.

Chapter 3 develops Galerkin-Legendre spectral schemes for time-space fractional partial differential equations with both smooth and nonsmooth initial data. It extends the analysis to include fractional Schrodinger-type equations in fractal and curved geometries, with rigorous proofs of stability and convergence.

Chapter 4 addresses the numerical simulation of variable-order fractional partial differential equations. These models reflect systems where the memory effect evolves over time or space, and are critical for accurately describing complex physical and biological processes. The chapter develops robust numerical strategies, including time-stepping methods with adaptive capabilities and appropriate function space formulations, ensuring precision and stability in scenarios involving continuously changing fractional orders. It proposes graded temporal meshes and convolution quadrature techniques that retain accuracy in long-time simulations and capture solution singularities.

Chapter 5 addresses inverse source problems for time-fractional differential equations. These problems involve determining unknown source terms from partial observations of the system's evolution, which is critical in applications such as medical imaging, geophysics, and environmental modeling. The chapter develops stable and accurate numerical schemes, including regularization strategies and iterative solvers, to handle the inherent ill-posedness of inverse problems in the fractional framework. It introduces stable numerical techniques that ac-

count for both hereditary effects and spatial-temporal coupling, with a focus on models arising in biological and ecological systems.

Chapter 6 presents a comprehensive numerical analysis of time-fractional nonlinear subdiffusion equations with delay. This chapter addresses the complexities arising from both the nonlinearity and the hereditary effects embedded in the model. Special focus is placed on developing and analyzing time-stepping schemes that maintain stability and accuracy over long time intervals. Additionally, the chapter leverages discrete fractional Gronwall-type inequalities to establish rigorous stability and convergence results for the proposed numerical methods. These inequalities are employed to rigorously justify the stability and convergence of the developed numerical schemes for time-fractional models.

Chapter 7 is dedicated to the development of grid-based numerical schemes for solving fractional integro-differential equations with weakly singular kernels. These problems are particularly relevant in the modeling of viscoelastic materials and anomalous transport processes. The chapter introduces a two-grid alternating direction implicit spline collocation method, which significantly reduces computational cost while maintaining high accuracy. Detailed theoretical analysis, including error estimates, is provided alongside comprehensive numerical experiments to validate the proposed schemes in multidimensional settings. for solving multidimensional fractional integro-differential equations with weakly singular kernels. The chapter offers detailed error estimates and performance validation.

Chapter 8 focuses on the development and integration of software packages tailored for solving fractional partial differential equations. It presents practical tools and frameworks implemented in environments such as Mathematica, designed to efficiently support the numerical schemes proposed throughout the thesis. The chapter also discusses implementation challenges, computational performance, user-interface design considerations, and the validation of software through benchmark problems and case studies. These software solutions play a key role in facilitating the practical application and dissemination of structure-preserving algorithms.

Chapter 9 Concludes by highlighting the main outcomes, situating them within the context of relevant literature, and indicating avenues for future investigation.

Each chapter is complemented by extensive numerical experiments designed to illustrate the accuracy, structure-preservation, and computational performance of the proposed methods. The dissertation as a whole constitutes a significant contribution to the field of computational fractional calculus, combining advanced numerical theory with practical algorithmic strategies tailored for complex nonlocal systems.

The thesis concludes with a comprehensive synthesis of the theoretical advancements and computational methodologies developed throughout the work. It highlights how each chapter collectively contributes to the overarching goal of constructing reliable, efficient, and structure-preserving schemes for fractional and memory-dependent partial differential equations. Special emphasis is placed on the theoretical innovations, such as the application of discrete fractional Gronwall-type inequalities, and the practical contributions, including the implementation of numerical schemes in versatile software environments.

Furthermore, the conclusion outlines several promising directions for future investigation. These include the extension of the developed techniques to stochastic fractional systems, their application to fractional optimal control problems, and the integration of machine learning strategies with structure-preserving numerical solvers for enhanced modeling of complex, data-driven phenomena. This forward-looking perspective reinforces the thesis's potential as a foundation for continued research in a wide range of scientific and engineering domains.