QGPINNs: A Physics-Informed Neural Network Framework for Nonlocal Differential Equations on Quantum Graphs

QGPINNs: A Physics-Informed Neural Network Framework for Nonlocal Differential Equations on Quantum Graphs Vaibhav Mehandirattaa,∗, Saket Ramchandraa a Department of Mathematics, Birla Institute of Technology and Science, Pilani, K K Birla Goa

arXiv:2608.28589v1 [cs.LG] 28 Aug 2026

Campus, Zuarinagar, Sancoale, Goa 403726, India

Abstract We propose QGPINNs, a physics-informed neural network framework developed in PyTorch [30] for the numerical solution of nonlocal differential equations on quantum graphs. The framework is designed as a general computational implementation in which the solution on each edge of the graph is approximated by a neural network, while a unified graph-based loss function enforces the governing equations together with initial, boundary, and vertex transmission conditions. In particular, the formulation incorporates standard continuity and Kirchhoff-Neumann vertex conditions and Dirichlet boundary conditions into the learning process to couple the local edge-wise neural approximations into a global solution on the graph. The framework is developed for two representative classes of nonlinear models: multi-order fractional elliptic problems and time-fractional evolution equations on quantum graphs. To improve accuracy and training stability, QGPINNs integrates several graph-adapted learning strategies, including soft and hard constraint enforcement, dynamic loss balancing, Fourier feature embeddings, and a learnable singularity-capturing feature for weakly singular solutions arising in the considered problems. The framework also extends naturally to inverse problems, including the identification of the orders of fractional operators and physical parameters from noisy observational data. We validate the accuracy, computational efficiency, and physical consistency of the proposed framework through numerical experiments on benchmark graph structures and real-world networks, including the IEEE 14-bus system and an open-channel agricultural drainage network. Keywords: Physics-Informed Neural Networks, Fractional Differential Equations, Quantum Graphs, Caputo Derivative, L1 Scheme, L2-1σ Scheme, Hard and Soft constraints, Dynamic Loss Function Weights, Fourier Embedding, Learnable Singularity Feature

1. Introduction In this paper, we establish a general framework based on the deep learning approach for solving a class of nonlocal differential equations defined on quantum graphs, and provide various examples and features of its usage and accuracy [32]. Over the past few decades, differential equation models involving non-local operators have gained significant interest in the fields of applied mathematics, physics, engineering, and computational science primarily due to their ability to capture long-range (i.e., nonlocal) interactions and anomalous diffusion processes [27], which cannot be adequately described by classical integer-order differential equations. Fractional differential operators, in particular, provide a powerful mathematical framework for modeling such nonlocal phenomena and have found applications in porous media flows, viscoelastic materials, biological systems transport in ∗ Corresponding author

Email addresses: vaibhavm@goa.bits-pilani.ac.in (Vaibhav Mehandiratta), saketrprofessional@gmail.com (Saket Ramchandra)

Figure 1: Quantum graph with the degree of each vertex labeled. Node colors indicate vertex degree.

heterogeneous media, and many other areas. However, many complex systems encountered in modern science and engineering do not evolve on standard Euclidean domains. Instead, they are naturally defined on network-like structures, such as connected urban traffic systems in a city, branched river networks, biological neural networks, gas and power transmission systems, etc. Consequently, in recent years, research focus has shifted significantly toward the mathematical modeling, statistical analysis, and numerical solution of the differential equations posed directly on network domains [2, 14, 25, 10, 11]. For such practical applications, the underlying spatial domains can be efficiently represented as compact metric graphs [6], where each edge is identified with an interval and is equipped with a metric structure. A quantum graph Γ is a metric graph endowed with additional structures, namely, assigning the differential operator, associated function spaces, and vertex conditions that describe the interaction of the edge-wise dynamics at the junctions of the network. More precisely, we define the graph Γ = (V, E) consisting of a finite |V| |E| set of nodes or vertices V = {vi }i=1 and a set of edges E = {ei }i=1 connecting the vertices, such that each edge of the graph is identified by an interval (0, ℓe ) with 0 < ℓe < ∞. That is, on each edge e, connecting two vertices vr and vs (vr → vs ), we impose a coordinate x that increases from 0 to ℓe in such a way that x = 0 corresponds to vertex vr , while x = ℓe corresponds to vertex vs . Moreover, in a more general setting, the edges of the underlying graph may not be straight connections between vertices, but can be induced by any regular curve r : [0, ℓe ] → Rd , d ∈ N, parametrized by arc-length. This geometric flexibility makes the quantum graph framework particularly suitable for modeling real-world networked structures, where the connections between vertices are naturally represented by continuous spatial curves rather than straight line segments. A visual representation of such a graph, the QGPINNs logo, is provided in Figure 1, along with the degree of each vertex, deg(v), defined as the number of edges connected to it. Several theoretical works have been investigated for classical and, more recently, for fractional models on quantum graphs, including spectral properties of differential operators, well-posedness of the associated differential problems, controllability, and qualitative behaviour of solutions; see in particular [5, 18, 3, 28, 26] and references therein. Nevertheless, the literature on numerical methods for solving differential equation models on such domains remains relatively sparse. In particular, for fractional differential equation (FDE) models on quantum graphs, traditional finite difference and finite element methods have been employed in recent years [23, 8, 15, 9]. However, these methods are restricted to discretization of each edge of the underlying graph together with the consistent treatment of the vertex conditions. Consequently, the traditional methods face computational bottlenecks and numerical instability, particularly for complex geometries and large-scale real networks such as those considered in this work. Moreover, for fractional or nonlocal models, 2

the computational difficulty further increases as the underlying operators involve memory effects and longrange interactions. These challenges motivate us to look beyond the existing traditional numerical methods for solving nonlocal differential equations on real network structures modeled by quantum graphs. Deep learning methods, particularly physics-informed neural networks (PINNs), provide an alternative computational framework for solving differential equations. In this approach, the solution is approximated by a neural network, and the governing equations, together with physical constraints such as initial and boundary conditions, are incorporated into a loss function. The key effectiveness of PINNs lies in employing automatic differentiation [4] to analytically derive the classical derivatives of neural network outputs. However, automatic differentiation cannot be directly employed for the fractional differential equation models due to the nonlocality of the fractional derivatives and the inapplicability of the classical chain rule in fractional calculus. To address this issue, fractional PINNs and nonlocal PINNs have been proposed [29, 13, 34], where one utilises conventional numerical techniques to discretize fractional operators, and these resulting approximations are then integrated into the residual loss. We refer to [36] for a nice comprehensive survey of PINNs with their different variants and applications. These existing nonlocal PINNs frameworks, however, are predominantly developed for problems defined on Euclidean domains and do not address the difficulties occur from graph topology and vertex conditions. On the other hand, recent PINN-based approaches for differential equations on quantum graphs have mainly focused on integer-order models [44, 16]. These works, in particular, shows that PINNs can be adapted to classical elliptic and evolution equations on graphs. However, to the best of our knowledge, a general PINN framework for nonlocal differential equations on quantum graphs has not yet been investigated. In particular, there is a need for a unified framework that simultaneously incorporates edge-wise nonlocal operators, graph topology, vertex transmission conditions, and neural-network-based training strategies. The main objective of the work is to describe and implement such a general computational framework. We propose QGPINNs, a unified physics-informed neural network framework for solving nonlocal differential equations on quantum graphs. The proposed framework is designed to handle both classical and nonlocal differential operators on quantum graphs while incorporating the vertex conditions that couple the edge-wise equations into a global graph problem. While we have striven to implement a general software framework for quantum graph computations using PINNs, the direction of development has been guided by two classes of nonlinear models, namely, fractional elliptic equations and time-fractional evolution equations on quantum graphs. Model 1: Multi-Order Fractional Elliptic Problem: The first class consists of elliptic-type fractional boundary value problems involving two spatial fractional operators of different orders. On each edge ei ∈ E, identified with the interval [0, ℓi ], we seek ui : [0, ℓi ] → R satisfying:   dui C γ C β + fi (x), x ∈ (0, ℓi ), i = 1, 2, . . . , E, (1.1) 0 Dx ui (x) = g ui (x), 0 Dx ui (x), dx where 1 < γ ≤ 2, 0 < β ≤ γ − 1, the function g : R3 → R encodes the possibly nonlinear reaction term, and fi is a prescribed forcing function on the edge ei . Model 2: Time-Fractional Evolution Equation: The second class consists of nonlocal evolution equations in which the classical first-order time derivative is replaced by a Caputo fractional derivative of order γ ∈ (0, 1). On each edge ei ∈ E, we seek ui : [0, ℓi ] × (0, T ] → R satisfying:   ∂ 2 ui ∂ui C γ (x, t), ui (x, t), (x, t) + fi (x, t), (x, t) ∈ (0, ℓi ) × (0, T ), i = 1, 2, . . . , E, (1.2) 0 Dt ui (x, t) = N ∂x ∂x2 subject to the initial condition

(0)

ui (x, 0) = ui (x),

x ∈ [0, ℓi ],

where N is a nonlinear differential operator and fi is the source term on edge ei . For a sufficiently smooth

3

γ function h, the nonlocal operator C 0 Dt of order γ > 0 is defined by C γ 0 Dt h(t) =

1 Γ(n − γ)

Z t

(t − τ )n−γ−1 h(n) (τ ) dτ,

(1.3)

t ∈ [0, T ],

0

where n − 1 < γ ≤ n, n ∈ N, and Γ(·) denotes the Euler gamma function. The spatial fractional operators C γ C β 0 Dx and 0 Dx in (1.1) are defined analogously. The above problems, containing both classical (local) and fractional (nonlocal) operators, are defined as a collection of equations on each edge of the graph. To define the global problem on the underlying graph Γ, we equipped the edge-wise equations with appropriate boundary conditions at the vertices, known as the vertex conditions in the context of quantum graphs. To this end, for a vertex v ∈ V, let Ev denote the set of all edges e incident to the vertex v. The vertex condition consists of Continuity condition: ui (v) = uj (v) ∀ei , ej ∈ Ev . (1.4) Kirchhoff–Neumann condition: X

(1.5)

∂ui (v) = 0,

ei ∈Ev

where ∂ui (v) denotes the directional derivative of ui at the vertex v taken in the direction pointing away from the vertex. That is, for an edge ei = [0, ℓi ], connecting the vertices vr → vs , one has ∂ui (vr ) = ∂ui (0) = u′i (0), and ui (vs ) = ∂ui (ℓi ) = −u′i (ℓi ). The above condition (1.5) is the natural extension of Neumann boundary conditions in the graph setting, as, in particular, for the boundary vertex, i.e., a vertex v with deg(v)=1, this reduces to the standard Neumann boundary condition. Moreover, this is the standard vertex condition for diffusion-type operators. For a general advection-diffusion or transport models, the derivative term in (1.5) may be replaced by the corresponding problem dependent flux which leads to a generalized Kirchhoff condition at the junction [7]. Such a flux-based junction condition has been used in the open-channel flow example in Section 5.3. Finally, the Dirichlet vertex condition is (1.6)

u(v) ≡ ui (v) = gv ,

which we impose only at boundary vertices as an alternative to the corresponding Neumann condition there. If one imposes the Dirichlet condition (1.6) for each vertex v ∈ V, then it would completely decouple the edgewise problems and essentially remove any connection between the edges adjacent to the vertex. Consequently, in this case, the topology of the graph becomes irrelevant, and the global problem on Γ reduces identically to a collection of independent problems posed locally on individual edges. v2

Neumann ∂u2 (v2 ) = 0

e2 = [0, ℓ2 ]

Interior Vertex Continuity:u1 (v0 ) = u2 (v0 ) = u3 (v0 )

v0

Kirchhoff-Neumann 3 X ∂ui (v0 ) = 0

e1 = [0, ℓ1 ]

v1 uDirichlet (v ) = 0 1

1

i=1

with ∂ui (v0 ) = u′i (0)

e3 = [0, ℓ3 ] v3

Neumann ∂u3 (v3 ) = 0

Figure 2: A star graph G with interior vertex v0 (filled) and pendant vertices v1 , v2 , v3 (shaded). At v0 , the solution must satisfy both the continuity condition (equal trace from all incident edges) and the Kirchhoff flux condition (outward derivatives sum to zero).

4

