1 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, [1]–[4]. 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 [5]. The main difficulties stem from the simultaneous presence of hypersingular nonlocal kernels, exterior Dirichlet constraints, and reduced boundary regularity of solutions [6]–[8]. These features make the design of accurate and robust numerical methods for fractional Laplacian problems on bounded domains especially demanding.

In this work, we first consider the fractional Poisson equation (fPE) with homogeneous Dirichlet exterior conditions: \[\label{eq:fra95lapla95problem} \left\{ \begin{align} (-\Delta)^{\alpha / 2} u(\boldsymbol{x}) &= f(\boldsymbol{x}), &&\boldsymbol{x} \in \Omega, \\ u(\boldsymbol{x}) &= 0, &&\boldsymbol{x} \in \Omega^c. \end{align} \right.\tag{1}\] Here \(\Omega\subset\mathbb{R}^d\) is either a hyperrectangle or the unit ball, and \(0<\alpha<2\). The fractional Laplacian operator \((-\Delta)^{\alpha / 2}\) admits several non-equivalent definitions in different settings [4]. In this paper, we adopt the Riesz definition with zero exterior extension: \[\label{eq:def95lap} (-\Delta)^{\alpha/ 2} u(\boldsymbol{x}) = C_{d,\alpha} \int_{\mathbb{R}^d} \frac{u(\boldsymbol{x}) - u(\boldsymbol{y})}{\| \boldsymbol{y} - \boldsymbol{x} \|^{d+\alpha}} d\boldsymbol{y} \quadwith\quad C_{d,\alpha}=\frac{2^{\alpha -1}\alpha\Gamma \left(\alpha/ 2+d/2 \right)}{\pi^{d/2} \Gamma (1-\alpha/ 2)}.\tag{2}\] We also consider the following time-dependent fractional partial differential equation (fPDE) posed on \(\Omega\) with zero exterior condition [9]: \[\label{eq:fPDE} \left\{ \begin{align} L_{t,\boldsymbol{x}}\{u(\boldsymbol{x},t)\}=\frac{\partial^\gamma u(\boldsymbol{x},t)}{\partial t^\gamma}+L_{\boldsymbol{x}}\{u(\boldsymbol{x},t)\} &=f(\boldsymbol{x},t), \quad \boldsymbol{x}\in\Omega,\;t\in(0,T],\\ u(\boldsymbol{x},t)&=0, \qquad\quad\:\: \boldsymbol{x}\in\Omega^c,\;t\in(0,T],\\ u(\boldsymbol{x},0)&=u_0(\boldsymbol{x}), \quad\:\: \boldsymbol{x}\in\Omega. \end{align} \right.\tag{3}\] Here \(L_{\boldsymbol{x}}\) denotes a linear spatial operator involving the fractional Laplacian and acting on \(u\). A typical example is \[\label{eq:space95L} L_{\boldsymbol{x}} = (-\Delta)^{\alpha/2}+\boldsymbol{v}\cdot\nabla,\tag{4}\] where \(\boldsymbol{v}\cdot\nabla\) is the advection term associated with a prescribed velocity field \(\boldsymbol{v}\). The time-fractional derivative in 3 is understood in the Caputo sense: \[\label{eq:Caputo} \frac{\partial^\gamma u(\boldsymbol{x},t)}{\partial t^\gamma} =\frac{1}{\Gamma(1-\gamma)} \int_0^t (t-\tau)^{-\gamma} \frac{\partial u(\boldsymbol{x},\tau)}{\partial \tau}\, d\tau, \qquad 0<\gamma\le 1.\tag{5}\] When \(\gamma=1\), the Caputo derivative reduces to the classical first-order time derivative. By choosing different values of \(\gamma\) and \(\alpha\), and by including or excluding the advection term, 3 recovers several important classes of fractional models. In particular, when \(0<\gamma<1\), \(\alpha=2\), and \(\boldsymbol{v}=0\), it reduces to a time-fractional diffusion equation [10], which describes subdiffusion with memory effects. When \(\gamma=1\) and \(\alpha\in(0,2)\) with \(\boldsymbol{v}=0\), it becomes a time-dependent space-fractional diffusion equation that captures nonlocal transport phenomena such as Lévy-flight-type diffusion. When \(0<\gamma<1\), \(\alpha\in(0,2)\), and \(\boldsymbol{v}\neq 0\), 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 [4].

Classical discretization methods, such as the finite difference method [11]–[13] and the finite element method [5], [14]–[16], 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 high-accuracy 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. [17], physics-informed neural networks (PINNs) [18] and their fractional extensions have been developed for a variety of forward and inverse problems. Representative examples include fPINNs [9], bi-orthogonal fPINNs [19], spectral-fPINNs [20], and general Monte Carlo PINN formulations on irregular domains [21]. 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. [8] 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. [22] proposed the so-called MC-fPINN method, in which the fractional Laplacian is approximated by Monte Carlo sampling after splitting the integral into integrals over a neighborhood \(B_{r_{0}}(\boldsymbol{x})\) of \(\boldsymbol{x}\) and its complement: \[(-\Delta)^{\alpha/ 2} u(\boldsymbol{x}) = C_{d,\alpha} \left( \int_{\boldsymbol{y} \in B_{r_{0}}(\boldsymbol{x})} \frac{u(\boldsymbol{x}) - u(\boldsymbol{y})}{\|\boldsymbol{x} - \boldsymbol{y}\|^{d+\alpha}} d\boldsymbol{y} + \int_{\boldsymbol{y} \notin B_{r_{0}}(\boldsymbol{x})} \frac{u(\boldsymbol{x}) - u(\boldsymbol{y})}{\|\boldsymbol{x} - \boldsymbol{y}\|^{d+\alpha}} d\boldsymbol{y} \right).\] 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. [23] 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 [24] 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 [25]–[29]. In particular, Lin et al. [10] proposed a tensor neural network subspace method to solve time-fractional partial integro-differential equations. The method combines a power-law temporal factor \(t^{\mu}\), 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 [22], Improved MC-fPINN [23], and especially the QE-MC-fPINN method [24], 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 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(\boldsymbol{x})^{\mu_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 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 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 4 to demonstrate the accuracy and efficiency of the proposed methods.

2 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)^{\beta_{1}}(1+x)^{\beta_{2}}\) using \(N\) nodes and weights (see [30], [31]): \[\label{Def95gauss95jacobi9511} \int_{-1}^1 (1-x)^{\beta_{1}}(1+x)^{\beta_{2}}f(x)dx = \sum_{i=1}^N \widehat{w}_{i}^{(\beta_{1},\beta_{2})} f(\widehat{x}_{i}^{(\beta_{1},\beta_{2})}) + R_N(f),\tag{6}\] with \(\beta_{1}> -1\), \(\beta_{2}> -1\), where \(R_N(f)\) is the integration error. The quadrature formula on a general interval \([a, b]\) is constructed as \[\label{Gauss95Jacobi952} \int_a^b (b-x)^{\beta_{1}}(x-a)^{\beta_{2}} f(x)\;dx \approx \sum_{k=1}^{N} w_{k}^{(\beta_{1},\beta_{2})} f(x_{k}^{(\beta_{1},\beta_{2})}),\tag{7}\] where the nodes and weights are obtained through affine transformation: \[x_{k}^{(\beta_{1},\beta_{2})} = \frac{b-a}{2} \widehat{x}_{k}^{(\beta_{1},\beta_{2})} + \frac{a+b}{2}, \quad w_{k}^{(\beta_{1},\beta_{2})} = \left( \frac{b-a}{2} \right)^{\beta_{1}+\beta_{2}+1} \widehat{w}_{k}^{(\beta_{1},\beta_{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 (2.3), it is also employed to evaluate the loss function when the residual exhibits boundary singularities.

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 [29].

The neural network is built with \(d+1\) subnetworks, and each subnetwork is a continuous mapping from a bounded closed set \(\Omega_{i}\subset\mathbb{R}\) (\(i=1, \cdots, d\)) to \(\mathbb{R}^p\), which can be expressed as \[\begin{align} \Phi_{i}(x_{i};\theta_{i})&=\big(\phi_{i,1}(x_{i};\theta_{i}), \phi_{i,2}(x_{i};\theta_{i}),\cdots,\phi_{i,p}(x_{i};\theta_{i})\big)^{\top},\\ \Phi_t(t;\theta_t)&=(\phi_{t,1}(t;\theta_t),\phi_{t,2}(t;\theta_t),\cdots ,\phi_{t,p}(t;\theta_t))^{\top}, \end{align}\] where each \(x_{i}\) denotes the one-dimensional input, \(\theta_{i}\) denotes the parameters of the \(i\)-th subnetwork, typically the weights and biases. As illustrated in 1, the neural network structure adopted in this paper is composed of \(d\) fully connected neural networks (FNNs) for the spatial basis functions \(\Phi_{i}(x_{i};\theta_{i}), i=1,2,\ldots,d\), and one FNN for the temporal basis function \(\Phi_t(t;\theta_t)\).

To improve numerical stability, we normalize each \(\phi_{i,j}(x_{i})\), \(\phi_{t,j}(t)\) and use the following normalized neural network structure: \[\label{def95NN95normed} \Psi (\boldsymbol{x},t;c,\theta )=\sum_{j=1}^p{c_j}\widehat{\phi }_{t,j}(t;\theta_t) \widehat{\varphi}_j(\boldsymbol{x};\theta_{\boldsymbol{x}}),\tag{8}\] where \[\widehat{\varphi}_j(\boldsymbol{x};\theta_{\boldsymbol{x}}) =\prod\limits_{i=1}^d{\widehat{\phi }_{i,j}(x_{i};\theta_{i})},\: \widehat{\phi }_{i,j}(x_{i};\theta _{i})=\frac{\phi _{i,j}(x_{i};\theta _{i})}{\left\| \phi _{i,j}(x_{i};\theta _{i}) \right\| _{L^{2}(\Omega _{i})}} , \: \widehat{\phi }_{t,j}(t;\theta _t)=\frac{\phi _{t,j}(t;\theta _t)}{\left\| \phi _{t,j}(t;\theta _t) \right\| _{L^{2}((0,T])}}.\] Here, the \(L^2\) norms are precomputed using high-order Gauss quadrature on each \(\Omega_i\) and on the time interval \((0,T]\).

We denote the neural network parameters as \(\theta=\{\theta_t,\theta_{\boldsymbol{x}}\}=\{\theta_t,\theta_1,\theta_2,...,\theta_d\}.\) To simplify the notation, we drop the hats hereafter and write \(\phi_{i,j}\), \(\phi_{t,j}\), and \(\varphi_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 \(L^2(\Omega\times(0,T])\) [29]. 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}(\Omega \times (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 \(\varphi_j(\boldsymbol{x})\) for \(\varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}})\) and \(\phi_{t,j}(t)\) for \(\phi_{t,j}(t;\theta_t)\).

2.3.1 Discrete scheme of the fractional Laplacian↩︎

We adopt a geometry-adaptive strategy proposed in [24] that splits the fractional Laplacian based on a spatially varying radius and directional distance-to-boundary information. For any \(\boldsymbol{x}\in\Omega\), define \(r_{0}(\boldsymbol{x})\) as the minimum distance from \(\boldsymbol{x}\) to the boundary \(\partial\Omega\). For each spatial basis function \(\varphi_j\), the fractional Laplacian is first divided into near-field and far-field parts: \[\frac{(-\Delta)^{\alpha/2} \varphi_j(\boldsymbol{x})}{C_{d,\alpha}} = \left(\int_{\|\boldsymbol{y}-\boldsymbol{x}\|_2 < r_{0}(\boldsymbol{x})} +\int_{\|\boldsymbol{y}-\boldsymbol{x}\|_2 \ge r_{0}(\boldsymbol{x})}\right) \frac{\varphi_j(\boldsymbol{x}) - \varphi_j(\boldsymbol{y})}{\|\boldsymbol{x} - \boldsymbol{y}\|_2^{d+\alpha}}\,d\boldsymbol{y} := I_{1,j}(\boldsymbol{x})+I_{2,j}(\boldsymbol{x}).\] It is derived in [24] that the near-field integral \[\begin{align} I_{1,j}(\boldsymbol{x}) &= r_{0}(\boldsymbol{x})^{-\alpha} \int_{\mathcal{S}_+^{d-1}} \int_0^1 \tau^{1-\alpha} \frac{F_{j}\left(\boldsymbol{x}, r_{0}(\boldsymbol{x})\tau, \boldsymbol{\xi}\right)}{\tau^2} \, d\tau \, J_{d}(\boldsymbol{\xi}) d\boldsymbol{\xi}, \end{align}\] with \[F_{j}(\boldsymbol{x}, r, \boldsymbol{\xi}) := 2\varphi_j(\boldsymbol{x}) - \varphi_j(\boldsymbol{x}+r\boldsymbol{\xi}) - \varphi_j(\boldsymbol{x}-r\boldsymbol{\xi}).\] Here \(\mathcal{S}_+^{d-1}\) denotes the upper hemisphere in \(\mathbb{R}^d\). To handle the singular behavior of \(I_{1,j}(\boldsymbol{x})\) around \(\boldsymbol{y}=\boldsymbol{x}\), the Gauss-Jacobi quadrature is used: \[\begin{align} I_{1,j}(\boldsymbol{x})\approx r_{0}(\boldsymbol{x})^{-\alpha} \sum_{\ell=1}^{N_0} \sum_{k=1}^{N} w_\ell J_{d}(\boldsymbol{\xi}_\ell) w_k^{(0,1-\alpha)} \frac{F_{j} \left(\boldsymbol{x}, r_{0}(\boldsymbol{x})\tau_k^{(0,1-\alpha)}, \boldsymbol{\xi}_\ell\right)}{\left(\tau_k^{(0,1-\alpha)}\right)^{2}}, \end{align}\] where \(\{\tau_k^{(0,1-\alpha)}\}_{k=1}^{N}\) and \(\{w_k^{(0,1-\alpha)}\}_{k=1}^{N}\) are Gauss-Jacobi points and weights on \([0,1]\). For \(d=2\), the angular integral over \(\mathcal{S}_+^{1}\) is discretized by an \(N_0\)-point Gauss-Legendre rule on \((0,\pi]\), producing the directional nodes \(\{\boldsymbol{\xi}_\ell\}_{\ell=1}^{N_0}\) and weights \(\{w_\ell\}_{\ell=1}^{N_0}\). For \(d=3\), a tensor-product Gauss-Legendre rule in the two spherical angles is used to generate the corresponding nodes and weights.

The distance from \(\boldsymbol{x}\) to the boundary along a unit direction \(\boldsymbol{\xi}\in\mathcal{S}^{d-1}\) was defined in [24] as \[d_{\boldsymbol{x}}(\boldsymbol{\xi}):=\min\{\,t>0 \mid \boldsymbol{x}+t\boldsymbol{\xi}\in\partial\Omega \,\}.\] Then the far-field integral \(I_{2,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 \[\begin{align} I_{2,j}(\boldsymbol{x}) &=\int_{\mathcal{S}^{d-1}} \left( \int_{r_{0}(\boldsymbol{x})}^{d_{\boldsymbol{x}}(\boldsymbol{\xi})} + \int_{d_{\boldsymbol{x}}(\boldsymbol{\xi})}^{\infty} \right) \frac{\varphi_j(\boldsymbol{x}) - \varphi_j(\boldsymbol{x}+r\boldsymbol{\xi})}{r^{1+\alpha}}\,dr\; J_{d}(\boldsymbol{\xi})d\boldsymbol{\xi}\\ &=\int_{\mathcal{S}^{d-1}} \left[\left(d_{\boldsymbol{x}}(\boldsymbol{\xi})-r_{0}(\boldsymbol{x})\right)\int_{0}^{1} \frac{\varphi_j(\boldsymbol{x}) - \varphi_j\left(\boldsymbol{x}+r(\boldsymbol{x},\boldsymbol{\xi},t)\boldsymbol{\xi}\right)}{r(\boldsymbol{x},\boldsymbol{\xi},t)^{1+\alpha}}dt + \frac{\varphi_j(\boldsymbol{x})}{\alpha\,[d_{\boldsymbol{x}}(\boldsymbol{\xi})]^{\alpha}} \right] J_{d}(\boldsymbol{\xi})d\boldsymbol{\xi}\\ &\approx\int_{\mathcal{S}^{d-1}} \left[\left(d_{\boldsymbol{x}}(\boldsymbol{\xi})-r_{0}(\boldsymbol{x})\right)\sum_{m=1}^{N'} w_{m}\frac{\varphi_j(\boldsymbol{x}) - \varphi_j\left(\boldsymbol{x}+r(\boldsymbol{x},\boldsymbol{\xi},t_{m})\boldsymbol{\xi}\right)}{r(\boldsymbol{x},\boldsymbol{\xi},t_{m})^{1+\alpha}} + \frac{\varphi_j(\boldsymbol{x})}{\alpha\,[d_{\boldsymbol{x}}(\boldsymbol{\xi})]^{\alpha}} \right] J_{d}(\boldsymbol{\xi})d\boldsymbol{\xi}\\ &:=\int_{\mathcal{S}^{d-1}} \left[Q_{1,j}(\boldsymbol{x},\boldsymbol{\xi})+Q_{2,j}(\boldsymbol{x},\boldsymbol{\xi})\right] J_{d}(\boldsymbol{\xi})d\boldsymbol{\xi}, \end{align}\] where \(\{t_{m}\}_{m=1}^{N'}\) and \(\{w_{m}\}_{m=1}^{N'}\) are the Gauss-Legendre nodes and weights on \([0,1]\). Here \(\mathcal{S}^{d-1}\) denotes the unit sphere in \(\mathbb{R}^d\). Then \(I_{2,j}\) is approximated by a discrete average over directions as follows: \[\begin{align} I_{2,j}(\boldsymbol{x})\approx \sum_{i=1}^{2}{\sum_{k=1}^{N_{i}}{\left(\pi/n_{i} \right) ^{d-1}J_{d}(\boldsymbol{\xi }_{ik})Q_{i,j}(\boldsymbol{x},\boldsymbol{\xi }_{ik})}}, \end{align}\] with \(N_{i} = 2n_{i}^{d-1}\) the number of distributed directions. For \(d\ge 2\), these directions are parameterized over the spherical coordinate domain \([0,2\pi] \times [0,\pi]^{d-2}\), while the associated weight for each direction is given by \(\left(\pi/n_{i} \right)^{d-1}J_{d}(\boldsymbol{\xi }_{ik})\). Direction vectors \(\boldsymbol{\xi}_{ik}\) in different dimensional cases were defined in [24].

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 2. The discretization scheme described here can, in principle, be generalized to higher dimensions. However, this work focuses on the high-accuracy 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 [22], Improved MC-fPINN [23], and QE-MC-fPINN method [24], 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 (\(\mu_1 = \alpha/2\)) to capture the operator-induced singularity. The present method further proposes the boundary and right-hand-side feature enhanced (BRFE) strategy: when the right-hand side (RHS) function is also singular, BRFE adjusts the exponent to \(\mu_1 = \alpha + 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 \ge 4\), while the present method trades higher per-dimension cost via deterministic directional quadrature for substantially superior accuracy, achieving relative \(L^2\) 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 low- to 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.

2.3.2.1 The near-field integral:

As derived in [24], \[\begin{align} \int_0^{r_0(\boldsymbol{x})}\frac{F_j(\boldsymbol{x},r,\boldsymbol{\xi})}{r^{1+\alpha}}\,dr =-\,\frac{r_0(\boldsymbol{x})^{2-\alpha}}{2-\alpha}\,\boldsymbol{\xi}^\top H(\varphi_j)(\boldsymbol{x})\boldsymbol{\xi} +O(r_0(\boldsymbol{x})^{4-\alpha}), \end{align}\] justifying the use of Gauss-Jacobi quadrature with parameters \((0,1-\alpha)\), which is tailored for integrands with this type of weak singularity. Since \(B_{r_0(\boldsymbol{x})}(\boldsymbol{x})\subset\Omega\), the local Taylor expansion is valid for \(0<r\le r_0(\boldsymbol{x})\) under standard regularity assumptions. Combined with a deterministic angular quadrature, the discretization yields high-accuracy and computationally efficient evaluation of \(I_{1,j}(\boldsymbol{x})\).

2.3.2.2 The far-field integral:

The angular integration over the unit sphere \(\mathcal{S}^{d-1}\) is discretized using a uniform directional grid determined by the angular resolution parameter \(n_{1}\), while the radial integration is treated separately for the terms \(Q_{1,j}\) and \(Q_{2,j}\). The term \(Q_{2,j}\) admits an analytical expression, which allows a relatively dense directional sampling to be used without additional computational cost. In contrast, \(Q_{1,j}\) is approximated by Gauss quadrature in the radial direction, and hence the principal error of \(I_{2,j}\) arises from the discretization of the spherical integral of \(Q_{1,j}\). For smooth integrands the uniform angular discretization yields an error decaying as \(O(n_1^{-2})\) (see, e.g., [32], [33]), so the accuracy of \(I_{2,j}\) can be improved by increasing \(n_1\).

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 [34], [35], solutions of fPDEs often exhibit an \(\alpha/2\)-order boundary singularity characterized by the asymptotic relation \[u(\boldsymbol{x}) \approx \operatorname{dist}(\boldsymbol{x},\partial\Omega)^{\alpha/2}u_{\mathrm{reg}}(\boldsymbol{x}), \quad \boldsymbol{x}\in\overline{\Omega},\] where \(u_{\mathrm{reg}}\) is comparatively smoother. Accurate treatment of this boundary behavior is therefore indispensable for the fractional problems considered here [24], [35].

To embed the homogeneous Dirichlet boundary condition directly into the ansatz, we introduce a prescribed boundary feature function \(b(\boldsymbol{x})\) satisfying \[\label{bd95fun} b(\boldsymbol{x})>0,\:\boldsymbol{x}\in\Omega,\quad and \quad b(\boldsymbol{x})=0, \:\boldsymbol{x}\in\Omega^{c},\tag{9}\] 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 \(\Psi(\boldsymbol{x}; c, \theta_{\boldsymbol{x}})\) which satisfies the homogeneous boundary conditions is then constructed as \[\label{eq:trai95fun} \Psi (\boldsymbol{x};c,\theta_{\boldsymbol{x}} ) =\sum_{j=1}^p{c_j\,b(\boldsymbol{x})^{\mu_{j}}\prod_{i=1}^{d}\phi _{i,j}(x_{i};\theta _{i})} :=\sum_{j=1}^p{c_j\,\varphi _j(\boldsymbol{x};\theta_{\boldsymbol{x}} ),}\tag{10}\] where \(\{\mu_{j}\}_{j=1}^{p}\) 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. To ensure the final output \(\Psi\) lies in the desired function space, we require that \(b(\boldsymbol{x})^{\mu_{j}}\in H_0^{\alpha/2}(\Omega)\). The values of the sequence \(\{\mu_{j}\}_{j=1}^p\) are determined by specifying only its first and last elements, \(\mu_{1}\) and \(\mu_{p}\). The remaining parameters are then defined via linear interpolation: \[\label{eq:def95muj} \mu_{j} = \mu_{1} + (j-1)(\mu_{p} - \mu_{1})/(p-1), \qquad j = 1, 2, \dots, p.\tag{11}\] The loss function is then defined as the \(L^{2}\) residual of the governing equation: \[\label{eq:loss95function} \mathcal{L}(c,\theta_{\boldsymbol{x}}) = \left\| L_{\boldsymbol{x}}\Psi(\boldsymbol{x};c,\theta_{\boldsymbol{x}}) - f(\boldsymbol{x}) \right\|_{L^{2}(\Omega)}.\tag{12}\] Minimizing this loss with respect to the linear coefficients \(c\) and the network parameters \(\theta_{\boldsymbol{x}}\) yields the approximate solution. For problems involving fractional Laplacian, however, the integrand in 12 may inherit a strongly singular boundary profile unless the dominant asymptotic factor is extracted explicitly. To this end, we define for \(\boldsymbol{x} \in \Omega\) \[\label{eq:s95condition} f_{\mathrm{reg}}(\boldsymbol{x}) := f(\boldsymbol{x})/b(\boldsymbol{x})^{s}, \quadwith \quad\limsup_{\substack{\boldsymbol{x} \to \partial \Omega \\ \boldsymbol{x} \in \Omega}} \frac{f(\boldsymbol{x})}{b(\boldsymbol{x})^{s}} = C \neq 0.\tag{13}\] Here the \(\limsup\) is taken over all possible interior paths approaching the boundary. With this choice, \(f_{\mathrm{reg}}\) 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 \(\mu_{1}\) in 11 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 \(\mu_{1}\) with the prevailing singularity mechanism of the problem. When \(s \geq 0\), we employ the boundary feature enhanced (BFE) strategy which has been used in [24], [35]. When \(-0.5<s < 0\), the boundary and RHS feature enhanced (BRFE) strategy is used. That is

1. BFE: As in [24], [35], we set \(\mu_{1}=\alpha/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 \(\Omega\).

2. BRFE: We set \(\mu_{1} = \alpha+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: \[\begin{align} \mathcal{L}(c,\theta_{\boldsymbol{x}})^{2} &=\left\|b(\boldsymbol{x})^{s}\left[b(\boldsymbol{x})^{-s}L_{\boldsymbol{x}}\Psi (\boldsymbol{x};c,\theta_{\boldsymbol{x}} ) -f_{\mathrm{reg}}(\boldsymbol{x}) \right] \right\|_{L^{2}(\Omega )}^{2}\\ &\approx\sum_{i_1=1}^{N_1}{\cdots \sum_{i_d=1}^{N_d}{\left( \prod_{k=1}^d{w_{i_k}^{(2s, 2s)}} \right) T^{2}\left( x_{1,i_1}^{(2s, 2s)},\dots ,x_{d,i_d}^{(2s, 2s)} \right)}}. \end{align}\]

Here \[T(\boldsymbol{y})=\sum_{j=1}^p{c_j\left[b(\boldsymbol{y})^{-s} L_{\boldsymbol{x}}\varphi_j(\boldsymbol{y};\theta_{\boldsymbol{x}} ) \right]}-f_{\mathrm{reg}}(\boldsymbol{y}),\] where \(\{x_{k,i_k}^{(2s, 2s)}\}\) and \(\{w_{i_k}^{(2s, 2s)}\}\) 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(\boldsymbol{x})\in H^1(\Omega)\) satisfying \[g(\boldsymbol{x})=u_b(\boldsymbol{x} ), \;\;\boldsymbol{x}\in \partial\Omega,\] is introduced. Then using the decomposition \(u(\boldsymbol{x}) = \widehat{u}(\boldsymbol{x}) + g(\boldsymbol{x} )\), the non-homogeneous boundary value problem is transformed to a homogeneous one for the solution \(\widehat{u}(\boldsymbol{x})\).

2.5 Trial spaces and loss function for time-space fractional PDEs↩︎

For the time-space fractional PDE, the trial function \(\Psi\) in 10 is then redefined as \[\label{eq:Psi95time} \Psi(\boldsymbol{x}, t; c, \theta) = \widehat{\Psi}(\boldsymbol{x}, t; c, \theta)+ u_0(\boldsymbol{x}),\quad \widehat{\Psi}(\boldsymbol{x}, t; c, \theta):=\sum_{j=1}^{p} c_j t^{\gamma} \phi_{t,j}(t) \varphi_j(\boldsymbol{x}).\tag{14}\] This design incorporates the factor \(t^{\gamma}\) into time basis functions \(\phi_{t,j}(t)\), following the approach proposed in [10]. 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 3 transforms the original problem into a homogeneous initial-boundary value problem for \(\widehat{\Psi}\) with RHS function \(\widehat{f}\). The corresponding loss function is then given by \[\label{eq:time95space95loss} \mathcal{L}(c,\theta )=\left\| \partial^\gamma \widehat{\Psi}(\boldsymbol{x},t;c,\theta)/\partial t^\gamma+L_ {\boldsymbol{x}}\widehat{\Psi}(\boldsymbol{x},t;c,\theta)-\widehat{f}(\boldsymbol{x},t) \right\|_{L^{2}(\Omega \times(0,T])}.\tag{15}\] 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 [10], [27].

3 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\ge 0\)) or the BRFE strategy (\(s<0\)) to choose the parameter \(\mu_1\). The remaining parameters \(\mu_j\) (\(j=2,\dots,p-1\)) are then obtained recursively from 11 , while \(\mu_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 14 into 3 , the original problem is then transformed into a homogeneous initial-boundary value problem for \(\widehat{\Psi}\), with the corresponding RHS function \(\widehat{f}\). To decompose the high-dimensional integrals appearing in the loss function into lower-dimensional temporal and spatial integrals, we express the time derivative of order \(\gamma\) of \(\widehat{\Psi}\), the operator \(L_{\boldsymbol{x}}\) applied to \(\widehat{\Psi}\), and \(\widehat{f}\) in separable forms as follows: \[\label{eq:separable} \left\{ \begin{align} \frac{\partial^{\gamma}\widehat{\Psi}(\boldsymbol{x},t)}{\partial t^{\gamma}} &=\sum_{j=1}^{p} c_j\frac{\partial^{\gamma}\left(t^{\gamma}\phi_{t,j}(t;\theta_t)\right)}{\partial t^{\gamma}} \varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}}) := \sum_{j=1}^{p} c_j L_{t,j}^{1}(t)L_{\boldsymbol{x},j}^{1}(\boldsymbol{x}),\\ L_{\boldsymbol{x}}\widehat{\Psi}(\boldsymbol{x},t) &=\sum_{j=1}^{p} c_jt^{\gamma}\phi_{t,j}(t;\theta_t)\, L_{\boldsymbol{x}}\varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}}) :=\sum_{j=1}^{p} c_j L_{t,j}^{2}(t)L_{\boldsymbol{x},j}^{2}(\boldsymbol{x}),\\ \widehat{f}(\boldsymbol{x},t) &:=\sum_{r=1}^{p_1} \widehat{f}_{t,r}(t)\widehat{f}_{\boldsymbol{x},r}(\boldsymbol{x}), \end{align} \right.\tag{16}\] where \(p_1\) denotes the number of terms in the spatiotemporally separable form of the source.

After that, we determine the boundary singularity index \(s\) from 13 . If \(s\ge 0\), we adopt the BFE strategy, setting \(\mu_1=\alpha/2\) and using Gauss quadrature for the trial function and the loss function. If \(s<0\), we adopt the BRFE strategy, setting \(\mu_1=\alpha+s\) and using Gauss-Jacobi quadrature with parameters \((2s,2s)\). After prescribing \(\mu_p\), we define the neural network subspace as \[V'_{p}(\theta) :=\operatorname{span}\left\{ t^{\gamma}\phi_{t,j}(t;\theta_t)\,\varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}}), \; j=1,\dots,p \right\},\: with \: \varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}})=b(\boldsymbol{x})^{\mu_j}\prod_{i=1}^d \phi_{i,j}(x_i;\theta_i).\]

The approximate solution is then sought in the form \(u_p=\widehat{u}_p+u_0\) with \(\widehat{u}_p\in V'_p\), by minimizing the residual in the least-squares sense. More precisely, we seek \(\widehat{u}_p\in V'_p\) such that \[\label{eq:least95squares} (L_{t,\boldsymbol{x}}\widehat{u}_p,L_{t,\boldsymbol{x}}\widehat{v}_p)=(\widehat{f},L_{t,\boldsymbol{x}}\widehat{v}_p), \quad \forall\,\widehat{v}_p\in V'_p,\tag{17}\] where \((\cdot,\cdot)\) denotes the \(L^2\) inner product over the space-time domain \(\Omega\times(0,T]\).

After the \(\ell\)-th training step, the trial function \(\Psi(\boldsymbol{x},t;c,\theta^{(\ell)})\) belongs to the subspace \[V'_{p}(\theta^{(\ell)}) := \mathrm{span} \left\{ t^\gamma \cdot \phi_{t,j}(t;\theta_t^{(\ell)})\cdot\varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}}^{(\ell)}), \; j=1,\dots,p \right\},\] where \(\phi_{t,j}(t;\theta_t^{(\ell)})\) and \(\varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}}^{(\ell)})\) are defined by 14 10 , respectively. Once the coefficients \(c\) are computed, the approximate solution \(\widehat{u}_p\) is obtained, and the subspace \(V'_p(\theta^{(\ell+1)})\) is updated through the loss optimization.

Assembling the discrete system on \(V'_p(\theta^{(\ell)})\) and exploiting the spatiotemporal separability of the basis and source representations yields the following stiffness matrix and load vector for the BFE (\(s\ge0\)) and BRFE (\(s<0\)) strategies: \[\begin{align} \tag{18} A_{m,n}^{(\ell)} &=\begin{cases} \displaystyle \sum_{j=1}^{2}\sum_{k=1}^{2} \left(L_{t,n}^{k},L_{t,m}^{j}\right)_t\cdot \left(L_{\boldsymbol{x},n}^{k},L_{\boldsymbol{x},m}^{j}\right)_{\boldsymbol{x}}, & s\ge 0,\\[1.2ex] \displaystyle \sum_{j=1}^{2}\sum_{k=1}^{2} \left(L_{t,n}^{k},L_{t,m}^{j}\right)_t\cdot \left(b^{-s}L_{\boldsymbol{x},n}^{k},b^{-s}L_{\boldsymbol{x},m}^{j}\right)_{\boldsymbol{x}, b^{2s}}, & s<0, \end{cases}\\ \tag{19} B_m^{(\ell)} &=\begin{cases} \displaystyle \sum_{r=1}^{p_1}\sum_{j=1}^{2} \left(\widehat{f}_{t,r},L_{t,m}^{j}\right)_t\cdot \left(\widehat{f}_{\boldsymbol{x},r},L_{\boldsymbol{x},m}^{j}\right)_{\boldsymbol{x}}, & s\ge 0,\\[1.2ex] \displaystyle \sum_{r=1}^{p_1}\sum_{j=1}^{2} \left(\widehat{f}_{t,r},L_{t,m}^{j}\right)_t\cdot \left(b^{-s}\widehat{f}_{\boldsymbol{x},r},\,b^{-s}L_{\boldsymbol{x},m}^{j}\right)_{\boldsymbol{x},\,b^{2s}}, & s<0, \end{cases} \end{align}\] for \(1\le m,n\le p\). Here, \((\cdot,\cdot)_t\) and \((\cdot,\cdot)_{\boldsymbol{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 \[\label{eq:inner95weight} (u,v)_{\boldsymbol{x},\,b^{2s}} :=\int_{\prod_{i=1}^{d}(a_i,b_i)} u(\boldsymbol{x})\,v(\boldsymbol{x})\,b(\boldsymbol{x})^{2s}\,d\boldsymbol{x}.\tag{20}\] The weighted integral in 20 is evaluated by Gauss-Jacobi quadrature with parameters \((2s,2s)\).

The parameters \(\theta^{(\ell)}\) determine the trial subspace \({V'_p}^{(\ell)}\), while the coefficient vector \(c^{(\ell)}\) determines the approximation within that subspace. Consequently, the optimization procedure can be naturally divided into two substeps.

First, with the neural network parameters \(\theta^{(\ell)}\) fixed, the optimal coefficient vector \(c^{(\ell+1)}\) is obtained in the least-squares sense from 17 by solving the linear system \[\label{linear95system95fpde} A^{(\ell)}c^{(\ell+1)}=B^{(\ell)}.\tag{21}\]

Second, with the coefficient vector \(c^{(\ell+1)}\) fixed, the neural network parameters \(\theta^{(\ell+1)}\) are updated by minimizing the loss function \[\label{eq:loss95ite} \mathcal{L}^{(\ell+1)}(c^{(\ell+1)},\theta^{(\ell)}) =\left(c^{(\ell+1)}\right)^{\top}A^{(\ell)}c^{(\ell+1)} -2\left(B^{(\ell)}\right)^{\top}c^{(\ell+1)} +(\widehat{f},\widehat{f})\tag{22}\] using gradient-based optimization methods such as Adam or L-BFGS. Thus, for fixed \(\theta^{(\ell)}\), minimizing the discrete loss is equivalent to solving the normal equations 21 .

The whole procedure is summarized in 3. This algorithm can be interpreted as an alternating optimization procedure for the unknown parameters \(c\) and \(\theta\) in the trial function 10 : the optimal coefficients \(c\) are determined in the least-squares sense, while the neural-network parameters are updated through the loss minimization.

3.2 Fractional Poisson equation↩︎

For the fPE defined in 1 , we introduce the neural-network trial subspace \[\label{eq:v95p} V_p(\theta_{\boldsymbol{x}}) :=\operatorname{span} \left\{ \varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}}), \; j=1,\dots,p \right\},\tag{23}\] where the approximate solution is found by minimizing the residual in the least-squares sense. More precisely, we seek \(u_p\in V_p(\theta_{\boldsymbol{x}})\) such that \[\label{eq:least95square95method} \left((-\Delta)^{\alpha/2}u_p,\,(-\Delta)^{\alpha/2}v_p\right)_{\boldsymbol{x}} =\left(f,\,(-\Delta)^{\alpha/2}v_p\right)_{\boldsymbol{x}}, \quad \forall\, v_p\in V_p(\theta_{\boldsymbol{x}}).\tag{24}\] After the \(\ell\)-th training step, the neural network \(\Psi(\boldsymbol{x};c,\theta_{\boldsymbol{x}}^{(\ell)})\) belongs to the subspace \[V_p^{(\ell)} :=\operatorname{span} \left\{\varphi_j(\boldsymbol{x};\theta_{\boldsymbol{x}}^{(\ell)}), \; j=1,\dots,p \right\}.\] By assembling the discrete system over \(V_p^{(\ell)}\), the stiffness matrix and load vector corresponding to the BFE and BRFE formulations can be written in the unified form \[\begin{align} \tag{25} A_{m,n}^{(\ell)} &=\begin{cases} \left((-\Delta)^{\alpha/2}\varphi_n^{(\ell)}, (-\Delta)^{\alpha/2}\varphi_m^{(\ell)}\right)_{\boldsymbol{x}}, & s\ge 0,\\[1ex] \left(b^{-s}(-\Delta)^{\alpha/2}(\varphi_n^{(\ell)}), b^{-s}(-\Delta)^{\alpha/2}(\varphi_m^{(\ell)})\right)_{\boldsymbol{x},\,b^{2s}}, & s<0, \end{cases}\\ \tag{26} B_m^{(\ell)} &=\begin{cases} \left(f,\,(-\Delta)^{\alpha/2}\varphi_m^{(\ell)}\right)_{\boldsymbol{x}}, & s\ge 0,\\[1ex] \left(b^{-s}f,\,b^{-s}(-\Delta)^{\alpha/2}(\varphi_m^{(\ell)})\right)_{\boldsymbol{x},\,b^{2s}}, & s<0, \end{cases} \end{align}\] for \(1\le m,n\le p\). Here, \((\cdot,\cdot)_{\boldsymbol{x},\,b^{2s}}\) is the weighted spatial inner product defined in 20 .

The parameters \(\theta_{\boldsymbol{x}}^{(\ell)}\) determine the trial subspace \(V_p^{(\ell)}\), while the coefficient vector \(c^{(\ell)}\) determines the approximation within this subspace. Accordingly, 4 can be interpreted as an alternating optimization procedure for fPEs. Since its implementation is analogous to 3, we omit the details here for brevity.

4 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.

The following two types of errors between the approximate solution \(\Psi (\boldsymbol{x},t;c,\theta)\) and the exact solution \(u\) are used to measure the convergence behavior and accuracy of the examples in this section:

• Relative \(L^{2}\) error \[e_{L^{2}}:=\frac{\|u-\Psi(\boldsymbol{x},t;c^{*},\theta^{*})\|_{L^{2}(\Omega\times(0,T])}}{\|u\|_{L^{2}(\Omega\times(0,T])}}.\] This norm is evaluated by high precision numerical quadrature.

• Relative \(L^{2}\) test error \[e_{\mathrm{test}}:=\frac{\sqrt{\sum_{k=1}^{K}\left(\Psi(\boldsymbol{x}^{k},t^{k};c,\theta)-u(\boldsymbol{x}^{k},t^{k})\right)^{2}}}{\sqrt{\sum_{k=1}^{K}\left(u(\boldsymbol{x}^{k},t^{k})\right)^{2}}},\] where the test points \(\{(\boldsymbol{x}^{k},t^{k})\}\) are placed on a uniform grid of \(300^{d+1}\) points over \(\Omega\times(0,T]\) for 1D and 2D, and on a uniform \(30^{4}\) grid for 3D.

For time-independent problems, the same definitions for the errors are used with the \(L^{2}\) norm taken over \(\Omega\) and test points \(\{\boldsymbol{x}^{k}\}\) uniformly distributed in \(\Omega\).

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 12 over the domain \(\prod_{i=1}^{d}[a_{i},b_{i}]\), we employ the quadrature point configurations specified in 1. For the rectangular domain, we define \[b(\boldsymbol{x}) = \max\left(\prod_{i=1}^d (x_i-a_i)(b_i-x_i), 0\right),\] which vanishes on the entire boundary, satisfying 9 .

In all experiments we set \(\mu_{p} = \mu_{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 4.

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 \(N_0 = 100\) Gauss quadrature points for the far-field radial integral.

We first consider the one-dimensional fPE on the interval \((0,1)\): \[\label{eq:1d95fra} (-\Delta)^{\alpha/2} u(x) = f(x), \quad x \in (0,1),\tag{27}\] 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)=\begin{cases} \left( \ell_{\alpha}(x,2)-2\ell_{\alpha}(x,3)+\ell_{\alpha}(x,4) \right)/(2\cos\left(\pi\alpha\right)),& u_{\text{exact}}=x^{2}(1-x)^{2},\\ \left( \ell_{\alpha}(x,3)-3\ell_{\alpha}(x,4)+3\ell_{\alpha}(x,5)-\ell_{\alpha}(x,6) \right)/(2\cos\left(\pi\alpha\right)), &u_{\text{exact}}=x^3(1-x)^3, \end{cases}\] where \[\ell_{\alpha}(x,n) = \frac{\Gamma(n+1)}{\Gamma(n+1-\alpha)} \left( x^{n-\alpha} + (1-x)^{n-\alpha} \right).\] Since \(s=2-\alpha>0\) and \(s=3-\alpha>0\) for the two cases, we use the BFE strategy. The results in 2 show fTNN provides a clear accuracy advantage over fPINN in this smooth setting.

We then consider a compactly supported function that is the sum of two bump functions with different exponents: \[\label{eq:compact95support95sum} u_{\text{exact}}(x) = (1 - x^{2})_+^{\beta_{1}} + (1 - x^{2})_+^{\beta_{2}}, \quad \text{for } x \in [-1,1],\quad \beta_{1} \leqslant \beta_{2}.\tag{28}\] The regularity of \(u_{\text{exact}}\) is determined by the smaller exponent \(\beta_{1}\), say \[u_{\text{exact}}\in\begin{cases} C^{\beta_{1} - 1,1}(\mathbb{R}) ,& \beta_{1} \in \mathbb{N},\\ C^{\lfloor \beta_{1} \rfloor, \beta_{1} - \lfloor \beta_{1} \rfloor}(\mathbb{R}), & \beta_{1} \notin \mathbb{N}. \end{cases}\]

For \(\beta_{i}\neq \alpha/2\), \[(-\Delta)^{\alpha/2}(1-x^{2})_{+}^{\beta_{i}} \sim C_{i}\,(1-x^{2})^{\beta_{i}-\alpha}, \quad |x|\to 1^{-}.\] Since \(\beta_{1}\le\beta_{2}\), the smallest exponent dominates, yielding \[f(x) \sim \text{const}\cdot (1-x^{2})^{\beta_{1}-\alpha}, \quad \beta_{1}<\beta_{2},\;\beta_{1}\neq\alpha/2.\] In the special case \(\beta_{1}=\beta_{2}= \alpha/2\), we have \[(-\Delta)^{\alpha/2} (1-x^{2})_{+}^{\alpha/2} = 2^{\alpha}\Gamma\left(\alpha/2+1/2\right)\,\Gamma\left(\alpha/2+1\right)/\sqrt{\pi},\] which is a non-zero constant, so \(f\) remains bounded and non-zero near the endpoints.

Consequently, the singularity index defined in 13 is \(s=\beta_{1}-\alpha\), except in the special case \(\beta_{1}=\beta_{2}=\alpha/2\) where \(s=0\). The appropriate strategy, BFE for \(s\ge 0\) and BRFE for \(s<0\), then follows as summarized in 3. Errors of fTNN are listed in 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\ge 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 5 and the error histories in 6 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 \(\alpha\) and exponent combinations \((\beta_{1},\beta_{2})\).

4.1.2 FPEs on the unit square/cube↩︎

We first consider the fPE 1 with an exact solution constructed as \[\label{eq:s95exact9523D} u_{\text{exact}}(\boldsymbol{x})=\begin{cases}\quad \sum\limits_{k=1}^{2}{\prod\limits_{i=1}^{2}{\left(x_{i}-x_{i}^{2} \right) ^{\beta _k}}}, &\boldsymbol{x}\in [0,1]^{2},\\ 100\sum\limits_{k=1}^{2}{\prod\limits_{i=1}^3{\left(x_{i}-x_{i}^{2} \right) ^{\beta _k}}}, &\boldsymbol{x}\in [0,1]^3.\\ \end{cases}\tag{29}\] where the exponents satisfy \(\beta_{1} \leqslant \beta_{2}\). Due to the complexity of deriving an explicit closed-form expression for the RHS function \(f(\boldsymbol{x}) = (-\Delta)^{\alpha / 2} u_{\text{exact}}(\boldsymbol{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 5.

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 = f_{1} + f_{2}\) with \[f_{k} = (-\Delta)^{ \alpha /2}\left(c_db(\boldsymbol{x})^{\beta_k}\right),\quad b(\boldsymbol{x}) = \prod\limits_{i=1}^d x_{i}(1-x_{i}),\] on \(\Omega = (0,1)^d\). Since \(\beta_{1} \le \beta_{2}\), the singularity of \(f\) is dominated by the term with the smaller exponent, i.e., \(f \sim f_{1}\) as \(\boldsymbol{x} \to \partial\Omega\). Near a face, say \(x_{i} = 0\), we have \(b(\boldsymbol{x}) \sim x_{i} B_{i}(\boldsymbol{x}')\), where \(\boldsymbol{x}'\) denotes the remaining coordinates and \(B_{i}(\boldsymbol{x}') = \prod_{j\neq i} x_{j}(1-x_{j})\) is smooth and non-zero on the face. Consequently, \(b(\boldsymbol{x})^{\beta_{1}} \sim x_{i}^{\beta_{1}} \left(B_{i}(\boldsymbol{x}')\right)^{\beta_{1}}\). A classical result for the fractional Laplacian (see [7]) yields the asymptotic expansion \[(-\Delta)^{ \alpha /2}\left(x_{i}^{\beta_{1}}\phi(\boldsymbol{x}')\right) \sim C(\beta_{1}, \alpha)\,\phi(\boldsymbol{x}')\,x_{i}^{\beta_{1}-\alpha}, \quad x_{i} \to 0^+,\] where \(C(\beta_{1},\alpha)\) is a constant. Applying this with \(\phi = B_{i}^{\beta_{1}}\) gives \[f(\boldsymbol{x}) \sim C(\beta_{1},\alpha)\,B_{i}(\boldsymbol{x}')^{\beta_{1}}\,x_{i}^{\beta_{1}-\alpha},\quad \beta_{1} \le \beta_{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(\boldsymbol{x})\) behaves like a constant (depending on the tangential coordinates) times \(\operatorname{dist}(\boldsymbol{x}, \partial\Omega)^{\beta_{1}-\alpha}\). To extract this maximal singular behavior, we require that the limit defined in 13 be finite and non-zero. Substituting the asymptotic forms shows that the exponent \(s=\beta_{1}-\alpha\). With this choice, the ratio tends to a non-zero smooth function on the faces, while it vanishes near edges and corners, hence the \(\limsup\) is indeed attained on the faces and equals a positive constant.

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 \(\beta_{1}=\alpha /2\) yields a constant RHS function for the term \((1-x^{2})^{ \alpha /2}\), in higher dimensions the product structure \(\prod_{i} (x_{i}-x_{i}^{2})^{\beta_k}\) does not produce a constant fractional Laplacian for any \(\beta_k\). This is because the cross-terms in the product lead to a genuine boundary singularity in \(f\) even when \(\beta_1 = \alpha/2\), justifying the exclusive use of \(\beta_1 - \alpha\) as the exponent \(s\) in the multidimensional setting. Consequently, the method choice is determined solely by the relative magnitudes of \(\beta_{1}\) and \(\alpha\):

• When \(\beta_{1} < \alpha\) (with any \(\beta_{2} \ge \beta_{1}\)), \(f(\boldsymbol{x})\) exhibits a genuine singularity \(\operatorname{dist}(\boldsymbol{x}, \partial\Omega)^{\beta_{1}-\alpha}\) near the boundary, necessitating the enriched approximation space of BRFE with the singular weight \(\mu_{1} = \beta_{1}\).

• When \(\beta_{1} \ge \alpha\) (with any \(\beta_{2} \ge \beta_{1}\)), \(f(\boldsymbol{x})\) remains bounded (or even smooth) on \(\overline{\Omega}\), and the BFE with \(\mu_{1} = \alpha /2\) suffices.

The relative errors \(e_{L^{2}}\) and \(e_{\mathrm{test}}\) of fTNN to solve 1 on the unit square/cube are listed in 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 \(n_{1}\). The average computational time per epoch of fTNN is drawn against \(n_{1}\) in 7, which shows a deterministic trend. Relative \(L^{2}\) errors in both two and three dimensional cases are shown in 7 8. The observed error reduction is consistent with the \(O(n_1^{-2})\) decay predicted by the spherical integration theory in 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 refining a deterministic angular discretization rather than from increasing the number of Monte Carlo samples as used in QE-MC-fPINN.

We then use fTNN to solve the fPE 1 with the RHS function \(f=1\) (clearly, \(s=0\)) and different \(\alpha\) in two dimensional case. Since the exact solutions are not known, we show the numerical solutions in 8 and corresponding loss curves in 9.

4.1.3 FPEs on the unit ball↩︎

We consider fPE 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 \[(-\Delta)^{\alpha / 2} \left(1-\|\boldsymbol{x}\|_2^{2}\right)^{1+\alpha/2} = 2^{\alpha}\,\Gamma(\alpha / 2+2)\,\Gamma(\alpha / 2+1)\,\left(1-(\alpha/2+1)\|\boldsymbol{x}\|_2^{2}\right).\] Using polar coordinates \(\boldsymbol{x}=(\rho\cos v,\rho\sin v)\) with \(\rho\in[0,1)\) and \(v\in[0,2\pi)\), the RHS function reads \[f(\rho,v)=2^{\alpha}\,\Gamma(\alpha / 2+2)\,\Gamma(\alpha / 2+1)\,\left(1-(\alpha / 2+1)\rho^{2}\right).\] Observe that \[\lim\limits_{\rho\to 1^{-}}f(\rho,v) = -2^{\alpha-1}\alpha\,\Gamma(\alpha / 2+2)\,\Gamma(\alpha / 2+1) \neq 0.\] Consequently, the limit \[\lim_{\rho\to 1^{-}}\frac{f(\rho,v)}{(1-\rho^{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: \[\Psi(\rho, v) = \sum\limits_{j=1}^p c_j \varphi_j(\rho, v).\] Then \[(-\Delta)^{\alpha /2} \varphi_j(\rho, v) = C_{2, \alpha} \left( I_{1,j}(\rho, v) + I_{2,j}(\rho, v) \right).\] The near-field integral is approximated as \[I_{1,j}(\rho,v) \approx r_{0}(\rho)^{-\alpha}\sum_{l=1}^{N_0} \sum_{k=1}^{N} w_lw_k^{(0,1-\alpha)} \frac{F_j(\rho, v, r_{0}(\rho)\tau_k^{(0,1-\alpha)}, \eta_l)}{(\tau_k^{(0,1-\alpha)})^{2}},\] where \(r_{0}(\rho)=1-\rho\), \(F_j(\rho, v, r, \eta) = 2\varphi_j(\rho, v) - \varphi_j(\rho', v') - \varphi_j(\rho'', v'')\), with \[\rho' = \sqrt{\rho^{2} + r^{2} + 2\rho r \cos(\eta - v)}, \quad \rho'' = \sqrt{\rho^{2} + r^{2} - 2\rho r \cos(\eta - v)},\] and \[v'=\arctan\left( \frac{\rho \sin v + r \sin\eta}{\rho \cos v + r \cos\eta} \right), \quad v'' = \arctan\left( \frac{\rho \sin v - r \sin\eta}{\rho \cos v - r \cos\eta} \right).\]

The far-field integral is then approximated as \[I_{2,j}(\rho,v) \approx \sum_{i=1}^{2} \sum_{k=1}^{2n_{i}}\frac{\pi}{n_{i}}Q_{i,j}(\rho,v, \eta_{ik}),\] with \(\eta_{ik}\) being discrete directional angles. \(Q_{1,j}\) and \(Q_{2,j}\) are defined as \[Q_{1,j}(\rho,v,\eta) \approx \left(d_{\boldsymbol{x}}(\eta)-r_{0}(\rho,v)\right) \sum_{m=1}^{N_0} w_m\frac{\varphi_j(\rho,v) - \varphi_j(\rho',v')}{\left(r(\eta,t_m)\right)^{1+\alpha}},\quad Q_{2,j}(\rho,v,\eta) = \frac{\varphi_j(\rho,v)}{\alpha \left(d_{\boldsymbol{x}}(\eta)\right)^{\alpha}},\] where \(r(\eta,t):=r_{0}(\rho,v)+\left(d_{\boldsymbol{x}}(\eta)-r_{0}(\rho,v)\right)t\) and the distance \(d_{\boldsymbol{x}}(\eta)\) from the point \((\rho,v)\) to the boundary along direction \(\eta\) is given by \[d_{\boldsymbol{x}}(\eta) = -\rho \cos(\eta - v) + \sqrt{1 - \rho^{2} \sin^{2}(\eta - v)}.\]

The above formulas provide the numerical approximation of fractional Laplacian in two-dimensional 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 [23]. The angular resolution \(n_1\) of fTNN is set to be \(32\) in 2D case and \(16\) in 3D case, while fPINN is configured as described in [9].

9 reports the relative \(L^2\) test errors \(e_{\mathrm{test}}\) of four different methods for solving two- and three-dimensional fPEs on the unit balls, with exact solutions \[\label{uesmooth} u_{\mathrm{exact}}=(1-\|\boldsymbol{x}\|^2)^{1+\alpha/2},\tag{30}\] and \[\label{uesingular} u_{\mathrm{exact}}=(1-\|\boldsymbol{x}\|^2)^{\alpha/2},\tag{31}\] for fractional orders \(\alpha=1.5\) and \(\alpha=1.9\).

For solution 30 with \(\alpha=1.5\), all five methods behave well, although fTNN achieves substantially lower errors than the other methods. For solution 30 with \(\alpha=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 31 in 2D case, the comparison results are similar to that for solution 30 with \(\alpha=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 10 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 11 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.

4.2 Fractional advection-diffusion equation↩︎

We next turn to time-dependent problems, for which 3 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 \(\boldsymbol{v} = (0.1,0.1)\) for \(d=2\) and \(\boldsymbol{v} = (0.1,0.1,0.1)\) for \(d=3\), which is a typical case of the fPDE defined in 3 . We employ the following quadrature scheme \[\begin{align} \Gamma(1-\gamma)\cdot{}_0^C D_t^{\gamma}\left( t^{\gamma}\phi_{t,j}(t) \right) = \int_0^1 (1-\tau)^{-\gamma} \tau^{\gamma-1} S_{1,j}(t,\tau) \mathrm{d}\tau \approx \sum_{k=1}^{N} w_k^{(-\gamma, \gamma-1)} S_{1,j}\left(t,\tau_k^{(-\gamma, \gamma-1)}\right), \end{align}\] for the time-fractional Caputo derivative, where \[S_{1,j}(t,\tau) = \gamma\phi_{t,j}(t\tau) + t\tau \phi_{t,j}'(t\tau),\] see [10] for further details.

For the spatial operator, we use the same deterministic discretization as in 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 \(\Psi\) in 10 as \[\label{ADE95trial} \Psi(\rho,\eta,t; c, \theta) = \widehat{\Psi}(\rho,\eta,t;c,\theta)+1-\rho^{2},\tag{32}\] with \[\widehat{\Psi}(\rho,\eta,t;c,\theta) = \sum\limits_{j=1}^{p}c_j t^{\gamma}\phi_{t,j}(t) (1-\rho^{2})^{\mu_{j}}\phi_{1,j}(\rho)\phi_{2,j}(\eta).\]

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 \[u_{\text{exact}}(\boldsymbol{x},t) = e^{-t}\left(1-\|\boldsymbol{x}\|_2^{2}\right)^{1+\alpha /2}.\] The corresponding source term \(f(\boldsymbol{x},t)\) for both two and three dimensions is provided in [9].

Substituting the decomposition 32 into 3 transforms the original problem into a homogeneous initial-boundary value problem for \(\widehat{\Psi}\), with the RHS function given by \[\begin{align} \widehat{f} = &2^{\alpha}\Gamma(\alpha/2+2)\Gamma(\alpha /2+1)(e^{-t}-1)\left(1-(\alpha /2+1)\rho^{2}\right) -(1-\rho^{2})^{1+\alpha /2}t^{1-\gamma}E_{1,2-\gamma}(-t) \\ &- 0.2\rho(\alpha /2+1)(e^{-t}-1)(1-\rho^{2})^{\alpha /2}(\cos\eta+\sin\eta), \end{align}\] where \(E_{\cdot, \cdot}(z)\) is the Mittag-Leffler function.

The relative \(L^{2}\) test errors of different methods are shown in 10. As seen in 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 12 reveal that fTNN exhibits smooth convergence histories. These results highlight the benefit of coupling deterministic spatial quadrature with the STSNN subspace formulation.

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 \[\label{uexact95long} u_{\text{exact}}(\boldsymbol{x},t)=(t^{\gamma_{2}}+t^{\gamma_{1}})(1-\|\boldsymbol{x}\|_2^{2})^{1+\alpha /2}.\tag{33}\] 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 11, the relative \(L^{2}\) 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.

We also compare the relative \(L^{2}\) test errors at \(T=100\) (for case A2 and case B1 in 11) of different methods in 12. The fPINN results are obtained using the DeepXDE implementation [36] with the \(L1\) scheme and Grünwald–Letnikov discretization configured as in [9]. The MC-fPINN, Improved MC-fPINN, and QE-MC-fPINN use the same hyperparameters as reported in [24]. All other solver configurations are kept the same as in their original references. The comparison result is similar to that in 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 long-time 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-MC-fPINN: 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 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.

The loss curves in 13 further indicate that the alternating STSNN subspace optimization is stable in the long-time regime: the relative \(L^{2}\) 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\,\mathrm{s}\) total runtime in 2D and about \(2100\,\mathrm{s}\) in 3D, which is acceptable given the difficulty of the nonlocal time-space operator.

5 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 systematic study of the angular resolution confirms that, after the radial singularity is removed, the error decays consistently with the \(O(n_1^{-2})\) 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 tensor-network architecture can simultaneously eliminate stochastic noise and mitigate ill-conditioning-two 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. School of Mathematics and Statistics & Hubei Key Laboratory of Mathematical Sciences, Central China Normal University, Wuhan 430079, China (maqingkui@mails.ccnu.edu.cn).↩︎

2. 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 (hhxie@lsec.cc.ac.cn).↩︎

3. School of Mathematics and Statistics & Key Laboratory of Nonlinear Analysis & Applications (Ministry of Education), Central China Normal University, Wuhan 430079, China (yinxb@ccnu.edu.cn).↩︎