# Accelerating Fractional PINNs using Operational Matrices of Derivative

Tayebeh Taheri Email: {ttaherii1401<EMAIL_ADDRESS>Department
of Computer and Data Sciences, Faculty of Mathematical Sciences, Shahid
Beheshti University, G.C. Tehran, Iran Alireza Afzal Aghaei††footnotemark:
Department of Computer and Data Sciences, Faculty of Mathematical Sciences,
Shahid Beheshti University, G.C. Tehran, Iran Kourosh Parand Email:
<EMAIL_ADDRESS>Corresponding author Department of Computer and Data
Sciences, Faculty of Mathematical Sciences, Shahid Beheshti University, G.C.
Tehran, Iran Department of Cognitive Modeling, Institute for Cognitive and
Brain Sciences, Shahid Beheshti University, G.C. Tehran, Iran  Department of
Statistics and Actuarial Science, University of Waterloo, Waterloo, Canada

###### Abstract

This paper presents a novel operational matrix method to accelerate the
training of fractional Physics-Informed Neural Networks (fPINNs). Our approach
involves a non-uniform discretization of the fractional Caputo operator,
facilitating swift computation of fractional derivatives within Caputo-type
fractional differential problems with $0<\alpha<1$. In this methodology, the
operational matrix is precomputed, and during the training phase, automatic
differentiation is replaced with a matrix-vector product. While our
methodology is compatible with any network, we particularly highlight its
successful implementation in PINNs, emphasizing the enhanced accuracy achieved
when utilizing the Legendre Neural Block (LNB) architecture. LNB incorporates
Legendre polynomials into the PINN structure, providing a significant boost in
accuracy. The effectiveness of our proposed method is validated across diverse
differential equations, including Delay Differential Equations (DDEs) and
Systems of Differential Algebraic Equations (DAEs). To demonstrate its
versatility, we extend the application of the method to systems of
differential equations, specifically addressing nonlinear Pantograph
fractional-order DDEs/DAEs. The results are supported by a comprehensive
analysis of numerical outcomes.

Keywords— Physics-informed neural networks, nonlinear differential equations,
Fractional derivative, Operational matrix

##  1 Introduction

Fractional calculus (FC) [1], a specialized field in mathematical analysis,
expands traditional differentiation and integration by embracing non-integer
orders. Fractional calculus has gained substantial prominence over four
decades. As an evolving arena, FC anticipates the introduction of numerous
models for real-world applications in science and engineering, particularly in
areas where nonlocality plays a crucial role. This burgeoning field
demonstrates profound applications across diverse scientific disciplines,
extending its reach into the dynamics of the complex real world, with new
ideas implemented and tested using real data. Here, some of these applications
are mentioned. In the domain of underwater sediment and biomedical
applications, fractional derivative models prove valuable for better
understanding wave propagation, providing insights into absorption mechanisms
grounded in the relaxation processes observed in materials such as polymers
[2]. In [3], a Continuous Time Random Walk (CTRW) optimal search framework is
introduced for specific targets with unknown locations. This approach employs
fractional calculus techniques to ascertain optimal distributions for both
search length and waiting time.

Fractional order systems, in contrast to memoryless integer-order dynamic
systems, feature long memory and the application of fractional calculus
improves and extends established control methods and strategies [4]. Its
applications in diverse scientific and engineering fields, including image
processing [5], underscore its recognition, providing valuable tools for
solving differential equations, integral differential equations, partial
differential equations (PDEs), differential-algebraic equations (DAEs), and
delay differential equations (DDEs). The utilization of discrete fractional
calculus emerges as a practical tool to address engineering challenges in
discrete time or space structures, particularly those entailing delays
inherent in delay differential equations. This approach provides stability
theory for fractional difference equations, facilitating long-term control and
mitigating potential errors linked to numerical discretization in continuous
fractional calculus [6].

DDEs represent a distinct class of differential equations where the unknown
functions depend on previous time states, setting them apart from conventional
Ordinary Differential Equations (ODEs). DDEs have diverse applications
spanning mathematics, biological systems, engineering, physics, economics and
finance, chemical reactions, ecology, communication networks, weather and
climate, and medicine [7, 8, 9]. For example, in biological systems, DDEs are
used to model the interactions and time delays in predator-prey systems,
disease spread, and ecological systems. Also, DDEs are employed to describe
neural oscillations and synchronization in the brain, considering the time
delays in signal transmission between neurons. In many of these applications,
DDEs provide a more accurate representation of real-world phenomena compared
to ODEs because they can capture the impact of time delays on system dynamics.
In addition to standard delay differential equations, there are several other
types of delay equations used in engineering to model systems with time
delays. Pantograph delay differential equations (PDDEs) introduce a more
intricate structure by incorporating both forward and backward time delays in
the equation. PDDEs have arisen in the modeling of pantograph control systems,
cell growth models [10], population growth [11], and electric locomotive [12].
Fractional differential equations [13] give significant meanings in physical
contexts, hence empowering researchers to capture more realistic and nuanced
behavior, as demonstrated in vibration theory [14], and biological systems
with memory [15].

Dealing with delay differential equations can be intricate due to the time
delays involved. The primary methods employed to solve these equations include
numerical, analytical, and semi-analytical approaches. One of the most widely
recognized analytical techniques for solving DDEs is the method of steps [16].
This involves breaking the DDE into small time intervals and solving it as a
sequence of ODEs at discrete time points. Laplace transform techniques, as
demonstrated in the work by [17], is another analytical method. Perturbation-
iteration methods, as discussed in the study by Bahsi [18], fall under the
category of semi-analytical approaches. The analytical and semi-analytical
technique enhances precision and stability while potentially incurring higher
computational costs.

To address the challenges posed by complex systems featuring time delays,
numerical methods for solving delay differential equations play a crucial
role. A variety of numerical techniques have been developed to efficiently
handle DDEs. Finite difference methods [19] constitute a category of numerical
techniques utilized for the solution of differential equations, involving the
approximation of derivatives through finite differences. Examples of finite
difference methods encompass Runge-Kutta techniques [20], Euler method [21],
and Linear multistep approaches, as presented in Shaalini’s work [22].

Spectral methods belong to a group of numerical techniques that make use of
spectral (or Fourier) representations of the solution. They prove to be
particularly effective when addressing DDEs that display periodic or
oscillatory characteristics. Illustrative instances of this approach for
solving DDEs involve collocation techniques [23], methods utilizing wavelets
in conjunction with collocation [24], the application of the Bernoulli
collocation method [25], and the utilization of the transferred Legendre
pseudospectral method [26]. The performance of spectral methods is profoundly
influenced by the selection of basis functions, including Fourier, Chebyshev,
or Legendre, as well as the initial values, introducing challenges in the
decision-making process. Finite element methods, such as the Chebyshev Tau
method [27] and the Chebyshev discretization scheme [28], along with meshless
methods like radial basis functions (RBFs) [29] and meshless local Petrov-
Galerkin techniques [30], represent promising numerical approaches for solving
DDEs. Meshless methods often rely on local approximations, such as radial
basis functions or moving least squares, leading to difficulties in capturing
global features of the solution. This limitation can result in reduced
accuracy, particularly when dealing with DDEs that exhibit complex and non-
local dynamics.

While DDEs capture systems with temporal dependencies, differential-algebraic
equations (DAEs) are a class of differential equations that describe dynamic
systems with algebraic constraints, and both play a crucial role in modeling
complex phenomena. Their applications are evident in a wide spectrum of
disciplines, from aerospace engineering and biomechanics to geophysics and
robotics. To effectively address the diverse applications of DAEs across
various domains, it is essential to employ specialized numerical methods for
their solution. Instances of these techniques comprise a variety of methods,
including the Adomian Decomposition Method (ADM) [31], Pade series [32],
homotopy perturbation [33], and homotopy analysis technique [34], among
others, all contributing to the effective solution of DAEs. Common issues and
limitations associated with numerical methods for DDEs and DAEs include
concerns related to stability, convergence, accuracy, and computational
demands.

Machine learning is gaining traction in solving differential equations,
offering flexibility and efficacy compared to traditional methods. Its
adaptability and ability to learn complex patterns result in more accurate
solutions and faster predictions with reduced computational costs, exhibiting
robust generalization across various data distributions and reduced
sensitivity to discretization choices. In specific contexts, such as solving
fractional-order differential equations, it has exhibited superior performance
through the utilization of techniques like least-squares support vector
regression [35]. In applying machine learning to solve differential equations,
a crucial aspect lies in the specialized use of artificial neural networks
(ANN), showcasing their effectiveness across a diverse range of differential
equation types. Highlighting instances of artificial neural network
applications in solving differential equations, Rahimkhani and Ordokhani [36]
introduced the orthonormal Bernoulli wavelets neural network methodology to
address the Lane–Emden equation. Sabir et al. [37] applied an optimization
approach for ANNs using a blend of genetic algorithms and sequential quadratic
programming to resolve multi-pantograph delay differential equations. Khan et
al. [38] devised intelligent networks with backpropagation, while Panghal et
al. [39] implemented feed-forward ANNs for solutions in delay differential
systems. Ye et al. [40] innovated with the shifted Legendre neural network
method, employing an extreme learning machine algorithm (SLeNN-ELM) for
solving fractional differential equations with constant and proportional
delays. Zhang et al. [41] suggested an approach utilizing neural networks to
solve delay differential equations and derive control policies from specific
regions of controlled systems. Ruthotto et al. [42] adopted a deep
convolutional neural network to solve partial differential equations.

Incorporating orthogonal functions as activation functions in neural networks
enhances their effectiveness in solving differential equations. Methods like
the orthonormal Bernoulli wavelets neural network and the shifted Legendre
neural network illustrate the synergy between spectral activation functions
and the networks’ capacity to capture complex mathematical relationships. In
diverse applications, orthogonal polynomials like Legendre, Chebyshev, and
Jacobi are effectively used as activation functions in neural networks. They
contribute to the accuracy of solving various integral and partial
differential equations, showcasing the versatility of employing specialized
functions in computational methodologies [43, 44].

Recently, a deep neural network framework known as Physics-Informed Neural
Networks (PINNs) has been introduced by Raissi et al. [45]. The PINNs
methodology combines the power of neural networks with the principles of
physics, offering a versatile and effective tool for addressing a wide range
of differential equations, including applications to solving stationary PDEs
[46], nonlinear integro-differential equations [47], heat transfer problems
[48], and dispersive PDEs [49]. Furthermore, PINNs have demonstrated success
in solving fractional differential equations, addressing challenges in areas
such as integro-differential equations of fractional order [50], fractional
differential equations revealing time-dependent parameters and data-driven
dynamics [51], and Hausdorff derivative Poisson equations [52]. In these
studies, the researchers employ discretization methods for fractional
derivatives. Unlike integer calculus, automatic differentiation cannot be
applied to fractional calculus due to the inapplicability of standard rules in
this context.

In this paper, we introduce an alternative non-uniform discretization approach
for Caputo-based fractional differential equations, specifically tailored for
integration into neural networks. Notably, we achieve the efficient
computation of the Caputo fractional derivative ($\alpha$) through the
utilization of operational matrices. Then, by utilizing a special extension of
PINNs, namely Legendre Deep Neural Network (LDNN) our method adeptly tackles
the complexities inherent in fractional differential equations. Overall, the
core aim of this article is to:

  * •

Proposed a physics-informed neural network for solving fractional-order
differential equations.

  * •

Proposed non-uniform discretization approach for Caputo-based FDEs.

  * •

Utilize hidden layers within this network, incorporating Legendre and
Chebyshev polynomials as activation functions.

  * •

Develop an efficient Artificial Neural Network for solving:

    * –

Pantograph delay differential equations.

    * –

Nonlinear delay differential equations.

    * –

Fractional delay differential equations.

    * –

Differential-algebraic equations.

    * –

Linear fractional differential-algebraic equations.

The rest of this paper is organized as follows: section 2 outlines techniques
for calculating both the fractional derivative and the operational matrix of
the Caputo fractional derivative. In the methodology section 3, we clarified
the proposed neural network method and its structure. To facilitate a
comprehensive understanding of the procedure, we introduce our network to
specifically address the resolution of differential equations with fractional
orders. We have assessed the efficacy of our method by employing it to
approximate various problems in the numerical examples section 4. Finally, an
overview of the proposed method is presented in the conclusion section 5.

##  2 Fractional Calculus

In this section, we present methods for computing the fractional derivative
and the operation matrix of the Caputo fractional derivative.

######  Definition 1.

Let $f$ be a continuous function on the interval $[a,b]$, the Reimann-
Liouville fractional derivative of arbitrary order $\alpha$ is defined as:

$\prescript{}{}{\mathcal{D}}^{\alpha}f(x)=\frac{1}{\Gamma(k-\beta)}\frac{d^{k}}{dx^{k}}\int_{a}^{x}(x-t)^{k-\alpha-1}f(t)dt,\quad
k-1\leq\alpha\leq k,$ (1)

where $\Gamma(.)$ is the Gamma function denoted by the following integral:

$\Gamma(z)=\int_{0}^{\infty}e^{-t}t^{z-1}dt,$ (2)

which is convergent on the right half of the complex plane $(Re(z)\geq 0)$.

######  Definition 2.

Caputo’s definition of the fractional order of $\alpha$ is expressed as
follows:

$\prescript{C}{}{\mathcal{D}}^{\alpha}f(x)=\frac{1}{\Gamma(n-\alpha)}\int_{a}^{x}\frac{f^{(n)}(t)dt}{(x-t)^{\alpha+1-n}},\quad
n-1<\alpha<n.$ (3)

######  Property 1.

Let’s denote the fractional order as $q$, where $q=n+\alpha$ and $n$ is the
integer part of $q$. The fractional derivative of order $q$ can be expressed
as:

$D^{q}[f(t)]=D^{n}[D^{1-\alpha}[f(t)]],$

here, $D^{n}$ represents the derivative of integer order $n$, and
$D^{1-\alpha}$ represents the fractional derivative of order $1-\alpha$.

######  Theorem 1.

Let $0<\alpha<1$ and the interval $[t_{0},t_{n}]$ is discretized to $n+1$
points, $0=t_{0}<t_{1}<\dots<t_{n}$. Then the following linear combination
approximates the Caputo fractional derivative of order $\alpha$:

${}_{0}^{C}\mathcal{D}_{t_{n}}^{\alpha}f{(t)}=\sum_{k=0}^{n}\omega_{k}f(t_{k}),$
(4)

where $f$ is the desired function and $\omega_{k}$ are real-valued weights.

###### Proof.

We start by splitting the integration into $n$ non-equidistant intervals

$\displaystyle{}_{0}^{C}\mathcal{D}_{t_{n}}^{\alpha}f{(t)}$
$\displaystyle=\frac{1}{\Gamma{(1-\alpha)}}\int_{0}^{t_{n}}\frac{f^{\prime}(x)}{(t_{n}-x)^{\alpha}}dx$
$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\sum_{k=0}^{n-1}\int_{t_{k}}^{t_{k+1}}\frac{1}{(t_{n}-x)^{\alpha}}f^{\prime}(x)dx$
$\displaystyle\approx\frac{1}{\Gamma(1-\alpha)}\sum_{k=0}^{n-1}\int_{t_{k}}^{t_{k+1}}\frac{1}{(t_{n}-x)^{k}}\frac{f(t_{k+1})-f(t_{k})}{t_{k+1}-t_{k}}dx.$

In each interval one may approximate the $f^{\prime}(x)$ using the forward
scheme finite difference approximation of the derivative:

$\displaystyle{}_{0}^{C}\mathcal{D}_{t_{n}}^{\alpha}f{(t)}$
$\displaystyle\approx\frac{1}{\Gamma(1-\alpha)}\sum_{k=0}^{n-1}\frac{f(t_{k+1})-f(t_{k})}{t_{k+1}-t_{k}}\int_{t_{k}}^{t_{k+1}}\frac{dx}{(t_{n}-x)^{\alpha}}.$

Then, the analytical integration can be utilized to simplify the formulation:

$\displaystyle{}_{0}^{C}\mathcal{D}_{t_{n}}^{\alpha}f{(t)}$
$\displaystyle\approx\frac{1}{\Gamma(1-\alpha)}\sum_{k=0}^{n-1}\frac{f(t_{k+1})-f(t_{k+1})}{t_{k+1}-t_{k}}\cdot\left[-\frac{(t_{n}-t_{k})^{1-\alpha}-(t_{n}-t_{k+1})^{1-\alpha}}{(\alpha-1)}\right]$
$\displaystyle\approx\frac{1}{\Gamma(2-\alpha)}\sum_{k=0}^{n-1}\left[\frac{(t_{n}-t_{k})^{1-\alpha}-t_{n-k+1})^{1-\alpha}}{t_{k+1}-t_{k}}\right]\left[f(t_{k})-f(t_{k+1})\right].$