The equations (1.1) and (1.2), together with the above-defined vertex conditions, define the fractional elliptic and the time-fractional evolution equations on quantum graphs, respectively, where, in particular, for the problem (1.2), the vertex conditions (1.4)-(1.6) are imposed for every t ∈ (0, T ). In recent years, these models, particularly time-fractional diffusion equations, have gained significant interest in modeling anomalous and memory-dependent processes on network structures and have been studied from both theoretical and computational aspects; see, for instance, [24, 23, 35, 12, 41, 40]. However, existing numerical methods for such models have been predominantly based on traditional discretization techniques, which, as mentioned above, are computationally costly for large-scale networks, as it requires the discretization of both the governing equations on each edge and the vertex conditions at the graph junctions. Therefore, these models require an efficient computational framework that is capable of handling local differentiation, nonlocal memory effects, and graph-based coupling in a unified manner. In the proposed QGPINNs framework, the solution on each edge of the graph is represented by a neural network, and the total loss function is then constructed via the residuals of the governing equations, initial conditions, boundary conditions, and vertex conditions. The local derivatives in the model are computed using automatic differentiation, whereas the nonlocal fractional derivatives are approximated using suitable numerical schemes, particularly the L1 and L2 − 1σ approximation schemes [21, 1], which are then integrated into the resulting residual loss derived for our framework. This strategy allows QGPINNs to preserve the flexibility of PINNs while accurately accounting for the nonlocal structure of fractional operators. In addition, the framework describes and implements various architectural and training features designed to improve accuracy and robustness. These include enforcement of hard and soft constraints, dynamic balancing of different loss components, Fourier feature embedding for the efficient representation of oscillatory solutions, and a singularity-capturing feature for the weak singular behaviour of the solutions to fractional-order models. These features provide substantial performance gains compared with standard PINNs, without introducing significant overhead in terms of runtime or memory. Moreover, in nonlocal models, the order of the operator is a key parameter because it determines the strength and range of nonlocal interactions. Therefore, our framework also incorporates parameter estimation, particularly the estimation of the order of the fractional operators, from noisy observed data. The main contributions of this paper are summarized as follows: • We propose QGPINNs, a unified physics-informed neural network framework for solving local and nonlocal differential equations on quantum graphs. • We formulate a unified loss function that incorporates edge-wise equation residuals together with Dirichlet and Neumann boundary conditions and continuity and Kirchhoff-Neumann vertex transmission conditions. • We introduce and implement several training and architectural enhancements, including hard and soft constraints, dynamic balancing of different loss components, Fourier feature embeddings, and learnable singularity features. • We demonstrate the applicability of the proposed framework to both forward and inverse problems, including identification of the fractional-orders and physical parameters from noisy observations. • We validate the accuracy, computational efficiency and physical consistency of the framework on several typical graph structures, including real-world network such as the IEEE 14-bus system and an openchannel agricultural drainage network. The rest of the paper is organised as follows. In section 2, we discuss important mathematical preliminaries including the approximation schemes used for the considered fractional operators. Section 3 introduces the QGPINNs framework, and , in particular, defines the loss function for both forward and inverse problems, the enforcement of vertex conditions, and various optimization strategies that significantly improve performance. Section 4 discusses the implementation details of the proposed framework. In Section 5, we demonstrate the accuracy and the effectiveness of the proposed framework using numerical examples on different quantum graph structures. Finally, the conclusion of the paper is provided in Section 6. 5

2. Preliminaries In this section, we briefly present the key mathematical tools used in the QGPINNs framework. We first describe the standard physics-informed neural network formulation for differential equation models and then provide the approximation of Caputo fractional derivatives, required for constructing the residual losses associated with nonlocal operators. 2.1. Physics-informed neural networks A neural network is a mathematical model consisting of different layers of neurons, including an input layer, several hidden layers, and an output layer, which are interconnected through nonlinear activation functions. If z 0 = z ∈ Rn0 denotes the input layer, then the hidden layers are defined recursively by
z k = σ W k z k−1 + bk , k = 1, 2, . . . , L − 1, and the output layer is given by

z L = W L z L−1 + bL ∈ RnL .

Here, W k ∈ Rnk ×nk−1 and bk ∈ Rnk denote the weight matrix and bias vector at the k-th layer, respectively, nk is the number of neurons in the k th layer, and σ is a nonlinear activation function. The collection of all trainable weights and biases is denoted by θ = {W k , bk }L k=1 and the neural network approximation of the unknown solution is then written as u b(z; θ) = z L . A loss function is associated with this neural network output, which essentially quantifies the discrepancies between the network prediction and the expected outcome. The underlying network is trained using backpropagation [33], together with a suitable optimization algorithm, to minimize the resulting loss function and update the trainable parameters θ. Physics-informed neural networks (PINNs) form a well-known class of neural network methods in which the underlying physics of the system is incorporated into the training process through the loss function. In particular, the neural network is trained not only from data but also from the governing differential equations and the associated initial and boundary conditions. To introduce the basic idea, consider the general differential problem ( L[u](z) = f (z), z ∈ Ω, B[u](z) = g(z), z ∈ ∂Ω, where Ω ⊂ Rn is the spatial domain, L[·] represents the governing differential operator, and B[·] denotes the boundary operator corresponding to the boundary conditions. For time-dependent problems, an initial condition I[u](z, 0) = u0 (z), z ∈ Ω is also prescribed. j Nb r Given the neural-network approximation u b(z; θ), let {zrj }N j=1 ⊂ Ω and {zb }j=1 ⊂ ∂Ω denote the residual and boundary collocation points, respectively. Then, the residual and boundary losses are defined by Lr (θ) =

Nr
2 1 X L[û](zrj ; θ) − f (zrj ) , Nr j=1

Lb (θ) =

Nb  2 1 X B[b u](zbj ; θ) − g(zbj ) , Nb j=1

and

6

respectively. For time-dependent problems, the initial condition is enforced through L0 (θ) =

N0  2 1 X u b(z0j , 0; θ) − u0 (z0j ) , N0 j=1

0 where {z0j }N j=1 ⊂ Ω denotes the set of collocation points corresponding to the initial conditions. The total loss function for PINN is then written as

L(θ) = λr Lr (θ) + λb Lb (θ) + λ0 L0 (θ), where λr , λb , λ0 > 0 are weights used to balance the different loss components. The trainable parameters are obtained by solving θ ∗ = arg min L(θ). θ

The trained function u b(·; θ ) is then considered as the PINN approximation of the solution. During the optimization process, the derivatives of u b appearing in the residual loss are evaluated using automatic differentiation for integer-order differential operator. However, for non-local operators L, suitable numerical schemes are employed to approximate such operators, and the resulting approximations are incorporated into the residual loss. In both cases, the key objective of PINN framework is to determine the optimal parameters that minimize a loss function containing the governing equations, boundary conditions, initial conditions, and any available observational data. In the proposed framework, the above standard PINN formulation has been employed edge-wise on a quantum graph. More precisely, for each edge ei ∈ E, identified with an interval [0, ℓi ], the underlying solution ui is approximated by a neural network ûi,θ , and the differential residual is evaluated separately on that edge. The edge losses are then coupled through additional vertex-condition loss functions associated with continuity, Kirchhoff-Neumann conditions, and boundary conditions at the vertices. The complete graph-based loss formulation is presented in Section 3. ∗

2.2. Approximation of non-local operators As mentioned above, the integer-order derivatives in the considered problems are computed via automatic differentiation. However, the non-local nature of the Caputo fractional derivative requires the use of suitable approximation methods to handle such operators within the neural network optimization. While PINNs inherently provide a generalized continuous approximate solution over the entire computational domain, the training of the model is executed on a finite set of collocation points. Therefore, the choice and distribution of these data points play an important role in the accuracy and convergence efficiency of the PINN framework. Furthermore, the considered nonlocal problems do not possess a smooth solution throughout the closed domain and often exhibit weak singularities [37]. Such singularities present a significant challenge for a PINN to approximate the solution when using a standard uniform discretization mesh. To tackle this, nonuniform meshes, in particular, graded meshes, are commonly employed. Specifically, let N be a positive integer, then the temporal graded mesh is defined by  n r tn = T , n = 0, 1, . . . , N, (2.1) N where r ≥ 1 is the mesh grading parameter chosen by the user. The graded mesh on each spatial edge ei = [0, ℓi ] is defined analogously. In particular, for r = 1, the graded mesh reduces to a uniform mesh. We shall employ two different discretizations, namely, the L1 scheme and L2 − 1σ of scheme, for the approximation for the Caputo fractional derivatives given in (1.1) and (1.2). In the following, we describe the approximation for the temporal Caputo fractional derivative; the corresponding spatial fractional derivative can be defined analogously on each edge.

7

uniform spacing

r=1 r=2 r=4

0

T t

dense near t = 0

Figure 3: Graded mesh with N = 12 points on time interval [0, T ]. Higher r clusters points near t = 0, providing finer resolution of the initial singularity region.

2.2.1. The L1 Scheme In view of the temporal graded mesh (2.1) defined over the interval (0, T ], the Caputo fractional derivative (1.3) of order γ ∈ (0, 1) at a mesh point t = tn can be written as C γ 0 Dt f (x, tn ) =

n−1 Z tk+1

X 1 Γ(1 − γ)

k=0

(tn − s)−γ

tk

∂f (x, s) ds ∂τ

Now, for n ≥ 1, the well-known L1 formula for approximating Caputo fractional derivative is expressed as [21] C γ 0 Dt f (x, tn ) ≈

=

n−1

Z tk+1

k=0

tk

X f (x, tk+1 ) − f (x, tk ) 1 Γ(1 − γ) τk+1 1 Γ(2 − γ)

n−1 X k=0

(tn − s)−γ ds

f (x, tk+1 ) − f (x, tk ) [(tn − tk )1−γ − (tn − tk+1 )1−γ ] τk+1

(2.2)

n−1

=

X 1 [wn,1 f (x, tn ) − wn,n f (x, 0) − f (x, tn−k )(wn,k − wn,k+1 )], Γ(2 − γ) k=1

where τk = tk − tk−1 and wn,k :=

(tn − tn−k )1−γ − (tn − tn−k+1 )1−γ , for k = 1, 2, . . . , n. τn−k+1

(2.3)

Based on the error analysis provided in [37], the above L1 approximation (2.2) on a graded mesh satisfies the convergence estimate O(τ min(rγ,2−γ) ), where τ = max τn . Consequently, the optimal grading mesh 1≤n≤N

parameter for achieving the maximum rate O(τ 2−γ ) is r = 2−γ γ . 2.2.2. The L2-1σ Scheme We now present a higher-order approximation scheme developed by Alikhanov [1], namely, the L2-1σ scheme on nonuniform mesh. The L2-1σ scheme improves on the L1 scheme by employing piecewise quadratic interpolation. To achieve higher accuracy, the fractional derivative is evaluated at the shifted point tn−σ = σtn−1 + (1 − σ)tn ,

σ =1−

First, we define the following local kernel components:

8

γ , 2

n = 1, 2, . . . , N.

An0 =

1 τn Γ(1 − γ)

Ann−k =

Z tn−σ

(tn−σ − s)−γ ds =

tn−1 Z tk

1 τk Γ(1 − γ)

(tn−σ − s)−γ ds =

tk−1

2

(tn−σ − tn−1 )1−γ (1 − σ)1−γ τn−γ = , τn Γ(2 − γ) Γ(2 − γ) (tn−σ − tk−1 )1−γ − (tn−σ − tk )1−γ , τk Γ(2 − γ)

Z tk

(tn−σ − s)−γ (s − tk−1/2 ) ds τk (tk+1 − tk−1 )Γ(1 − γ) tk−1   τk 2 (tn−σ − tk−1 ))2−γ − (tn−σ − tk )2−γ − ((tn−σ − tk−1 )1−γ + (tn−σ − tk )1−γ ) + = τk (τk + τk+1 )Γ(2 − γ) 2 2−γ

n Bn−k =

where tk−1/2 = (tk−1 + tk )/2 and k = 1, 2, . . . , n − 1. Based on the above components, we construct the (n,γ)

discrete kernels Ck

(1,γ)

. For n = 1, we set C0 = A10 and for 2 ≤ n ≤ N , we define  n n  k = 1, An−1 − Bn−1 , (n,γ) n n n Cn−k = An−k + ρk−1 Bn−k+1 − Bn−k , 2 ≤ k ≤ n − 1,   n A0 + ρn−1 B1n , k = n,

(2.4)

where ρk = τk /τk+1 . Therefore, the final discrete approximation for the Caputo fractional derivative of order γ ∈ (0, 1) at the shifted point tn−σ is: (n,γ) C γ f (x, tn ) − 0 Dt f (x, tn−σ ) ≈ C0

n−1 X

(n,γ) (n,γ)
(n,γ) Cn−k−1 − Cn−k f (x, tk ) − Cn−1 f (x, t0 ),

1 ≤ n ≤ N.

(2.5)

k=1

The above L2 − 1σ scheme on nonuniform mesh achieves the convergence rate O(τ min(rγ,2) ) [20], under suitable regularity assumptions. Therefore, the optimal grading exponent for obtaining the rate O(N −2 ) is r = 2/α. α Remark 2.1. In view of Definition (1.3), the Caputo fractional derivative C 0 Dx f of order 1 < α ≤ 2 can be expressed as
C α C α−1 f ′ (x) . 0 Dx f (x) = 0 Dx

Thus, in view of above relation, the approximation of the fractional derivative in (1.1) can be obtained by first evaluating the classical derivative f ′ (x) through automatic differentiation, and then applying the L1 or L2 − 1σ approximation given in (2.2) and (2.5), respectively, to the resulting quantity with fractional order γ = α − 1. 2.2.3. Unified Matrix Representation Both the L1 and L2-1σ discretizations of the Caputo fractional derivative can be efficiently implemented  T as a matrix-vector product. To this end, let uθ = uθ (ζ0 ), uθ (ζ1 ), . . . , uθ (ζN ) denote the vector of neural network outputs evaluated at the chosen mesh or collocation points ζi , i = 0, 1, . . . , N . Then, the discrete Caputo operator based on the L1 scheme can be written in the operator form C γ 0 Dζ uθ ≈ Wγ uθ ,

where Wγ ∈ R(N +1)×(N +1) is the lower triangular matrix given by

9



0 w1,1 w2,2 − w2,1 w3,3 − w3,2 .. .

