fTNN: a tensor neural network for fractional PDEs

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

fTNN: a tensor neural network for fractional PDEs

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

fTNN: a tensor neural network for fractional PDEs Qingkui Ma∗ Hehu Xie† and Xiaobo Yin‡

arXiv:2606.27140v1 [cs.LG] 25 Jun 2026

Abstract We develop the fTNN, a deterministic tensor neural network subspace method for problems involving the fractional Laplacian on bounded domains, taking the fractional Poisson equation and time-dependent fractional advection-diffusion equation as typical representatives. The work employs a geometry-adapted integration split featuring a spatially dependent near-field radius, which decomposes the fractional Laplacian into three contributions: a singular near field, a regular interior far field, and an analytical exterior far field. Then the singular radial integrals are treated by Gauss-Jacobi quadrature, the regular radial integrals by Gauss quadrature, and the angular variables by deterministic angular quadrature, yielding a fully deterministic integration framework of the fractional Laplacian operator. To accurately resolve low-regularity solutions and the associated loss functional, we construct boundary-singularity-aware trial functions enriched with explicit boundary features, and propose two strategies for automatically selecting the leading exponent and evaluating the loss function from the singularity structure induced by the fractional operator, or jointly by the fractional operator and the source term. For time-dependent fractional PDEs, we design a spatiotemporally separable neural network that factorizes the time-space residual into a sum of low-dimensional temporal and spatial integrals, and we integrate this representation with an alternating neural network subspace optimization strategy for efficient training. Numerical experiments show that the proposed framework attains high accuracy on the tested benchmarks and improves substantially over existing fPINN and Monte Carlo baselines, particularly for problems with strong boundary singularities and long-time simulations. Keywords. tensor neural network, deterministic integration framework, fractional Laplacian, boundary singularity, fractional advection-diffusion.

Introduction

Fractional partial differential equations involving the fractional Laplacian have attracted sustained attention in analysis, scientific computing, and applications because they arise naturally in anomalous transport, long-range interaction, and nonlocal diffusion models; see, for example, [3, 6, 7, 19]. Among these, models involving the spatial fractional Laplacian are substantially more challenging to treat numerically than those involving only time-fractional derivatives, especially on bounded domains [30]. The main difficulties stem from the simultaneous presence of hypersingular nonlocal kernels, exterior Dirichlet constraints, and reduced boundary regularity of solutions [24, 27, 29]. These features make the design of accurate and robust numerical methods for fractional Laplacian problems on bounded domains especially demanding. ∗

School of Mathematics and Statistics & Hubei Key Laboratory of Mathematical Sciences, Central China Normal University, Wuhan 430079, China ([email protected]). † SKLMS, NCMIS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, No.55, Zhongguancun Donglu, Beijing 100190, China, and School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing, 100049, China ([email protected]). ‡ School of Mathematics and Statistics & Key Laboratory of Nonlinear Analysis & Applications (Ministry of Education), Central China Normal University, Wuhan 430079, China ([email protected]).

In this work, we first consider the fractional Poisson equation (fPE) with homogeneous Dirichlet exterior conditions: (−∆)α/2 u(x) = f (x), x ∈ Ω,