Defining $\mu_{k}$ as the weight part of the summation, and a simple
reformulation of it yields:

$\displaystyle{}_{0}^{C}\mathcal{D}_{t_{n}}^{\alpha}f{(t)}$
$\displaystyle\approx\frac{1}{\Gamma{(2-\alpha)}}\sum_{k=0}^{n-1}\mu_{k}\left[f{(t_{k})}-f{(t_{k+1})}\right]$
$\displaystyle\approx\frac{1}{\Gamma(2-\alpha)}\sum_{k=0}^{n-1}\left[\mu_{k}-\mu_{k-1}\right]f(t_{k})$
$\displaystyle\approx\sum_{k=0}^{n}\omega_{k}f(t_{k}),$

where
$\omega_{k}=\nicefrac{{\left[\mu_{k}-\mu_{k-1}\right]}}{{\Gamma(2-\alpha)}}$
and:

$\mu_{k}=\begin{cases}\displaystyle\frac{(t_{n}-t_{k})^{1-\alpha}-t_{n-k+1})^{1-\alpha}}{t_{k+1}-t_{k}}&0\leq
k<n\\\ 0&otherwise.\end{cases}$

∎

######  Corollary 1.

for the equidistant data points $t_{k}=kh$ for step size $h$, the formulation
of $\mu_{k}$ can be simplified as:

$\mu_{k}=\begin{cases}\displaystyle\frac{((n-k)h)^{1-\alpha}-{(n-k-1)h})^{1-\alpha}}{h}&0\leq
k<n\\\ 0&otherwise.\end{cases}$

###### Proof.

The proof is straightforward. ∎

######  Corollary 2.

For computing the fractional derivative of Caputo type for $\alpha>1$, one can
calculate the derivative of integer order, followed by the computation of the
derivative of order $1-\alpha$ where $\alpha$ represents the fractional part
of the specified order. (see property 1)

######  Theorem 2.

The operation matrix of the Caputo fractional derivative can be obtained using
the lower triangular matrix $\mathcal{A}$:

$\mathcal{A}=\begin{bmatrix}0\\\ \omega_{0}^{(1)}&\omega_{1}^{(1)}\\\
\omega_{0}^{(2)}&\omega_{1}^{(2)}&\omega_{2}^{(2)}\\\ &\vdots\\\
\omega_{0}^{(n-2)}&\omega_{1}^{(n-2)}&\omega_{2}^{(n-2)}&\omega_{3}^{(n-2)}&\cdots&\omega_{n-2}^{(n-2)}\\\
\omega_{0}^{(n-1)}&\omega_{1}^{(n-1)}&\omega_{2}^{(n-1)}&\omega_{3}^{(n-1)}&\cdots&\omega_{n-2}^{(n-1)}&\omega_{n-1}^{(n-1)}\end{bmatrix}$

Hence, for the vector-valued function $\mathbf{f}_{i}=f(x_{i})$ and arbitrary
nodes $x_{i}$ the Caputo fractional derivative of order $\alpha$ can be
efficiently computed using the operational matrix of derivative:

${}^{C}\mathbf{f}^{\alpha}\approx\mathcal{A}\mathbf{f}.$

###### Proof.

The proof can be completed if we approximate the $i$-th element of
${}^{C}\mathbf{f}^{\alpha}$ using theorem 1:

${}^{C}\mathbf{f}^{\alpha}_{i}=\sum_{k=0}^{i}\omega_{k}^{(i)}\mathbf{f}_{i},$

where
$\omega_{k}^{(i)}=\nicefrac{{\left[\mu_{k}^{(i)}-\mu_{k-1}^{(i)}\right]}}{{\Gamma(2-\alpha)}}$
and:

$\mu_{k}^{(i)}=\begin{cases}\displaystyle\frac{(t_{i}-t_{k})^{1-\alpha}-t_{i-k+1})^{1-\alpha}}{t_{k+1}-t_{k}}&0\leq
k<n\\\ 0&otherwise.\end{cases}$

∎

##  3 Methodology

In this section, we introduce neural network architecture designed for solving
fractional delay differential equations. Importantly, we highlight the
inherent adaptability of this approach, indicating its potential extension to
address differential-algebraic equations. First, we consider a general form of
fractional delay differential equation presented below:

$\phi^{(\alpha)}(\tau)=\chi(\tau,\phi(\tau),\phi^{\prime}(\tau),\phi^{\prime\prime}(\tau),\dots,\phi^{(\lfloor\alpha\rfloor)}(\tau),\phi(\delta(\tau))),$
(5)

where $\delta$ represents the non-negative delay term, and $\alpha$ denotes
the order of the derivative. Now, in the vector form the residual function can
be defined as follows:

$\mathcal{R}(\boldsymbol{\tau})=\mathcal{A}\tilde{\phi}(\boldsymbol{\tau})-\chi(\boldsymbol{\tau},\tilde{\phi}(\boldsymbol{\tau}),\tilde{\phi}^{\prime}(\boldsymbol{\tau}),\tilde{\phi}^{\prime\prime}(\boldsymbol{\tau}),\dots,\tilde{\phi}^{(\lfloor\alpha\rfloor)}(\boldsymbol{\tau}),\tilde{\phi}(\delta(\boldsymbol{\tau}))),$
(6)