0 0 w2,1 w3,2 − w3,1 .. .

0 0 0 w3,1 .. .

··· ··· ··· ··· .. .

0 0 0 0 .. .

wN,N − wN,N −1

wN,N −1 − wN,N −2

···

wN,2 − wN,1

wN,1

0  −w1,1   −w2,2 1  Wγ =  Γ(2 − γ)  −w3,3  ..  . −wN,N

     .   

(2.6)

Here, the weights wn,k , k = 1, 2, . . . , n, are defined analogously to (2.3) with ti replaced by ζi . Similarly, for the L2 − 1σ scheme (2.5), we have C γ 0 Dζ uθ ≈ Aγ uθ , where 0  −C (1,γ) 0   −C (2,γ)  1 Aγ =   −C2(3,γ)  ..   . 

(N,γ)

−CN −1

0 (1,γ) C0 (2,γ) (2,γ) C1 − C0 (3,γ) (3,γ) C2 − C1

.. . (N,γ) (N,γ) CN −1 − CN −2

0 0 (2,γ) C0 (3,γ) (3,γ) C1 − C0

.. . (N,γ) (N,γ) CN −2 − CN −3

··· ··· ··· ··· .. .

0 0 0 (3,γ)

C0

.. . ···

(N,γ)

C1

0 0 0 0 .. . (N,γ)

− C0

(N,γ)

     ,    

(2.7)

C0

(n,γ)

with coefficients Ck are given by (2.4), again interpreted with respect to the mesh points ζi . For a fixed graded mesh parameter, the above matrices can be precomputed once outside the PINN architecture. Therefore, during the training process, the matrix-vector products Wγ uθ or Aγ uθ are directly used to approximate the corresponding temporal or spatial fractional derivatives. This allows the nonlocal fractional operators to be incorporated efficiently into the residual loss function. 3. Loss Formulations and Learning Framework In this section, we describe the framework used to obtain the approximate solutions of the proposed nonlocal differential equations on quantum graphs. The framework consists of loss formulations on the underlying graph together with suitable training strategies for enforcing governing equations, vertex transmission conditions, boundary conditions, and initial conditions. Moreover, additional features such as hard and soft constraint enforcement, Fourier feature embeddings, and learnable singularity features are incorporated to improve the accuracy and stability of the proposed PINN solver. 3.1. Unified Loss Formulation As described in Section 2.1, the primary distinction between a conventional neural network and a PhysicsInformed Neural Network lies in the formulation of the loss function. In this section, we derive a rigorous, unified residual loss formulation for the proposed problems on a quantum graph Γ = (V, E). To this end, let u bi,θ denote the neural network approximation of the solution on an edge ei ∈ E, where θ represents the trainable parameters of the neural network. For elliptic problems (1.1), the approximation is of the form u bi,θ : [0, ℓi ] → R, whereas for the time-dependent problems (1.2), it maps onto the spatiotemporal domain: u bi,θ : [0, ℓi ] × [0, T ] → R.. The total loss function Ltotal,θ over the quantum graph Γ is constructed as the sum of edge-wise residual loss and the loss associated with the vertex conditions. That is, X Ltotal,θ := Lei ,θ + λLnode,θ , (3.1) ei ∈E

where λ > 0 is an adaptive parameter that balances the relative weight of the different residuals. The term Lei ,θ and Lnode,θ in (3.1) corresponds to the edge loss and node loss, respectively, which essentially 10

compute the discrepancies between the target data and PINN predictions along the edges and the nodes of the quantum graph. We shall now provide the mathematical expression for each term. x For the training process, we define spatial collocation points {xi,j }N j=1 ⊂ (0, ℓi ) on each edge ei ∈ E, t i = 1, 2, . . . , |E|, and temporal points {tk }N k=1 ⊂ (0, T ]. The edge loss function associated to elliptic and evolutionary problem on each edge is formulated as N

Lelliptic = ei ,θ Levolution = ei ,θ

x 1 X 2 N [ûi,θ ](xi,j ) − fi (xi,j ) , Nx j=1

Nx Nt X Nx
2 1 X 1 X 2 (0) N [ûi,θ ](xi,j , tk ) − fi (xi,j , tk ) + ûi (xi,j , 0) − ui (xi,j ) , i = 1, 2, . . . , |E|, Nx Nt N x j=1 j=1 k=1

respectively, where fi represents the source terms on edge ei and N denotes the governing differential operator in the considered problems containing both local and non-local derivatives of the approximate solution ûi,θ . Now, to define the loss function for the nodes or vertices, we partition the set of all vertices V into a set of boundary (degree-one) vertices Ω ⊂ V, and a set of interior junction vertices V \ Ω. The node loss term aggregates the continuity and Kirchhoff-Neumann conditions in junction (interior) nodes, and the conditions on boundary nodes: X X Lelliptic Lbc,v + (Lcont,v + LKN,v ) . (3.2) node,θ = v∈Ω

v∈V\Ω

Here, for boundary vertices v ∈ Ω subject to a Dirichlet condition gv or a Neumann condition hv , the boundary penalty is
2
2 Lbc,v = ûi,θ (xv ) − gv or û′i,θ (xv ) − hv , where the index i corresponds to the single edge ei , incident on v, and based on the convention adopted xv ∈ {0, ℓi } is the local coordinate for the vertex. At an interior junction v ∈ V \ Ω of degree dv , let {ej1 , ej2 , . . . , ejdv } denotes the set of all edges incident to v. Then, the continuity condition is enforced by penalizing pairwise differences: X
2 Lcont,v = ûja ,θ (xv ) − ûjb ,θ (xv ) . 1≤a<b≤dv

The loss function corresponds to the Kirchhoff-Neumann condition at the junction node v reads LKN,v =

dv X

!2 sk û′jk ,θ (xv )

,

k=1

where sk ∈ {−1, +1} depending on whether the direction of edge ejk is incoming or outgoing with respect to the junction v. For evolutionary problems, we define the node loss function by taking the temporal average of the above vertex residuals over the discrete time levels tk :   Nt X X X 1  (3.3) Levolution = Lbc,v (tk ) + (Lcont,v (tk ) + LKN,v (tk )) , node,θ Nt k=1

v∈Ω

v∈V\Ω

11

where each constituent is explicitly defined as 
2
2  or û′i,θ (xv , tk ) − hv (tk ) , for v ∈ Ω, Lbc,v (tk ) = ûi,θ (xv , tk ) − gv (tk ) 
2 P Lcont,v (tk ) = 1≤a<b≤dv ûja ,θ (xv , tk ) − ûjb ,θ (xv , tk ) , for v ∈ V \ Ω,    2  LKN,v (tk ) = Pdv sm û′ , for v ∈ V \ Ω. jm ,θ (xv , tk ) m=1 Finally, we provide the loss formulation for inverse problems where the objective is to estimate governing system parameters from sparse and noisy observations. This is achieved by augmenting the above-defined loss function with an explicit empirical data residual term, defined as Ldata,i =

Nd 1 X 2 ûi,θ (xi,m ) − udata (xi,m ) , i = 1, 2, . . . , |E|, Nd m=1

(3.4)

Nd where {xi,m }m=1 denotes the set of discrete observation data points on edge ei and udata denotes the observed values. The total per-edge loss is then directly modified by including the data-driven loss term (3.4), which allows the total loss function for the inverse problems to retain the structure of the forward problem formulation, that is X Linverse := Li + λLnode,θ total i∈E

= λd

X

Ldata,i +

i∈E

X

Lei ,θ + λLnode,θ .

(3.5)

i∈E

Note that λd represents the adaptive weight of the data-driven loss term and its computation has been done identically to the adaptive node weight λ. During the training phase, the unknown physical parameters are treated as learnable variables and are optimized along with the neural network parameters using back-propagation. Notably, for estimating the fractional order of the considered nonlocal operators, the corresponding operational matrices derived in (2.6) and (2.7) must be continuously differentiable with respect to the fractional order parameter. This requirement is addressed by initializing the discrete matrix coefficients as differentiable analytical tensors[29]. 3.2. Hard and Soft Constraint Enforcement The proposed framework supports two modes for enforcing conditions at the graph vertices. The first mode corresponds to soft constraints, in which the residuals associated with the vertex conditions are incorporated into the loss function and penalized during training. The second mode is hard constraints, introduced in the context of neural-network-based differential equation solvers by Lagaris et al. [17], which modify the network architecture such that the prescribed boundary conditions are satisfied by construction. Therefore, in the latter case, the corresponding boundary residuals need not be included in the loss functional. Soft constraints utilize the raw neural network prediction ûi,θ without any modification and treat it directly as the approximate solution on each edge ei . The boundary conditions at the degree-one vertices, together with the continuity and Kirchhoff-Neumann conditions at the interior vertices, are imposed through the node loss term as described in (3.2) and (3.3). This approach is the default operational mode in our framework since it applies uniformly to different types of graph constraints. The hard constraint mode is formulated only for Dirichlet conditions imposed at boundary vertices. Thus, rather than using raw network predictions, the trial solution is modified by means of suitable lifting functions to ensure that prescribed boundary values are satisfied identically. Consequently, the corresponding Dirichlet residuals are excluded from Lbc,v in this mode. The modification for hard constraint mode in our framework is obtained by adapting the boundary-distance concept of Sukumar and Srivastava [38] to the localized graph setting through a onesided lifting. To this end, let Fi,θ denote the raw neural network output on an edge ei , and let xv ∈ {0, ℓi } denotes the local coordinate for a boundary vertex v ∈ Ω incident with edge ei . Since the other vertex of the same 12

edge connects to an interior node, the hard constraint enforcement in the graph setting is one-sided, that is, it leaves the other endpoint of the edge unconstrained. For an elliptic problem, let gv be the prescribed Dirichlet value at a boundary vertex xv . Then, the constrained network output is defined as     |x − (ℓi − xv )| |x − xv | ûi,θ (x) = gv + Fi,θ (x), i = 1, 2, . . . , |E|. (3.6) ℓi ℓi Upon evaluating the above expression, it is evident that at the boundary vertex x = xv , the ansatz gives ûi,θ (xv ) = gv , which ensures that Dirichlet conditions are satisfied by construction. Moreover, at the other vertex, i.e., the interior junction where |x − xv | = ℓi , the output ûi,θ reduces to the raw neural network prediction Fi,θ , which remains unconstrained. In case of evolutionary problems with time-varying Dirichlet boundaries gv (t), the construction must accommodate the initial profile u0i (x) and boundary data simultaneously. The exact mathematical ansatz chosen depends on the compatibility of the initial and boundary conditions at t = 0. If gv (0) = u0i (xv ), then the spatiotemporal extension naturally takes the form:     |x − xv | |x − (ℓi − xv )| gv (t) + Fi,θ (x, t). (3.7) ûi,θ (x, t) = ℓi ℓi When the conditions are incompatible, i.e., gv (0) ̸= u0i (xv ), then the initial and boundary data cannot be jointly satisfied by a single continuous hard constraint lifting function. In this case, we employ a timedecaying weight function: w(t) = e−κt ,

κ = 20,

which provides the corresponding spatiotemporal ansatz     
|x − (ℓi − xv )| |x − xv | ûi,θ (x, t) = w(t)u0i (x) + 1 − w(t) gv (t) + Fi,θ (x, t) . ℓi ℓi Under this formulation, at t = 0 the weight w(0) = 1 recovers the initial condition u0i (x) exactly. As t increases, w(t) decays to zero and the ansatz approaches the boundary-constrained lifting (3.7). Thus, in this case, the above construction provides a transition from the initial data to the boundary constrained representation, rather than enforcing both constraints exactly at t = 0. Remark 3.1. In our framework, the hard constraint formulation is restricted only to Dirichlet conditions at boundary vertices. The Neumann conditions at boundary vertices, along with continuity and KirchhoffNeumann conditions at junction nodes, are more naturally imposed as soft constraints. Although hardconstrained formulation for Neumann data at boundary vertices can be constructed using higher-order lifting functions in the spirit of (3.6) and (3.7), such constructions are less natural on metric graphs because they may affect the behavior of the function at the opposite endpoint, i.e., at an interior junction. This can interfere with the simultaneous enforcement of continuity and Kirchhoff–Neumann vertex conditions across incident edges. Therefore, these conditions are retained as soft constraints to maintain a flexible and stable optimization formulation. The advantage of hard constraints is that they enforce prescribed boundary data exactly and may provide better accuracy. However, their performance depends on the choice of lifting function and is highly sensitive to the distribution of collocation points, especially for nonlocal fractional models where graded meshes are used to resolve weak singularities. Soft constraints, on the other hand, provide a more flexible and highly robust formulation, since the same loss formulation can be used for different types of boundary conditions, graph topologies, and mesh distributions. We compare the hard and soft constraints and demonstrate the limitations of hard constraints empirically through an example. To quantify the error in the estimation of the approximate solution provided by our

13

framework, we define the relative discrete L2 error for elliptic problems on the graph Γ as ∥b uθ − u∥L2 (Γ) ≈ ∥u∥L2 (Γ)

P|E| PNi −1

2

i=1 m=0 |ûi,θ (xi,m ) − ui (xi,m )| P|E| PNi −1 2 i=1 m=0 |ui (xi,m )|

!1/2 ,

(3.8)