(

u(x) = 0,

(1.1)

x ∈ Ωc .

Here Ω ⊂ Rd is either a hyperrectangle or the unit ball, and 0 < α < 2. The fractional Laplacian operator (−∆)α/2 admits several non-equivalent definitions in different settings [19]. In this paper, we adopt the Riesz definition with zero exterior extension: (−∆)

α/2

u(x) = Cd,α

u(x) − u(y) dy d+α d R ∥y − x∥

with

Cd,α =

2α−1 αΓ (α/2 + d/2) . π d/2 Γ(1 − α/2)

(1.2)

We also consider the following time-dependent fractional partial differential equation (fPDE) posed on Ω with zero exterior condition [25]:  ∂ γ u(x, t)   + Lx {u(x, t)} = f (x, t), L {u(x, t)} =  t,x  ∂tγ    

x ∈ Ω, t ∈ (0, T ],

u(x, t) = 0,

x ∈ Ωc , t ∈ (0, T ],

u(x, 0) = u0 (x),

x ∈ Ω.

(1.3)

Here Lx denotes a linear spatial operator involving the fractional Laplacian and acting on u. A typical example is Lx = (−∆)α/2 + v · ∇, (1.4) where v · ∇ is the advection term associated with a prescribed velocity field v. The timefractional derivative in (1.3) is understood in the Caputo sense: ∂ γ u(x, t) 1 = ∂tγ Γ(1 − γ)

Z t 0

(t − τ )−γ

∂u(x, τ ) dτ, ∂τ

0 < γ ≤ 1.

(1.5)

When γ = 1, the Caputo derivative reduces to the classical first-order time derivative. By choosing different values of γ and α, and by including or excluding the advection term, (1.3) recovers several important classes of fractional models. In particular, when 0 < γ < 1, α = 2, and v = 0, it reduces to a time-fractional diffusion equation [18], which describes subdiffusion with memory effects. When γ = 1 and α ∈ (0, 2) with v = 0, it becomes a time-dependent space-fractional diffusion equation that captures nonlocal transport phenomena such as Lévyflight-type diffusion. When 0 < γ < 1, α ∈ (0, 2), and v ̸= 0, (1.3) becomes a time-space fractional advection-diffusion equation, incorporating both temporal memory and spatial nonlocality. Such models are widely used to describe anomalous transport phenomena [19]. Classical discretization methods, such as the finite difference method [14, 23, 31] and the finite element method [1, 2, 8, 30], remain important tools for solving fractional PDEs. However, the nonlocal character of fractional operators typically leads to dense matrices or globally coupled discrete systems, which are expensive to assemble, store, and solve. In addition, boundary singularities reduce regularity and complicate mesh design, quadrature, and error control, particularly in three dimensions. These issues motivate the development of alternative highaccuracy mesh-free methods. In recent years, neural network-based solvers have provided a complementary framework for PDEs. After the pioneering work of Lagaris et al. [15], physics-informed neural networks (PINNs) [26] and their fractional extensions have been developed for a variety of forward and inverse problems. Representative examples include fPINNs [25], bi-orthogonal fPINNs [21], spectral-fPINNs [36], and general Monte Carlo PINN formulations on irregular domains [32]. 2

For high-dimensional fPDEs, Monte Carlo discretization has become an important direction because it avoids the explicit construction of dense nonlocal matrices. From a probabilistic perspective, Sheng et al. [29] proposed an efficient Monte Carlo solver for fractional PDEs. Their approach extends the walk-on-spheres strategy by using a Feynman-Kac representation based on the Green’s function for the unit ball, enabling computations in complex domains and high dimensions. Within the PINN framework, Guo et al. [11] proposed the so-called MCfPINN method, in which the fractional Laplacian is approximated by Monte Carlo sampling after splitting the integral into integrals over a neighborhood Br0 (x) of x and its complement: (−∆)

α/2

u(x) = Cd,α

u(x) − u(y) dy + d+α y∈Br0 (x) ∥x − y∥

u(x) − u(y) dy . d+α y ∈B / r0 (x) ∥x − y∥

Splitting the integral into a singular near-field part and a regular far-field part enables tailored numerical treatment. Direct Monte Carlo sampling of the singular part, however, tends to introduce substantial variance and sampling error, limiting attainable accuracy and slowing convergence. Moreover, the performance of such sampling-based methods is often sensitive to the choice of splitting radius and to the way the singular component is truncated. To reduce sampling variance and improve radial integration accuracy, Hu et al. [13] replaced the Monte Carlo estimate of the one-dimensional singular radial integral by Gauss-Jacobi quadrature for the near-field part. Extending this strategy, we recently proposed the Quadrature-Enhanced Monte Carlo fPINN (QE-MC-fPINN) method [22] which adopts a geometry-adaptive decomposition of the fractional Laplacian based on a spatially varying near-field radius. Furthermore, the method combines deterministic radial quadrature with Monte Carlo angular sampling, and embeds the resulting operator approximation into a feature-enhanced neural network trial space tailored to low-regularity solutions. Although it improves upon representative state-of-the-art methods, QE-MC-fPINN also reveals a critical bottleneck: the dominant error is no longer the radial singularity, but the stochastic noise from the angular Monte Carlo sampling. Eliminating this noise while preserving the geometry-adaptive decomposition is the primary motivation to develop a fully deterministic approach in this paper. A second, complementary ingredient of our framework comes from recent neural network subspace and tensor neural network methods [16, 17, 34, 35, 33]. In particular, Lin et al. [18] proposed a tensor neural network subspace method to solve time-fractional partial integrodifferential equations. The method combines a power-law temporal factor tµ , Gauss-Jacobi quadrature, and an alternating optimization strategy. However, the multidimensional nonlocality of the fractional Laplacian poses significant challenges for reducing high-dimensional integrals to one-dimensional representations. To address this, we propose a spatiotemporally separable neural network (STSNN). The STSNN constructs a structured solution subspace such that the PDE residual is decomposed into decoupled temporal and spatial terms, enabling efficient representation and optimization for time-dependent nonlocal problems. Building on the MC-fPINN [11], Improved MC-fPINN [13], and especially the QE-MCfPINN method [22], the present work develops a deterministic and subspace-based framework for fractional PDEs on bounded domains. We retain the geometry-adaptive near-/far-field philosophy of the QE-MC-fPINN, but replace the remaining angular Monte Carlo sampling by deterministic angular quadrature, incorporate boundary-singularity-aware trial spaces through adaptive exponent selection, and combine these ideas with STSNN and alternating subspace optimization for time-space fractional PDEs. The main contributions of this paper are summarized as follows: 1. A fully deterministic quadrature framework for the fractional Laplacian is established, which employs a spatially adaptive near-field radius to decompose the operator into three 3

parts. The radial integrals are evaluated via Gauss-Jacobi or Gauss quadrature, and the angular integrals for all three parts are approximated using deterministic quadrature rules. 2. We construct boundary-singularity-aware neural-network trial spaces with explicit boundary features b(x)µj , and propose two strategies, BFE and BRFE, to determine the leading exponent according to the singularity structure induced by the fractional operator, or jointly by the fractional operator and the source term, respectively. These two strategies supply effective ways to treat low-regularity behaviors of the exact solutions. 3. For time-space fractional PDEs, we propose a STSNN subspace method. The resulting time-space residual factorizes into sums of products of lower-dimensional temporal and spatial integrals, and this structure is combined naturally with alternating subspace optimization. Consequently, the proposed method handles both short- and long-time simulations efficiently: the separable structure allows dense temporal Gauss quadrature without memory explosion, which is particularly advantageous for the nonlocal memory of the Caputo derivative. The remainder of this paper is organized as follows. In Section 2, we recall Gauss-Jacobi quadrature, introduce the STSNN architecture, detail the deterministic discretization of the fractional Laplacian, and construct the corresponding loss functions. In Section 3, we introduce the neural-network subspace procedure used to determine the trial basis and linear coefficients. Numerical experiments in one, two, and three dimensional cases are reported in Section 4 to demonstrate the accuracy and efficiency of the proposed methods.

Machine learning methods for fPDEs

This section presents the analytical and computational components of the fully deterministic framework developed in this work: Gauss-Jacobi quadrature, the STSNN architecture, the deterministic discretization of the fractional Laplacian, and the associated loss constructions.

2.1

Gauss-Jacobi quadrature

Gauss-Jacobi quadrature is an efficient high-precision method for approximating definite integrals. It computes integrals over [−1, 1] with the weight function (1 − x)β1 (1 + x)β2 using N nodes and weights (see [5, 28]): Z 1 −1

(1 − x)β1 (1 + x)β2 f (x)dx =

N X (β ,β ) (β ,β ) bi 1 2 f (x bi 1 2 ) + RN (f ), w

(2.1)

i=1

with β1 > −1, β2 > −1, where RN (f ) is the integration error. The quadrature formula on a general interval [a, b] is constructed as Z b a

(b − x)β1 (x − a)β2 f (x) dx ≈

N X

(β ,β )

wk 1 2 f (xk 1 2 ),

(2.2)

k=1

where the nodes and weights are obtained through affine transformation: b − a (β1 ,β2 ) a + b b − a β1 +β2 +1 (β1 ,β2 ) (β ,β ) bk bk x + , wk 1 2 = w . 2 2 2 This quadrature rule will be repeatedly used for the singular radial integrals arising in the fractional Laplacian and for the weakly singular temporal integrals in the Caputo derivative. Moreover, in the BRFE strategy (Section 2.3), it is also employed to evaluate the loss function when the residual exhibits boundary singularities.

(β ,β )

xk 1 2 =

2.2

Spatiotemporally separable neural network

This subsection introduces the STSNN architecture, which constitutes the main methodological extension beyond QE-MC-fPINN in the time-dependent setting. Its approximation properties and the computational complexity of related integral evaluations are closely connected to the tensor neural network in [33]. Sum of each product

Element-wise product of each output

……

…

……

…

Outputs of subnetworks

…

Hidden layers of subnetworks

Inputs of subnetworks

Figure 1: Architecture of the spatiotemporally separable neural network (STSNN). Black arrows denote linear (or affine) transformations. Each blue arrow indicates that the ending node is the product of all starting nodes of the same color. The final output is obtained by summing the contributions from the red arrows. The neural network is built with d + 1 subnetworks, and each subnetwork is a continuous mapping from a bounded closed set Ωi ⊂ R (i = 1, · · · , d) to Rp , which can be expressed as Φi (xi ; θi ) = ϕi,1 (xi ; θi ), ϕi,2 (xi ; θi ), · · · , ϕi,p (xi ; θi )

⊤

Φt (t; θt ) = (ϕt,1 (t; θt ), ϕt,2 (t; θt ), · · · , ϕt,p (t; θt ))⊤ , where each xi denotes the one-dimensional input, θi denotes the parameters of the i-th subnetwork, typically the weights and biases. As illustrated in Figure 1, the neural network structure adopted in this paper is composed of d fully connected neural networks (FNNs) for the spatial basis functions Φi (xi ; θi ), i = 1, 2, . . . , d, and one FNN for the temporal basis function Φt (t; θt ). To improve numerical stability, we normalize each ϕi,j (xi ), ϕt,j (t) and use the following normalized neural network structure: Ψ(x, t; c, θ) =

p X

cj ϕbt,j (t; θt )φbj (x; θx ),

(2.3)

j=1

where φbj (x; θx ) =

d Y i=1

ϕbi,j (xi ; θi ), ϕbi,j (xi ; θi ) =

ϕi,j (xi ; θi ) ϕt,j (t; θt ) , ϕbt,j (t; θt ) = . ∥ϕi,j (xi ; θi )∥L2 (Ωi ) ∥ϕt,j (t; θt )∥L2 ((0,T ])

Here, the L2 norms are precomputed using high-order Gauss quadrature on each Ωi and on the time interval (0, T ]. 5

We denote the neural network parameters as θ = {θt , θx } = {θt , θ1 , θ2 , ..., θd }. To simplify the notation, we drop the hats hereafter and write ϕi,j , ϕt,j , and φj for the normalized basis functions unless otherwise stated. Owing to its separable structure, the proposed network can be viewed as a CANDECOMP/PARAFAC decomposition in L2 (Ω × (0, T ]) [33]. For the present paper, this structure is crucial because it allows the time-space residual to be factorized into lower-dimensional temporal and spatial integrals, so dense temporal quadrature becomes feasible without sacrificing the deterministic spatial treatment of the fractional Laplacian. To demonstrate the effectiveness of solving fPDEs using the proposed neural network method, we introduce the following approximation result for functions in the space H m (Ω × (0, T ]). Theorem 2.1. [33] Assume that each Ωi is a bounded closed interval in R for i = 1, · · · , d, Ω = Ω1 × · · · × Ωd , and the function f (x, t) ∈ H m (Ω × (0, T ]) with a non-negative integer m. Then for any tolerance ε > 0, there exists a positive integer p and the corresponding neural network basis functions defined by (2.3) such that the following approximation property holds ∥f (x, t) − Ψ(x, t; c, θ)∥H m (Ω×(0,T ]) < ε. This result guarantees the approximation of smooth functions. Since the possible boundary singularities are handled by the explicit boundary feature factor in the trial functions, the remaining smoother residual is then well approximated by the neural network components.

2.3

Computation of the fractional Laplacian

This subsection describes the spatial operator discretization. Compared with the QE-MC-fPINN method, this paper uses the same geometry-adaptive splitting, while the main difference lies in that the angular Monte Carlo sampling is replaced by deterministic directional quadrature. For simplicity of notation, we write φj (x) for φj (x; θx ) and ϕt,j (t) for ϕt,j (t; θt ). 2.3.1

Discrete scheme of the fractional Laplacian

We adopt a geometry-adaptive strategy proposed in [22] that splits the fractional Laplacian based on a spatially varying radius and directional distance-to-boundary information. For any x ∈ Ω, define r0 (x) as the minimum distance from x to the boundary ∂Ω. For each spatial basis function φj , the fractional Laplacian is first divided into near-field and far-field parts: (−∆)α/2 φj (x) = Cd,α

Z ∥y−x∥2 <r0 (x)

Z ∥y−x∥2 ≥r0 (x)

φj (x) − φj (y) dy := I1,j (x) + I2,j (x). ∥x − y∥d+α 2

It is derived in [22] that the near-field integral I1,j (x) = r0 (x)−α with

Z 1

Z d−1 S+

τ 1−α

Fj (x, r0 (x)τ, ξ) dτ Jd (ξ)dξ, τ2

Fj (x, r, ξ) := 2φj (x) − φj (x + rξ) − φj (x − rξ).

d−1 Here S+ denotes the upper hemisphere in Rd . To handle the singular behavior of I1,j (x) around y = x, the Gauss-Jacobi quadrature is used:

I1,j (x) ≈ r0 (x)−α

N0 X N X

(0,1−α)

F x, r0 (x)τk (0,1−α) j

wℓ Jd (ξ ℓ )wk

ℓ=1 k=1

(0,1−α) 2

τk

, ξℓ

(0,1−α)

where {τk }N }N k=1 and {wk k=1 are Gauss-Jacobi points and weights on [0, 1]. For d = 1 2, the angular integral over S+ is discretized by an N0 -point Gauss-Legendre rule on (0, π], N0 0 producing the directional nodes {ξℓ }ℓ=1 and weights {wℓ }N ℓ=1 . For d = 3, a tensor-product Gauss-Legendre rule in the two spherical angles is used to generate the corresponding nodes and weights. Square domain

Circular domain

Figure 2: Discretization of the fractional Laplacian on 2D domains. The black curve denotes the boundary ∂Ω. For each interior point x (colored dots), a local ball Br0 (x) (x) is defined. R2 is split into three regions: (1) near-field from x to ∂Br0 (x) (x), using N0 = 10 symmetric Gaussian directions; (2) interior far-field from ∂Br0 (x) (x) to ∂Ω, with N1 = 32 uniform directions; and (3) exterior far-field from ∂Ω to infinity, with N2 = 125 uniform directions. Line segments show integration paths, and colors distinguish different evaluation points. The distance from x to the boundary along a unit direction ξ ∈ S d−1 was defined in [22] as dx (ξ) := min{ t > 0 | x + tξ ∈ ∂Ω }. Then the far-field integral I2,j is split into the part inside the domain and the part over the exterior domain. Using the quadrature rule in the radial direction of the inside part, we have I2,j (x) =

Z dx (ξ)

S d−1

r0 (x)

Z ∞ ! dx (ξ)

φj (x) − φj (x + rξ) dr Jd (ξ)dξ r1+α φj (x) − φj (x + r(x, ξ, t)ξ) φj (x) dt + Jd (ξ)dξ r(x, ξ, t)1+α α [dx (ξ)]α

Z 1



N′ X

(dx (ξ) − r0 (x))

S d−1



φj (x) − φj (x + r(x, ξ, tm )ξ) φj (x)  (dx (ξ) − r0 (x)) ≈ wm + Jd (ξ)dξ 1+α r(x, ξ, tm ) α [dx (ξ)]α S d−1 m=1 Z

Z S d−1 ′

[Q1,j (x, ξ) + Q2,j (x, ξ)] Jd (ξ)dξ, ′

N d−1 where {tm }N m=1 and {wm }m=1 are the Gauss-Legendre nodes and weights on [0, 1]. Here S denotes the unit sphere in Rd . Then I2,j is approximated by a discrete average over directions as follows:

I2,j (x) ≈

Ni 2 X X

(π/ni )d−1 Jd (ξ ik )Qi,j (x, ξik ),

i=1 k=1

with Ni = 2nd−1 the number of distributed directions. For d ≥ 2, these directions are paramei terized over the spherical coordinate domain [0, 2π] × [0, π]d−2 , while the associated weight for 7

each direction is given by (π/ni )d−1 Jd (ξ ik ). Direction vectors ξik in different dimensional cases were defined in [22]. So the fractional Laplacian is ultimately decomposed into three contributions: a singular near-field, a regular interior far-field, and an analytical exterior far-field. A schematic illustration for a two-dimensional domain is provided in Figure 2. The discretization scheme described here can, in principle, be generalized to higher dimensions. However, this work focuses on the highaccuracy neural network method for real-world physical problems with dimensions up to three. Readers interested in neural network methods for higher dimensional fPDEs are referred to MC-fPINN [11], Improved MC-fPINN [13], and QE-MC-fPINN method [22], etc. Although the present method adopts the same geometry-adaptive decomposition of the fractional Laplacian as QE-MC-fPINN, the two approaches differ fundamentally in network architecture and discretization, boundary treatment, accuracy characteristics, and optimization strategy. These differences will be expounded one by one in the subsequent text. First, in terms of network architecture and discretization, QE-MC-fPINN adopts a standard fully-connected PINN and employs Monte Carlo sampling for both the directional discretization of the fractional Laplacian and the loss evaluation, whereas the present method employs a tensor neural network that approximates a multivariate function as sums of products of univariate subnetworks, thereby enabling fully deterministic quadrature framework for both the angular discretization and the loss computation. Thus, QE-MC-fPINN is inherently stochastic and the present method is fully deterministic. For time-space fractional PDEs, the present method further employs a STSNN that factorizes the residual into independent temporal and spatial components, a structural capability unavailable in conventional PINNs. Second, for boundary treatment, both methods adopt the boundary feature enhanced (BFE) strategy (µ1 = α/2) to capture the operator-induced singularity. The present method further proposes the boundary and right-hand-side feature enhanced (BRFE) strategy: when the righthand side (RHS) function is also singular, BRFE adjusts the exponent to µ1 = α + s and employs a weighted Gauss-Jacobi quadrature, thereby handling the joint singularity mechanism with more stable convergence and higher accuracy. Third, with respect to dimensional applicability and accuracy, QE-MC-fPINN leverages the dimension-favorable scaling of Monte Carlo sampling and is well-suited to d ≥ 4, while the present method trades higher per-dimension cost via deterministic directional quadrature for substantially superior accuracy, achieving relative L2 errors one to three orders of magnitude lower on d = 1, 2, 3 benchmarks. Fourth, regarding optimization, QE-MC-fPINN optimizes all parameters jointly via gradient descent, whereas the present method decouples the optimization, that is, linear coefficients are solved by least squares and only network weights are updated by gradient descent-improving conditioning and accelerating convergence. In essence, the present method is a fully deterministic, high-accuracy solver tailored for lowto moderate-dimensional fractional PDEs, with adaptive boundary singularity treatment and subspace-based optimization that jointly distinguish it from the stochastic, Monte Carlo-based QE-MC-fPINN paradigm. 2.3.2

Analysis of the discretization scheme

In this subsection, we provide a detailed analysis of the discretization of fractional Laplacian to demonstrate the robustness and accuracy of the proposed method.

The near-field integral: Z r0 (x) 0

As derived in [22],

r0 (x)2−α ⊤ Fj (x, r, ξ) dr = − ξ H(φj )(x)ξ + O(r0 (x)4−α ), r1+α 2−α

justifying the use of Gauss-Jacobi quadrature with parameters (0, 1 − α), which is tailored for integrands with this type of weak singularity. Since Br0 (x) (x) ⊂ Ω, the local Taylor expansion is valid for 0 < r ≤ r0 (x) under standard regularity assumptions. Combined with a deterministic angular quadrature, the discretization yields high-accuracy and computationally efficient evaluation of I1,j (x). The far-field integral: The angular integration over the unit sphere S d−1 is discretized using a uniform directional grid determined by the angular resolution parameter n1 , while the radial integration is treated separately for the terms Q1,j and Q2,j . The term Q2,j admits an analytical expression, which allows a relatively dense directional sampling to be used without additional computational cost. In contrast, Q1,j is approximated by Gauss quadrature in the radial direction, and hence the principal error of I2,j arises from the discretization of the spherical integral of Q1,j . For smooth integrands the uniform angular discretization yields an error decaying as O(n−2 1 ) (see, e.g., [4, 9]), so the accuracy of I2,j can be improved by increasing n1 .

2.4

Trial spaces and loss function for fractional Laplacian

This subsection addresses the second bottleneck left open after improving operator discretization, namely the construction of trial spaces that remain accurate for boundary-singular solutions. For fractional Laplacian problems defined on bounded domains, a generic PINN ansatz is often inefficient because the residual loss and the boundary behavior are strongly coupled. As demonstrated in [10, 12], solutions of fPDEs often exhibit an α/2-order boundary singularity characterized by the asymptotic relation u(x) ≈ dist(x, ∂Ω)α/2 ureg (x),

x ∈ Ω,

where ureg is comparatively smoother. Accurate treatment of this boundary behavior is therefore indispensable for the fractional problems considered here [12, 22]. To embed the homogeneous Dirichlet boundary condition directly into the ansatz, we introduce a prescribed boundary feature function b(x) satisfying b(x) > 0, x ∈ Ω,

and

b(x) = 0, x ∈ Ωc ,

(2.4)

whose role is not only to impose the boundary condition, but also to provide an explicit mechanism for encoding the dominant boundary singularity inside the trial space. The trial function Ψ(x; c, θx ) which satisfies the homogeneous boundary conditions is then constructed as Ψ(x; c, θx ) =

p X j=1

µj

cj b(x)

d Y

ϕi,j (xi ; θi ) :=

i=1

p X

cj φj (x; θx ),

(2.5)

j=1

where {µj }pj=1 denotes the exponents of the boundary feature functions. These exponents are the crucial parameters of the boundary-singularity-aware trial family. We use them to construct the approximation space according to the asymptotic boundary structure of the target solution. α/2 To ensure the final output Ψ lies in the desired function space, we require that b(x)µj ∈ H0 (Ω). p The values of the sequence {µj }j=1 are determined by specifying only its first and last elements, µ1 and µp . The remaining parameters are then defined via linear interpolation: µj = µ1 + (j − 1)(µp − µ1 )/(p − 1), 9

j = 1, 2, . . . , p.

(2.6)

The loss function is then defined as the L2 residual of the governing equation: L(c, θx ) = ∥Lx Ψ(x; c, θx ) − f (x)∥L2 (Ω) .

(2.7)

Minimizing this loss with respect to the linear coefficients c and the network parameters θx yields the approximate solution. For problems involving fractional Laplacian, however, the integrand in (2.7) may inherit a strongly singular boundary profile unless the dominant asymptotic factor is extracted explicitly. To this end, we define for x ∈ Ω freg (x) := f (x)/b(x)s ,

with

f (x) = C ̸= 0. s x→∂Ω b(x)

lim sup

(2.8)

x∈Ω

Here the lim sup is taken over all possible interior paths approaching the boundary. With this choice, freg captures the remainder after the leading power-law boundary behavior has been removed. We assume s > −0.5 which guarantees that the resulting residual is square-integrable and hence that the weighted least-squares formulation is numerically well posed. The first exponent µ1 in (2.6) determines the dominant boundary behavior of the trial space, and therefore has a strong influence on the approximation quality. In this paper, we formulate a rule to align µ1 with the prevailing singularity mechanism of the problem. When s ≥ 0, we employ the boundary feature enhanced (BFE) strategy which has been used in [12, 22]. When −0.5 < s < 0, the boundary and RHS feature enhanced (BRFE) strategy is used. That is (a) BFE: As in [12, 22], we set µ1 = α/2 to match the leading-order boundary profile of bounded-domain problems. This operator-driven choice applies when the RHS function introduces no additional boundary singularity. In this case, the loss is evaluated using standard Gauss-Legendre quadrature over Ω. (b) BRFE: We set µ1 = α + s. This choice incorporates the singularity index s from the RHS function, aligning the trial space with the composite boundary asymptotics. This right-hand-side-aware selection applies when s < 0. Then, the loss is approximated by a tensor-product Gauss-Jacobi rule as follows: 2

L(c, θx )2 = b(x)s b(x)−s Lx Ψ(x; c, θx ) − freg (x)

≈

N1 X

···

i1 =1

Here T (y) =

p X

Nd X

d Y

id =1

k=1

! (2s,2s)

wik

(2s,2s)

T 2 x1,i1

L2 (Ω) (2s,2s)

, . . . , xd,id

cj b(y)−s Lx φj (y; θx ) − freg (y),

j=1 (2s,2s)

(2s,2s)

where {xk,ik } and {wik } are the nodes and weights of the Gauss-Jacobi quadrature. This weighted formulation is especially important when s < 0, because the deterministic quadrature can resolve the singular residual accurately and makes the BRFE trial space numerically stable in the regime where a standard residual formulation becomes unreliable. For non-homogeneous boundary conditions, a function g(x) ∈ H 1 (Ω) satisfying g(x) = ub (x), x ∈ ∂Ω, b(x) + g(x), the non-homogeneous boundis introduced. Then using the decomposition u(x) = u b(x). ary value problem is transformed to a homogeneous one for the solution u

2.5

Trial spaces and loss function for time-space fractional PDEs

For the time-space fractional PDE, the trial function Ψ in (2.5) is then redefined as b Ψ(x, t; c, θ) = Ψ(x, t; c, θ) + u0 (x),

b Ψ(x, t; c, θ) :=

p X

cj tγ ϕt,j (t)φj (x).

(2.9)

j=1

This design incorporates the factor tγ into time basis functions ϕt,j (t), following the approach proposed in [18]. Combined with the boundary-singularity-aware spatial features, this yields the time-space ansatz on which the STSNN subspace method is constructed. Substituting this decomposition into (1.3) transforms the original problem into a homogeneous initial-boundary b with RHS function fb. The corresponding loss function is then given by value problem for Ψ b b L(c, θ) = ∂ γ Ψ(x, t; c, θ)/∂tγ + Lx Ψ(x, t; c, θ) − fb(x, t)

L2 (Ω×(0,T ])

(2.10)

The BFE and BRFE can be extended to this setting in the same manner as that for fPE. We omit the details for brevity. For more information about the treatment for non-homogeneous initial/boundary value conditions, please refer to [18, 34].

fTNN: a tensor neural network subspace method for fPDEs

Based on the neural network architecture, the discretization of the fractional Laplacian, and the associated loss constructions, we propose in this section a neural network subspace method for fPDEs. We first determine the dominant boundary singularity index s based on the RHS function. According to the sign of s, we adopt either the BFE strategy (s ≥ 0) or the BRFE strategy (s < 0) to choose the parameter µ1 . The remaining parameters µj (j = 2, . . . , p − 1) are then obtained recursively from (2.6), while µp is prescribed as a fixed constant. Using these parameters, we construct a neural-network trial subspace spanned by the output and feature functions, and compute the approximate solution in this subspace in the least-squares sense. The subspace is iteratively updated during training to improve the approximation accuracy. We first present the STSNN subspace method for time-space fractional PDEs and then describe the corresponding neural network subspace method for the fPEs. Since we use tensor neural network to solve the fractional PDEs, we name this method fTNN.

3.1

Time-space fractional PDEs

Substituting (2.9) into (1.3), the original problem is then transformed into a homogeneous b with the corresponding RHS function fb. To decompose initial-boundary value problem for Ψ, the high-dimensional integrals appearing in the loss function into lower-dimensional temporal b the operator Lx applied and spatial integrals, we express the time derivative of order γ of Ψ, b b to Ψ, and f in separable forms as follows:  p p b X  ∂ γ Ψ(x, t) X ∂ γ (tγ ϕt,j (t; θt ))   = c φ (x; θ ) := cj L1t,j (t)L1x,j (x),  j j x  γ  ∂tγ ∂t  j=1 j=1     p p  X X b Lx Ψ(x, t) = cj tγ ϕt,j (t; θt ) Lx φj (x; θx ) := cj L2t,j (t)L2x,j (x),   j=1 j=1     p 1  X     fb(x, t) := fbt,r (t)fbx,r (x),  r=1

(3.1)

where p1 denotes the number of terms in the spatiotemporally separable form of the source. After that, we determine the boundary singularity index s from (2.8). If s ≥ 0, we adopt the BFE strategy, setting µ1 = α/2 and using Gauss quadrature for the trial function and the loss function. If s < 0, we adopt the BRFE strategy, setting µ1 = α + s and using GaussJacobi quadrature with parameters (2s, 2s). After prescribing µp , we define the neural network subspace as Vp′ (θ) := span {tγ ϕt,j (t; θt ) φj (x; θx ), j = 1, . . . , p} , with φj (x; θx ) = b(x)µj

d Y

ϕi,j (xi ; θi ).

i=1

bp + u0 with u bp ∈ Vp′ , by The approximate solution is then sought in the form up = u bp ∈ Vp′ such that minimizing the residual in the least-squares sense. More precisely, we seek u bp , Lt,x vbp ) = (fb, Lt,x vbp ), (Lt,x u

∀ vbp ∈ Vp′ ,

(3.2)

where (·, ·) denotes the L2 inner product over the space-time domain Ω × (0, T ]. After the ℓ-th training step, the trial function Ψ(x, t; c, θ(ℓ) ) belongs to the subspace n

(ℓ)

Vp′ (θ(ℓ) ) := span tγ · ϕt,j (t; θt ) · φj (x; θx(ℓ) ), j = 1, . . . , p , (ℓ)

(ℓ)

where ϕt,j (t; θt ) and φj (x; θx ) are defined by (2.9) and (2.5), respectively. Once the coeffibp is obtained, and the subspace Vp′ (θ (ℓ+1) ) is cients c are computed, the approximate solution u updated through the loss optimization. Assembling the discrete system on Vp′ (θ(ℓ) ) and exploiting the spatiotemporal separability of the basis and source representations yields the following stiffness matrix and load vector for the BFE (s ≥ 0) and BRFE (s < 0) strategies:  2 2 X X   j k j k   , · L , L L , L x,n x,m t,n t,m   x t

s ≥ 0,

j=1 k=1 A(ℓ) m,n = X 2 X 2

(ℓ) Bm =

(3.3)

 j k −s k −s j   · b L , b L , s < 0, L , L  t,n x,n x,m t,m  t x,b2s j=1 k=1 p 2 1 X X     fbt,r , Ljt,m · fbx,r , Ljx,m , s ≥ 0,   t x r=1 j=1

p1 X 2  X  −s b −s j  bt,r , Lj  f · b f , b L , x,r  x,m t,m  x, b2s t

(3.4)

s < 0,

r=1 j=1

for 1 ≤ m, n ≤ p. Here, (·, ·)t and (·, ·)x denote the temporal and spatial inner products, respectively, both approximated by Gauss quadrature. For the BRFE strategy, the boundary singular factor is absorbed into the spatial terms via the factor b−s , which naturally leads to the weighted spatial inner product Z

(u, v)x, b2s := Qd

u(x) v(x) b(x)2s dx.

(a ,b ) i=1 i i

(3.5)

The weighted integral in (3.5) is evaluated by Gauss-Jacobi quadrature with parameters (2s, 2s). The parameters θ(ℓ) determine the trial subspace Vp′ (ℓ) , while the coefficient vector c(ℓ) determines the approximation within that subspace. Consequently, the optimization procedure can be naturally divided into two substeps. 12

First, with the neural network parameters θ(ℓ) fixed, the optimal coefficient vector c(ℓ+1) is obtained in the least-squares sense from (3.2) by solving the linear system A(ℓ) c(ℓ+1) = B (ℓ) .

(3.6)

Second, with the coefficient vector c(ℓ+1) fixed, the neural network parameters θ(ℓ+1) are updated by minimizing the loss function

L(ℓ+1) (c(ℓ+1) , θ(ℓ) ) = c(ℓ+1)

⊤

A(ℓ) c(ℓ+1) − 2 B (ℓ)

⊤

c(ℓ+1) + (fb, fb)

(3.7)

using gradient-based optimization methods such as Adam or L-BFGS. Thus, for fixed θ(ℓ) , minimizing the discrete loss is equivalent to solving the normal equations (3.6). Algorithm 1: fTNN method for the time-space fractional PDE b using (2.9), with RHS fb. 1. Transform (1.3) into a homogeneous problem for Ψ b Lx Ψ, b and fb in spatiotemporally separable form, see (3.1). 2. Express ∂tγ Ψ,

3. Select the strategy. Determine the boundary singularity index s from (2.8). • If s ≥ 0 (BFE), set µ1 = α/2 and use Gauss quadrature points. • If s < 0 (BRFE), set µ1 = α + s and use Gauss-Jacobi quadrature with (2s, 2s). 4. Initialize. Construct the initial trial function Ψ(x, t; c(0) , θ(0) ) defined in (2.5) with the chosen µ1 and µp . Set the maximum number of iterations M and initialize ℓ = 0. 5. Repeat until convergence or ℓ = M : (a) Using µ1 and µp to construct the trial subspace n

(ℓ)

Vp′ (θ(ℓ) ) := span tγ ϕt,j (t; θt ) φj (x; θx(ℓ) ), j = 1, . . . , p . (b) Assemble the stiffness matrix A(ℓ) and load vector B (ℓ) as in (3.3) and (3.4). (c) Solve the linear system A(ℓ) c = B (ℓ) and set c(ℓ+1) = c. (d) Update θ(ℓ+1) by minimizing the loss function (3.7), with c(ℓ+1) fixed. (e) Set ℓ ← ℓ + 1.

The whole procedure is summarized in Algorithm 1. This algorithm can be interpreted as an alternating optimization procedure for the unknown parameters c and θ in the trial function (2.5): the optimal coefficients c are determined in the least-squares sense, while the neuralnetwork parameters are updated through the loss minimization.

3.2

Fractional Poisson equation

For the fPE defined in (1.1), we introduce the neural-network trial subspace Vp (θx ) := span {φj (x; θx ), j = 1, . . . , p} ,

(3.8)

where the approximate solution is found by minimizing the residual in the least-squares sense. More precisely, we seek up ∈ Vp (θx ) such that

(−∆)α/2 up , (−∆)α/2 vp

= f, (−∆)α/2 vp 13

∀ vp ∈ Vp (θx ).

(3.9)

Algorithm 2: fTNN method for the fractional Poisson equation 1. Select the strategy. Determine the boundary singularity index s from (2.8). • If s ≥ 0 (BFE), set µ1 = α/2 and use Gauss quadrature points. • If s < 0 (BRFE), set µ1 = α + s and use Gauss-Jacobi quadrature with (2s, 2s). 2. Initialize. Construct the initial trial function Ψ(x; c(0) , θ(0) ) defined in (2.5) with the chosen µ1 and µp . Set the maximum number of iterations M and initialize ℓ = 0. 3. Repeat until convergence or ℓ = M : (a) Define the trial subspace (3.8) using µ1 and µp n

Vp(ℓ) := span φj (x; θx(ℓ) ), j = 1, . . . , p . (b) Assemble the stiffness matrix A(ℓ) and load vector B (ℓ) using the BFE or BRFE strategy, see (3.10) and (3.11). (c) Solve the linear system A(ℓ) c = B (ℓ) and set c(ℓ+1) = c. (ℓ+1)

(d) Update θx

(ℓ)

by minimizing the loss function L(ℓ+1) (c(ℓ+1) , θx ) with c(ℓ+1) fixed.

(e) Set ℓ ← ℓ + 1.

(ℓ)

After the ℓ-th training step, the neural network Ψ(x; c, θx ) belongs to the subspace n

Vp(ℓ) := span φj (x; θx(ℓ) ), j = 1, . . . , p . (ℓ)

By assembling the discrete system over Vp , the stiffness matrix and load vector corresponding to the BFE and BRFE formulations can be written in the unified form  (ℓ) (ℓ)   (−∆)α/2 φn , (−∆)α/2 φm , x A(ℓ) m,n =  −s (ℓ) (ℓ)  b (−∆)α/2 (φn ), b−s (−∆)α/2 (φm ) , x, b2s  (ℓ)   f, (−∆)α/2 φm , s ≥ 0, x (ℓ) Bm = (ℓ)   b−s f, b−s (−∆)α/2 (φm ) , s < 0, 2s

s ≥ 0, s < 0,

(3.10)

(3.11)

x, b

for 1 ≤ m, n ≤ p. Here, (·, ·)x, b2s is the weighted spatial inner product defined in (3.5). (ℓ)

(ℓ)

The parameters θx determine the trial subspace Vp , while the coefficient vector c(ℓ) determines the approximation within this subspace. Accordingly, Algorithm 2 can be interpreted as an alternating optimization procedure for fPEs. Since its implementation is analogous to Algorithm 1, we omit the details here for brevity.

Numerical examples

This section check the performance of the proposed method on a sequence of one-, two-, and three-dimensional benchmarks for which the deterministic spatial quadrature is computationally meaningful. All the experiments are conducted on NVIDIA GeForce RTX 4090 D GPUs. 14

The following two types of errors between the approximate solution Ψ(x, t; c, θ) and the exact solution u are used to measure the convergence behavior and accuracy of the examples in this section: • Relative L2 error eL2 :=

∥u − Ψ(x, t; c∗ , θ∗ )∥L2 (Ω×(0,T ]) . ∥u∥L2 (Ω×(0,T ])

This norm is evaluated by high precision numerical quadrature. • Relative L2 test error q

k k k k 2 k=1 (Ψ(x , t ; c, θ) − u(x , t ))

etest :=

q PK

k k k=1 (u(x , t ))

where the test points {(xk , tk )} are placed on a uniform grid of 300d+1 points over Ω×(0, T ] for 1D and 2D, and on a uniform 304 grid for 3D. For time-independent problems, the same definitions for the errors are used with the L2 norm taken over Ω and test points {xk } uniformly distributed in Ω. Table 1: Neural network architectures and quadrature point configurations for different dimensions using BFE (Gauss points) and BRFE (Gauss-Jacobi points) d 1 2 3

Each FNN [1, 50, 50, 50] [1, 50, 50, 10] [1, 10, 10, 5 ]

Num of quadrature points BFE

BRFE

128 [16, 16] [10, 10, 10]

48 [16, 16] [10, 10, 10]

In the following numerical experiments, all FNNs consist of two hidden layers with the Tanh activation function applied to each layer. To accurately compute the integral in the loss function Q (2.7) over the domain di=1 [ai , bi ], we employ the quadrature point configurations specified in Table 1. For the rectangular domain, we define b(x) = max

d Y

(xi − ai )(bi − xi ), 0 ,

i=1

which vanishes on the entire boundary, satisfying (2.4). Remark 4.1. Near a face of the rectangular domain, e.g., x1 = a1 , b(x) ∼ (x1 − a1 ) C(x′ ), where C(x′ ) is smooth and strictly positive on that face. Near edges or corners, however, b(x) decays faster than the minimal distance from boundary: near a corner where xi = ai , Q iq i = i1 , i2 , · · · , iq , it holds that b(x) ∼ i=i (xi − ai ). Consequently, on lower-dimensional 1 boundary strata, b(x)α/2 = o dist(x, ∂Ω)α/2 . This property makes b(x) a natural factor in the trial functions: it faithfully reproduces the leading singular behavior of solutions near faces, while automatically providing higher-order decay near edges and corners, thereby preventing over-singularity on lower-dimensional boundary strata. The remaining singular or complex features of the solution are left to be learned by the neural network components in the trial space. 15

In all experiments we set µp = µ1 + 0.5, which was found to provide a sufficiently rich approximation space for the tested parameter ranges. The number of basis functions p is set to the output dimension of each subnetwork (e.g., 50 for 1D, 10 for 2D, 5 for 3D). During training, the neural network (for the fPE) or the STSNN (for the fPDE) is optimized using Adam with learning rate 0.003 for 1000 epochs, then fine-tuned by L-BFGS with initial step size 0.1 for 200 iterations (stationary) or 400 iterations (time-dependent).

4.1

Fractional Poisson equations

In this subsection, the fTNN refers to Algorithm 2. 4.1.1

One dimensional case

For the numerical discretization of the fractional Laplacian, we employ N = 10 Gauss-Jacobi quadrature points for the near-field radial integral and N0 = 100 Gauss quadrature points for the far-field radial integral. We first consider the one-dimensional fPE on the interval (0, 1): (−∆)α/2 u(x) = f (x),

x ∈ (0, 1),

(4.1)

subject to homogeneous Dirichlet boundary conditions u(0) = u(1) = 0. The forcing terms with the fabricated smooth solutions are given respectively by (

f (x) =

(ℓα (x, 2) − 2ℓα (x, 3) + ℓα (x, 4)) /(2 cos (πα)), (ℓα (x, 3) − 3ℓα (x, 4) + 3ℓα (x, 5) − ℓα (x, 6)) /(2 cos (πα)),

where ℓα (x, n) =

uexact = x2 (1 − x)2 , uexact = x3 (1 − x)3 ,

Γ(n + 1) xn−α + (1 − x)n−α . Γ(n + 1 − α)

Since s = 2 − α > 0 and s = 3 − α > 0 for the two cases, we use the BFE strategy. The results in Table 2 show fTNN provides a clear accuracy advantage over fPINN in this smooth setting. Table 2: Relative L2 test error etest for (4.1) with smooth solutions. u = x2 (1 − x)2

u = x3 (1 − x)3

fPINN

fTNN

fPINN

fTNN

1.5 1.9

3.79e-2 2.98e-3

2.92e-7 5.70e-7

1.67e-02 2.83e-04

6.94e-7 1.38e-6

We then consider a compactly supported function that is the sum of two bump functions with different exponents: uexact (x) = (1 − x2 )β+1 + (1 − x2 )β+2 ,

for x ∈ [−1, 1],

β1 ⩽ β2 .

The regularity of uexact is determined by the smaller exponent β1 , say C β1 −1,1 (R), C ⌊β1 ⌋,β1 −⌊β1 ⌋ (R),

(

uexact ∈ For βi ̸= α/2,

β1 ∈ N, β1 ∈ / N.

(−∆)α/2 (1 − x2 )β+i ∼ Ci (1 − x2 )βi −α , 16

|x| → 1− .

(4.2)

Table 3: Summary of parameter regimes and strategy selection. Case

Conditions

Strategy

µ1

Case I Case II Case III Case IV Case V

β1 = α/2, β2 = α/2 or β2 ≥ α β1 = α/2, α/2 < β2 < α β1 < α/2 (any β2 ≥ β1 ) α/2 < β1 < α (any β2 ≥ β1 ) β1 ≥ α (any β2 ≥ β1 )

BFE BRFE BRFE BRFE BFE

α/2 β2 β1 β1 α/2

Since β1 ≤ β2 , the smallest exponent dominates, yielding f (x) ∼ const · (1 − x2 )β1 −α ,

β1 < β2 , β1 ̸= α/2.

In the special case β1 = β2 = α/2, we have √ α/2 (−∆)α/2 (1 − x2 )+ = 2α Γ (α/2 + 1/2) Γ (α/2 + 1) / π, which is a non-zero constant, so f remains bounded and non-zero near the endpoints. Table 4: Errors of fTNN to solve (4.1) with exact solution (4.2). BRFE (s < 0)

BFE (s ⩾ 0)

Case: (α, β1 , β2 )

e L2

etest

Case: (α, β1 , β2 )

eL2

etest

II: (0.40, 0.20, 0.30) II: (1.60, 0.80, 1.20) II: (1.10, 0.55, 0.88) II: (1.90, 0.95, 1.60) III: (0.50, 0.05, 0.10) III: (0.40, 0.08, 0.30) III: (0.20, 0.08, 0.30) III: (0.70, 0.30, 1.30) IV: (1.30, 0.92, 1.30) IV: (1.42, 1.20, 3.30) IV: (1.90, 1.50, 1.80) IV: (1.90, 1.80, 2.20)

6.266e-5 3.492e-5 1.922e-5 2.932e-4 8.875e-5 5.980e-5 5.465e-5 1.804e-5 6.899e-6 8.611e-7 2.135e-7 2.182e-7

6.595e-5 3.657e-5 1.902e-5 7.300e-6 3.074e-3 3.122e-3 5.391e-5 5.391e-5 6.436e-6 8.577e-7 2.158e-7 2.180e-7

I: (0.20, 0.10, 0.10) I: (1.80, 0.90, 0.90) I: (1.60, 0.80, 1.90) I: (0.90, 0.45, 0.90) V: (0.10, 0.10, 0.20) V: (0.20, 0.20, 0.30) V: (0.40, 0.41, 1.20) V: (0.40, 1.20, 1.90) V: (1.10, 1.15, 1.61) V: (1.80, 1.90, 2.30) V: (1.70, 1.99, 2.90) V: (1.90, 2.30, 3.20)

3.156e-5 8.615e-7 2.859e-6 3.854e-5 1.501e-5 2.889e-5 1.885e-5 7.622e-7 1.006e-6 1.596e-7 4.110e-7 6.370e-7

3.188e-3 8.389e-7 2.852e-6 3.854e-5 2.087e-3 4.785e-5 1.910e-5 4.270e-6 1.009e-6 1.623e-7 4.027e-7 6.348e-7

Consequently, the singularity index defined in (2.8) is s = β1 − α, except in the special case β1 = β2 = α/2 where s = 0. The appropriate strategy, BFE for s ≥ 0 and BRFE for s < 0, then follows as summarized in Table 3. Errors of fTNN are listed in Table 4 for five regimes with different parameter settings. As is shown, when s < 0, BRFE keeps the errors uniformly small by aligning the leading basis exponent with the composite singularity induced by the operator and the forcing term; when s ≥ 0, BFE recovers the canonical operator-driven boundary profile and achieves a comparable level of accuracy. In this sense, fTNN with BFE/BRFE achieves a further performance improvement compared with QE-MC-fPINN. The solution profiles in Figure 3 and the error histories in Figure 4 show that fTNN with two strategies remains stable throughout training, with the largest errors concentrated near the endpoints where the singularity is the strongest. These results validate the robustness and accuracy of fTNN across the full range of fractional orders α and exponent combinations (β1 , β2 ).

2e-05

6e-02

1.5

1e-05

4e-02

1.0

0.5 0.0

0.0

0.2

Exact Prediction 0.4 0.6 0.8

1.0

2e-02

0.5

0e+00

0.0

0.2

0.4

0.6

0.8

1.0

(a) I: (α, β1 , β2 ) = (0.20, 0.10, 0.10), BFE. 2e-01

2.0

Absolute Error

1.0 0.5 0.0

0.0

0.2

Exact Prediction 0.4 0.6 0.8

1.0

8e-02

1.0

5e-02

0.5

3e-02 0.2

0.4

0.6

0.8

4e-04 3e-04 3e-04 3e-04 2e-04 2e-04 1e-04 5e-05

0.5 0.0

0.0

0.2

Exact Prediction 0.4 0.6 0.8

1.0

0.0

0.2

Exact Prediction 0.4 0.6 0.8

2.0

1.0 0.5

0.2

0.4

0.6

0.8

1.0

0.8

(e) V: (α, β1 , β2 ) = (0.20, 0.20, 0.30), BFE.

2e-06 1e-06 5e-07 0.0

0.2

( , 1, 2) = (1.9, 1.8, 2.2) 0.4 0.6 0.8

1.0

4e-07

1.5

0.0

0.4

2e-06

0e+00

1.0

0.2

(d) IV: (α, β1 , β2 ) = (1.90, 1.80, 2.20), BRFE.

( , 1, 2) = (0.2, 0.2, 0.3)

Absolute Error

1.5

0.0

1.0

0.0

3e-06

1.5

0.0

5e-06

3e-06

2.0

1e-01

0e+00

1e-05

0e+00

1.0

1.5

0.2

( , 1, 2) = (1.9, 0.95, 1.6)

(b) II: (α, β1 , β2 ) = (1.90, 0.95, 1.60), BRFE.

( , 1, 2) = (0.5, 0.05, 0.1)

1e-01

0.0

Exact Prediction 0.4 0.6 0.8

Absolute Error

1.0

Absolute Error

1.5

( , 1, 2) = (0.2, 0.1, 0.1)

Absolute Error

2.0

0.0

1.0

0.0

0.2

Exact Prediction 0.4 0.6 0.8

2e-07 1e-07

0e+00

1.0

3e-07

0.0

0.2

( , 1, 2) = (1.8, 1.9, 2.3) 0.4 0.6 0.8 1.0

(f) V: (α, β1 , β2 ) = (1.80, 1.90, 2.30), BFE.

1e+00

1e-01

1e-02

Relative L2 error (eL2 )

Figure 3: Numerical results of fTNN to solve (4.1) with exact solution (4.2).

1e-02 1e-03

Adam

L-BFGS

1e-04 1e-05 II : (α, β1 , β2 ) = (0.40, 0.20, 0.30) III: (α, β1 , β2 ) = (0.50, 0.05, 0.10)

1e-06

1e-03 1e-04

Adam

V: (α, β1 , β2 ) = (0.20, 0.20, 0.30)

1e-06

V: (α, β1 , β2 ) = (1.70, 1.99, 2.90)

IV: (α, β1 , β2 ) = (1.90, 1.50, 1.80)

1e-07

L-BFGS

1e-05

I : (α, β1 , β2 ) = (1.80, 0.90, 0.90)

1e-07

Training steps (102 )

Figure 4: Relative L2 errors of fTNN with BFE (left) and BRFE (right) strategies to solve (4.1) with exact solution (4.2).

4.1.2

FPEs on the unit square/cube

We first consider the fPE (1.1) with an exact solution constructed as

uexact (x) =

2 Q 2 P

   

k=1 i=1 2 Q 3 P

  100

xi − x2i

k=1 i=1

β k

xi − x2i

x ∈ [0, 1]2 ,

β k

(4.3)

x ∈ [0, 1]3 .

where the exponents satisfy β1 ⩽ β2 . Due to the complexity of deriving an explicit closed-form expression for the RHS function f (x) = (−∆)α/2 uexact (x), we resort to a numerical approximation to obtain a high-accuracy RHS function. To this end, we employ dense quadrature rules whose parameters are specified in Table 5. Table 5: Parameters setting of fTNN to solve (1.1) with f ≈ (−∆)α/2 uexact (x). (−∆)α/2 uexact

RHS f

d−1 Angular directions on S+ Gauss-Jacobi radial points

10 10

120 8

50 200

300 32

Angular resolution for Q1,j Angular resolution for Q2,j Gauss radial points for Q1,j

32 250 100

16 35 16

600 1000 200

35 35 40

Integral

Parameters

Description

I1,j (x)

N0 N

I2,j (x)

n1 n2 N′

Although an explicit expression for the RHS function is unavailable, its singular behavior near the domain boundary can be characterized theoretically. This characterization is crucial for integration of the loss, as it enables the incorporation of the appropriate singular weight into the approximation space. By linearity, f = f1 + f2 with

fk = (−∆)α/2 cd b(x)βk ,

b(x) =

d Y

xi (1 − xi ),

i=1

on Ω = (0, 1)d . Since β1 ≤ β2 , the singularity of f is dominated by the term with the smaller exponent, i.e., f ∼ f1 as x → ∂Ω. Near a face, say xi = 0, we have b(x) ∼ xi Bi (x′ ), where x′ Q denotes the remaining coordinates and Bi (x′ ) = j̸=i xj (1 − xj ) is smooth and non-zero on the face. Consequently, b(x)β1 ∼ xβi 1 (Bi (x′ ))β1 . A classical result for the fractional Laplacian (see [27]) yields the asymptotic expansion

(−∆)α/2 xβi 1 ϕ(x′ ) ∼ C(β1 , α) ϕ(x′ ) xβi 1 −α ,

xi → 0+ ,

where C(β1 , α) is a constant. Applying this with ϕ = Biβ1 gives f (x) ∼ C(β1 , α) Bi (x′ )β1 xβi 1 −α ,

β1 ≤ β2 .

Analogous expansions hold near other faces and edges, with the distance to the boundary raised to the appropriate power. Thus, near any face, f (x) behaves like a constant (depending on the tangential coordinates) times dist(x, ∂Ω)β1 −α . To extract this maximal singular behavior, we require that the limit defined in (2.8) be finite and non-zero. Substituting the asymptotic forms shows that the exponent s = β1 − α. With this choice, the ratio tends to a non-zero smooth function on the faces, while it vanishes near edges and corners, hence the lim sup is indeed attained on the faces and equals a positive constant. 19

Based on the above asymptotic analysis, the selection of the strategy (BFE or BRFE) in the multi-dimensional setting follows a simpler criterion compared to the one-dimensional case. Unlike in 1D, where the special case β1 = α/2 yields a constant RHS function for the Q term (1 − x2 )α/2 , in higher dimensions the product structure i (xi − x2i )βk does not produce a constant fractional Laplacian for any βk . This is because the cross-terms in the product lead to a genuine boundary singularity in f even when β1 = α/2, justifying the exclusive use of β1 − α as the exponent s in the multidimensional setting. Consequently, the method choice is determined solely by the relative magnitudes of β1 and α: • When β1 < α (with any β2 ≥ β1 ), f (x) exhibits a genuine singularity dist(x, ∂Ω)β1 −α near the boundary, necessitating the enriched approximation space of BRFE with the singular weight µ1 = β1 . • When β1 ≥ α (with any β2 ≥ β1 ), f (x) remains bounded (or even smooth) on Ω, and the BFE with µ1 = α/2 suffices. Table 6: Errors of fTNN to solve (1.1) on the unit square/cube. 2D: unit square BRFE (s < 0), N1 = 64 (n1 = 32) BFE (s ⩾ 0), N1 = 64 (n1 = 32) α

(β1 , β2 )

eL2

etest

(β1 , β2 )

eL2

etest

0.50 0.60 1.00 1.20 1.60 1.90 1.96

(0.08,0.20) (0.15,2.20) (0.55,0.80) (0.80,1.20) (1.50,1.79) (1.65,1.85) (1.45,2.30)

3.030e-5 3.598e-5 6.233e-5 7.202e-5 3.160e-5 1.070e-5 4.088e-5

6.353e-3 1.237e-3 6.226e-5 7.203e-5 3.160e-5 1.070e-5 4.088e-5

0.20 0.10 0.80 1.00 1.20 1.80 1.76

(0.20,0.40) (0.51,0.80) (0.90,1.20) (1.21,1.53) (1.80,1.90) (1.90,1.95) (1.70,2.20)

3.497e-5 5.246e-5 7.308e-5 8.097e-5 7.726e-5 2.446e-5 3.470e-5

2.771e-4 5.306e-5 7.150e-5 7.850e-5 7.631e-5 2.362e-5 3.454e-5

3D: unit cube BRFE (s < 0), N1 = 512 (n1 = 16) BFE (s ⩾ 0), N1 = 512 (n1 = 16) α

(β1 , β2 )

eL2

etest

(β1 , β2 )

eL2

etest

0.40 1.00 1.10 1.20 1.60 1.90 1.96

(0.08,0.20) (0.55,0.80) (0.65,2.30) (0.90,1.70) (1.20,1.75) (1.65,1.85) (1.65,2.45)

9.335e-4 4.838e-3 5.813e-3 9.087e-4 3.369e-3 2.333e-4 1.495e-4

2.936e-2 3.389e-3 4.128e-3 7.887e-4 3.174e-3 2.306e-4 1.488e-4

0.20 0.80 1.10 1.20 1.40 1.70 1.90

(0.20,0.30) (0.80,1.80) (1.20,1.99) (1.30,1.70) (1.50,1.85) (1.86,2.88) (1.90,2.80)

7.941e-4 6.025e-4 5.747e-4 5.687e-4 4.874e-4 3.635e-4 1.005e-4

1.118e-3 5.931e-4 5.689e-4 5.646e-4 4.858e-4 3.637e-4 1.006e-4

The relative errors eL2 and etest of fTNN to solve (1.1) on the unit square/cube are listed in Table 6. Both types of errors lie below 10−4 for most cases, confirming the high accuracy of fTNN in two and three dimensions. We also examine how the accuracy-cost trade-off depends on the directional resolution n1 . The average computational time per epoch of fTNN is drawn against n1 in Figure 5, which shows a deterministic trend. Relative L2 errors in both two and three dimensional cases are shown in Tables 7 and 8. The observed error reduction is consistent with the O(n−2 1 ) decay predicted by the spherical integration theory in Section 2.3.2, confirming that the angular resolution dominates accuracy once the radial singularities are properly resolved. This trend is precisely the one targeted by the present paper: after the geometry-adaptive decomposition has removed the difficulty caused by radial singularity, the further improvement comes from 20

0.25

0.20

0.15

0.10

0.05

1 8

0.5

Adam of per epoch (s)

L-BFGS of per epoch (s)

Adam of per epoch (s)

Adam L-BFGS

0.30

Adam L-BFGS

0.4 6

0.3

0.2

0.1

2 2

L-BFGS of per epoch (s)

0.35

Figure 5: The average computational time per epoch of fTNN for solving (1.1) with direction resolution n1 on the unit square (left) and the unit cube (right). refining a deterministic angular discretization rather than from increasing the number of Monte Carlo samples as used in QE-MC-fPINN. Table 7: Relative L2 errors and time of fTNN solving the 2D fPE with different n1 . (α, β1 , β2 )

n1 = 4

n1 = 8

n1 = 16

n1 = 32

n1 = 64

eL2

(1.80, 1.90, 1.95) (0.80, 0.90, 1.20) (1.60, 1.50, 1.79)

5.130e-4 3.545e-3 3.545e-3

6.639e-5 3.696e-4 3.696e-4

2.908e-5 1.273e-4 1.273e-4

2.446e-5 7.308e-5 7.308e-5

1.257e-5 2.297e-5 2.297e-5

Time (s)

t1 t2

0.033 1.111

0.052 1.111

0.093 1.915

0.173 3.554

0.339 6.888

t1 /t2 is the average time per epoch of Adam/L-BFGS.

Table 8: Relative L2 errors and time of fTNN solving the 3D fPE with different n1 . n1 = 2

n1 = 4

n1 = 8

n1 = 16

eL2

(1.60, 1.20, 1.75) (1.96, 1.65, 2.45) (1.70, 1.86, 2.88)

7.970e-2 1.523e-2 4.431e-2

9.383e-3 1.317e-3 7.411e-3

4.268e-3 3.590e-4 1.850e-3

3.369e-3 2.333e-4 3.635e-4

Time (s)

t1 t2

0.054 1.122

0.066 1.370

0.119 2.426

0.329 6.735

1.0

0.7

0.8

0.3

0.4

0.20

0.6

0.2

1.0

0.25

0.4

0.6

0.4 0.3

Predicted u(x1, x2) (f = 1 and = 1.2)

0.8

0.5

0.6

1.0 0.5

0.6

0.8

Predicted u(x1, x2) (f = 1 and = 0.6)

0.15

Predicted u(x1, x2) (f = 1 and = 1.8)

0.10

0.8

0.08

0.6

0.06

Predicted u(x1, x2) (f = 1 and = 0.4)

1.0

(α, β1 , β2 )

0.4

0.10

0.4

0.04

0.2

0.05

0.2

0.02

0.00

0.0 0.0

0.2

0.2 0.0 0.0

0.1

0.2

0.4

(a) α = 0.40

0.6

0.8

1.0

0.2 0.0 0.0

0.1

0.2

0.4

(b) α = 0.60

0.6

0.8

0.0 0.0

1.0

0.2

0.4

0.6

0.8

1.0

0.2

0.4

0.6

0.8

1.0

0.00

(d) α = 1.80

Figure 6: Solutions of fTNN to solve the 2D fPE with RHS function f = 1. We then use fTNN to solve the fPE (1.1) with the RHS function f = 1 (clearly, s = 0) and different α in two dimensional case. Since the exact solutions are not known, we show the numerical solutions in Figure 6 and corresponding loss curves in Figure 7. 21

3e-02

1e-01 α = 0.40 α = 0.60

α = 1.20 α = 1.80

Loss

1e-02

1e-02 Adam

L-BFGS

1e-03 Adam

L-BFGS

2e-03

4e-04

Training steps (102 )

Figure 7: Loss curves of fTNN to solve the 2D fPE with RHS function f = 1. 4.1.3

FPEs on the unit ball

We consider fPE (1.1) on the unit ball, for which the corresponding RHS function could be analytically derived for a given solution. For brevity, we present the two-dimensional case, while the extension to three dimensions is analogous. In two dimensional case, it holds that

(−∆)α/2 1 − ∥x∥22

1+α/2

= 2α Γ(α/2 + 2) Γ(α/2 + 1) 1 − (α/2 + 1)∥x∥22 .

Using polar coordinates x = (ρ cos v, ρ sin v) with ρ ∈ [0, 1) and v ∈ [0, 2π), the RHS function reads f (ρ, v) = 2α Γ(α/2 + 2) Γ(α/2 + 1) 1 − (α/2 + 1)ρ2 . Observe that

lim f (ρ, v) = −2α−1 α Γ(α/2 + 2) Γ(α/2 + 1) ̸= 0.

ρ→1−

Consequently, the limit

f (ρ, v)

lim

ρ→1− (1 − ρ2 )s

is finite and nonzero if and only if the denominator tends to a constant, i.e., s = 0. Hence, we employ the BFE strategy for this problem. The approximate solution is expanded in terms of neural network basis functions: Ψ(ρ, v) =

p X

cj φj (ρ, v).

j=1

Then

(−∆)α/2 φj (ρ, v) = C2,α (I1,j (ρ, v) + I2,j (ρ, v)) .

The near-field integral is approximated as I1,j (ρ, v) ≈ r0 (ρ)−α

N0 X N X

(0,1−α) , ηl ) (0,1−α) Fj (ρ, v, r0 (ρ)τk , (0,1−α) 2 (τk )

wl wk

l=1 k=1

where r0 (ρ) = 1 − ρ, Fj (ρ, v, r, η) = 2φj (ρ, v) − φj (ρ′ , v ′ ) − φj (ρ′′ , v ′′ ), with ρ′ =

ρ2 + r2 + 2ρr cos(η − v),

and v ′ = arctan

ρ sin v + r sin η , ρ cos v + r cos η

ρ′′ =

ρ2 + r2 − 2ρr cos(η − v),

v ′′ = arctan 22

ρ sin v − r sin η . ρ cos v − r cos η

The far-field integral is then approximated as I2,j (ρ, v) ≈

2ni 2 X X π i=1 k=1

Qi,j (ρ, v, ηik ),

with ηik being discrete directional angles. Q1,j and Q2,j are defined as Q1,j (ρ, v, η) ≈ (dx (η) − r0 (ρ, v))

N0 X

m=1

φj (ρ, v) − φj (ρ′ , v ′ ) , (r(η, tm ))1+α

Q2,j (ρ, v, η) =

φj (ρ, v) , α (dx (η))α

where r(η, t) := r0 (ρ, v) + (dx (η) − r0 (ρ, v)) t and the distance dx (η) from the point (ρ, v) to the boundary along direction η is given by dx (η) = −ρ cos(η − v) +

1 − ρ2 sin2 (η − v).

The above formulas provide the numerical approximation of fractional Laplacian in twodimensional polar coordinates. All computations are performed in polar coordinates, with the neural network basis functions taking polar coordinates as inputs. The experiments on the unit ball provide a well-controlled benchmark for evaluating fTNN against several state-of-the-art methods based on neural network. The availability of analytic solutions exhibiting canonical boundary singularities, combined with the spherical geometry that admits highly accurate deterministic integration, ensures a meaningful comparison. In the following experiments, MC-fPINN, Improved MC-fPINN, and QE-MC-fPINN all employ the same PINN backbone and training protocol. The specific hyperparameter settings and implementation details are adopted from [13]. The angular resolution n1 of fTNN is set to be 32 in 2D case and 16 in 3D case, while fPINN is configured as described in [25]. Table 9: Relative L2 test errors using different methods to solve the fPEs on the unit balls. uexact

d=2

Method

d=3

α = 1.5

α = 1.9

α = 1.5

α = 1.9

(1 − ∥x∥2 )1+α/2

MC-fPINN Improved MC-fPINN fPINN QE-MC-fPINN fTNN

1.78e-2 4.88e-3 3.61e-3 2.00e-3 2.20e-6

5.31e-1 4.96e-1 1.12e-3 1.33e-3 1.70e-5

1.37e-2 7.66e-3 5.02e-3 1.74e-3 3.97e-4

8.11e-1 4.35e-1 5.59e-3 1.18e-3 3.56e-4

(1 − ∥x∥2 )α/2

MC-fPINN Improved MC-fPINN fPINN QE-MC-fPINN fTNN

7.01e-1 4.24e-1 1.36e-2 3.41e-2 2.17e-5

7.20e-1 5.12e-1 1.34e-2 6.05e-3 1.72e-5

4.43e-1 2.65e-1 4.95e-1 3.07e-2 3.74e-4

8.57e-1 1.88e-1 1.30e-0 3.00e-2 1.58e-4

Table 9 reports the relative L2 test errors etest of four different methods for solving twoand three-dimensional fPEs on the unit balls, with exact solutions

and

uexact = (1 − ∥x∥2 )1+α/2 ,

(4.4)

uexact = (1 − ∥x∥2 )α/2 ,

(4.5)

for fractional orders α = 1.5 and α = 1.9. 23

1e+00

Relative L2 error (etest)

1e-01 2D, MC fPINN 3D, MC fPINN 2D, Improved MC-fPINN 3D, Improved MC-fPINN 2D, QE-MC-fPINN 3D, QE-MC-fPINN

1e-01

1e-02

1e-03 5e-04

1e-03 Adam

1e-04

L-BFGS

1e-05

2D, BFE (n1 = 32) 3D, BFE (n1 = 16)

Adam

1e-02

8 4

2e-06

Training steps (102 )

Training steps (10 )

Figure 8: Relative L2 test error curves of fPINN-based methods (left column) and fTNN (right column) for solving fPE with exact solution (4.4), α = 1.5. For solution (4.4) with α = 1.5, all five methods behave well, although fTNN achieves substantially lower errors than the other methods. For solution (4.4) with α = 1.9, there are more pronounced differences among the errors of different methods. To be specific, the relative errors of MC-fPINN and Improved MC-fPINN are the greatest, fPINN and QE-MC-fPINN reduce them by about two orders of magnitude, while fTNN achieves a further reduction. For the more singular solution (4.5) in 2D case, the comparison results are similar to that for solution (4.4) with α = 1.9. While for that in 3D case, fPINN exhibits unsatisfactory performance. Notably, fTNN achieves substantially lower errors than the other methods across all test cases. The relative error curves in Figure 8 demonstrate that fTNN converges smoothly and attains markedly lower errors compared with MC-fPINN, Improved MC-fPINN, and QE-MC-fPINN. We plot the solutions in Figure 9 which indicates that fTNN accurately captures the singular structure near boundary without the erratic oscillations and variance inherent in Monte Carlo-based methods. Overall, these results highlight that the combination of deterministic integration, boundary-aware trial functions, and adaptive singularity treatment confers substantial improvements in both accuracy and stability for low- and moderate-dimensional fPEs.

0.5 1.0

1.0

0.5

0.0

0.5

x1 Exact u(x1, x2)

1.0

0.5 0.0 0.5 1.0

1.0

0.5

0.0

0.5

1.0

0.96 0.84 0.72 0.60 0.48 0.36 0.24 0.12 0.00

0.0 0.5 1.0

1.0

0.5

0.0

0.5

x1 Predicted u(x1, x2)

1.0

0.5 0.0 0.5 1.0

1.0

0.5

0.0

0.5

1.0

0.96 0.84 0.72 0.60 0.48 0.36 0.24 0.12 0.00

Absolute Error

1.0

1e 6

0.5

0.0

Predicted u(x1, x2)

1.0

0.5

0.96 0.84 0.72 0.60 0.48 0.36 0.24 0.12 0.00

0.0 0.5 1.0

1.0

0.5

0.0

x1 Absolute Error

0.5

1.0

1e 5

0.5

Exact u(x1, x2)

1.0

0.0 0.5 1.0

1.0

0.5

0.0

0.5

1.0

8 7 6 5 4 3 2 1 0

2.7 2.4 2.1 1.8 1.5 1.2 0.9 0.6 0.3 0.0

Figure 9: Plots of fTNN to solve fPE (1.1) on the unit ball. Top: exact solution (4.4) with α = 1.5, numerical solution, and the absolute error. Bottom: exact solution (4.5) with α = 1.9, numerical solution, and the absolute error.

4.2

Fractional advection-diffusion equation

We next turn to time-dependent problems, for which Algorithm 1 is used (This is what the fTNN refers to in this subsection). We consider the time-space fractional advection-diffusion equation (ADE) on the unit ball, with convection coefficient v = (0.1, 0.1) for d = 2 and v = (0.1, 0.1, 0.1) for d = 3, which is a typical case of the fPDE defined in (1.3). We employ the following quadrature scheme γ γ Γ(1 − γ) · C 0 Dt (t ϕt,j (t)) =

Z 1 0

(1 − τ )−γ τ γ−1 S1,j (t, τ )dτ ≈

N X

(−γ,γ−1)

S1,j t, τk

k=1

for the time-fractional Caputo derivative, where S1,j (t, τ ) = γϕt,j (tτ ) + tτ ϕ′t,j (tτ ), see [18] for further details. For the spatial operator, we use the same deterministic discretization as in Section 4.1.3; the new feature in this subsection is that this spatial solver is coupled with the spatiotemporally separable subspace representation and alternating optimization. In the two-dimensional case, we define the trial function Ψ in (2.5) as b Ψ(ρ, η, t; c, θ) = Ψ(ρ, η, t; c, θ) + 1 − ρ2 ,

with b Ψ(ρ, η, t; c, θ) =

p X

(4.6)

cj tγ ϕt,j (t)(1 − ρ2 )µj ϕ1,j (ρ)ϕ2,j (η).

j=1

We first examine the performance of the new method in short-time simulations. The time interval [0, 1] is divided into 4 subintervals, with 10 Gauss quadrature points selected within each subinterval to ensure sufficient accuracy for both two- and three-dimensional time-space fractional ADEs. The exact solution is taken as

uexact (x, t) = e−t 1 − ∥x∥22

1+α/2

The corresponding source term f (x, t) for both two and three dimensions is provided in [25]. Table 10: Relative L2 test errors of the time-space fractional ADEs using fPINN and fTNN. Problem

fPINN

2D, γ = 1.0 2D, γ = 0.5 3D, γ = 1.0 3D, γ = 0.5

1.066e-3 1.241e-3 2.359e-3 2.758e-3

fTNN n1 = 2

n1 = 4

n1 = 8

n1 = 16

2.051e-4 1.041e-3 2.294e-2 2.302e-2

1.803e-5 6.891e-5 3.402e-3 3.444e-3

6.233e-6 6.130e-5 4.110e-4 7.942e-4

7.217e-6 7.408e-5 3.158e-4 2.102e-4

Substituting the decomposition (4.6) into (1.3) transforms the original problem into a hob with the RHS function given by mogeneous initial-boundary value problem for Ψ,

fb =2α Γ(α/2 + 2)Γ(α/2 + 1)(e−t − 1) 1 − (α/2 + 1)ρ2 − (1 − ρ2 )1+α/2 t1−γ E1,2−γ (−t) − 0.2ρ(α/2 + 1)(e−t − 1)(1 − ρ2 )α/2 (cos η + sin η), where E·,· (z) is the Mittag-Leffler function. 25

1e-01 2D: n1 = 2 2D: n1 = 4 2D: n1 = 8 2D: n1 = 16

6e-03

Relative L2 error (etest)

2e-02

1e-04

3D: n1 = 2 3D: n1 = 4 3D: n1 = 8 3D: n1 = 16

1e-02

1e-03 Adam

Adam

3e-05

L-BFGS

2e-04

Training steps (102 )

Figure 10: Relative L2 test errors for solving time-space fractional ADE (γ = 0.5, α = 1.50). The relative L2 test errors of different methods are shown in Table 10. As seen in Table 9, fPINN achieves lower error than MC-fPINN and Improved MC-fPINN methods, here we list just fPINN among those fPINN-based methods to prevent redundancy. Relative to fPINN, fTNN reduces the errors to the 10−5 level in 2D and to the 10−4 level in 3D. The relative error curves in Figure 10 reveal that fTNN exhibits smooth convergence histories. These results highlight the benefit of coupling deterministic spatial quadrature with the STSNN subspace formulation. Table 11: Relative L2 test errors of fTNN for solving time-space fractional ADEs and average per-epoch times (seconds) of Adam and L-BFGS for long-time simulations. Case

(α, γ, γ1 , γ2 )

T = 50

T = 100

T = 150

Adam

L-BFGS

2D, case A1 2D, case A2 3D, case B1 3D, case B2

(1.90, 0.85, 0.85, 1.70) (1.30, 0.45, 0.75, 0.80) (1.90, 0.85, 0.85, 1.00) (1.30, 0.45, 0.75, 0.80)

3.363e-5 2.558e-4 9.140e-5 4.258e-4

9.074e-5 1.889e-4 9.486e-5 4.376e-4

3.405e-5 1.353e-4 8.933e-5 4.411e-4

0.119s 0.121s 0.413s 0.419s

0.744s 0.758s 8.445s 8.544s

We then proceed to discuss the case of long-time simulations. Since long-time simulation of fPDEs faces far more severe theoretical and computational obstacles than short-time counterpart, the performance of fTNN deserves focused and in-depth discussion. To this end, the exact solution is set as uexact (x, t) = (tγ2 + tγ1 )(1 − ∥x∥22 )1+α/2 . (4.7) During the computation of the loss function, the time interval [0, T ] is divided into 25 subintervals, with 16 Gauss quadrature points selected within each subinterval. As shown in Table 11, the relative L2 test errors at final times T = 50, 100, 150 remain in the 10−4 -10−5 range, indicating that the chosen temporal quadrature is sufficiently dense to resolve the memory term even over long horizons. Table 12: Relative L2 test errors using different methods (long-time simulation, T = 100). Case

MC-fPINN

Improved MC-fPINN

fPINN

QE-MC-fPINN

fTNN

2D, case A2 3D, case B1

1.001e-01 9.851e-01

2.321e-01 3.717e-01

4.066e-02 3.371e-03

6.208e-03 1.222e-03

1.889e-4 9.486e-5

We also compare the relative L2 test errors at T = 100 (for case A2 and case B1 in Table 11) of different methods in Table 12. The fPINN results are obtained using the DeepXDE implementation [20] with the L1 scheme and Grünwald–Letnikov discretization configured as in [25]. The MC-fPINN, Improved MC-fPINN, and QE-MC-fPINN use the same hyperparameters as 26

reported in [22]. All other solver configurations are kept the same as in their original references. The comparison result is similar to that in Section 4.1.3: the relative errors of MC-fPINN and Improved MC-fPINN are the greatest, fPINN and QE-MC-fPINN reduce them by one or two orders of magnitude, while fTNN achieves a further reduction in error. These results further verify the advantages of STSNN subspace formulation for temporally nonlocal dynamics, which is the third contribution of this paper. We may further expound upon the function that STSNN serves in fTNN for solving longtime fractional PDEs. The STSNN decouples the spatial and temporal dimensions, reducing the evaluation of high-dimensional time-space fractional integrals to lower-dimensional temporal and spatial fractional integrals. For long-time simulations, this allows the use of a dense set of Gauss quadrature points in the time direction without causing prohibitive memory growth. This is precisely the aspect that distinguishes the present work most clearly from QE-MCfPINN: beyond improving the spatial operator evaluation, it furnishes a practical subspace mechanism for accurate long-time simulation of time-space fractional ADEs on bounded domains. The per-epoch costs reported in Table 11 remain moderate in 2D and still manageable in 3D, which shows that the deterministic/subspace reformulation improves accuracy without destroying computational feasibility on structured domains. 5e-02

Relative L2 error (etest)

1e+04

1e+03

Loss

1e+02

1e+01

Adam 2D, case A1 2D, case A2 3D, case B2 3D, case B1

1e+00

1e-01

Training steps (102 )

Adam

1e-03

1e-04 6e-05

L-BFGS

1e-02

2D, case A1 2D, case A2 3D, case B2 3D, case B1

L-BFGS

Training steps (102 )

Figure 11: The loss curve and the relative L2 test error of fTNN for solving time-space fractional ADE (long-time simulation, T = 150). The loss curves in Figure 11 further indicate that the alternating STSNN subspace optimization is stable in the long-time regime: the relative L2 test errors reach the 10−4 level within the first 100 epochs and then settle after roughly 200 epochs. At the final horizon T = 150, this corresponds to about 430 s total runtime in 2D and about 2100 s in 3D, which is acceptable given the difficulty of the nonlocal time-space operator.

Conclusion

We have developed a fully deterministic, subspace-based neural framework for fractional PDEs on bounded domains with structured geometries up to three dimensions. It rests on three mutually reinforcing components: (i) geometry-adapted quadrature that resolves singular radial integrals with Gauss-Jacobi rules and replaces Monte Carlo sampling by deterministic angular quadrature; (ii) boundary-singularity-aware trial spaces whose leading exponents are adaptively selected (BFE/BRFE) to match the asymptotic structure induced by the operator and the source term; and (iii) a spatiotemporally separable neural network with alternating subspace optimization that factorizes the time-space residual into lower-dimensional integrals. Numerical experiments from one to three dimensions show that the method yields accurate and stable approximations, with the largest gains in regimes of strong boundary singularities and long-time simulations-precisely where existing neural network solvers struggle most. A 27

systematic study of the angular resolution confirms that, after the radial singularity is removed, the error decays consistently with the O(n−2 1 ) rate predicted by spherical integration theory, making deterministic angular quadrature the primary accuracy lever without Monte Carlo noise. More broadly, this framework demonstrates that on structured domains, a careful synthesis of problem-adapted singular quadrature, explicit boundary enrichment, and separable tensornetwork architecture can simultaneously eliminate stochastic noise and mitigate ill-conditioningtwo longstanding obstacles in neural network-based fractional PDE solvers. Future work will extend the approach to more general geometries, nonlinear fractional models, and hybrid deterministic-stochastic strategies that preserve the present accuracy while improving scalability to higher dimensions.

References [1] Gabriel Acosta, Francisco M Bersetche, and Juan Pablo Borthagaray. A short FE implementation for a 2d homogeneous Dirichlet problem of a fractional Laplacian. Computers & Mathematics with Applications, 74(4):784–816, 2017. [2] Mark Ainsworth and Christian Glusa. Towards an efficient finite element method for the integral fractional Laplacian on polygonal domains. In Contemporary Computational Mathematics-A Celebration of the 80th Birthday of Ian Sloan, pages 17–57. Springer, 2018. [3] Andrea Bonito, Juan Pablo Borthagaray, Ricardo H Nochetto, Enrique Otárola, and Abner J Salgado. Numerical methods for fractional diffusion. Computing and Visualization in Science, 19(5):19–46, 2018. [4] Johann Brauchart, E Saff, I Sloan, and R Womersley. QMC designs: optimal order quasi Monte Carlo integration schemes on the sphere. Mathematics of Computation, 83(290):2821–2851, 2014. [5] Dariusz W Brzezinski. Computation of Gauss-Jacobi quadrature nodes and weights with arbitrary precision. In 2018 Federated Conference on Computer Science and Information Systems (FedCSIS), pages 297–306. IEEE, 2018. [6] Luis Caffarelli and Luis Silvestre. An extension problem related to the fractional Laplacian. Communications in Partial Differential Equations, 32(8):1245–1260, 2007. [7] Marta D’Elia, Qiang Du, Christian Glusa, Max Gunzburger, Xiaochuan Tian, and Zhi Zhou. Numerical methods for nonlocal and fractional models. Acta Numerica, 29:1–124, 2020. [8] Jing Gao, Meng Zhao, Ning Du, Xu Guo, Hong Wang, and Jiwei Zhang. A finite element method for space–time directional fractional diffusion partial differential equations in the plane and its error analysis. Journal of Computational and Applied Mathematics, 362:354– 365, 2019. [9] Peter J Grabner and Tetiana A Stepanyuk. Upper and lower estimates for numerical integration errors on spheres of arbitrary dimension. Journal of Complexity, 53:113–132, 2019. [10] Gerd Grubb. Fractional Laplacians on domains, a development of Hörmander’s theory of µ-transmission pseudodifferential operators. Advances in Mathematics, 268:478–528, 2015.

[11] Ling Guo, Hao Wu, Xiaochen Yu, and Tao 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:115523, 2022. [12] Yixiao Guo and Pingbing Ming. A deep learning method for computing eigenvalues of the fractional Schrödinger operator. Journal of Systems Science and Complexity, 37(2):391– 412, 2024. [13] Zheyuan Hu, Kenji Kawaguchi, Zhongqiang Zhang, and George Em Karniadakis. Tackling the curse of dimensionality in fractional and tempered fractional PDEs with physicsinformed neural networks. Computer Methods in Applied Mechanics and Engineering, 432:117448, 2024. [14] Weizhang Huang and Jinye Shen. A grid-overlay finite difference method for the fractional Laplacian on arbitrary bounded domains. SIAM Journal on Scientific Computing, 46(2):A744–A769, 2024. [15] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks, 9(5):987–1000, 1998. [16] Yongxin Li, Zhongshuo Lin, Yifan Wang, and Hehu Xie. Tensor neural network interpolation and its applications. arXiv preprint arXiv:2404.07805, 2024. [17] Yangfei Liao, Zhongshuo Lin, Jianghao Liu, Qingyuan Sun, Yifan Wang, Teng Wu, and Hehu Xie. Solving Schrödinger equation using tensor neural network. arXiv preprint arXiv:2209.12572, 2022. [18] Zhongshuo Lin, Qingkui Ma, Hehu Xie, and Xiaobo Yin. Solving time-fractional partial integro-differential equations using tensor neural network. SIAM Journal on Scientific Computing, 48(1):C164–C189, 2026. [19] Anna Lischke, Guofei Pang, Mamikon Gulian, Fangying Song, Christian Glusa, Xiaoning Zheng, Zhiping Mao, Wei Cai, Mark M Meerschaert, Mark Ainsworth, et al. What is the fractional Laplacian? A comparative review with new results. Journal of Computational Physics, 404:109009, 2020. [20] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis. DeepXDE: a deep learning library for solving differential equations. SIAM Review, 63(1):208–228, 2021. [21] Lei Ma, Fanhai Zeng, Ling Guo, George Em Karniadakis, et al. Bi-Orthogonal fPINN: a physics-informed neural network method for solving time-dependent stochastic fractional PDEs. Communications in Computational Physics, 34(4):1133–1176, 2023. [22] Qingkui Ma, Hehu Xie, and Xiaobo Yin. Quadrature-Enhanced Monte Carlo fPINN method for high-dimensional fractional PDEs. arXiv preprint arXiv:2604.19601, 2026. [23] Mark M Meerschaert, Hans-Peter Scheffler, and Charles Tadjeran. Finite difference methods for two-dimensional fractional dispersion equation. Journal of Computational Physics, 211(1):249–261, 2006. [24] Guofei Pang, Wen Chen, and Zhuojia Fu. Space-fractional advection–dispersion equations by the Kansa method. Journal of Computational Physics, 293:280–296, 2015. 29

[25] Guofei Pang, Lu Lu, and George Em Karniadakis. fPINNs: fractional physics-informed neural networks. SIAM Journal on Scientific Computing, 41(4):A2603–A2626, 2019. [26] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019. [27] Xavier Ros-Oton and Joaquim Serra. Boundary regularity for fully nonlinear integrodifferential equations. Duke Mathematical Journal, 165(11):2079–2154, 2016. [28] Jie Shen, Tao Tang, and Li-Lian Wang. Spectral methods: algorithms, analysis and applications, volume 41. Springer Science & Business Media, 2011. [29] Changtao Sheng, Bihao Su, and Chenglong Xu. Efficient Monte Carlo method for integral fractional Laplacian in multiple dimensions. SIAM Journal on Numerical Analysis, 61(5):2035–2061, 2023. [30] Changtao Sheng, Li-Lian Wang, Hongbin Chen, and Huiyuan Li. Fast implementation of FEM for integral fractional Laplacian on rectangular meshes. Communications in Computational Physics, 36(3):673–710, 2024. [31] Charles Tadjeran and Mark M Meerschaert. A second-order accurate numerical method for the two-dimensional fractional diffusion equation. Journal of Computational Physics, 220(2):813–823, 2007. [32] Shupeng Wang and George Em Karniadakis. GMC-PINNs: a new general Monte Carlo PINNs method for solving fractional partial differential equations on irregular domains. Computer Methods in Applied Mechanics and Engineering, 429:117189, 2024. [33] Yifan Wang, Pengzhan Jin, and Hehu Xie. Tensor neural network and its numerical integration. Journal of Computational Mathematics, 42(6):1714–1742, 2024. [34] Yifan Wang, Zhongshuo Lin, Yangfei Liao, Haochen Liu, and Hehu Xie. Solving highdimensional partial differential equations using tensor neural network and a posteriori error estimators. Journal of Scientific Computing, 101(3):67, 2024. [35] Yifan Wang and Hehu Xie. Computing multi-eigenpairs of high-dimensional eigenvalue problems using tensor neural networks. Journal of Computational Physics, 506:112928, 2024. [36] Tianxin Zhang, Dazhi Zhang, Shengzhu Shi, and Zhichang Guo. Spectral-fPINNs: spectral method based fractional physics-informed neural networks for solving fractional partial differential equations. Nonlinear Dynamics, 113(11):12565, 2025.

Record · ID 310825 · SHA-256 0ee6a89f4f1c1738
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.