where $\tilde{\phi}(\boldsymbol{\tau})$ are the predicted values by the neural
network for the input batch $\boldsymbol{\tau}$ and $\mathcal{A}$ is the
operation matrix to calculate the fractional derivative, that is defined in
the equation (2).

In crafting the loss function, we draw inspiration from the Physics-Informed
Neural Network, integrating both the initial values of the differential
equations and the residual values. It is crucial to take into consideration
the initial conditions of the delay differential equation, as these conditions
are inherently guided by the physical significance of the problem. Below, we
have explicitly outlined these conditions.

If $\phi^{(p)}(a)=k_{p},p=0,1,...,\lceil\alpha\rceil$ are the initial values
of the differential equation, we define the boundary conditions as follows:

$\mathcal{B}_{p}=\tilde{\phi}^{(p)}(a)-k_{p}.$ (7)

As a result, the loss function is defined as:

$loss(\boldsymbol{\tau})=\lambda\|\mathcal{R}({\boldsymbol{\tau}})\|_{2}+\sum_{p=0}^{\lceil\alpha\rceil}\mathcal{B}_{p}^{2},$
(8)

where $\|.\|$ is the $L_{2}$ norm. To enhance the effectiveness of the
residual function in the loss function, we multiplied its value by the weight
$\lambda$.

Extending the proposed neural network method initially devised for solving
fractional delay differential equations, we now apply this approach to the
general case fractional-order system of differential-algebraic equations,
broadening the scope of its applicability. The general form of a DAE can be
expressed as:

$\displaystyle\chi\left(\tau,\phi(\tau),\phi^{\prime}(\tau),\phi^{\prime\prime}(\tau),\dots,\phi^{(\alpha)}(\tau)\right)$
$\displaystyle=0,$ $\displaystyle\zeta\left(\tau,\phi(\tau)\right)$
$\displaystyle=0,$

where

  * •

$\tau$ is the independent variable,

  * •

$\phi\in\mathbb{R}^{M}$ is the vector of state variables,

  * •

$\chi$ is a vector-valued function that describes the dynamic behavior of the
system,

  * •

$\zeta$ is a vector-valued function representing algebraic constraints.

Now, the residuals for each $j$ are formulated as follows:

$\mathcal{R}_{j}(\boldsymbol{\tau})=\chi_{j}\left(\boldsymbol{\tau},\tilde{\phi}(\boldsymbol{\tau}),\tilde{\phi}^{\prime}(\boldsymbol{\tau}),\tilde{\phi}^{\prime\prime}(\boldsymbol{\tau}),\dots,\mathcal{A}\tilde{\phi}(\boldsymbol{\tau})\right)-\zeta_{j}\left(\boldsymbol{\tau},\tilde{\phi}(\boldsymbol{\tau})\right),$
(9)

where $m$ is the number of equations in a DAE and $\boldsymbol{\tau}$ is the
input batch. Each equation is addressed with a dedicated neural network
architecture, akin to the previously introduced models for solving DDEs,
ensuring uniformity across the networks. This design integrates activation
functions consistent with those employed in our prior frameworks. In
addressing the system of DAEs, the loss function is formulated, drawing
inspiration from the previously discussed neural network architectures for
DDEs.

If $\phi_{j}(a)=k_{j}$ are the initial values of a DAE, the boundary
conditions for each $j$ are defined as follows:

$\mathcal{B}_{j}=\tilde{\phi_{j}}(a)-k_{j}.$ (10)

The loss function is defined as:

$loss(\boldsymbol{\tau})=\lambda\frac{1}{m}\sum_{i=1}^{m}\|\mathcal{R}_{i}(\tau)\|_{2}+\left(\sum_{j=1}^{m}(\mathcal{B}_{j}^{2})\right),$
(11)

where $n$ is the number of discrete points chosen for solving the equation and
$\tilde{\phi_{j}}(a)$ are the predicted values by the neural network. To
amplify the impact of the residual function within the loss function, its
value is weighted by the parameter $\lambda$.

The architecture of the network in the proposed approach can be any neural
network type. For example a pure MLP network [45], an Recurrent Neural Network
(RNN) [53] or a Legendre block [44] can be employed. Here we chose the
Legendre block which is an extension of the MLP network with better accuracy.

Fig. 1 visually depicts the graphical representation of neural block.
Formally, this block is defined by

$\mathbf{y}_{i}=v_{i}\left(\tanh(\sum_{j=1}^{k}\mathbf{w}_{j}\mathbf{x}_{j}+b)\right),\quad
i=1,2,\dots,n,$

where $v_{n}(x)$ represents $n$-th Legendre polynomial, serving as an
orthogonal activation function, defined by:

$v_{0}(x)=1$, $v_{1}(x)=x$,

$(n+1)v_{n+1}(x)=(2n+1)xv_{n}(x)-nv_{n-1}(x)$, $n\geq 1$.

The choice of $tanh(x)$ is due to the domain restriction of Legendre
polynomials. This integration is motivated by the unique mathematical
properties of Legendre polynomials, synergizing with the nonlinear
characteristics of $tanh(x)$. Crucially, the computational complexity
associated with these blocks is relatively low. Leveraging operational
matrices during the back-propagation phase facilitates the computation of the
gradient of the cost function concerning network weights. This operational
matrix is defined by:

$\frac{\mathrm{d}}{\mathrm{d}x}\mathcal{V}(x)=A\mathcal{V}(x)\;,$ (12)

where

$A=(a_{i,j})=\begin{cases}2j+1&j=i-k\\\ 0&o.w.\end{cases},$ (13)

where

$k=\begin{cases}1,3,...,m&\;if\;m\;is\;Odd\\\ 1,3,...,m-1&\;if\;m\;is\;Even\\\
\end{cases}.$ (14)

Consequently, the overall complexity is minimized, with potential intricacies
in the calculation of Legendre polynomials being negligible due to the use of
a small $n$ in the blocks. These blocks exhibit flexibility, allowing
placement anywhere in the network. Their incorporation enables the embedding
of polynomial spaces within the neural network, contributing significantly to
the network’s expressive capacity.

Figure 1: Legendre Neural Block. Figure 2: physics-informed neural networks.

In Fig. 2, we depict the architecture of the suggested physics-informed neural
network. This network consists of residual terms derived from differential
equations and initial conditions. The network processes its inputs to generate
corresponding outputs which is the dynamics of the problem. The subsequent
component is the physics-informed network, where the output $\phi$ is utilized
to calculate derivatives based on provided equations, alongside the evaluation
of boundary and initial conditions. The ultimate phase involves a feedback
mechanism, wherein these residuals are minimized through optimization with
respect to a learning rate, resulting in the adjustment of the neural
network’s parameters. To enhance model performance by parameter adjustments
and minimizing the loss function, two common optimization techniques are
employed: Adam and L-BFGS. Adam is a well-known first-order optimizer with an
adaptive learning rate, while the latter is a second-order algorithm.

######  Theorem 3.

The Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimization
algorithm updates the current value $\mathbf{x}_{k}$ using the formula:

$\mathbf{x}_{k+1}=\mathbf{x}_{k}-\alpha_{k}H_{k}\nabla f(\mathbf{x}_{k}),$

where:

  * •

$\alpha_{k}$ is the step size,

  * •

$H_{k}$ is the approximation to the inverse Hessian matrix, and

  * •

$\nabla f(\mathbf{x}_{k})$ is the gradient of the objective function at
$\mathbf{x}_{k}$.