where xi,m represents the discrete spatial collocation points along edge ei , and Ni represents the total number of spatial validation points sampled per edge, ûi,θ is the predicted solution by QGPINNs on edge ei and ui represents the exact solution on that edge. To ensure edges of varying lengths are considered equally, the number of validation points per unit length is uniform across all edges. This ensures that longer edges contain proportionally more validation points. Moreover, the Ni validation points are randomly distributed throughout the spatial domain ei = [0, ℓi ] to ensure that the L2 errors are not biased towards the points used in network training. Example 3.1. We consider the fractional elliptic problem C α C β 0 Dx ui (x) + 0 Dx ui (x) = fi ,

i = 1, 2, 3,

on a metric star graph with 3 edges (outward pointing) as given in Figure 2 with vertex conditions (1.4)(1.5) at node v0 and Dirichlet vertex conditions (1.6) at boundary nodes vi , i = 1, 2, 3. We take α = 1.5 and β = 0.5, and the source terms fi are chosen such that u1 (x) = x3 − 2x, u2 (x) = x1.2 + x and u3 (x) = x2 + x forms the exact solution of the problem. Table 1: Effect of the mesh grading factor r on the relative discrete L2 error (3.8) for soft and hard constraints on a metric star graph. Mesh Grading Factor r=1 r=2 r=3

Soft Constraints

Hard Constraints

1.78 × 10−3

1.76 × 10−3 3.35 × 10−3 1.01 × 10−3

2.17 × 10−3 1.66 × 10−4

Table 1 compares the global relative L2 error (3.8) obtained using hard and soft constraint enforcement for different mesh grading factors. For the uniform mesh r = 1, both approaches produce nearly identical errors. However, as the grading factor increases, the accuracy of the hard constraints framework deteriorates. In particular, for high grading factor r = 3, the soft-constraint formulation gives much smaller error than the hard-constraint formulation. This indicates that, although hard constraints can be effective for simple geometries and weakly graded meshes, their accuracy may deteriorate when the mesh is strongly redistributed to capture weak singular behaviour in nonlocal models. The soft-constraint formulation is therefore more stable with respect to mesh grading in this example and is better suited as the default enforcement strategy in the proposed framework. We also emphasise that the above comparison can only be interpreted as an empirical rather than a theoretical guarantee. In general, hard constraints may still be advantageous when the boundary data are simple and the lifting function is well adapted to the structure of the graph. However, for nonlocal problems on metric graphs, where boundary, continuity, and Kirchhoff–Neumann conditions must be handled simultaneously, the soft formulation provides a more uniform and robust implementation. 3.3. Dynamic Balancing of Loss Function Weight The choice of the weighting parameter λ in the total loss function (3.1) plays an important role in balancing the loss residuals and significantly affects the efficacy of the PINN. For small λ, the neural network may reduce the governing equation residual while violating the vertex constraints, and for large λ, the 14

optimization process may overemphasize the residuals associated with boundary and vertex conditions at the expense of the governing equation. Therefore, a fixed chosen value of the weight parameter λ, in general, may not be equally effective during the training process. To mitigate the task of manually tuning the parameter, our framework implements methods to dynamically update the value of λ during training. We describe two distinct adaptive strategies: boundary data matching strategy and gradient pathology balancing, which are then integrated into a mixed dual update rule. Let λ(n) denotes the value of the weight parameter at the current training epoch n, where n = 0, 1, . . . , E, and E is the total number of training epochs, such that the normalized training progress is denoted as (n) n . Each adaptive strategy produces a target weight λtarget . However, a noisy loss floor may cause pn = E fluctuations in the target weight; thus, to avoid abrupt changes in λ, we use a log-space exponential moving average to define a smoothed target weight   (n) (n−1)
(n)
λtarget = exp β ln λtarget + (1 − β) ln λtarget where β ∈ (0, 1) is a user-defined smoothing coefficient with a default value of 0.95. Given the smoothed target weight, we introduce a user-defined update rate ρ (default set to ρ = 0.1). To prevent unbounded growth of the loss weight λ, we apply strict bounds and express the updated value of λ as   
(n)
λ(n+1) = Π[λmin ,λmax ] exp (1 − ρ) ln λ(n) + ρ ln λtarget where Π[a,b] (x) = max(a, min(x, b)) is the standard clipping operator mapped to [λmin , λmax ]. The above mechanism ensures that the relative adaptive weights of the loss function remain stable. We shall now (n) provide the details to compute the target weight λtarget from both the strategies. Boundary Data Matching Mechanism (BDMM) The boundary data matching strategy balances the objective function components by driving the loss associated with graph vertices towards a self calibrating target function. The aim is to prevent the vertex constraint contribution from becoming either negligible or excessively dominant during training. Let X X X Ledge,θ = Lei ,θ , Lnode,θ = Lbc,v + (Lcont,v + LKN,v ) ei ∈E

v∈Ω

v∈V\Ω

denote the total edge residual loss, and the total node residual loss, respectively, as defined in Section 3.1. Then, at epoch n, we define the instantaneous loss ratio as (n)

(n)

λtarget =

Lnode,θ (n)

,

max(Ledge,θ , ε)

where ε > 0 is a small stabilization parameter which has been used to avoid division by zero. The target weight is then applied to the smooth update rule described above. Gradient Pathology Balancing The second strategy is adapted from the gradient balancing methodology proposed by Wang et al. [42]. In PINNs, different components of the loss may generate gradients of substantially different magnitudes, which may lead to slower convergence or can cause the optimization to prioritize certain components of the loss over the others. To reduce this effect, we dynamically adapt λ using the relative magnitudes of the gradients associated with the edge residual and the node residual across the network layers. Let θ the set of all trainable parameters in the neural network. At epoch n, we define the gradient-based target weight by

15

r

(n)

P

p∈θ

(n)

λtarget =

∇p Ledge,θ

r max

P

p∈θ

(n) ∇p Lnode,θ

2 2

!

2

,ε 2

where ∥·∥2 denotes the standard L2 norm. The numerator captures the global gradient contribution from the edge residual, whereas the denominator provides the scale for the gradients generated by the node residual. This target weight is again used in the same smoothed update rule for λ. Mixed Dual Strategy The above-defined adaptive strategies exhibit complementary regularizing characteristics. The BDMM monitors the relative size of the boundary residual, whereas the gradient pathology balancing accounts for the relative magnitude of the optimization gradients governed by different loss components. By leveraging their complementary properties, we introduce a four-phase mixed dual strategy that will be used in the proposed framework. To this end, let E be the total number of training epochs and let τ ∈ (0, 1) denote the epoch fraction after which the adaptive update is frozen. 1. Phase 0 (0 ≤ n < 0.02E): The value of the weight parameter λ is kept fixed at its initial value, and no adaptive update is applied during this stage. 2. Phase 1 (0.02E ≤ n < 0.1E): The BDMM controller is active, and the target weight is chosen as the BDMM target. 3. Phase 2 (0.1E ≤ n < τ E): The gradient pathology controller is active, and the target weight is chosen as the gradient-based target. 4. Phase 3 (τ E ≤ n ≤ E): Both controllers are inactive, and therefore λ remains fixed during the final part of the first-order optimization stage. Notably, after the completion of the first-order Adam optimization stage, the balancing parameter λ is fixed to the median of its values over the final Nmedian Adam iterations. This value is then remains static throughout the subsequent L-BFGS refinement stage. In our implementation, we choose Nmedian = 1000. We now demonstrate the comparative analysis of these strategies through a numerical example. For evolutionary problems, the relative discrete L2 error on the graph Γ is defined as ∥b uθ − u∥L2 (0,T ;Γ) ≈ ∥u∥L2 (0,T ;Γ)

P|E| PNi −1 PNt −1 i=1

2 m=0 k=0 |ûi,θ (xi,m , tk ) − ui (xi,m , tk )| P|E| PNi −1 PNt −1 2 k=0 |ui (xi,m , tk )| i=1 m=0

!1/2 ,

(3.9)

where tk denotes the k-th temporal validation point distributed along the temporal domain, and Nt represents the total number of temporal validation points sampled across the temporal domain [0, T ]. The other notations are the same as in (3.8). Example 3.2. We consider the following time-fractional evolution equation on the tree-like metric graph given in Figure 4: ∂ui ∂ C γ (x, t) − ν 0 Dt ui (x, t) + ui (x, t) ∂x

2

ui (x, t) + cui (x, t) = fi (x, t), ∂x2

x ∈ (0, 1),

t ∈ (0, 1],

(3.10)

where the fractional order is chosen as γ = 0.5, the kinematic viscosity (diffusion coefficient) is ν = 0.5, and the linear reaction coefficient is c = 5.0. The initial condition is ui (x, 0) = 0, for x ∈ [0, 1] and homogeneous

16

Dirichlet boundary conditions are applied at the boundary nodes v0 , v2 , v4 . The source terms fi (x, t) are chosen such that exact solution on the edge i is: ui (x, t) = Ai sin(kx)t2 ,

(3.11)

with k = 4π, and A = [2.0, 1.0, 1.0, 1.0],

where the coefficients Ai have been obtained by enforcing the continuity (1.4) and Kirchhoff-Neumann (1.5) constraints at the node v1 and v3 . Strategy

v2 e1 v0

e0

λ=1 λ = 30 Gradient Pathology BDMM Dual

v1

e2

v3

e3

v4

Relative L2 Error 1.26 × 10−2 5.40 × 10−3 4.86 × 10−3 1.19 × 10−2 4.36 × 10−3

Table 2: Relative L2 error (3.9) for Example 3.2 obtained using different weighting strategies.

Figure 4: A four-edge cascading tree metric graph.

For this comparison, we utilize a configuration of 4 hidden layers with 128 neurons each. The Adam optimization phase is run for 10000 epochs followed by a L-BFGS refinement phase. All factional derivatives were computed using the L2 − 1σ scheme. Table 2 compares the relative L2 error (3.9) obtained using a fixed value of λ and the proposed adaptive balancing strategies for fractional order γ = 0.5. The adaptive strategies outperformed the fixed λ = 1 and λ = 30 variants without requiring any problem-specific calibration. Figure 5 shows the evolution of the weight λ during training for the adaptive strategies. These results support the use of dynamic balancing in the proposed framework, since it reduces sensitivity to a fixed choice of the loss-weighting parameter.

Figure 5: Dynamic evolution of the balancing parameter λ during training for different adaptive strategies.

Therefore, the above numerical results show that dynamic balancing of λ can improve the performance of the proposed PINN framework as compared to a fixed manually chosen loss weight, where the proposed dual strategy provided the most accurate result for this example. Nevertheless, the relative performance of the individual adaptive strategies may depend on the structure of the problem, the type of boundary and vertex conditions, and the relative strength of the edge and node residuals. For more problems and discussion, we refer to [32]. 17

3.4. Spectral Bias and Fourier Feature Embedding Rather than utilizing only a conventional multi-layer perceptron (MLP) architecture, as described in Section 2.1, our framework also incorporates random Fourier feature embeddings. The underlying feature is implemented to mitigate the phenomenon of spectral bias, which is inherent in standard neural networks, as described by Rahaman et al. [31]. Spectral bias describes the tendency of MLPs to preferentially learn low-frequency components of a target function over the high-frequency components. Consequently, standard PINNs may fail to accurately resolve oscillatory solutions or sharp variations during the early stages of training. To address this issue, we project the raw spatial or spatio-temporal input coordinates into a higherdimensional feature space using a Fourier feature mapping adapted from Tancik et al. [39]. Let x ∈ Rd denote the input coordinate vector, which for elliptic problems represents the spatial coordinate, and for evolutionary problems includes both spatial and temporal coordinates. The Fourier feature map is defined by  ⊤ γ(x) = cos(2πBx), sin(2πBx) where B ∈ Rm×d represents a static frequency matrix whose entries are sampled once prior to training from a Gaussian distribution N (0, σ 2 I). The parameter σ > 0 controls the range of frequencies introduced into the input representation. In the Fourier feature PINN framework, the embedded coordinate γ(x) is passed through the MLP instead of the raw coordinate x. This allows the network to represent oscillatory components more effectively and reduces the low-frequency preference of standard MLP-based PINNs. From the neural tangent kernel (NTK) perspective, Fourier feature mappings modify the effective kernel associated with the MLP. As shown by Tancik et al. [39], this transformation can convert the effective NTK into a stationary kernel with a tunable bandwidth, which improves the learning of high-frequency components. This motivates the use of Fourier embeddings in our framework for nonlocal problems on metric graphs, where the solution may exhibit oscillatory behaviour, weak singularities, or multiscale spatial structures. Standard PINN

Fourier feature PINN

Input x

Input x

raw coordinates

raw coordinates

MLP layers

Fourier embedding γ(x)

tanh activations

[cos(2πBx), sin(2πBx)]⊤

where B ∼ N (0, σ 2 I) fixed, sampled once before training

Spectral bias weak high-frequency learning

MLP layers linear combination of features

ûθ (x) misses high frequencies

Tunable NTK bandwidth improved high-frequency learning

ûθ (x) better resolves high frequencies Figure 6: Framework comparison between a standard PINN and a Fourier feature PINN. The Fourier embedding γ(x) modifies the input representation by introducing oscillatory features and thus, mitigating spectral bias and improving the learning of high-frequency components.

To quantify the performance enhancement driven by this Fourier mapping under non-local dynamics, we consider the following time-fractional diffusion equation on the metric star graph given in Figure 2: 18

∂ C γ 0 Dt ui (x, t) =

2

ui (x, t) + fi (x, t), for i = 1, 2, 3, ∂x2 subject to the homogeneous initial condition ui (x, 0) = 0 for all i ∈ {1, 2, 3}. On boundary nodes v1 and v2 , ∂u2 2 2 1 the Neumann conditions are imposed as ∂u ∂x (1, t) = 24πt and ∂x (1, t) = 0.3π cos(0.3π)t , while on v3 , the Dirichlet condition u3 (1, t) = −5.4t2 has been imposed. The source terms fi has been chosen such that ui (x, t) = ci t2 sin(ωi x),

for i = 1, 2, 3,

forms the exact solution. The values of ci are chosen to satisfy continuity and Kirchhoff-Neumann conditions at the junction node v0 . The edge-wise spatial frequencies are: ω1 = 6π; ω2 = 0.3π; ω3 = 4.5π with coefficients c1 = 4.0; c2 = 1.0; c3 = −5.4. We apply our optimization strategies on the fractional order γ = 0.5 and use a configuration of 4 hidden layers with 128 neurons each. The Fourier embedding variant has 128 dimensions with σ = 1. The Fourier dimensions represents the number of waves sampled, hence 128 Fourier dimensions translate to 64 sine and 64 cosine waves. We use the dual scheme and train the Adam phase for 10000 iterations. All fractional derivatives are computed using the L2 − 1σ scheme. Table 3: Performance comparison of an embedded Fourier feature PINN with a standard PINN. Strategy

Relative L2 Error 1.90 × 10−1 5.92 × 10−3

Fourier Disabled Fourier Enabled

Figure 7: Performance of the network with Fourier embedding enabled and disabled. a) represents the edge residual loss, b) shows the node residual loss.

The numerical results in Table 3 demonstrates how Fourier feature embedded significantly reduce the relative L2 approximation error between the exact and predicted solution. This confirms that the the standard MLP-based PINN struggles to resolve the oscillatory components of the solution, while the Fourier feature architecture provides a more suitable input parametrization for high-frequency behavior. Further, Figure 7 shows this improvement by comparing the training behavior with and without Fourier embedding. The Fourier-enabled model significantly reduces the edge residual loss and node loss, and thus, support the use of Fourier feature embeddings as a default component in the proposed framework. 19

3.5. Singularity-Capturing Feature While non-uniform graded meshes are employed to minimize truncation errors in traditional discrete schemes, their direct implementation within PINNs may create severe optimization bottlenecks. Graded meshes heavily concentrate collocation points near the singular point of the domain, such as x = 0 or t = 0, which may cause the coefficients of the discrete operators defined in (2.3) and (2.4) to become large and highly uneven. Moreover, the smooth structure of the activation functions in standard MLPs makes it difficult to accurately capture singular derivatives near such points. To mitigate this behaviour, we incorporate a singularity-capturing feature in the spirit of Li and Tan [19]. Rather than forcing the network to approximate singular derivatives at the origin using only raw spatial or spatiotemporal coordinates, we introduce an explicit, non-smooth auxiliary singularity feature defined by Z(c) = cξ , where ξ > 0 is the singularity-capturing parameter, and c corresponds to the singular coordinate. In particular, c = t for time-fractional evolution problems, while c = x may be used for space-fractional elliptic problems on an edge parametrized from the singular endpoint. Although graded meshes increase the resolution near the singular point, they do not by themselves change the approximation space of the neural network. The singularity-capturing feature complements mesh grading by enriching the input space with a non-smooth coordinate Z(c) = cξ , which allows the network to represent fractional-type singular behaviour more naturally. In our graph-based framework, particularly for time-dependent problems, the neural network on each edge is evaluated using the augmented input u bi,θ (x, t) = Fi,θ (x, t, Z(t)). The mathematical utility of this feature becomes evident from the classical chain rule. For instance, ∂Fi,θ ∂Fi,θ dZ ∂Fi,θ ∂Fi,θ ξ−1 ∂b ui,θ = + = + ξt . ∂t ∂t ∂Z dt ∂t ∂Z which is singular at t = 0 whenever 0 < ξ < 1. Therefore, the augmented input allows the neural network to represent solution components whose derivatives exhibit fractional-type singularities near the initial time. The exponent ξ may be prescribed when the expected singularity order is known. However, in practical computational settings, especially with inverse parameter estimation, the analytical order of the initial singularity is unknown. Therefore, in our implementation, ξ is treated as a trainable parameter and optimized together with the network parameters θ. This allows the network to adapt the strength of the singular feature during training. We initialize ξ using the corresponding fractional order and update it dynamically via backpropagation. We compare the performance of the proposed singularity-capturing PINN with a standard PINN by considering the following problem defined on the tree graph in Figure 4. ∂ui ∂ C γ (x, t) − ν 0 Dt ui (x, t) + ui (x, t) ∂x

2

ui (x, t) + cui (x, t) = fi (x, t), ∂x2

x ∈ (0, 1),

t ∈ (0, 1],

(3.12)

where the fractional order is chosen as γ = 0.5, the kinematic viscosity (diffusion coefficient) is ν = 0.5, and the linear reaction coefficient is c = 5.0. The initial condition is ui (x, 0) = 0, for x ∈ [0, 1] and zero Dirichlet boundary conditions are applied at the boundary nodes v0 , v2 , v4 such that u(x, t) = 0, for t ∈ (0, 1]. The source terms fi (x, t) are chosen such that exact solution on the edge i is: ui (x, t) = Ai sin(kx)tγ ,

with k = 4π, and A = [2.0, 1.0, 1.0, 1.0].

(3.13)

For this comparison we utilize a configuration of 4 hidden layers with 128 neurons each. The Adam optimization phase is run for 10000 epochs followed by a L-BFGS refinement phase. We choose the dual adaptive weighting and Fourier feature embedding with 128 dimensions. All fractional derivatives were computed using the L2 − 1σ scheme. 20

Table 4: Effect of mesh grading and singularity-capturing feature on the relative L2 error Mesh grading factor r=1 r=2 r=3 r=4

Z(t) disabled

Z(t) enabled

2.1 × 10−3

9.71 × 10−4 7.65 × 10−4 1.31 × 10−3 1.44 × 10−3

2.00 × 10−3 2.38 × 10−3 2.78 × 10−3

Figure 8: Comparison of PINN approximations with and without the singularity-capturing feature Z(t) = tξ .

Table 4 highlights the benefits of the singularity capturing feature when the governing system has an initial singularity. The ability of the PINN to represent the solution as a function of tξ has resulted in a significant improvement in accuracy and lower PDE and boundary residuals, as seen in figure 8. The above ablation comparing the engine at various grading parameters highlights an interesting behaviour of our collocation methodology. Our network uses a unified mesh approach where the graded mesh (2.1) used for the approximation of the Caputo fractional derivative directly serves as the neural network training points. Graded meshes concentrate points towards the initial portion of the domain. While this reduces the truncation error of the approximation schemes (minimum at 2/α for L2 − 1σ ), it also starves later portions of the domain of collocation points. Therefore, for such implementations, the choice of the grading parameter r is a trade-off between minimizing the truncation error and ensuring that the entire domain is well represented with collocation points. 4. Framework Architecture and Implementation Blueprint To bridge the mathematical methods established in Section 3 with an efficient computational framework, we present a modular implementation architecture for QGPINNs. The implementation is designed to separate the graph topology, the nonlocal physical operators, and the neural network optimization components. This separation allows the framework to handle different graph configurations without requiring structural modifications to the underlying deep learning architecture. 4.1. Neural Network Structure on a Metric Graph In case of Euclidean domains, a standard neural network typically takes the coordinates as input through a single MLP. However, in the context of metric graphs, such a global configuration may become inefficient as physical coefficients, edge lengths, and solution gradients can vary significantly across different edges 21

of a complex network. Consequently, a single set of weights have to represent substantially different local behaviour across the graph, which can affect the convergence. This effect can be further amplified by the nonlocality of the fractional differential operators. Therefore, in our framework, we assign an independent edge-wise sub-network Ni (·; θi ) to each edge ei ∈ E, where θi denotes the trainable parameters of the subnetwork on the edge ei . Although the subnetworks are parametrized edge-wise, they are trained through a global graph loss. Thus, the approximation on each edge are coupled during optimization through the continuity, Kirchhoff– Neumann, and boundary residuals defined at the graph vertices. The user specifies the graph topology as an edge list |E|

Glist = {(ui , vi , ℓi )}i=1 ,

(4.1)

where ui , vi ∈ V denote the vertices of the edge ei with ℓi ∈ R+ denotes its length. During initialization, the computational engine builds a networkx graph object and allocates an independent sub-network to each edge. We illustrate a minimal code snippet to create a graph object below. class StarGraph : nodes = [0 , 1 , 2 , 3] edges = [(0 , 1 , L ) , (0 , 2 , L ) , (0 , 3 , L ) ]

Algorithm 1 describes the graph construction and edge-wise subnetwork allocation process. Algorithm 1 Graph Construction and Sub-Network Allocation |E|

Require: Edge list Glist = {(ui , vi , ℓi )}i=1 1: G ← BuildGraph(Glist ) ▷ networkx graph object 2: for each edge ei = (ui , vi , ℓi ) ∈ E do 3: Ni ← InitializeMLP(hidden_layers, hidden_dim) 4: if use_fourier_features then 5: Ni ← PrependFourierEmbedding(Ni , embed_dim, σ) 6: end if 7: AssignSubNetwork(ei , Ni ) ▷ edge-wise parameters coupled through the global graph loss 8: end for |E| 9: return G, {Ni }i=1

4.2. Separation of the Physics Layer from the Solver Engine Frameworks that define the network architecture and governing physics in a single script are difficult to generalize across different graph topologies and differential operators. We therefore separate the implementation into two components: a solver engine and a run-script. The solver engine handles all neural network architecture, while in run-script the user specifies only the graph topology and problem-specific information, such as the governing equation on each edge. In particular, the engine handles the creation of the fractional differentiation matrices and loss function definitions internally based on the parameters of the run script. The following code snippet describes the creation of the physics object. # The physics module for a time - fractional Burgers equation on edges def F ( self , x , t , u , u_x , u_xx , dt_alpha_u , f ) : return dt_alpha_u + u * u_x - self . nu * u_xx - f def get_ic ( self , x ) : # initial condition if necessary return torch . sin ( math . pi * x )

22

The above separation allows the same solver engine to be reused for different equations and graph topologies. Once the graph topology and governing equations are defined, the user configures the solver object using object attributes. Table 5 summarizes the main configuration attributes and their roles. Table 5: Important solver attributes. The framework also provides a secondary set of attributes to allow granular control over network architecture

Stage

Role

set_constraints() set_frac_scheme() set_mesh() set_architecture() set_singularity_capture() compile() train()

Specify boundary condition type (soft/hard) and values per vertex Select fractional discretization (L1 or L2-1σ) Choose spatial/temporal collocation density and grading Configure sub-network width, depth, and Fourier features Configure the singularity capturing feature. Initialise the network architecture Choose the weight balancing strategy, training iterations and begin optimization

|E|

Glist = {(ui , vi , ℓi )}i=1 (user-defined graph topology)

Global structure (topology)

Local structure (edge physics)

networkx Graph Object vertices V, edges E, lengths ℓi

Physics Definition (run script) F(x,t,u,u_x,u_xx,dt_alpha_u) get_ic(x), boundary data

Per-Edge Sub-Network Allocation edge-wise MLP Ni (x, t; θi ) per edge (+ optional Fourier feature embedding)

Fractional Discretization Caputo ∂tγ u via L1 / L2 − 1σ scheme non-local assembly matrix Wγ /Aγ

Solver Instantiation ParabolicPINNSolver(graph, physics_list) /

EllipticPINNSolver(...)

Configuration of neural network parameters and various mathematical schemes.

Training dual-phase optimization (Adam → L-BFGS) singularity-capturing feature, adaptive λ balancing continuity and Kirchhoff-Neumann coupling enforced at internal junction nodes Trained QGPINN Solution edge-wise approximations u bi,θ on ei globally coupled through vertex conditions

Figure 9: System architecture and pipeline flow of QGPINNs.

Boundary conditions are specified at boundary vertices by assigning the condition type (dirichlet/neumann) and the corresponding value or time-dependent function. the vertex transmission conditions at junction nodes are handled separately by the engine from the graph topology and enforced through the continuity and Kirchhoff–Neumann losses defined in (1.4) and (1.5). The following code snippet describes the syntax used to define the boundary condition. bc_types = {1: ’ neumann ’ , 2: ’ neumann ’ , 3: ’ dirichlet ’}

23

bc_values = { 1: lambda t : - math . pi * 5* torch . pow ( t . clamp ( min =1 e -10) , ALPHA ) , 2: 0.0 , 3: lambda t : 5* torch . pow ( t . clamp ( min =1 e -10) , ALPHA ) } solver . set_constraints ( ’ soft ’ , bc_types = bc_types , bc_values = bc_values )