The step size $\alpha_{k}$ is typically determined through line search
methods. The L-BFGS algorithm maintains a limited memory of past gradients and
parameter changes to approximate the inverse Hessian matrix efficiently. The
update formula for $H_{k}$ involves vectors $s_{k}$ and $y_{k}$ representing
the change in parameters and gradient, respectively:

$H_{k+1}=(I-\rho_{k}s_{k}y_{k}^{T})H_{k}(I-\rho_{k}y_{k}s_{k}^{T})+\rho_{k}s_{k}s_{k}^{T},$

where $I$ is the identity matrix and $\rho_{k}=\frac{1}{y^{T}_{k}s_{k}}$.

##  4 Numerical Examples

In order to demonstrate the efficacy of our proposed method, this section
presents a comprehensive numerical solution to a diverse set of well-known
delay differential equations and differential-algebraic equations. This
exemplary set includes linear and nonlinear, equations involving fractional
derivatives, pantograph DDEs, fractional pantograph delay differential
equations, and fractional differential-algebraic equations. To gauge the
precision of our calculations, we employ a suite of error metrics, including
Mean Absolute Error $(MAE)$, $L_{1}$ norm, $L_{2}$ norm, and $L_{\infty}$
norm, each defined as illustrated in Tab. 1. Here, $y_{i}$ represents the
exact solution, and $\tilde{y_{i}}$ signifies the predicted output generated
by our algorithm. In addition to presenting the numerical results and error
metrics, we provide graphical representations of the exact and computed
solutions. The experiments were conducted using Python programming language,
Google Colab for computational resources, PyTorch framework for neural network
development, and employing a learning rate of $0.01$ with both Adam and L-BFGS
optimizers.

Table 1: Error formula Metrics |  Formula
---|---
$L_{1}$ Norm |  $\sum_{i=1}^{n}\left|\hat{y}_{i}-y_{i}\right|$
$L_{2}$ Norm |  $\sqrt{\sum_{i=1}^{n}\left|\hat{y}_{i}-y_{i}\right|^{2}}$
$L_{\infty}$ Norm |  $\max\left|\hat{y}_{i}-y_{i}\right|,\;\;\;\;i=1,...,n$
Relative $\;L_{2}$ |  $\frac{\left\|y_{i}-\hat{y}_{i}\right\|_{2}}{\left\|y_{i}\right\|_{2}}$
MAE |  $\frac{1}{n}\sum_{i=1}^{n}\left|\hat{y}_{i}-y_{i}\right|$

######  Example 1.

We consider the following nonlinear delay differential equations:

$y^{\prime}(x)+e^{x}y^{\prime}(y-\sin(x^{2}))+\cos(x)y(x-\sin(x))=-e^{-x}-e^{\sin(x)^{2}}+\cos(x)e^{\sin(x)-x},\quad
x\in[0,1]$ (15)

with initial condition $y(0)=1$ and the exact solution $y(x)=e^{-x}$. Tab. 2
presents the computed error values for this equation, and in Fig. 3, the
difference between the actual equation solution and the predicted value, as
well as the residual value, have been plotted. The architecture comprises six
sequentially arranged trainable layers. It begins with the first layer as a
Legendre Block with 16 nodes, followed by three fully connected layers with
$\tanh(x)$ activations: the second layer with 32 nodes, the third layer with
64 nodes, and the fourth layer with 32 nodes. The final layer is another
Legendre Block with 5 nodes.

Table 2: Approximate Solutions and Error Values for Ex. 1

Name | Values
---|---
$L_{1}$ Norm | $1.88e-02$
$L_{2}$ Norm | $1.25e-03$
$L_{\infty}$ Norm | $1.43e-04$
Relative $L_{2}$ | $1.09e-04$
Mean Absolute Error | $6.27e-05$

x | y
---|---
$0.0000000000000000$ | $0.9999691812777030$
$0.1000000000000000$ | $0.9048038580658174$
$0.2000000000000000$ | $0.8187027715677497$
$0.5000000000000000$ | $0.6064693827773249$
$1.0000000000000000$ | $0.3677365537432892$

Figure 3: (a) Comparison of Results (b) Absolute Error of Ex.1.

######  Example 2.

Consider the following nonlinear delay differential equation:

$\displaystyle
y^{\prime}(x)+\sqrt{\cos(x)}y^{\prime}(\sqrt{x})+(\sin(\sqrt{x})+e^{x})y(\sin(x))=e^{x}+\sqrt{\cos(x)}e^{\sqrt{x}}+(\sin(\sqrt{x})+e^{x})e^{\sin(x)},\quad
x\in[0,1]$ (16)

with an initial condition of $y(0)=1$ and the exact solution being
$y(x)=e^{x}.$ In Tab. 3, approximate values and error values have been
obtained. Furthermore, in Fig. 4, it is evident that the residual error values
have been depicted with precision to five decimal places. The architecture is
identical to that in example 1.

Table 3: Approximate Solutions and Error Values for Ex. 2

Name | Values
---|---
$L_{1}$ Norm | $2.96e-02$
$L_{2}$ Norm | $1.94e-03$
$L_{\infty}$ Norm | $1.84e-04$
Relative $L_{2}$ | $6.26e-05$
Mean Absolute Error | $9.87e-05$

x | y
---|---
$0.0000000000000000$ | $1.0001576081160477$
$0.1000000000000000$ | $1.1053201283380465$
$0.2000000000000000$ | $1.2214592190076856$
$0.5000000000000000$ | $1.6489053832574243$
$1.0000000000000000$ | $2.7183227015316742$

Figure 4: (a) Comparison of Results (b) Absolute Error of Ex.2.

######  Example 3.

We have the following delay differential equation, known as a Pantograph
equation:

$y^{\prime}(x)=\frac{1}{2}y(qx)-y(x)-\frac{1}{2}e^{-qx},\quad
x\in[0,1],y(0)=1.$ (17)

The exact solution to this Pantograph equation is given by $y(x)=e^{-x}$ for
$0<q<1$. In Figure 5, the exact and estimated solutions are plotted for
$q=0.5$. The architecture consists of six trainable layers arranged
sequentially, commencing with the Legendre Block as the initial layer.
Subsequently, fully connected layers are incorporated, employing $\tanh(x)$
activation functions. The node configuration for each layer is specified as
follows: 16 nodes in the Legendre Block layer, 32 nodes in both the second and
fourth layers, 64 nodes in the third layer, and 5 nodes in the final layer.
(see Tab. 4)

Table 4: Approximate Solutions and Error Values for Ex. 3

Name | Values
---|---
$L_{1}$ Norm | $7.16e-04$
$L_{2}$ Norm | $4.75e-05$
$L_{\infty}$ Norm | $4.59e-06$
Relative $L_{2}$ | $4.17e-06$
Mean Absolute Error | $2.39e-06$

x | y
---|---
$0.0000000000000000$ | $1.0000022668994222$
$0.1000000000000000$ | $0.9048372111735443$
$0.2000000000000000$ | $0.8187349625675631$
$0.5000000000000000$ | $0.6065352015875470$
$1.0000000000000000$ | $0.3678792914460537$

Figure 5: (a) Comparison of Results (b) Absolute Error of Ex.3.

######  Example 4.

Consider the following delay differential equation, which is a Pantograph
equation:

$y^{\prime}(x)=\frac{1}{2}y(x)+\frac{1}{2}e^{\frac{x}{2}}y\left(\frac{x}{2}\right),\quad
0\leq x\leq R,\quad y(0)=1.$ (18)

The exact solution to this Pantograph equation is $y(x)=e^{x}$, where we take
$R=1$. Fig. 6 and Tab. 5 present a comparative analysis of results, showcasing
absolute errors and approximate solutions. The architecture is the same as the
one in example 3.

Table 5: Approximate Solutions and Error Values for Ex. 4

Name | Values
---|---
$L_{1}$ Norm | $6.51e-03$
$L_{2}$ Norm | $4.12e-04$
$L_{\infty}$ Norm | $4.07e-05$
Relative $L_{2}$ | $1.33e-05$
Mean Absolute Error | $2.17e-05$

x | y
---|---
$0.0000000000000000$ | $0.9999892555021601$
$0.1000000000000000$ | $1.1051505603136731$
$0.2000000000000000$ | $1.2213964499059211$
$0.5000000000000000$ | $1.6487125644152036$
$1.0000000000000000$ | $2.7182473537536378$

Figure 6:  (a) Comparison of Results (b) Absolute Error of Ex. 4.

######  Example 5.

We have the following fractional delay differential equation:

$D^{0.3}y(x)=y(x-1)-y(x)+1-3x+3x^{2}+\frac{2000x^{2.7}}{1071\Gamma(0.7)},\qquad
y(0)=0,$ (19)

where the analytical solution is $y(x)=x^{3}$ for $x\geq 0$. In Fig. 7, a
comparison of outcomes is presented, illustrating absolute errors.
Correspondingly, Tab. 6 provides the details of approximate solutions and
absolute errors.

Table 6: Approximate Solutions and Error Values for Ex. 5

Name | Values
---|---
$L_{1}$ Norm | $5.02e-01$
$L_{2}$ Norm | $1.96e-02$
$L_{\infty}$ Norm | $1.09e-03$
Relative $L_{2}$ | $1.63e-03$
Mean Absolute Error | $5.02e-04$

x | y
---|---
$0.0000000000000000$ | $-0.0000532769776176$
$0.1000000000000000$ | $0.0014034548414760$
$0.2000000000000000$ | $0.0083881827921961$
$0.5000000000000000$ | $0.1237044632150650$
$1.0000000000000000$ | $0.9992348206216326$

Figure 7: (a) Comparison of Results (b) Absolute Error of Ex. 5.

######  Example 6.

Consider the following fractional-order delay differential equation:

$D^{\frac{1}{2}}y(x)=y(\frac{1}{3}x)+y^{2}(x)+g(x),\quad y(0)=0,$ (20)

where
$g(x)=\frac{8}{3\sqrt{\pi}}x^{\frac{3}{2}}-\frac{2}{\sqrt{\pi}}x^{\frac{1}{2}}-x^{4}+2x^{3}-\frac{10}{9}x^{2}+\frac{1}{3}x$,
and the equation has an analytical solution given by $y(x)=x^{2}-x$. Fig. 8
displays a comparative analysis of outcomes, highlighting absolute errors.
Additionally, Tab. 7 details the corresponding approximate solutions and
absolute errors. The layer architecture in examples 5 and 6 is identical to
that in example 3.

Table 7: Approximate Solutions and Error Values for Ex. 6

Name | Values
---|---
$L_{1}$ Norm | $3.09e-02$
$L_{2}$ Norm | $1.11e-03$
$L_{\infty}$ Norm | $7.50e-05$
Relative $L_{2}$ | $1.92e-04$
Mean Absolute Error | $3.09e-05$

x | y
---|---
$0.0000000000000000$ | $0.0000065602927927$
$0.1000000000000000$ | $-0.0899926607745801$
$0.2000000000000000$ | $-0.1599556891727947$
$0.5000000000000000$ | $-0.2499623639493245$
$1.0000000000000000$ | $-0.0000014611903909$

Figure 8: (a) Comparison of Results (b) Absolute Error of Ex. 6.

In the following examples, we address a differential-algebraic equation and a
fractional-order differential-algebraic equation.

######  Example 7.

consider the following DAE:

$\begin{cases}y^{\prime}_{1}(x)=y_{1}(x)-y_{3}(x)y_{2}(x)+sin(x)+xcos(x),\\\
y^{\prime}_{2}(x)=xy_{3}(x)+y^{2}_{1}(x)+sec^{2}(x)-x^{2}(cos(x)+sin^{2}(x)),\\\
0=y_{1}(x)-y_{3}(x)+x(cos(x)-sin(x)),\quad x\in[0,1],\\\ \end{cases}$ (21)

the exact solutions are given with the initial conditions
$y_{1}(0)=y_{2}(0)=y_{3}(0)=0$: $y_{1}(x)=x\sin(x)$, $y_{2}(x)=\tan(x)$, and
$y_{3}(x)=x\cos(x)$. In Tab. 8, we present the error for this case.
Additionally, in Fig. 9, we have plotted the absolute error and the comparison
between exact and calculated solutions for $y_{1}(x)$, $y_{2}(x)$, and
$y_{3}(x)$. Three multilayer perceptron models with identical layer
architectures are utilized. Each model starts with a Legendre Block having 40
nodes, followed by a series of fully connected layers: a layer with 5 nodes
and $\tanh(x)$ activation, a layer with 32 nodes and $\tanh(x)$ activation,
another layer with 5 nodes and $\tanh(x)$ activation, followed by linear
layers with 10 and output nodes.

Table 8: Error Values for Ex. 7 Name |  $y_{1}$ Error Values |  $y_{2}$ Error Values |  $y_{3}$ Error Values
---|---|---|---
$L_{1}$ Norm | $7.45e-02$ | $2.48e-02$ | $1.84e-01$
$L_{2}$ Norm | $2.65e-03$ | $1.27e-03$ | $8.01e-03$
$L_{\infty}$ Norm | $1.09e-04$ | $1.07e-04$ | $8.15e-04$
Relative $L_{2}$ | $2.37e-04$ | $6.03e-05$ | $6.75e-04$
Mean Absolute Error | $9.31e-05$ | $3.10e-05$ | $2.30e-04$
Figure 9: Absolute error and comparison of results for Ex. 7.

######  Example 8.

In this example, we compute the solution for a linear fractional differential-
algebraic equation

$\begin{cases}D^{\frac{1}{2}}y_{1}(x)+2y_{1}(x)-\frac{\Gamma(\frac{7}{2})}{\Gamma(3)}y_{2}(x)+y_{3}(x)=2x^{\frac{5}{2}}+sin(x),\\\
D^{\frac{1}{2}}y_{2}(x)+y_{2}(x)+y_{3}(x)=\frac{\Gamma(3)}{\Gamma(\frac{5}{2})}x^{\frac{3}{2}}+x^{2}+sin(x),\\\
2y_{1}(x)+y_{2}(x)-y_{3}(x)=2x^{\frac{5}{2}}+x^{2}-sin(x),\quad x\in[0,1]\\\
\end{cases}$ (22)

with zero initial conditions, the exact solutions are as follows:
$y_{1}(x)=x^{5/2}$, $y_{2}(x)=x^{2}$, and $y_{3}(x)=sin(x)$. Within Tab. 9, we
showcase the error corresponding to this scenario. Furthermore, Fig. 10
illustrates both the absolute error and the contrast between the exact and
computed solutions for $y_{1}(x)$, $y_{2}(x)$, and $y_{3}(x)$. Three PyTorch
double-precision multilayer perceptron models exhibit a consistent layer
configuration: Legendre Block (40 nodes), fully connected layers (5 nodes with
$\tanh(x)$, 32 nodes with $\tanh(x)$, 5 nodes with $\tanh(x)$), Chebyshev
Block (10 nodes), and an output layer.