After graph construction and solver configuration, the training procedure assembles the total loss (3.1). Depending on the selected configuration, the framework may also incorporate Fourier feature embeddings, hard constraint enforcement for Dirichlet boundary data, dynamic loss balancing, and singularity-capturing inputs. The optimization-phase follows a dual-stage schedule. First, the model is trained using the Adam optimizer together with the dynamic balancing strategy to update the loss weight λ. After the first-order stage, the value of λ is frozen, and the network parameters are further refined using L-BFGS. The trained model provides edge-wise approximations u bi,θ on all edges, which together define the approximate solution on the quantum graph. Figure 9 provides a complete overview of the implementation pipeline of QGPINNs, from graph input and physics definition to solver configuration, training, and output generation. The above modular and structured design of QGPINNs allows the framework to handle different graph topologies, fractional operators, boundary conditions, and training configurations without changing the core solver architecture. The next section demonstrates these capabilities through numerical experiments on various graphs, including real network structures. 5. Numerical Experiments This section evaluates the proposed framework against standard implementations and physical topologies to assess its accuracy and computational performance. Problem 1 is defined as a nonlinear fractional elliptic system on a tadpole graph. Tadpole graphs are well-studied in the quantum graph literature and are considered to be building blocks for more sophisticated graphs. We utilize the problem to provide insight into the efficacy of the various optimization schemes introduced in Section 3. For problem 2, we consider a time-fractional diffusion equation on the same tadpole graph. We compare the relative discrete L2 error (3.9) based on two different discretizations, namely, an auxiliary graded mesh with distributed training points [19], and a unified mesh, where the graded mesh 2.1, used for the approximation of the Caputo fractional derivative, directly serves as the neural network training points. Problem 3 describes a time fractional Burger’s equation defined on an open channel drainage system. The presence of heavy vegetation along the banks may induce non-local wave propagation behaviour, for which the considered Caputo derivative is suitable for the modeling. Problem 4 describes a time-fractional telegraph equation on the well-known bus network system. High-frequency voltage events are susceptible to memory effects, this implies that the usage of the Caputo derivative is better suited to describe such a phenomenon. Problems 3 and 4 do not contain an exact solution, thus, we employ secondary physical constraints to test the approximation’s validity. To mitigate initialization bias and ensure reproducibility, each problem is benchmarked across three fixed random seeds for 500 epochs. The optimal seed initialization is then selected for complete network training. Examples with an exact solution are tested for numerical accuracy based on the relative discrete L2 error given in (3.8) and (3.9). All experiments were executed on NVIDIA T4 and P100 GPUs. 5.1. Non-linear elliptic fractional differential equation on a Tadpole graph This problem evaluates the accuracy and computational overhead of the optimization schemes described in Section 3. The optimization schemes are successively enabled on a problem to identify their performance benefits. Consider the following elliptic differential equation defined on the tadpole graph given in Figure 10:  
2 dui β C α − C = fi (x), x ∈ (0, Li ), i ∈ {0, 1, 2, 3}, D u (x) + V u (x) i 0 Dx ui (x) 0 x i dx 24

where α = 1.5, and β = 0.5 and we choose the advection velocity coefficient V = 1.0. The topology has Dirichlet Constraints at boundary vertices v3 and v4 such that u3 (v3 ) = u3 (1) = 4 and u4 (v4 ) = u4 (1) = 1.6 · 20.6 . The forcing terms fi (x) are chosen to allow the following solution on each edge: for x ∈ e1 (Loop)

u1 (x) = x1.6 − 20.6 x, 1.6

u2 (x) = x

for x ∈ e2 (Inner Tail)

+ x,

1.6

+ A3 x + B3 x + 2,

for x ∈ e3 (Outer Tail)

1.6

+ B4 x,

for x ∈ e4 (Branch)

u3 (x) = x u4 (x) = x

2

The parameters A3 , B3 , and B4 have been obtained by enforcing the continuity (1.4) and Kirchhoff-Neumann (1.5) constraints at the nodes v1 and v2 . v4 e4 e2

e1 v1

e3 v2

v3

Figure 10: Tadpole graph representation. We choose the length of edges as L1 = 2.0; L2 = 1.0; L3 = 1.0 and L4 = 1.0, where Li is the length of the ith edge.

We compare the performance and accuracy of the proposed method by considering several variants with different features. The neural network parameters for each variant are defined in Table 6. Moreover, to quantify the computational load of each variant, we compare the peak GPU memory (MB) and a relative utilization metric defined as   X Utili GPU-seconds = ∆ti 100 i Table 6: Training and Mesh Parameters

Parameter

Variant 1

Variant 2

Variant 3

Variant 4

Spatial point density (Nx /Li ) Spatial Grading (r) Training Iterations Fractional scheme Neurons per layer Hidden Layers Constraint type Adaptive Weight Strategy Singularity Capture Fourier Dimension Fourier σ

100 2 10000 L2 − 1σ 128 4 Soft Fixed Disabled 0 -

100 2 10000 L2 − 1σ 128 4 Soft Dual Disabled 0 -

100 2 10000 L2 − 1σ 128 4 Soft Dual Enabled 0 -

100 2 10000 L2 − 1σ 128 4 Soft Dual Enabled 128 1

Table 7 shows the relative performance and computational load of each variant. The enabling of the singularity capturing feature in variant 3 improved accuracy significantly, pointing towards the capability of such a feature to efficiently and accurately capture singularities. Enabling Fourier feature embedding further in Variant 4 improved performance and reduced the global L2 error to the order of 10−5 . Figure 11(a) illustrates the computational load of each variant during training. The computational cost of each variant highlights the relative efficiency of these methods. While such numerical accuracy trends can also be obtained by increasing point densities and neural network sizes, the above optimization strategies are significantly more efficient. We observe that Variant 4 required an additional 141 MB of memory, however it improves 25

the accuracy significantly as compared to other variants. Figure 11(b) depicts the predicted solution on the tadpole graph for Variant 1, where one can observe that the predicted solution satisfies the continuity constraint (1.4) at the junction nodes v1 and v2 . Table 7: Performance comparison of variants. Variant

Global Relative L2 Error

GPU Memory (MB)

GPU-S

Variant 1 Variant 2 Variant 3 Variant 4

2.62 × 10−2

270 270 270 348

250 256 340 391

2.51 × 10−2 1.33 × 10−4 8.50 × 10−5

(a)

(b)

Figure 11: (a) Computational demand of each variant. (b) A 3-Dimensional plot of the predicted solution.

A major benefit of PINNs over traditional solvers is the ability to estimate system parameters using sparse and noisy data. Our framework supports the estimation of parameters using the set_inverse() method. For this problem, we choose 200 data points for an edge. Since the above system has an exact solution, these points are chosen at random from each edge. For realistic settings, we systemically add Gaussian noise to the data set. For our experiment, we utilize the network parameters of Variant 4 and run the model for 6000 iterations with an initial guesses of α̂ = 1.2, β̂ = 0.2 and V̂ = 0.6. Table 8 highlights the ability of the proposed framework to estimate system parameters and Figure 12 illustrates the convergence of the parameters from their initial guess to their final values. Table 8: Simultaneous Parameter Recovery Profile under Synthetic Noise

Noise (σ) 0% 1% 5%

Estimated α̂

Estimated β̂

Estimated V̂

1.4954 1.4929 1.4947

0.5032 0.4997 0.4961

1.0014 0.9983 0.9938

26

Figure 12: Convergence of Parameters for each noise level

5.2. Time-fractional diffusion equation on a Tadpole graph We consider the following time-fractional diffusion equation on a tadpole graph: ∂ C γ 0 Dt ui (x, t) −

2

ui = fi (x, t), ∂x2

x ∈ (0, Li ),

i ∈ {1, 2, 3, 4}, t ∈ [0, 1].

The graph topology employs Neumann constraints at the boundary nodes v3 (x = 1) and v4 (x = 1) such ∂u3 γ γ 2 that ∂u ∂x (1, t) = πt and ∂x (1, t) = πt , along with the continuity and Kirchhoff-Neumann constraints at node v1 and v2 . The initial condition across all edges is set to ui (x, 0) = 0. The forcing terms fi (x, t) has been chosen such that for x ∈ e1 (Loop)

u1 (x, t) = tγ (sin(πx) + cos(πx)) , γ

for x ∈ e2 (Inner Tail)

γ

for x ∈ e3 (Outer Leaf)

γ

for x ∈ e4 (Branch Leaf)

u2 (x, t) = t (sin(πx) + cos(πx)) , u3 (x, t) = t (− sin(πx) − cos(πx)) , u4 (x, t) = t (− sin(πx) + cos(πx)) ,

forms the exact solution of the problem. In this problem, we compare two different collocation strategies for approximating the Caputo fractional derivative. The first strategy follows the approach of Li et al. [19], in which each training point is associated with an auxiliary graded mesh which is particularly used for approximating the fractional derivative. The primary training points and the auxiliary discretization are therefore decoupled, which allows the Caputo derivative to be computed using a locally constructed graded mesh at every collocation point. The second strategy is the propose unified mesh where a single graded mesh is used for both the neural network collocation points and the discretization of the fractional derivative. Consequently, the corresponding Caputo differentiation matrix is assembled only once during the initialization stage and reused throughout the optimization process. Table 9: Training and Mesh Parameters

27

Parameter

Value

Spatial point density (Nx /Li ) Temporal points (Nt ) Temporal Grading (r) Training Iterations Fractional scheme Neurons per layer Hidden Layers Constraint type Adaptive Weight Strategy Singularity Capture Fourier Dimension Fourier σ

80 50 2 10000 L2 − 1σ 128 4 Soft Dual Enabled 128 2

The usage of auxiliary graded meshes poses significant challenges when adapting the approach to governing systems defined on quantum graphs. The computation of the Caputo derivative matrix at every training point introduces significant computational bottlenecks, especially when defining problems on complex realworld topologies. Therefore, in our framework we prioritize the usage of a single unified graded mesh where the Caputo fractional matrix is computed once during initialization as described in Section 2.2. We compare these two strategies to assess their influence on the accuracy and computational performance of the proposed framework. The neural network parameters for both collocation methodologies are given in Table 9. As the auxiliary mesh implementation requires each training point to be supplemented with a unique auxiliary mesh, we choose 2 variations. Variant 1 supplements each training point with 50 auxiliary points, variant 2 supplements each training point with 25 auxiliary points. Each auxiliary mesh is graded with a grading factor of r = 4. Table 10: Performance comparison of variants for γ = 0.5 Variant Unified Mesh Auxiliary Mesh (50 points) Auxiliary Mesh (25 points)

Global Relative L2 Error

GPU Memory (MB)

GPU-S

7.73 × 10−4

1146 6814 3932

473.6 2909.3 1694.0

5.83 × 10−4 7.48 × 10−4

Table 11: Loss values at Epoch 10000

Mesh Type 50 Auxiliary Points 25 Auxiliary Points Unified Mesh

PDE Loss 7.63 × 10−5 8.07 × 10−5 3.36 × 10−4

BC Loss 5.49 × 10−5 8.52 × 10−5 4.23 × 10−5

Table 10 provides the relative L2 accuracy and computational cost of each collocation strategy. The results highlight the drawbacks of using auxiliary meshes in our framework. They provide very minimal performance benefits while requiring significantly more computational resources. Figure 13 highlights the computational resource use of each variant. Their inability to provide superior accuracy is due to mesh decoupling. Computing the Caputo derivative using localized auxiliary meshes results in the network attempting to overfit to the local derivative computation. This is evident from the loss data shown in 11 where the auxiliary mesh with 50 points had a significantly lower PDE loss but slightly higher boundary condition loss. However, this lower PDE loss did not translate to an effective reduction in the approximation error. This is a consequence of the derivative overfitting as described above. This becomes more evident when comparing the auxiliary mesh with 25 points and the unified mesh approach. We therefore promote the use of a single unified mesh which serves as an ideal tradeoff between approximation accuracy and computational efficiency. Figure 14 illustrates the comparison of exact and predicted solution on each edge of the graph at t = 1. It is evident that the predicted solution obtained from the proposed framework matches closely with the exact solution of the considered problem. 28

(a)

(b)

Figure 13: Computational cost of each collocation methodology

Figure 14: Predicted vs exact solution on each edge at t = 1s.

Furthermore, the exact solution allows us to test the framework’s ability to estimate parameters from x noisy and sparse data. We choose a spatial point density as N Li = 100 and Nt = 200 temporal points. We incorporate 200 data points from each edge into the loss function, with initial guesses of γ̂ = 0.35 and V̂ = 0.6. Similar to the elliptic problem 5.1, we simulate real world data for time-fractional diffusion problem by adding varying levels of Gaussian noise. Using the unified mesh strategy, we compare the accuracy of the L1 and L2 − 1σ fractional schemes, while all other neural network parameters remain identical to the forward problem. Table 12 demonstrates the framework’s ability to extract the governing parameters from noisy data. The neural network using the higher-order L2 − 1σ scheme achieved superior accuracy. Figure 15 illustrates the convergence of the estimated parameters to their final values during training for each fractional scheme. 5.3. Time-Fractional Burgers’ Type Equation for Open-Channel Flow Routing We simulate sub-diffusive water wave propagation in an open channel drainage system in Japan, following the network topology structured by Yoshioka et al. [43]. The system models a real world drainage system present in Japan as shown in Figure 16. The presence of heavy vegetation along the banks of the drainage 29

Table 12: Simultaneous Parameter Recovery Profile under Synthetic Noise

Fractional Scheme Used L2 − 1σ

L1 Noise (σ) 0% 1% 5%

Estimated γ̂

Estimated V̂

Estimated γ̂

Estimated V̂

0.4920 0.4918 0.4866

1.0003 0.9855 0.9611

0.4974 0.4979 0.4871

1.0010 0.9925 0.9680

(a) Estimated parameter value vs. iterations (L1 scheme)

(b) Estimated parameter value vs. iterations (L2 − 1σ scheme)

Figure 15: Parameter convergence history.

system introduces anomalous wave propagation and memory effects that the Caputo derivative is suited to handle. The problem is described by embedding the real world topology on a connected metric graph 17 consisting of six boundary vertices and five edges (reaches), normalized over a physical length Li = 20.0 m and a temporal window Tmax = 10.0 s. The governing differential equation for the depth fluctuation ui (x, t) = hi (x, t) − h0 over a steady equilibrium depth h0 = 5.0 m on edge ei is given by: C γ 0 Dt ui (x, t) +

i ∂ h ∂ 2 ui V (ui )ui − D 2 = 0, ∂x ∂x

x ∈ (0, Li ),

t ∈ (0, Tmax ],

(5.1)

where γ = 0.85 models sub-diffusive porous interactions along heavily vegetated banks, D = 25.0 m2 /s is the physical diffusion coefficient, and g = 9.81 m/s2 . The advection parameter V (hi ) is defined by [43]: V (hi ) =

i 3 3 2√ h g (hi + h0 ) 2 − h02 . 3

To map the spatial and temporal coordinates into unit intervals [0, 1], we introduce the scaling variables x̂ = x/Li ∈ [0, 1] and t̂ = t/Tmax ∈ [0, 1], where Li = 20.0 m and Tmax = 10.0 s. Substituting these transformation variables into Eq. (5.1) and multiplying through by (Tmax )γ yields the non-dimensional residual used for training: h i ∂ 2 ui V (ui )ui − Λdiff 2 = 0, ∂ x̂ ∂ x̂

∂ C γ 0 Dt̂ ui (x̂, t̂) + Λadv

x̂ ∈ [0, 1],

t̂ ∈ (0, 1],

where the advection and diffusion coefficients Λadv and Λdiff , respectively, preserve the physical scale: Λadv =

(Tmax )γ , Li

Λdiff =

30

(Tmax )γ D . L2i

(5.2)

B

A

E0

E3

E

E2

E4

D

F E1

C Figure 17: Agricultural drainage network topology from Japan. Inlets are located at vertices A (Sakai River), B (Drainage Canal), and C (Eigo River). The network terminates at downstream outlet D, which discharges into Lake Biwa.

Figure 16: Satellite imagery of the topology

At t̂ = 0, the channel network is initially at rest with zero water depth fluctuation, ui (x̂, 0) = 0. Dynamic high-energy boundary surges are applied at t̂ > 0 via time-dependent Dirichlet conditions at all inlets (A, B, C), that is
u3 (0, t̂) = u0 (0, t̂) = u1 (0, t̂) = 2.0 sin2 π t̂ . (5.3) A Neumann boundary condition ( ∂u ∂ x̂ = 0) is enforced at the downstream outlet vertex D. As our open channel drainage system is governed by advection, general flux conservation does not capture the true physical system. Therefore, we define a custom conservative flux operator Q̂i (x̂, t̂) that combines both the non linear advection transport and the diffusive wave motion along edge i: ∂ui (5.4) ∂ x̂ The system satisfies the Kirchhoff-Neumann junction conditions using this complete flux operator. At any internal junction node v, the total flux in the incoming edges Ein (v) must be equal to the flux in the outgoing edges Eout (v). That is, X X Q̂i (1, t̂) = Q̂j (0, t̂) (5.5) Q̂i (x̂, t̂) = Λadv V (ui )ui − Λdiff

i∈Ein (v)

j∈Eout (v)

To define such a custom flux, our framework uses the inbuilt get_flux method. This allows the user to define a flux operator. The neural network parameters for this problem are given in table 13 Table 13: Training and mesh parameters for time-fractional Burgers’ equation.

31

Parameter

Value

Network topology x Spatial point density( N Li ) Temporal points(Nt ) Temporal grading (r) Training iterations Fractional scheme Neurons per layer Hidden layers Constraint type Adaptive weight strategy Singularity capture Fourier dimension

6 nodes, 5 edges 200 200 2 15000 L2 − 1σ 128 4 Soft Dual Enabled 128

Since the system does not have an exact solution, we verify the numerical solution using secondary physical constraints. In classical (integer-order) systems where the instantaneous flux is defined by Qnet = dM dT , the following conservation law holds: Z t

M (t) − M (0) =

Q(τ )dτ 0

However, the considered system is sub-diffusive, thus, the rate of mass accumulation is not instantaneous due to the presence of the Caputo derivative. This models the phenomenon of the vegetative banks absorbing, storing and gradually releasing water over time. Assuming water density and channel width are constant with time, the above quantity can serve as a proxy for global mass fluctuation, and hence satisfies the global fractional integral conservation balance: 1 M̂ (t̂) − M̂ (0) = It̂γ (Q̂net ) = Γ(γ)

Z t̂

(t̂ − τ )γ−1 Q̂net (τ )dτ,

(5.6)

0

where Q̂net (t̂) = Q̂0 (0, t̂) + Q̂1 (0, t̂) + Q̂3 (0, t̂) − Q̂4 (1, t̂) represents the net boundary flux into the system. The global stored volume fluctuation per unit width M (t) is defined as M̂ (t̂) =

4 Z 1 X i=0

ui (x̂, t̂) dx̂.

0

Table 14: Open-System Mass Conservation, γ = 0.85

Time t̂ = 0.20 t̂ = 0.50 t̂ = 0.80 t̂ = 1.00

M̂ (t̂) − M̂ (0) −1

5.9106 × 10 8.6204 × 100 7.4173 × 100 4.0244 × 100

I γ [Q̂net ](t̂) −1

5.9061 × 10 8.6190 × 100 7.3766 × 100 4.0561 × 100

Rel. Error (%) 0.077 0.016 0.549 0.790

The derivatives and integrals involved in the above expressions are evaluated using automatic differentiation and scipy.integrate respectively. Table 14 describes the mass conservation metrics. Since deviations from these conservation laws are not explicitly penalized in the QGPINNs framework during training, the small observed errors indicate that the framework is able to recover physically consistent behaviour from the governing residuals and graph constraints. Moreover, the right-hand side of the expression (5.6) requires fractional integration of order γ, and therefore provides an additional diagnostic for verifying fractional-order dynamics captured by the model. In Table 15, we examine the Kirchhoff-Neumann and continuity constraints 32

Table 15: Junction Continuity Residuals, γ = 0.85

Junction E E E E F F F

Time t̂ = 0.20 t̂ = 0.50 t̂ = 0.80 t̂ = 1.00 t̂ = 0.50 t̂ = 0.80 t̂ = 1.00

Q̂in

Q̂out −2

Rel. Residual (%) −2

2.0283 × 10 1.6337 × 101 7.8457 × 100 2.0642 × 100 1.1534 × 101 1.4285 × 101 5.1354 × 100

2.1648 × 10 1.6341 × 101 7.8419 × 100 2.0688 × 100 1.1523 × 101 1.4284 × 101 5.1481 × 100

6.308 0.026 0.048 0.224 0.093 0.006 0.248