Table 9: Error Values for Ex. 8 Name |  $y_{1}$ Error Values |  $y_{2}$ Error Values |  $y_{3}$ Error Values
---|---|---|---
$L_{1}$ Norm | $2.16e-02$ | $1.57e-01$ | $1.33e-01$
$L_{2}$ Norm | $8.85e-04$ | $5.32e-03$ | $4.92e-03$
$L_{\infty}$ Norm | $1.47e-04$ | $4.00e-04$ | $3.17e-04$
Relative $L_{2}$ | $6.85e-05$ | $3.76e-04$ | $2.98e-04$
Mean Absolute Error | $2.16e-05$ | $1.57e-04$ | $1.33e-04$
Figure 10: Absolute error and comparison of results for Ex. 8. Figure 11: Time
complexity of training fPINN using the operational matrix of Caputo derivative
using 1000 Adam epochs.

##  5 Conclusion

While automatic differentiation is employed for derivative operations of
integer order, we proposed a non-uniform discretization of the fractional
Caputo operator to handle differential equations of fractional order.

In our proposed framework, distinguished by non-uniform discretization, we
observe that it does not impact the computational speed in the computation of
fractional orders. When compared to automatic differentiation for integral
orders, there is no significant decrease in speed as shown in Fig. 11. This is
because the operational matrix $\mathcal{A}$ can be computed before the
training process. In addition, computing fractional derivatives in the
training process can be done by matrix-vector multiplication.

Successfully applied to Delay Differential Equations, pantograph Delay
Differential Equations, and Differential-algebraic Equations characterized by
fractional orders, the method demonstrates efficacy in accurately solving
these equations.

Implemented on the Legendre Neural Block architecture, which integrates
Legendre polynomials into the artificial neural network structure, our
approach capitalizes on properties such as nonlinearity, computational
efficiency, and the ability to address vanishing/exploding gradient issues,
ensuring correct backpropagation through continuous differentiability.
Additionally, matrix operations streamline the computation of derivatives for
Legendre polynomials in terms of themselves.

However, the fPINN exhibits limitations, notably the absence of guaranteed
convergence in optimization error and relatively lower speed and accuracy
compared to other machine learning methods, such as LS-SVM. Future endeavors
may extend this method to address diverse types of fractional differential
equations, encompassing partial and integral differential equations.
Furthermore, exploring the approach for fractional orders greater than two
holds promise for advancing research and applications.