at the intermediate junctions E and F by measuring flow residuals across various temporal levels. The significant deviation at junction E at t̂ = 0.2 can be observed as the wave has not propagated to b out (0.2) − Q b in (0.2)) = 1.36 × 10−3 appears to be a the junction, therefore, the small absolute residual (Q substantial large when measured relative to the total flux. Once the wave propagates to the junction and reaches non-negligible levels, the underlying residual is reduced to below 0.3%.

(a)

(b)

Figure 16: (a) Water level profile along the path B → E → F → D at time t after reversing the normalization. (b) Heatmap representing the water height fluctuation at time t = 5s.

Figure 16(a) illustrates the water level fluctuation along the path B → E → F → D. The normalization is reversed in order to interpret the solution on the original real-world topology. The behaviour of the predicted solution is consistent with the imposed junction constraints. The sharp upward kinks, for example at junction E when t = 5s, are the consequence of the periodic influx of water from multiple inlet branches, which causes a change in slope across the junction. The same qualitative behaviour is observed at junction F, where the two incoming branches combine to supply the downstream reach. This problem demonstrates the capability of the framework to simulate complex real-world network dynamics in the absence of an exact solution. Figure 16(b) shows the water height fluctuation at t = 5s, where the qualitative effect of the Kirchhoff–Neumann flux balance can be observed across the network.

33

5.4. Fractional Telegraph Model for High-Frequency Surge Propagation on a Transmission Network This problem aims to test the robustness of our framework on complex real-world topologies. We simulate the propagation of a high frequency, finite speed wave across an electric transmission network. Such phenomena are characteristic of voltage surges due to lightning strikes and other such events. Such systems undergo frequency attenuation, where the high frequency characteristics dampen over time. Standard integer-order wave models fail to capture these power-law dissipation dynamics. We consider the well-known IEEE-14 bus network as the metric graph benchmark shown in Figure 17. The governing differential equation is a time-fractional telegraph equation that models the surge voltage ui (x, t) on edge ei . Each edge present in the IEEE-14 system has a corresponding wave ci > 0. 11 1.1

1.2

1.2

4

1.1

5

12

1.0

1.0

0

0.9

10

0.8

0.

1.1

1

3

1.

1.3

0.9

6 1.4

1.0

9

8

1.

0.7

1.0

0

3

2

8 0.

1.0

9

13

7 Figure 17: IEEE 14-Bus System Topology (0-Indexed)

2 ∂ 2 ui C γ 2 ∂ ui (x, t) + η D u (x, t) − c (x, t) = 0, x ∈ (0, Li ), t ∈ (0, T ]. (5.7) i i 0 t i ∂t2 ∂x2 We choose the fractional order γ = 0.5 and similar to problem considered in Section 5.3, the system is normalized before training. The spatial coordinate is handled by the engine’s per-edge input normalization x̂ = x/Li while time is mapped to t̂ = t/T ∈ [0, 1], with T = 5.0 s. Substituting into Eq. (5.7) and multiplying through by T 2 gives the residual used in training,

∂ 2 ui ∂ 2 ui γ + Λdamp,i C = 0, 0 Dt̂ ui − Λwave,i 2 ∂x2 ∂ t̂ where

Λdamp,i = ηi T 2−γ ,

(5.8)

Λwave,i = c2i T 2 .

To simulate a lightning strike, we impose a spatial initial surge profile centered at vertex v0 (node 0) such that for all edges ek incident to v0 we specify "   # 2 x ∂ui uk (x, 0) = A exp − , (x, 0) = 0, (5.9) λi ∂t where A = 1 is the peak surge amplitude and λi = 0.20 controls the spatial penetration length along edge ei . The value of these constants are chosen to ensure that the pulse quickly decays to machine precision within the spatial domain of the initial edge. Therefore, edges ej not incident to v0 are initialized with homogeneous ∂u initial conditions, that is, uj (x, 0) = 0 and ∂tj (x, 0) = 0. At all the junction nodes, Kirchhoff-Neumann and continuity conditions have been imposed to ensure that the solution is mathematically valid. The governing 34

system contains a classical second-order time derivative. While our engine does not natively support such derivatives, the addition of a run script helper class, which computes and passes the required derivative using torch.autograd(), served as enough of a modification, which also highlights the flexibility of the given framework. However, we emphasise that the singularity capturing feature describe in Section 3.5 is incompatible with the considered problem. The term z(t) = t0.5 , when passed through the automatic differentiator in time, 3 results in t− 2 . This resulting term blows up to incredibly high values near t = 0, creating artificial gradients. Therefore, specifically for this problem, we disable the singularity capturing feature. The neural network parameters for the problem are given in Table 16. Table 16: Training and mesh parameters for the fractional Telegraph problem.

Parameter

Value

Network topology x Spatial point density ( N Li ) Temporal points (Nt ) Temporal grading Training iterations Fractional scheme Neurons per layer Hidden layers Constraint type Adaptive weight strategy Singularity capture Fourier dimension Fourier σ

IEEE 14-bus, 20 edges 60 80 2 15000 L2 − 1σ 128 4 Soft Dual Disabled 128 1

Figure 18: Plot of energy vs time in Log-Log axes

To verify physical consistency, we track the total mechanical-type network energy functional: " 2  2 # 20 Z 1 X Li ∂ui ∂ui 2 E(t) = (x, t) + ci (x, t) dx. 2 i=1 0 ∂t ∂x 35

(5.10)

Table 17: Comprehensive Physical-Consistency and Energy Diagnostics Across Fractional Orders γ.

γ = 0.50

γ = 0.75

γ = 1.00

Initial-Boundary Conditions Initial Pulse Relative RMS Error Initial Pulse L2 Error Initial Velocity L2 Error

0.388% 4.206 × 10−4 1.516 × 10−3

0.558% 6.038 × 10−4 2.488 × 10−3

0.467% 5.057 × 10−4 2.324 × 10−3

Network Transmission Constraints Junction Continuity L2 Error Physical Kirchhoff Flux L2 Residual Grounded Terminal (Node 0) L2 Residual

3.707 × 10−3 9.375 × 10−4 1.598 × 10−3

3.872 × 10−3 1.494 × 10−3 1.804 × 10−3

2.493 × 10−3 1.109 × 10−3 1.300 × 10−3

Global Energy Physics Initial Network Energy E(0) Final Network Energy E(Tmax ) Energy Fraction Remaining (Efinal /Einitial ) Energy Monotonicity Fraction Late-Time Log-Log Slope

3.728 0.376 10.09% 98.74% −1.564

3.726 0.082 2.21% 98.11% −2.599

3.737 0.020 0.53% 98.74% −3.523

Pulse Propagation Peak Surge Amplitude (t = 0.0 s) Peak Surge Amplitude (t = 0.05 s) Peak Surge Amplitude (t = 0.10 s) Peak Surge Amplitude (t = 0.20 s) Peak Surge Amplitude (t = 1.00 s)

1.001 0.492 0.463 0.412 0.116

1.001 0.487 0.446 0.373 0.068

1.002 0.483 0.427 0.335 0.043

Metric / Quantity

The energy is computed after reversing the normalization, and the spatial integrals are computed numerically using scipy.integrate. Since the underlying system is dissipative, the total network energy is expected to be non-increasing in time due to frequency attenuation. We therefore plot E(t) as a diagnostic quantity. In addition to the energy behaviour, we also verify junction continuity and Kirchhoff-Neumann flux conditions over the temporal validation grid. We also compute several global physical quantities and report their deviations in Table 17. Similar to problem 5.3, the majority of these quantities were not explicitly minimized as part of the loss function during training. Hence, their small deviations depicted in Table 17 shows that the learned solution from QGPINNs is physically consistent. The framework captures the initial driving pulse with high accuracy, while maintaining the continuity and Kirchhoff–Neumann constraints across the large graph topology. A particularly important attribute is the global energy decay. As shown in Figure 18, the expected non-increasing behaviour of E(t) is observed over 98.74% of the temporal domain. The same figure also compares the late-time log–log decay slopes for the fractional orders 0.5 and 0.75 and the integer-order case 1.0. The slower energy dissipation observed for the fractional-order models indicates that QGPINNs captures an algebraic power-law decay rather than the standard exponential decay associated with integer-order dynamics [22]. We also emphasise that GPU memory limitations restricted the spatial and temporal point densities on such large topologies. Neverthless, the observed loss trends suggest that increasing the number of validation and training points may further improve the accuracy of the QGPINNs. Figure 19 illustrates the propagation of the voltage spike over time. The spatial arrangement in Figure 19 differs slightly from the graph shown in Figure 17 due to automatic graph layout generated by networkx, which dynamically optimize node distribution and canvas usage. Despite this difference in layout, the movement of the pulse along the graph edges and the vertices can be clearly observed.

36

Figure 19: Propogation of the voltage surge at time t

6. Conclusion The proposed QGPINN framework provides a physics-informed deep learning approach for computing the numerical solutions to nonlocal differential equations on quantum graphs. While fractional differential equations have been widely studied on Euclidean domains, their extension to quantum graphs has primarily been addressed through traditional numerical methods. The proposed framework provides a physics-informed alternative by combining edge-wise neural approximations with a unified graph-based loss formulation. The addition of the various optimization schemes and architectural strategies described in section 3 allows the framework to provide superior accuracy when compared to standard PINN frameworks without substantially increasing the computational cost. The numerical experiments in Section 5 demonstrates the robustness of the method on a range of graph structures. The QGPINNs framework successfully computed the numerical solution to fractional differential equations defined on real world complex topologies. Although the present implementation focuses predominantly on elliptic and time-fractional evolution equations given in (1.1) and (1.2), respectively, the modular design of the framework remain compatible to a broad class of problem structures and graph topologies, as illustrated in Problem 5.4. This modularity also allows users to solve such systems without requiring extensive familiarity with the internal details of PyTorch. Although the design of the QGPINNs reduces the need for problem specific calibration, the choice of the appropriate adaptive weighting strategy, the number of Fourier dimensions, and the size of the neural network may still depend on the particular application. The framework initializes an edge-wise neural network to each edge instead of a single neural network over the entire graph, which allows it to capture the edgedependent behaviour and the memory effects associated with nonlocal Caputo derivatives. However, such design may not be computationally viable on extremely large networks. Future work in this direction will focus on incorporating fewer problem specific parameters and also improving the computational efficiency of the framework to allow it to handle large and more complex topologies.

37

Acknowledgments The first author acknowledges the Anusandhan National Research Foundation (ANRF), Government of India, for supporting this work through the Prime Minister Early Career Research Grant (PM-ECRG) under Grant No. ANRF/ECRG/2024/002135/PMS. Declaration of generative AI and AI-assisted technologies in the manuscript preparation process During the preparation of this work the authors used Gemini 3.1 Pro and Claude Sonnet 5 in order to assist with code implementation and debugging, organization of the manuscript, and language and grammar checks of the manuscript. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article. References [1] A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, Journal of Computational Physics, 280 (2015), pp. 424–438. [2] M. Arioli and M. Benzi, A finite element method for quantum graphs, IMA J. Numer. Anal., 38 (2018), pp. 1119–1163, https://doi.org/10.1093/imanum/drx029. [3] R. Band and G. Lévy, Quantum graphs which optimize the spectral gap, in Annales Henri Poincaré, vol. 18, Springer, 2017, pp. 3269–3323. [4] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of machine learning research, 18 (2018), pp. 1–43. [5] J. V. Below, Classical solvability of linear parabolic equations on networks, Journal of differential equations, 72 (1988), pp. 316–337. [6] G. Berkolaiko and P. Kuchment, Introduction to quantum graphs, vol. 186 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 2013, https://doi.org/10.1090/ surv/186. [7] G. Berkolaiko and P. Kuchment, Introduction to Quantum Graphs, vol. 186, American Mathematical Society, Providence, RI, 2013. [8] D. Bolin, M. Kovács, V. Kumar, and A. B. Simas, Regularity and numerical approximation of fractional elliptic differential equations on compact metric graphs, Mathematics of Computation, 93 (2024), pp. 2439–2472. [9] D. Bolin, L. Riera-Segura, and A. B. Simas, Numerical approximation of fractional diffusion equations on metric graphs, arXiv preprint arXiv:2608.01932, (2026). [10] D. Bolin, A. B. Simas, and J. Wallin, Gaussian Whittle–Matérn fields on metric graphs, Bernoulli, 30 (2024), pp. 1611–1639. [11] D. Bolin, A. B. Simas, and J. Wallin, Statistical inference for Gaussian Whittle–Matérn fields on metric graphs, Journal of the Royal Statistical Society Series B: Statistical Methodology, (2026). [12] M. Faheem and A. Khan, A collocation method for time-fractional diffusion equation on a metric star graph with η edges, Mathematical Methods in the Applied Sciences, 46 (2023), pp. 8895–8914.

38

[13] L. Guo, H. Wu, X. Yu, and T. Zhou, Monte Carlo fPINNs: Deep learning method for forward and inverse problems involving high dimensional fractional partial differential equations, Computer Methods in Applied Mechanics and Engineering, 400 (2022), p. 115523. [14] M. Kovács and E. Sikolya, Stochastic reaction–diffusion equations on networks, Journal of Evolution Equations, 21 (2021), pp. 4213–4260. [15] S. Kumari, M. Mehra, and V. Mehandiratta, Finite difference approximation of time-fractional advection-diffusion equation on a metric star graph, International Journal of Computer Mathematics, 102 (2025), pp. 2032–2050. [16] C. L. Laczkó, M. A. Vághy, and M. Kovács, A transferable PINN-based method for quantum graphs with unseen structure, IFAC-PapersOnLine, 59 (2025), pp. 67–72. [17] I. E. Lagaris, A. Likas, and D. I. Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE transactions on neural networks, 9 (1998), pp. 987–1000. [18] J. E. Lagnese, G. Leugering, and E. Schmidt, Control of planar networks of Timoshenko beams, SIAM journal on control and optimization, 31 (1993), pp. 780–811. [19] H. Li and Z. Tan, A novel singularity-and discontinuity-capturing pinn for time-fractional diffusion equations involving initial singularities and interfaces on complex curved surfaces, Applied Mathematics and Computation, 519 (2026), p. 129917. [20] H. L. Liao, W. McLean, and J. Zhang, A second-order scheme with nonuniform time steps for a linear reaction-sudiffusion problem, Communications in Computational Physics, 30 (2021), pp. 567–601. [21] Y. Lin and C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, Journal of computational physics, 225 (2007), pp. 1533–1552. [22] F. Mainardi, Fractional calculus and waves in linear viscoelasticity, World Scientific Publishing Company, 2022. [23] V. Mehandiratta and M. Mehra, A difference scheme for the time-fractional diffusion equation on a metric star graph, Applied Numerical Mathematics, 158 (2020), pp. 152–163. [24] V. Mehandiratta, M. Mehra, and G. Leugering, Existence and uniqueness results for a nonlinear Caputo fractional boundary value problem on a star graph, Journal of Mathematical Analysis and Applications, 477 (2019), pp. 1243–1264. [25] V. Mehandiratta, M. Mehra, and G. Leugering, Optimal control problems driven by timefractional diffusion equations on metric graphs: optimality system and finite difference approximation, SIAM Journal on Control and Optimization, 59 (2021), pp. 4216–4242. [26] V. Mehandiratta, M. Mehra, and G. Leugering, Well-posedness, optimal control and discretization for time-fractional parabolic equations with time-dependent coefficients on metric graphs, Asian Journal of Control, 25 (2023), pp. 2360–2377. [27] R. H. Nochetto, E. Otárola, and A. J. Salgado, A PDE approach to fractional diffusion in general domains: a priori error analysis, Foundations of Computational Mathematics, 15 (2015), pp. 733– 791. [28] A. Odžak and L. Šćeta, On the Weyl law for quantum graphs, Bulletin of the Malaysian Mathematical Sciences Society, 42 (2019), pp. 119–131. [29] G. Pang, L. Lu, and G. E. Karniadakis, fPINNs: Fractional physics-informed neural networks, SIAM Journal on Scientific Computing, 41 (2019), pp. A2603–A2626. 39

[30] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al., Pytorch: An imperative style, high-performance deep learning library, Advances in neural information processing systems, 32 (2019). [31] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. Hamprecht, Y. Bengio, and A. Courville, On the spectral bias of neural networks, in International conference on machine learning, PMLR, 2019, pp. 5301–5310. [32] S. Ramchandra and V. Mehandiratta, QGPINNs-framework-to-solve-fdes-on-quantum-graphs, 2026, https://github.com/Saket2006/QGPINNs-Framework-to-solve-FDEs-on-Quantum-Graphs. GitHub repository. [33] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Learning representations by back-propagating errors, nature, 323 (1986), pp. 533–536. [34] A. K. Singh, M. Mehra, and R. Pulch, Non-local physics informed neural networks for forward and inverse problems containing non-local operators, Neural Comput. Appl, 37 (2024), pp. 1–22. [35] Z. Sobirov, K. U. Rakhimov, and R. E. Ergashov, Green’s function method for time-fractional diffusion equation on the star graph with equal bonds, Nanosystems: Physics, Chemistry, Mathematics, 12 (2021), pp. 271–278. [36] A. A. Sophiya, A. K. Nair, S. Maleki, and S. K. Krishnababu, A comprehensive analysis of PINNs: Variants, applications, and challenges, arXiv preprint arXiv:2505.22761, (2025). [37] M. Stynes, E. O’Riordan, and J. L. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM Journal on Numerical Analysis, 55 (2017), pp. 1057–1079. [38] N. Sukumar and A. Srivastava, Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks, Computer Methods in Applied Mechanics and Engineering, 389 (2022), p. 114333. [39] M. Tancik, P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. Barron, and R. Ng, Fourier features let networks learn high frequency functions in low dimensional domains, Advances in neural information processing systems, 33 (2020), pp. 7537–7547. [40] Y. Vats, M. Mehra, and D. Oelz, Time-fractional cable equation for anomalous electrodiffusion on a metric star graph with existence, uniqueness and numerical solution, Mathematics and Computers in Simulation, (2026). [41] J. Wang, Y. Yu, J. Hou, and X. Meng, Local error analysis of L1 scheme for time-fractional diffusion equation on a star-shaped pipe network, Physica Scripta, 99 (2024), p. 125274. [42] S. Wang, Y. Teng, and P. Perdikaris, Understanding and mitigating gradient flow pathologies in physics-informed neural networks, SIAM Journal on Scientific Computing, 43 (2021), pp. A3055–A3081. [43] H. Yoshioka, K. Unami, and M. Fujihara, Burgers type equation models on connected graphs and their application to open channel hydraulics, RIMS Kôkyûroku, 1890 (2014), pp. 160–171. [44] Y. Zhao and M. L. Pasini, A deep learning approach to solve forward differential problems on graphs, arXiv preprint arXiv:2210.03746, (2022).

40

Related documents

• Physics-Informed Neural Networks for Quantum Eigenvalue Problems#19800
• QGPINNs: A Physics-Informed Neural Network Framework for Nonlocal Differential Equations on Quantum Graphs#792581
• QGPINNs: A Physics-Informed Neural Network Framework for Nonlocal Differential Equations on Quantum Graphs#788581
• Physics-Informed Classical and Quantum Neural Networks for One-Dimensional Schrodinger Eigenvalue Problems#1036872
• Continuous-Time Quantum Walks based Graph Neural Network#659292