## References

  * [1] Bruce J West “Fractional calculus view of complexity: tomorrow’s science” CRC Press, 2016
  * [2] Vikash Pandey and Sverre Holm “Linking the fractional derivative and the Lomnitz creep law to non-Newtonian time-varying viscosity” In _Physical Review E_ 94.3 APS, 2016, pp. 032606
  * [3] Caibin Zeng and YangQuan Chen “Optimal random search, fractional dynamics and fractional calculus” In _Fractional Calculus and Applied Analysis_ 17.2 Versita, 2014, pp. 321–332
  * [4] Concepción A Monje et al. “Fractional-order systems and controls: fundamentals and applications” Springer Science & Business Media, 2010
  * [5] Jun Zhang, Zhihui Wei and Liang Xiao “Adaptive fractional-order multi-scale method for image denoising” In _Journal of Mathematical Imaging and Vision_ 43 Springer, 2012, pp. 39–49
  * [6] Guo–Cheng Wu, Dumitru Baleanu and Sheng–Da Zeng “Finite-time stability of discrete fractional delay systems: Gronwall inequality and stability criterion” In _Communications in Nonlinear Science and Numerical Simulation_ 57 Elsevier, 2018, pp. 299–308
  * [7] Hal L Smith “An introduction to delay differential equations with applications to the life sciences” springer New York, 2011
  * [8] Rodney David Driver “Ordinary and delay differential equations” Springer Science & Business Media, 2012
  * [9] Fathalla A Rihan “Delay differential equations and applications to biology” Springer, 2021
  * [10] B Van Brunt and GC Wake “A Mellin transform solution to a second-order pantograph equation with linear dispersion arising in a cell growth model” In _European Journal of Applied Mathematics_ 22.2 Cambridge University Press, 2011, pp. 151–168
  * [11] Walter G Aiello, HI Freedman and Jianhong Wu “Analysis of a model representing stage-structured population growth with state-dependent time delay” In _SIAM Journal on Applied Mathematics_ 52.3 SIAM, 1992, pp. 855–869
  * [12] John Richard Ockendon and Alan B Tayler “The dynamics of a current collection system for an electric locomotive” In _Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences_ 322.1551 The Royal Society London, 1971, pp. 447–468
  * [13] Andrzej Dzieliński et al. “Fractional Variable-Order Derivative and Difference Operators and Their Applications to Dynamical Systems Modelling” In _Fractional Dynamical Systems: Methods, Algorithms and Applications_ Springer, 2022, pp. 107–133
  * [14] Ismail T Huseynov and Nazim I Mahmudov “A class of Langevin time-delay differential equations with general fractional orders and their applications to vibration theory” In _Journal of King Saud University-Science_ 33.8 Elsevier, 2021, pp. 101596
  * [15] Laécio Carvalho de Barros et al. “The memory effect on fractional calculus: an application in the spread of COVID-19” In _Computational and Applied Mathematics_ 40 Springer, 2021, pp. 1–21
  * [16] Robert J Bauer, Gary Mo and Wojciech Krzyzanski “Solving delay differential equations in S-ADAPT by method of steps” In _Computer methods and programs in biomedicine_ 111.3 Elsevier, 2013, pp. 715–734
  * [17] Zhanwen Yang, Qi Li and Zichen Yao “A stability analysis for multi-term fractional delay differential equations with higher order” In _Chaos, Solitons & Fractals_ 167 Elsevier, 2023, pp. 112997
  * [18] M Mustafa Bahşi and Mehmet Çevik “Numerical solution of pantograph-type delay differential equations using perturbation-iteration algorithms” In _Journal of Applied Mathematics_ 2015 Hindawi, 2015
  * [19] Andi Johnson “Investigation of Network Models Finite difference Method” In _Eurasian Journal of Chemical, Medicinal and Petroleum Research_ 2.1 Amir Samimi Company (ASC), 2023, pp. 1–9
  * [20] N Senu, KC Lee, A Ahmadian and SNI Ibrahim “Numerical solution of delay differential equation using two-derivative Runge-Kutta type method with Newton interpolation” In _Alexandria Engineering Journal_ 61.8 Elsevier, 2022, pp. 5819–5835
  * [21] Pham Huu Anh Ngoc and Ky Quan Tran “On stability of solutions of stochastic delay differential equations” In _Systems & Control Letters_ 169 Elsevier, 2022, pp. 105384
  * [22] VJ Shaalini and SE Fadugba “A new multi-step method for solving delay differential equations using Lagrange interpolation” In _Journal of the Nigerian Society of Physical Sciences_ , 2021, pp. 159–164
  * [23] Nikhil Sriwastav and Amit K Barnwal “Numerical solution of Lane-Emden pantograph delay differential equation: stability and convergence analysis” In _International Journal of Mathematical Modelling and Numerical Optimisation_ 13.1 Inderscience Publishers (IEL), 2023, pp. 64–83
  * [24] Mo Faheem, Akmal Raza and Arshad Khan “Wavelet collocation methods for solving neutral delay differential equations” In _International Journal of Nonlinear Sciences and Numerical Simulation_ 23.7-8 De Gruyter, 2022, pp. 1129–1156
  * [25] Waleed Adel and Zulqurnain Sabir “Solving a new design of nonlinear second-order Lane–Emden pantograph delay differential model via Bernoulli collocation method” In _The European Physical Journal Plus_ 135.5 Springer, 2020, pp. 1–12
  * [26] H Jafari, M Mahmoudi and MH Noori Skandari “A new numerical method to solve pantograph delay differential equations with convergence analysis” In _Advances in Difference Equations_ 2021.1 SpringerOpen, 2021, pp. 1–12
  * [27] Waleed M Abd-Elhameed, José A Tenreiro Machado and Youssri H Youssri “Hypergeometric fractional derivatives formula of shifted Chebyshev polynomials: tau algorithm for a type of fractional delay differential equations” In _International Journal of Nonlinear Sciences and Numerical Simulation_ 23.7-8 De Gruyter, 2022, pp. 1253–1268
  * [28] Muyang Liu, Ioannis Dassios and Federico Milano “On the stability analysis of systems of neutral delay differential equations” In _Circuits, Systems, and Signal Processing_ 38 Springer, 2019, pp. 1639–1653
  * [29] Mostafa Abbaszadeh and Mehdi Dehghan “A Galerkin meshless reproducing kernel particle method for numerical solution of neutral delay time-space distributed-order fractional damped diffusion-wave equation” In _Applied Numerical Mathematics_ 169 Elsevier, 2021, pp. 44–63
  * [30] Fariba Takhtabnoos and Ahmad Shirzadi “A local meshless method based on the finite collocation and local integral equations method for delay PDEs” In _Engineering Analysis with Boundary Elements_ 83 Elsevier, 2017, pp. 67–73
  * [31] Ercan Celik, Mustafa Bayram and Turgut Yeloglu “Solution of differential-algebraic equations (DAEs) by adomian decomposition method” In _International Journal Pure & Applied Mathematical Sciences_ 3.1, 2006, pp. 93–100
  * [32] Ercan Çelik and Mustafa Bayram “On the numerical solution of differential-algebraic equations by Padé series” In _Applied mathematics and computation_ 137.1 Elsevier, 2003, pp. 151–160
  * [33] Fahimeh Soltanian, Mehdi Dehghan and Seyed-Mehdi Karbassi “Solution of the differential algebraic equations via homotopy perturbation method and their engineering applications” In _International Journal of Computer Mathematics_ 87.9 Taylor & Francis, 2010, pp. 1950–1974
  * [34] Mohammad Zurigat, Shaher Momani and Ahmad Alawneh “Analytical approximate solutions of systems of fractional algebraic–differential equations by homotopy analysis method” In _Computers & Mathematics with Applications_ 59.3 Elsevier, 2010, pp. 1227–1235
  * [35] Tayebeh Taheri, Alireza Afzal Aghaei and Kourosh Parand “Bridging machine learning and weighted residual methods for delay differential equations of fractional order” In _Applied Soft Computing_ Elsevier, 2023, pp. 110936
  * [36] Parisa Rahimkhani and Yadollah Ordokhani “Orthonormal Bernoulli wavelets neural network method and its application in astrophysics” In _Computational and Applied Mathematics_ 40 Springer, 2021, pp. 1–24
  * [37] Zulqurnain Sabir et al. “Intelligent computing technique for solving singular multi-pantograph delay differential equation” In _Soft Computing_ 26.14 Springer, 2022, pp. 6701–6713
  * [38] Imtiaz Khan et al. “Design of backpropagated intelligent networks for nonlinear second-order Lane–Emden pantograph delay differential systems” In _Arabian Journal for Science and Engineering_ 47.2 Springer, 2022, pp. 1197–1210
  * [39] Shagun Panghal and Manoj Kumar “Neural network method: delay and system of delay differential equations” In _Engineering with Computers_ 38.Suppl 3 Springer, 2022, pp. 2423–2432
  * [40] Yinlin Ye et al. “SLeNN-ELM: A shifted Legendre neural network method for fractional delay differential equations based on extreme learning machine” In _NHM_ 18.1, 2023, pp. 494–512
  * [41] Jingdong Zhang, Qunxi Zhu, Wei Yang and Wei Lin “SYNC: Safety-aware neural control for stabilizing stochastic delay-differential equations” In _The Eleventh International Conference on Learning Representations_ , 2023
  * [42] Lars Ruthotto and Eldad Haber “Deep neural networks motivated by partial differential equations” In _Journal of Mathematical Imaging and Vision_ 62 Springer, 2020, pp. 352–364
  * [43] Hajimohammadi, Zeinab and Baharifard, Fatemeh and Ghodsi, Ali and Parand, Kourosh “Fractional Chebyshev deep neural network (FCDNN) for solving differential models” In _Chaos, Solitons & Fractals_ 153 Elsevier, 2021, pp. 111530
  * [44] K Parand et al. “A neural network approach for solving nonlinear differential equations of Lane–Emden type” In _Engineering with Computers_ Springer, 2023, pp. 1–17
  * [45] Maziar Raissi, Paris Perdikaris and George E Karniadakis “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations” In _Journal of Computational physics_ 378 Elsevier, 2019, pp. 686–707
  * [46] Pai Peng, Jiangong Pan, Hui Xu and Xinlong Feng “RPINNs: Rectified-physics informed neural networks for solving stationary partial differential equations” In _Computers & Fluids_ 245 Elsevier, 2022, pp. 105583
  * [47] Lei Yuan, Yi-Qing Ni, Xiang-Yun Deng and Shuo Hao “A-PINN: Auxiliary physics informed neural networks for forward and inverse problems of nonlinear integro-differential equations” In _Journal of Computational Physics_ 462 Elsevier, 2022, pp. 111260
  * [48] Shengze Cai et al. “Physics-informed neural networks for heat transfer problems” In _Journal of Heat Transfer_ 143.6 American Society of Mechanical Engineers, 2021, pp. 060801
  * [49] Genming Bai, Ujjwal Koley, Siddhartha Mishra and Roberto Molinaro “Physics informed neural networks (PINNs) for approximating nonlinear dispersive PDEs” In _arXiv preprint arXiv:2104.05584_ , 2021
  * [50] Guofei Pang, Lu Lu and George Em Karniadakis “fPINNs: Fractional physics-informed neural networks” In _SIAM Journal on Scientific Computing_ 41.4 SIAM, 2019, pp. A2603–A2626
  * [51] Ehsan Kharazmi et al. “Identifiability and predictability of integer-and fractional-order epidemiological models using physics-informed neural networks” In _Nature Computational Science_ 1.11 Nature Publishing Group US New York, 2021, pp. 744–753
  * [52] Guozheng Wu, Fajie Wang and Lin Qiu “Physics-informed neural network for solving Hausdorff derivative Poisson equations” In _Fractals_ World Scientific, 2023, pp. 2340103
  * [53] Tobias Hagge, Panos Stinis, Enoch Yeung and Alexandre M Tartakovsky “Solving differential equations with unknown constitutive relations as recurrent neural networks” In _arXiv preprint arXiv:1710.02242_ , 2017