Modeling Unknown Nonlocal PDE Systems via Flow Map Learning

# Modeling Unknown Nonlocal PDE Systems via Flow Map Learning
Source: [https://arxiv.org/html/2608.00400](https://arxiv.org/html/2608.00400)
\\AddToHook

package/hyperref/before

\\headers

Flow map learning of nonlocal PDEsZ\.Xu, Y\.Li, Y\.Zhang, and D\.Xiu

Zhongshu Xu\. Funding: This work was partially supported by AFOSR FA9550\-24\-1\-0237Ying Li11footnotemark:1Yanzhi ZhangDepartment of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409, USA\. Email:Dongbin Xiu11footnotemark:1

###### Abstract

Nonlocal partial differential equations arise in many applications but are often difficult to model and learn because of the presence of nonlocal operators\. We present a flow\-map learning \(FML\) framework for modeling unknown nonlocal PDEs directly from solution data\. Rather than learning or approximating the underlying nonlocal operators, the proposed approach learns the finite\-time evolution operator in either modal or nodal space\. Two complementary formulations are developed for spectral and grid\-based solution representations\. Numerical experiments on one\- and two\-dimensional fractional diffusion and wave equations demonstrate accurate and stable long\-time prediction using only short observation windows\. The proposed approach provides an effective data\-driven framework for learning unknown nonlocal dynamics without explicit evaluation of nonlocal operators\.

###### keywords:

Data\-driven modeling, flow\-map learning, nonlocal PDEs

## 1Introduction

Nonlocal partial differential equations \(PDEs\) have attracted considerable attention because of their ability to model long\-range interactions and anomalous transport phenomena that cannot be adequately described by classical local PDEs\. Representative applications include anomalous diffusion, peridynamics in solid mechanics, porous media flow, and other multiscale systems governed by nonlocal interactions\[du2012analysis,kirkpatrick2016,metzler2000random,benson2000application,silling2000reformulation,caffarelli2011nonlinear\]\. Despite substantial progress in the development of nonlocal mathematical models, deriving predictive governing equations from first principles remains difficult in many practical settings because of complex multiscale interactions, incomplete physical knowledge, or limited observational access\. Moreover, even when the governing equations are known, the associated nonlocal operators often introduce considerable analytical and computational complexity\[delia2020numerical,Duo2018,Duo2019,Zhou2024\]\. These challenges have motivated increasing interest in data\-driven approaches for modeling nonlocal dynamics directly from observations\.

Data\-driven methods for PDEs have received considerable attention in recent years and pursue several distinct objectives, One class seeks to discover the governing equations in explicit analytical form\. Such methods typically employ sparse regression over a prescribed dictionary of candidate differential operators and identify the terms that best fit the available data\[Rudy2017,Schaeffer2017,Long2019\]\. A second class focuses on approximating PDE solutions when the governing equations are known\. Physics\-informed neural networks \(PINNs\) and related methods incorporate differential\-equation residuals, boundary conditions, and other physical constraints into the training objective and have been applied to a broad range of forward and inverse PDE problems\[Raissi2019,Zhu2019,Sun2020,Geneva2020\]\.

Extensions of physics\-informed methods to fractional and nonlocal PDEs include fractional and nonlocal PINNs, high\-dimensional variants, and peridynamic formulations\[Pang2019,Pang2020,guo2022monte,Haghighat2021\]\. Other approaches reduce the burden of directly evaluating nonlocal residuals through variational formulations, auxiliary extensions, or alternative numerical representations of known fractional operators\[gu2022deep,hao2023neural\]\. Neural\-operator approaches have also been developed for learning nonlocal constitutive relations from paired field data, often with the nonlocal structure incorporated into the operator architecture\[jafarzadeh2024peridynamic\]\. Nevertheless, these methods generally require some representation of the underlying nonlocal mechanism during training, either through an explicit residual, a known reformulation, a prescribed architecture, or paired input\-output data generated by the governing system\. Since automatic differentiation does not directly apply to fractional and nonlocal operators\[lischke2020fractional,Nezza2012\], the training procedure may still require repeated evaluation or approximation of singular kernels, integral operators, or auxiliary discretizations\.

In this work, we pursue a different strategy based on Flow Map Learning \(FML\)\[ChurchillXiu\_FML2023\]\. Rather than identifying the governing equation or approximating its nonlocal operators, FML learns the finite\-time evolution operator, or flow map, directly from solution data\. Given two solution states separated by a fixed time interval, the method constructs a data\-driven mapping between them\. Once trained, the learned flow map can be applied recursively to generate long\-time predictions\. FML has been developed for autonomous and non\-autonomous dynamical systems\[qin2019data,Qin2021a\], unknown PDEs\[WU2020109307,CHEN2022110782,churchill2023dnn\], stochastic systems\[chen2024learning,xu2024modeling\], partially observed systems\[Fu2020,churchill2023dnn\], and fluid simulation and active control\[chen2026targeted,liu2026numerical\]\.

The objective of this paper is to extend FML to unknown nonlocal PDE systems and to investigate its effectiveness for long\-time prediction\. We develop two complementary formulations\. The first operates in modal space, where the solution is represented through coefficients in a finite\-dimensional basis\. This formulation is particularly natural for periodic problems involving fractional operators, since the fractional Laplacian admits a simple representation in Fourier space\. The second formulation operates directly in nodal space and learns the evolution of solution values on a spatial grid, making it applicable when a convenient modal representation is unavailable\. In both cases, the learned FML model uses only solution snapshots and does not require explicit knowledge of the governing nonlocal operator\.

The proposed methods are evaluated on one\- and two\-dimensional nonlocal diffusion and wave equations\. The numerical results demonstrate that both modal\- and nodal\-space FML can provide accurate and stable long\-time predictions from short training trajectories\. Additional parameter\-recovery experiments show that the predicted modal dynamics retain sufficient information to recover the underlying fractional orders over appropriate time intervals\. These results suggest that FML provides a practical data\-driven framework for modeling unknown nonlocal PDE systems without explicit approximation of their nonlocal operators\.

The remainder of the paper is organized as follows\. Section[2](https://arxiv.org/html/2608.00400#S2)introduces the problem setup and the observational data used for training\. Section[3](https://arxiv.org/html/2608.00400#S3)presents the modal\-space FML formulation, and Section[4](https://arxiv.org/html/2608.00400#S4)develops the nodal\-space FML formulation\. Section[5](https://arxiv.org/html/2608.00400#S5)contains numerical experiments for representative nonlocal diffusion and wave equations\. Conclusions are presented in Section[6](https://arxiv.org/html/2608.00400#S6)\.

## 2Problem Setup and Data

We consider an autonomous nonlocal partial differential equation of the form

\(1\)∂u∂t=ℒ​\(u\),\(x,t\)∈Ω×ℝ\+,\\frac\{\\partial u\}\{\\partial t\}=\\mathcal\{L\}\(u\),\\qquad\(x,t\)\\in\\Omega\\times\\mathbb\{R\}^\{\+\},subject to appropriate boundary and initial conditions

u​\(x,0\)=u0​\(x\),x∈Ω\.u\(x,0\)=u\_\{0\}\(x\),\\qquad x\\in\\Omega\.HereΩ⊂ℝd\\Omega\\subset\\mathbb\{R\}^\{d\}, ford=1,2,3d=1,2,3, is the spatial domain,u=u​\(x,t\)u=u\(x,t\)is the unknown solution, andℒ\\mathcal\{L\}denotes a possibly nonlinear nonlocal operator\. Throughout this paper, the explicit form ofℒ\\mathcal\{L\}is assumed to be unknown\.

### 2\.1Data

We shall assume that solution data are available at discrete time instances over a constant time stepΔ​t\>0\\Delta t\>0,

\(2\)0=t0<t1<⋯,tℓ\+1−tℓ=Δ​t,0=t\_\{0\}<t\_\{1\}<\\cdots,\\qquad t\_\{\\ell\+1\}\-t\_\{\\ell\}=\\Delta t,and a set of spatial grid points

\(3\)XN=\{x1,…,xN\}⊂Ω\.X\_\{N\}=\\\{x\_\{1\},\\ldots,x\_\{N\}\\\}\\subset\\Omega\.We use vector notation to represent the solution

\(4\)𝐮​\(t\)=u​\(x,t\)\|XN=\(u​\(x1,t\),…,u​\(xN,t\)\)T∈ℝN\.\\mathbf\{u\}\(t\)=u\(x,t\)\|\_\{X\_\{N\}\}=\\left\(u\(x\_\{1\},t\),\\ldots,u\(x\_\{N\},t\)\\right\)^\{T\}\\in\\mathbb\{R\}^\{N\}\.The available solution data consist ofNtraj≥1N\_\{\\rm traj\}\\geq 1solution trajectories

\(5\)𝐮\(m\)​\(tℓ\),ℓ=0,…,L,m=1,…,Ntraj,\\mathbf\{u\}^\{\(m\)\}\(t\_\{\\ell\}\),\\qquad\\ell=0,\\ldots,L,\\qquad m=1,\\ldots,N\_\{\\rm traj\},generated from different initial conditions\. HereL≥1L\\geq 1is the length of each trajectory and assumed to be a constant for notational convenience\. The collection of these trajectory data is our raw data set

\(6\)𝒮=\{𝐮0\(m\),𝐮1\(m\),…,𝐮L\(m\)\}m=1Ntraj,\\mathcal\{S\}=\\left\\\{\\mathbf\{u\}^\{\(m\)\}\_\{0\},\\mathbf\{u\}^\{\(m\)\}\_\{1\},\\dots,\\mathbf\{u\}^\{\(m\)\}\_\{L\}\\right\\\}\_\{m=1\}^\{N\_\{\\rm traj\}\},where the subscript0≤ℓ≤L0\\leq\\ell\\leq Lprovides the time index\. Consequently, the time stampstℓt\_\{\\ell\}are not required\.

### 2\.2Learning Objective

For the unknown nonlinear problem \([1](https://arxiv.org/html/2608.00400#S2.E1)\), our objective is not to recover or approximateℒ\\mathcal\{L\}itself, but rather to learn the evolution of the solution from the observational data \([6](https://arxiv.org/html/2608.00400#S2.E6)\)\.

LetEΔ​tE\_\{\\Delta t\}denote the exact flow map of \([1](https://arxiv.org/html/2608.00400#S2.E1)\) over a time intervalΔ​t\>0\\Delta t\>0:

\(7\)u​\(⋅,t\+Δ​t\)=EΔ​t​\(u​\(⋅,t\)\)\.u\(\\cdot,t\+\\Delta t\)=E\_\{\\Delta t\}\\left\(u\(\\cdot,t\)\\right\)\.FML seeks an approximation of this evolution operator directly from data\. Once such an approximation is learned, long\-time prediction can be performed by recursively applying the learned flow map operator, without requiring knowledge of the governing equation \([1](https://arxiv.org/html/2608.00400#S2.E1)\)\.

In the following two sections, we present two FML formulations, one in modal space and the other in nodal space\. The modal formulation naturally aligns with the Fourier representation of fractional operators\. The nodal formulation works directly with solution values on spatial grids\.

## 3Learning Nonlocal PDEs in Modal Space

### 3\.1Motivation

The fractional Laplacian plays a significant role in modeling anomalous diffusion and has become a fundamental nonlocal operator for describing long\-range interactions in a wide range of scientific applications\. It admits the following pseudo\-differential representation\[kilbas1993fractional,Nezza2012\]:

\(8\)\(−Δ\)α/2​u=ℱ−1​\[\|ξ\|α​ℱ​\[u\]\],α\>0,\(\-\\Delta\)^\{\\alpha/2\}u=\\mathcal\{F\}^\{\-1\}\\left\[\|\\xi\|^\{\\alpha\}\\mathcal\{F\}\[u\]\\right\],\\qquad\\alpha\>0,whereℱ\\mathcal\{F\}denotes the Fourier transform andξ\\xiis the frequency variable\. This representation reveals an important feature of the fractional Laplacian: Although it is nonlocal in physical space, it becomes a multiplication operator in frequency space\. Modal coefficients therefore provide a natural representation for the dynamics of many nonlocal PDEs\.

To further illustrate this, consider the fractional diffusion equation

\(9\)∂tu​\(x,t\)=−\(−Δ\)α/2​u,x∈ℝd\.\\partial\_\{t\}u\(x,t\)=\-\(\-\\Delta\)^\{\\alpha/2\}u,\\qquad x\\in\{\\mathbb\{R\}\}^\{d\}\.Expanding the solution in Fourier modes

\(10\)u​\(x,t\)=∑k∈ℤdu^k​\(t\)​ei​k⋅xu\(x,t\)=\\sum\_\{k\\in\\mathbb\{Z\}^\{d\}\}\\widehat\{u\}\_\{k\}\(t\)e^\{\\mathrm\{i\}k\\cdot x\}withi=−1\\mathrm\{i\}=\\sqrt\{\-1\}being the imaginary unit and using the definition \([8](https://arxiv.org/html/2608.00400#S3.E8)\), we obtain

\(11\)d​u^k​\(t\)d​t=−\|k\|α​u^k,k∈ℤd\.\\frac\{d\\widehat\{u\}\_\{k\}\(t\)\}\{dt\}=\-\|k\|^\{\\alpha\}\\widehat\{u\}\_\{k\},\\qquad k\\in\{\\mathbb\{Z\}\}^\{d\}\.Thus, in modal coordinates, the PDE is transformed into a dynamical system for the Fourier coefficients, analogous to the modal evolution obtained for classical local PDEs\. In fact, whenα=2\\alpha=2, the definition \([8](https://arxiv.org/html/2608.00400#S3.E8)\) reduces to the spectral representation of the classical negative Laplacian and \([9](https://arxiv.org/html/2608.00400#S3.E9)\) corresponds to the standard diffusion equation\. The exact one\-step evolution is

\(12\)u^k​\(t\+Δ​t\)=exp⁡\(−\|k\|α​Δ​t\)​u^k​\(t\),α\>0\.\\widehat\{u\}\_\{k\}\(t\+\\Delta t\)=\\exp\\left\(\-\|k\|^\{\\alpha\}\\Delta t\\right\)\\widehat\{u\}\_\{k\}\(t\),\\qquad\\alpha\>0\.This simple example demonstrates that, despite its nonlocal nature in physical space, modal\-space learning can still be natural for nonlocal diffusion\-type problems\. For more general nonlocal PDEs, the modal dynamics can be far more complex and unknown, but the evolution map can still be learned from data\.

### 3\.2Finite\-dimensional modal representation

LetVVbe a function space containing the solution of \([1](https://arxiv.org/html/2608.00400#S2.E1)\), and let

\(13\)Vn=span​\{ϕ1,…,ϕn\}⊂VV\_\{n\}=\{\\rm span\}\\\{\\phi\_\{1\},\\ldots,\\phi\_\{n\}\\\}\\subset Vbe annn\-dimensional approximation space, where\{ϕj\}j=1n\\\{\\phi\_\{j\}\\\}\_\{j=1\}^\{n\}are the basis functions\. We approximate the solutionu∈Vu\\in Vbyun∈Vnu\_\{n\}\\in V\_\{n\}through a projection operatorPn:V→VnP\_\{n\}:V\\rightarrow V\_\{n\}\. The approximate solutionunu\_\{n\}can be expressed as

\(14\)un​\(x,t\)=Pn​u​\(x,t\)=∑j=1nvj​\(t\)​ϕj​\(x\),u\_\{n\}\(x,t\)=P\_\{n\}u\(x,t\)=\\sum\_\{j=1\}^\{n\}v\_\{j\}\(t\)\\phi\_\{j\}\(x\),where

\(15\)𝐯​\(t\)=\(v1​\(t\),…,vn​\(t\)\)T∈ℝn\\mathbf\{v\}\(t\)=\(v\_\{1\}\(t\),\\ldots,v\_\{n\}\(t\)\)^\{T\}\\in\\mathbb\{R\}^\{n\}is the vector of modal coefficients\. The modal coefficients depend on the choice of the approximation spaceVnV\_\{n\}and the projection operatorPnP\_\{n\}\.

Define the coefficient\-to\-function map

\(16\)Π:ℝn→Vn,Π​𝐯=∑j=1nvj​ϕj\.\\Pi:\\mathbb\{R\}^\{n\}\\rightarrow V\_\{n\},\\qquad\\Pi\\mathbf\{v\}=\\sum\_\{j=1\}^\{n\}v\_\{j\}\\phi\_\{j\}\.When the basis functions are linearly independent,Π\\Piis a one\-to\-one correspondence between coefficient vectors inℝn\\mathbb\{R\}^\{n\}and functions inVnV\_\{n\}\.

Consider the infinite\-dimensional flow mapEΔ​tE\_\{\\Delta t\}\([7](https://arxiv.org/html/2608.00400#S2.E7)\):

\(17\)EΔ​t:V→V,u​\(⋅,t\+Δ​t\)=EΔ​t​\(u​\(⋅,t\)\)\.E\_\{\\Delta t\}:V\\to V,\\qquad u\(\\cdot,t\+\\Delta t\)=E\_\{\\Delta t\}\(u\(\\cdot,t\)\)\.Upon projecting the solutions via \([14](https://arxiv.org/html/2608.00400#S3.E14)\) onto the finite\-dimensional spaceVnV\_\{n\}, i\.e\.,

\(18\)un​\(⋅,t\)=Pn​u​\(⋅,t\),un​\(⋅,t\+Δ​t\)=Pn​u​\(⋅,t\+Δ​t\),u\_\{n\}\(\\cdot,t\)=P\_\{n\}u\(\\cdot,t\),\\qquad u\_\{n\}\(\\cdot,t\+\\Delta t\)=P\_\{n\}u\(\\cdot,t\+\\Delta t\),one can define a finite dimensional evolution operatorE~Δ​t\\widetilde\{E\}\_\{\\Delta t\}:

\(19\)E~Δ​t:Vn→Vn,un​\(⋅,t\+Δ​t\)=E~Δ​t​\(un​\(⋅,t\)\)\.\\widetilde\{E\}\_\{\\Delta t\}:V\_\{n\}\\to V\_\{n\},\\qquad u\_\{n\}\(\\cdot,t\+\\Delta t\)=\\widetilde\{E\}\_\{\\Delta t\}\(u\_\{n\}\(\\cdot,t\)\)\.Subsequently, in modal coordinates, this finite\-dimensional evolution is represented by a map

\(20\)MΔ​t:ℝn→ℝn,𝐯​\(t\+Δ​t\)=MΔ​t​\(𝐯​\(t\)\),M\_\{\\Delta t\}:\\mathbb\{R\}^\{n\}\\rightarrow\\mathbb\{R\}^\{n\},\\qquad\\mathbf\{v\}\(t\+\\Delta t\)=M\_\{\\Delta t\}\(\\mathbf\{v\}\(t\)\),which satisfies

\(21\)MΔ​t=Π−1​E~Δ​t​Π\.M\_\{\\Delta t\}=\\Pi^\{\-1\}\\widetilde\{E\}\_\{\\Delta t\}\\Pi\.

### 3\.3Training of modal\-space FML

The purpose of modal\-space FML is to approximateMΔ​tM\_\{\\Delta t\}\([21](https://arxiv.org/html/2608.00400#S3.E21)\) directly from the data \([6](https://arxiv.org/html/2608.00400#S2.E6)\)\. The derivation \([18](https://arxiv.org/html/2608.00400#S3.E18)\)–\([20](https://arxiv.org/html/2608.00400#S3.E20)\) readily gives us the following learning procedure\.

#### 3\.3\.1Training Data Set

Given the training data set of the solution snapshots𝐮\\mathbf\{u\}\([6](https://arxiv.org/html/2608.00400#S2.E6)\), our first task is to convert it into a training data set for the modal coefficients𝐯\\mathbf\{v\}\. This is accomplished by projecting each snapshot𝐮\\mathbf\{u\}ontoVnV\_\{n\}via the projectionPnP\_\{n\}\([14](https://arxiv.org/html/2608.00400#S3.E14)\)\. However, since𝐮∈ℝN\\mathbf\{u\}\\in\\mathbb\{R\}^\{N\}is finite dimensional, we need to define a discrete projectionP^n:ℝN→Vn\\hat\{P\}\_\{n\}:\\mathbb\{R\}^\{N\}\\to V\_\{n\}such that

\(22\)u^n​\(x,t\)=P^n​𝐮​\(t\)=∑j=1nv^j​\(t\)​ϕj​\(x\),\\hat\{u\}\_\{n\}\(x,t\)=\\hat\{P\}\_\{n\}\\mathbf\{u\}\(t\)=\\sum\_\{j=1\}^\{n\}\\hat\{v\}\_\{j\}\(t\)\\phi\_\{j\}\(x\),where

\(23\)𝐯^​\(t\)=\(v^1​\(t\),…,v^n​\(t\)\)T∈ℝn\\hat\{\\mathbf\{v\}\}\(t\)=\(\\hat\{v\}\_\{1\}\(t\),\\ldots,\\hat\{v\}\_\{n\}\(t\)\)^\{T\}\\in\\mathbb\{R\}^\{n\}is the vector of modal coefficients\. These coefficients are determined by the choice of the discrete projection operatorP^n\\hat\{P\}\_\{n\}, which is essentially a numerical approximation ofPnP\_\{n\}via the finite gridXNX\_\{N\}\. The resulting difference betweenunu\_\{n\}andu^n\\hat\{u\}\_\{n\}, resp\. between𝐯\\mathbf\{v\}and𝐯^\\hat\{\\mathbf\{v\}\}, is the so\-called aliasing error\. Upon applying the discrete projectionP^n\\hat\{P\}\_\{n\}to each solution snapshot in the training data set \([6](https://arxiv.org/html/2608.00400#S2.E6)\) and thenΠ−1\\Pi^\{\-1\}\([16](https://arxiv.org/html/2608.00400#S3.E16)\) to retract the \(approximate\) modal coefficients, we obtain the following data set

\(24\)𝒟=\{𝐯^0\(m\),𝐯^1\(m\),…,𝐯^L\(m\)\}m=1Ntraj,\\mathcal\{D\}=\\left\\\{\\hat\{\\mathbf\{v\}\}^\{\(m\)\}\_\{0\},\\hat\{\\mathbf\{v\}\}^\{\(m\)\}\_\{1\},\\ldots,\\hat\{\\mathbf\{v\}\}^\{\(m\)\}\_\{L\}\\right\\\}\_\{m=1\}^\{N\_\{\\rm traj\}\},which will be used as our training data set for model\-space learning\.

#### 3\.3\.2FML Model Construction and Prediction

In modal\-space FML, we seek an operatorM~Δ​t:ℝn→ℝn\\widetilde\{M\}\_\{\\Delta t\}:\\mathbb\{R\}^\{n\}\\to\\mathbb\{R\}^\{n\}as an approximation to the true unknown modal flow mapMΔ​tM\_\{\\Delta t\}\([20](https://arxiv.org/html/2608.00400#S3.E20)\)\. More specifically, we use a residual neural network \(ResNet\) to represent it, i\.e\.,

\(25\)M~Δ​t=𝐈n\+𝒩​\(⋅;Θ\),\\widetilde\{M\}\_\{\\Delta t\}=\\mathbf\{I\}\_\{n\}\+\\mathcal\{N\}\(\\cdot;\\Theta\),where𝐈n\\mathbf\{I\}\_\{n\}is the identity matrix of sizen×nn\\times nand𝒩:ℝn→ℝn\\mathcal\{N\}:\\mathbb\{R\}^\{n\}\\to\\mathbb\{R\}^\{n\}is the mapping operator of a fully connected DNN withnninput nodes,nnoutput nodes, and hyperparametersΘ\\Theta\.

To learn this operator, we use the training data set \([24](https://arxiv.org/html/2608.00400#S3.E24)\)\. We randomly sampleNpair\>1N\_\{\\rm pair\}\>1pairs of consecutive entries from \([24](https://arxiv.org/html/2608.00400#S3.E24)\)

\(26\)\(𝐯^ik\(mk\),𝐯^ik\+1\(mk\)\),k=1,…,Npair\.\\left\(\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\},\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+1\}\\right\),\\qquad k=1,\\dots,N\_\{\\rm pair\}\.Typically, such random sampling is conducted via uniform distribution over the subscript indices of \([24](https://arxiv.org/html/2608.00400#S3.E24)\)\. We then seek the operator \([25](https://arxiv.org/html/2608.00400#S3.E25)\) by minimizing the following mean\-squared loss

\(27\)ℒ1​\(Θ\)=1Npair​∑k=1Npair‖M~Δ​t​\(𝐯^ik\(mk\);Θ\)−𝐯^ik\+1\(mk\)‖22\.\\mathcal\{L\}\_\{1\}\(\\Theta\)=\\frac\{1\}\{N\_\{\\rm pair\}\}\\sum\_\{k=1\}^\{N\_\{\\rm pair\}\}\\left\\\|\\widetilde\{M\}\_\{\\Delta t\}\\left\(\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\};\\Theta\\right\)\-\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+1\}\\right\\\|\_\{2\}^\{2\}\.In practice, multi\-step loss \(or, roll\-over loss\) is preferred because it enhances the numerical stability and accuracy of the trained model \(\[ChurchillXiu\_FML2023\]\)\. Letnr≥1n\_\{r\}\\geq 1be the number of multi\-step loss\. We then randomly sample, from the data set \([24](https://arxiv.org/html/2608.00400#S3.E24)\),Nseq\>1N\_\{\\rm seq\}\>1sequences of consecutive entries,

\(28\)\(𝐯^ik\(mk\),𝐯^ik\+1\(mk\),⋯,𝐯^ik\+nr\(mk\)\),k=1,…,Nseq\.\\left\(\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\},\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+1\},\\cdots,\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+n\_\{r\}\}\\right\),\\qquad k=1,\\dots,N\_\{\\rm seq\}\.The multi\-step loss function is defined as

\(29\)ℒnr​\(Θ\)=1Nseq​∑k=1Nseq1nr​∑j=0nr−1‖M~Δ​t​\(𝐯~ik\+j\(mk\);Θ\)−𝐯^ik\+j\+1\(mk\)‖22,\\mathcal\{L\}\_\{n\_\{r\}\}\(\\Theta\)=\\frac\{1\}\{N\_\{\\rm seq\}\}\\sum\_\{k=1\}^\{N\_\{\\rm seq\}\}\\frac\{1\}\{n\_\{r\}\}\\sum\_\{j=0\}^\{n\_\{r\}\-1\}\\left\\\|\\widetilde\{M\}\_\{\\Delta t\}\\left\(\\widetilde\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+j\};\\Theta\\right\)\-\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+j\+1\}\\right\\\|\_\{2\}^\{2\},where

\(30\)𝐯~ik\(mk\)=𝐯^ik\(mk\),M~Δ​t​\(𝐯~ik\+j\(mk\);Θ\)=𝐯~ik\+j\+1\(mk\),j=0,…,nr−1\.\\widetilde\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\}=\\hat\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\},\\qquad\\widetilde\{M\}\_\{\\Delta t\}\\left\(\\widetilde\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+j\};\\Theta\\right\)=\\widetilde\{\\mathbf\{v\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+j\+1\},\\quad j=0,\\dots,n\_\{r\}\-1\.Whennr=1n\_\{r\}=1, this loss becomes the one\-step loss \([27](https://arxiv.org/html/2608.00400#S3.E27)\)\.

Upon minimizing the loss to a satisfactorily low level, the FML model training is considered finished, and the hyperparametersΘ\\Thetabecome fixed\. We then obtain an FML model for system prediction\. Given a new initial condition𝐯0\\mathbf\{v\}\_\{0\}, the FML prediction follows

\(31\)\{𝐯~0=𝐯0,𝐯~k\+1=M~Δ​t​\(𝐯~k\),k=0,1,…\\left\\\{\\begin\{aligned\} &\\widetilde\{\\mathbf\{v\}\}\_\{0\}=\\mathbf\{v\}\_\{0\},\\\\ &\\widetilde\{\\mathbf\{v\}\}\_\{k\+1\}=\\widetilde\{M\}\_\{\\Delta t\}\(\\widetilde\{\\mathbf\{v\}\}\_\{k\}\),\\quad\\qquad k=0,1,\\dots\\end\{aligned\}\\right\.The solution at any time instance can be obtained by applying \([16](https://arxiv.org/html/2608.00400#S3.E16)\), i\.e\.,u~n​\(x,tk\)=Π​𝐯~k\\widetilde\{u\}\_\{n\}\(x,t\_\{k\}\)=\\Pi\\widetilde\{\\mathbf\{v\}\}\_\{k\},∀k\\forall k\.

## 4Learning Nonlocal PDEs in Nodal Space

### 4\.1Motivation

The modal\-space formulation is natural when an appropriate basis is available and the solution data can be readily projected using the basis\. In many applications with complex geometry, modal\-space representation can be cumbersome, if not impossible\. It is therefore useful to construct a learning framework directly in physical space\.

Consider the solution vector𝐮​\(t\)\\mathbf\{u\}\(t\)\([4](https://arxiv.org/html/2608.00400#S2.E4)\) over the spatial gridXNX\_\{N\}\([3](https://arxiv.org/html/2608.00400#S2.E3)\), its semi\-discrete evolution follows

\(32\)d​𝐮​\(t\)d​t=𝐅​\(𝐮\),\\frac\{d\\mathbf\{u\}\(t\)\}\{dt\}=\\mathbf\{F\}\(\\mathbf\{u\}\),where𝐅:ℝN→ℝN\\mathbf\{F\}:\\mathbb\{R\}^\{N\}\\rightarrow\\mathbb\{R\}^\{N\}is an unknown evolution operator containing discrete representations of various operators in the unknown governing PDE equation \([1](https://arxiv.org/html/2608.00400#S2.E1)\)\. For nonlocal PDEs involving nonlocal differential operators, their discretization can be expressed in terms of matrix\-vector products\. For example, the fractional Laplacian may be approximated in nodal space as\[Duo2018,Duo2019,Zhou2024\]

\(33\)\(−Δ\)α/2​u≈𝐃α​𝐮,\(\-\\Delta\)^\{\\alpha/2\}u\\approx\\mathbf\{D\}\_\{\\alpha\}\\mathbf\{u\},where𝐃α∈ℝN×N\\mathbf\{D\}\_\{\\alpha\}\\in\\mathbb\{R\}^\{N\\times N\}is a dense differentiation matrix\.

Integrating the semi\-discrete system \([32](https://arxiv.org/html/2608.00400#S4.E32)\) over one time step\[tn,tn\+1\]\[t\_\{n\},t\_\{n\+1\}\]gives the exact identity

\(34\)𝐮​\(tn\+1\)=𝐮​\(tn\)\+∫tntn\+1𝐅​\(𝐮​\(t\),𝐃α1​𝐮​\(t\),𝐃α2​𝐮​\(t\),⋯\)​𝑑t,\\mathbf\{u\}\(t\_\{n\+1\}\)=\\mathbf\{u\}\(t\_\{n\}\)\+\\int\_\{t\_\{n\}\}^\{t\_\{n\+1\}\}\\mathbf\{F\}\(\\mathbf\{u\}\(t\),\\mathbf\{D\}\_\{\\alpha\_\{1\}\}\\mathbf\{u\}\(t\),\\mathbf\{D\}\_\{\\alpha\_\{2\}\}\\mathbf\{u\}\(t\),\\cdots\)\\,dt,where𝐃αj\\mathbf\{D\}\_\{\\alpha\_\{j\}\},j=1,…j=1,\\dots, are the differentiation matrices for various differential operators\. The formulation is Markovian and effectively defines a mapping

\(35\)𝐆Δ​t:ℝN→ℝN,𝐮​\(tn\+1\)=𝐆Δ​t​\(𝐮​\(tn\)\)\.\\mathbf\{G\}\_\{\\Delta t\}:\\mathbb\{R\}^\{N\}\\to\\mathbb\{R\}^\{N\},\\qquad\\mathbf\{u\}\(t\_\{n\+1\}\)=\\mathbf\{G\}\_\{\\Delta t\}\(\\mathbf\{u\}\(t\_\{n\}\)\)\.We remark again that this is an exact formulation and not an Euler forward approximation of \([32](https://arxiv.org/html/2608.00400#S4.E32)\)\. Hereafter, we shall refer to𝐆Δ​t\\mathbf\{G\}\_\{\\Delta t\}as the nodal flow map\.

### 4\.2Nodal\-space FML model

We approximate the nodal flow map𝐆Δ​t\\mathbf\{G\}\_\{\\Delta t\}by a residual neural network

\(36\)𝐆~Δ​t=𝐈N\+𝒩​\(⋅;Θ\),\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}=\\mathbf\{I\}\_\{N\}\+\\mathcal\{N\}\(\\cdot;\\Theta\),where𝐈N\\mathbf\{I\}\_\{N\}is the identity matrix of sizeN×NN\\times Nand𝒩:ℝN→ℝN\\mathcal\{N\}:\\mathbb\{R\}^\{N\}\\rightarrow\\mathbb\{R\}^\{N\}is the mapping operator of a DNN withNNinput nodes,NNoutput nodes, and parametersΘ\\Theta\.

Following\[CHEN2022110782\], we take𝒩\\mathcal\{N\}in \([36](https://arxiv.org/html/2608.00400#S4.E36)\) to be a disassembly–assembly network, whose structure is illustrated in Figure[1](https://arxiv.org/html/2608.00400#S4.F1)\. It consists of a disassembly block, an assembly layer, and an output layer\.

The disassembly block comprisesJ≥1J\\geq 1fully connected feedforward channels operating in parallel, each withndn\_\{d\}hidden layers ofnwn\_\{w\}neurons\. Let

\(37\)ℱi:ℝN→ℝnw,i=1,…,J,\\mathcal\{F\}\_\{i\}:\\mathbb\{R\}^\{N\}\\rightarrow\\mathbb\{R\}^\{n\_\{w\}\},\\qquad i=1,\\ldots,J,denote the operator defined by theiith channel\. Collecting theJJchannels, the disassembly block defines the operator

\(38\)𝒟=⨂1≤i≤Jℱi:ℝN→ℝnw×J\.\\mathcal\{D\}=\\bigotimes\_\{1\\leq i\\leq J\}\\mathcal\{F\}\_\{i\}:\\mathbb\{R\}^\{N\}\\rightarrow\\mathbb\{R\}^\{n\_\{w\}\\times J\}\.The assembly layer is a fully connected feedforward network withnan\_\{a\}hidden layers ofJJneurons\. It operates componentwise on the output of the disassembly block, taking aJJ\-dimensional input to a scalar,

\(39\)𝒜:ℝJ→ℝ,\\mathcal\{A\}:\\mathbb\{R\}^\{J\}\\rightarrow\\mathbb\{R\},so that thenw×Jn\_\{w\}\\times Joutput of the disassembly block is reduced to annwn\_\{w\}\-vector\. The output layer

\(40\)𝒪:ℝnw→ℝN\\mathcal\{O\}:\\mathbb\{R\}^\{n\_\{w\}\}\\rightarrow\\mathbb\{R\}^\{N\}is a fully connected feedforward network that maps the assemblednwn\_\{w\}\-vector back to theNN\-dimensional nodal space\. With these operators, the network in \([36](https://arxiv.org/html/2608.00400#S4.E36)\) reads

\(41\)𝐆~Δ​t=𝐈N\+𝒩=𝐈N\+𝒪∘𝒜∘𝒟,\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}=\\mathbf\{I\}\_\{N\}\+\\mathcal\{N\}=\\mathbf\{I\}\_\{N\}\+\\mathcal\{O\}\\circ\\mathcal\{A\}\\circ\\mathcal\{D\},where the residual term𝒪∘𝒜∘𝒟\\mathcal\{O\}\\circ\\mathcal\{A\}\\circ\\mathcal\{D\}gives the one\-step increment of the nodal flow map\. The construction of the disassembly\-assembly structure is motivated by the representation of the various unknown differentiation matrices in \([34](https://arxiv.org/html/2608.00400#S4.E34)\)\. For the mathematical details, see\[CHEN2022110782\]\.

![Refer to caption](https://arxiv.org/html/2608.00400v1/x1.png)Figure 1:Structure of the disassembly–assembly network used for nodal\-space FML\.An essential feature of the proposed DNN is that there is no need to explicitly construct a fractional differentiation matrix or evaluate a singular integral operator\. All the nonlocal dependence of the unknown system is embedded in the nodal flow map𝐆Δ​t\\mathbf\{G\}\_\{\\Delta t\}, to be learned from the data \([6](https://arxiv.org/html/2608.00400#S2.E6)\)\.

### 4\.3Training of nodal\-space FML

The training of the nodal FML operator \([36](https://arxiv.org/html/2608.00400#S4.E36)\) can be accomplished by directly using the training data set \([6](https://arxiv.org/html/2608.00400#S2.E6)\)\. Again, letnr≥1n\_\{r\}\\geq 1be the number of multi\-step loss\. We randomly sampleNseq\>1N\_\{\\rm seq\}\>1sequences of consecutive entries from \([6](https://arxiv.org/html/2608.00400#S2.E6)\),

\(42\)\(𝐮ik\(mk\),𝐮ik\+1\(mk\),⋯,𝐮ik\+nr\(mk\)\),k=1,…,Nseq\.\\left\(\{\\mathbf\{u\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\},\{\\mathbf\{u\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+1\},\\cdots,\{\\mathbf\{u\}\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+n\_\{r\}\}\\right\),\\qquad k=1,\\dots,N\_\{\\rm seq\}\.Let

𝐆~Δ​tℓ=𝐆~Δ​t∘𝐆~Δ​t∘⋯∘𝐆~Δ​t⏟ℓ​times,ℓ≥1,\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}^\{\\ell\}=\\underbrace\{\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}\\circ\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}\\circ\\cdots\\circ\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}\}\_\{\\ell\\text\{ times\}\},\\qquad\\ell\\geq 1,beℓ\\elltimes recurrent operations of𝐆~Δ​t\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}\. We then seek to minimize the following loss

\(43\)ℒnr​\(Θ\)=1Nseq​∑k=1Nseq1nr​∑ℓ=1nr‖𝐆~Δ​tℓ​\(𝐮ik\(mk\);Θ\)−𝐮ik\+ℓ\(mk\)‖22\.\\mathcal\{L\}\_\{n\_\{r\}\}\(\\Theta\)=\\frac\{1\}\{N\_\{\\rm seq\}\}\\sum\_\{k=1\}^\{N\_\{\\rm seq\}\}\\frac\{1\}\{n\_\{r\}\}\\sum\_\{\\ell=1\}^\{n\_\{r\}\}\\left\\\|\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}^\{\\ell\}\\left\(\\mathbf\{u\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\};\\Theta\\right\)\-\\mathbf\{u\}^\{\(m\_\{k\}\)\}\_\{i\_\{k\}\+\\ell\}\\right\\\|\_\{2\}^\{2\}\.Upon satisfactory training, the parametersΘ\\Thetaare fixed, and we obtain a nodal FML predictive model\. For any given solution snapshot𝐮0\\mathbf\{u\}\_\{0\}as the initial condition, the FML model prediction follows

\(44\)\{𝐮~0=𝐮0,𝐮~k\+1=𝐆~Δ​t​\(𝐮~k\),k=0,1,…\\left\\\{\\begin\{aligned\} &\\widetilde\{\\mathbf\{u\}\}\_\{0\}=\\mathbf\{u\}\_\{0\},\\\\ &\\widetilde\{\\mathbf\{u\}\}\_\{k\+1\}=\\widetilde\{\\mathbf\{G\}\}\_\{\\Delta t\}\(\\widetilde\{\\mathbf\{u\}\}\_\{k\}\),\\quad\\qquad k=0,1,\\dots\\end\{aligned\}\\right\.
Compared with the methods that require evaluating nonlocal operators during training, the nodal\-space FML formulation uses only solution snapshots\. This feature is particularly attractive for unknown nonlocal PDEs, where the operator form, the kernel, or the fractional order are often not available\.

## 5Numerical Examples

This section evaluates the nodal\- and modal\-space FML methods on three nonlocal PDEs\. The governing equations are used only to generate synthetic solution data, either analytically or with high\-resolution numerical solvers, for training and validation of the FML models\. In all error plots, the errors are computed from the discrete solution vectors on the corresponding computational grid\. Unless otherwise specified, the reported relative error is defined as

‖𝐮pred−𝐮ref‖2‖𝐮ref‖2,\\frac\{\\\|\\mathbf\{u\}\_\{\\rm pred\}\-\\mathbf\{u\}\_\{\\rm ref\}\\\|\_\{2\}\}\{\\\|\\mathbf\{u\}\_\{\\rm ref\}\\\|\_\{2\}\},where𝐮pred\\mathbf\{u\}\_\{\\rm pred\}is the FML model prediction and𝐮ref\\mathbf\{u\}\_\{\\rm ref\}the reference solution, over the same grid\.

### 5\.1Anomalous Diffusion Equation

We first consider the 1D anomalous diffusion equation\[delia2020numerical\]:

\(45\)\{∂u​\(x,t\)∂t=−\(−Δ\)α2​u​\(x,t\),x∈\(0,2​π\),t\>0,u​\(x,0\)=u0​\(x\),x∈\[0,2​π\]\.\\left\\\{\\begin\{aligned\} &\\frac\{\\partial u\(x,t\)\}\{\\partial t\}=\-\(\-\\Delta\)^\{\\frac\{\\alpha\}\{2\}\}u\(x,t\),\\quad x\\in\(0,2\\pi\),\\ \\ t\>0,\\\\ &u\(x,0\)=u\_\{0\}\(x\),\\quad x\\in\[0,2\\pi\]\.\\end\{aligned\}\\right\.with periodic boundary conditions\. We setα=1\.5\\alpha=1\.5in this experiment\.

The training data are the numerical solution of \([45](https://arxiv.org/html/2608.00400#S5.E45)\), obtained using the Fourier pseudo\-spectral method for spatial discretization and exact time integration\[Duo2016\], withNx=51N\_\{x\}=51spatial points\. We generateNtrain=105N\_\{\\rm train\}=10^\{5\}training sequences from random initial conditions of the form

\(46\)u​\(x,0\)=a0\+∑n=1Nc\(an​cos⁡\(n​x\)\+bn​sin⁡\(n​x\)\),x∈\[0,2​π\],u\(x,0\)=a\_\{0\}\+\\sum\_\{n=1\}^\{N\_\{c\}\}\\big\(a\_\{n\}\\cos\(nx\)\+b\_\{n\}\\sin\(nx\)\\big\),\\quad x\\in\[0,2\\pi\],We sampleNc∼U​\{0,1,…,7\}N\_\{c\}\\sim U\\\{0,1,\\ldots,7\\\},a0∼U​\[−2,2\]a\_\{0\}\\sim U\[\-2,2\], andan,bn∼U​\[−1/n,1/n\]a\_\{n\},b\_\{n\}\\sim U\[\-1/n,1/n\]for1≤n≤Nc1\\leq n\\leq N\_\{c\}, whereUUdenotes the uniform distribution\. For each initial condition, we evolve \([45](https://arxiv.org/html/2608.00400#S5.E45)\) for 30 steps withΔ​t=0\.05\\Delta t=0\.05\(up toT=1\.5T=1\.5\)\. We then select six consecutive snapshots, corresponding to five recurrent prediction steps, to form

\{𝐮\(i\)​\(ti\),𝐮\(i\)​\(ti\+Δ​t\),…,𝐮\(i\)​\(ti\+5​Δ​t\)\}i=1Ntrain,\\left\\\{\\mathbf\{u\}^\{\(i\)\}\(t\_\{i\}\),\\mathbf\{u\}^\{\(i\)\}\(t\_\{i\}\+\\Delta t\),\\ldots,\\mathbf\{u\}^\{\(i\)\}\(t\_\{i\}\+5\\Delta t\)\\right\\\}\_\{i=1\}^\{N\_\{\\rm train\}\},whereti∼U​\{0,Δ​t,…,25​Δ​t\}t\_\{i\}\\sim U\\\{0,\\Delta t,\\ldots,25\\Delta t\\\}\. Thus,NtrainN\_\{\\rm train\}denotes the number of training sequences rather than the number of snapshots\. For nodal FML, we use a disassembly–assembly network withJ=3J=3,nd=na=1n\_\{d\}=n\_\{a\}=1, andnw=51n\_\{w\}=51\. For modal FML, we use a one\-block ResNet with six hidden layers of 50 neurons each\. Both models use a five\-step recurrent loss and a cyclic learning rate that ranges from10−710^\{\-7\}to10−310^\{\-3\}with exponential decay factor0\.99999970\.9999997\. We train each network for 10,000 epochs with a batch size of 50\.

For validation, we generate 100 test samples independently from the same distribution and evaluate 10 independently trained models over 500 time steps \(up toT=25T=25\)\. Figures[2](https://arxiv.org/html/2608.00400#S5.F2)and[3](https://arxiv.org/html/2608.00400#S5.F3)show nodal and modal predictions, respectively, for one test sample\.

![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_1d/traj_RNN5_id2.png)Figure 2:Reference and nodal FML predictions for the solution of \([45](https://arxiv.org/html/2608.00400#S5.E45)\)\.![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_1d_modal/id0_long.png)Figure 3:Reference and modal FML predictions for the solution of \([45](https://arxiv.org/html/2608.00400#S5.E45)\)\.Both predictions follow the reference trajectory up toT=25T=25\. Figure[4](https://arxiv.org/html/2608.00400#S5.F4)shows the meanℓ2\\ell\_\{2\}errors over the 100 test samples: approximately10−410^\{\-4\}for nodal FML and10−510^\{\-5\}for modal FML\. This demonstrates that the learned FML model remains accurate under recursive application for several hundred prediction steps, despite being trained using only short trajectory segments\.

\(a\)![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_1d/error_RNN5_c10_0.png)\(b\)![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_1d_modal/l2_error.png)

Figure 4:Averageℓ2\\ell\_\{2\}error over 100 test samples from nodal \(a\) and modal \(b\) FML\.
### 5\.2Nonlocal Wave Equation

We consider the 1D nonlocal wave equation of the form:

\(47\)\{ut​t​\(x,t\)=−D​\(−Δ\)α2​u​\(x,t\),for​x∈Ω,t\>0,u​\(x,0\)=u0​\(x\),ut​\(x,0\)=v0​\(x\),for​x∈Ω¯\\left\\\{\\begin\{aligned\} &u\_\{tt\}\(x,t\)=\-D\(\-\\Delta\)^\{\\frac\{\\alpha\}\{2\}\}u\(x,t\),\\quad\\ \\text\{for\}\\ \\,x\\in\\Omega,\\ \\ t\>0,\\\\ &u\(x,0\)=u\_\{0\}\(x\),\\quad u\_\{t\}\(x,0\)=v\_\{0\}\(x\),\\quad\\ \\text\{for\}\\ \\,x\\in\\bar\{\\Omega\}\\end\{aligned\}\\right\.with periodic boundary conditions\. In this experiment,Ω=\(0,2​π\)\\Omega=\(0,2\\pi\),α=0\.5\\alpha=0\.5, andD=1D=1\. The numerical solutions are generated with the Fourier pseudo\-spectral method in space withNx=51N\_\{x\}=51and exact time integration\.

For the training set, we generateNtrain=106N\_\{\\rm train\}=10^\{6\}sequences from random initial displacementu0u\_\{0\}and velocityv0v\_\{0\}\. Their nodal vectors are sampled from periodic Gaussian vectors as follows\. Set the grid points to beθj=2​π​jNx−1\\theta\_\{j\}=\\frac\{2\\pi j\}\{N\_\{x\}\-1\}forj=0,…,Nx−1j=0,\\ldots,N\_\{x\}\-1, where the periodic grid includes the duplicated endpoint\. We generate𝐬=\(s0,…,sNx−1\)\\mathbf\{s\}=\(s\_\{0\},\\dots,s\_\{N\_\{x\}\-1\}\)such thatCov⁡\(sj,sk\)=exp⁡\(−\(θj−θk\)2/b2\)\\operatorname\{Cov\}\(s\_\{j\},s\_\{k\}\)=\\exp\\big\(\-\\left\(\\theta\_\{j\}\-\\theta\_\{k\}\\right\)^\{2\}/b^\{2\}\\big\)withbbthe correlation length by

sj\\displaystyle s\_\{j\}=∑n=−KGKGRn​e−i​n​θj,j=0,…,Nx−1,\\displaystyle=\\sum\_\{n=\-K\_\{G\}\}^\{K\_\{G\}\}R\_\{n\}e^\{\-\\mathrm\{i\}n\\theta\_\{j\}\},\\quad j=0,\\ldots,N\_\{x\}\-1,Rn\\displaystyle R\_\{n\}∼𝒩​\(0,Cn\),Cn=∫02​πe−θ2b2​cos⁡\(n​θ\)​𝑑θ\.\\displaystyle\\sim\\mathcal\{N\}\(0,C\_\{n\}\),\\qquad C\_\{n\}=\\int\_\{0\}^\{2\\pi\}e^\{\-\\frac\{\\theta^\{2\}\}\{b^\{2\}\}\}\\cos\(n\\theta\)d\\theta\.The sampled displacement and velocity vectors are𝐮0=𝐬\+mu\\mathbf\{u\}\_\{0\}=\\mathbf\{s\}\+m\_\{u\}and𝐯0=𝐬\+mv\\mathbf\{v\}\_\{0\}=\\mathbf\{s\}\+m\_\{v\}, wheremu∼U​\[−1,1\]m\_\{u\}\\sim U\[\-1,1\],mv∼U​\[−0\.1,0\.1\]m\_\{v\}\\sim U\[\-0\.1,0\.1\],b∼U​\[0\.1,1\.1\]b\\sim U\[0\.1,1\.1\], andKG=10K\_\{G\}=10\. We represent the discrete wave state by

𝐳​\(t\)=\(𝐮​\(t\)𝐮˙​\(t\)\)∈ℝ2​Nx\.\\mathbf\{z\}\(t\)=\\begin\{pmatrix\}\\mathbf\{u\}\(t\)\\\\ \\dot\{\\mathbf\{u\}\}\(t\)\\end\{pmatrix\}\\in\\mathbb\{R\}^\{2N\_\{x\}\}\.For each initial condition, we solve \([47](https://arxiv.org/html/2608.00400#S5.E47)\) for 20 steps withΔ​t=0\.05\\Delta t=0\.05and form one\-step training pairs for nodal FML:

\{𝐳\(i\)​\(ti\),𝐳\(i\)​\(ti\+Δ​t\)\}i=1Ntrain,\\Big\\\{\\mathbf\{z\}^\{\(i\)\}\(t\_\{i\}\),\\mathbf\{z\}^\{\(i\)\}\(t\_\{i\}\+\\Delta t\)\\Big\\\}\_\{i=1\}^\{N\_\{\\rm train\}\},whereti∼U​\{0,Δ​t,…,19​Δ​t\}t\_\{i\}\\sim U\\\{0,\\Delta t,\\ldots,19\\Delta t\\\}\.

For nodal FML, we use a disassembly–assembly network with input dimensionN=2​Nx=102N=2N\_\{x\}=102,J=3J=3,nd=na=1n\_\{d\}=n\_\{a\}=1, andnw=51n\_\{w\}=51\. For modal FML, we use the sameNtrain=106N\_\{\\rm train\}=10^\{6\}trajectories and select six consecutive snapshots, corresponding to five recurrent prediction steps:

\{𝐳\(i\)​\(ti\),𝐳\(i\)​\(ti\+Δ​t\),…,𝐳\(i\)​\(ti\+5​Δ​t\)\}i=1Ntrain,\\left\\\{\\mathbf\{z\}^\{\(i\)\}\(t\_\{i\}\),\\mathbf\{z\}^\{\(i\)\}\(t\_\{i\}\+\\Delta t\),\\ldots,\\mathbf\{z\}^\{\(i\)\}\(t\_\{i\}\+5\\Delta t\)\\right\\\}\_\{i=1\}^\{N\_\{\\rm train\}\},whereti∼U​\{0,Δ​t,…,15​Δ​t\}t\_\{i\}\\sim U\\\{0,\\Delta t,\\ldots,15\\Delta t\\\}\. The modal model is a two\-block ResNet with three hidden layers per block and 40 neurons per layer, trained with a five\-step recurrent loss\.

For validation, we independently generate 100 test samples from the training distribution and evaluate 10 independently trained models for 500 time steps \(up toT=25T=25\)\.

![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/wave_1d/RNN5_id2_long.png)Figure 5:Reference and nodal FML predictions for the solution of \([47](https://arxiv.org/html/2608.00400#S5.E47)\)\.![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/wave_1d_modal/id67_long.png)Figure 6:Reference and modal FML predictions for the solution of \([47](https://arxiv.org/html/2608.00400#S5.E47)\)\.Figures[5](https://arxiv.org/html/2608.00400#S5.F5)and[6](https://arxiv.org/html/2608.00400#S5.F6)show that both methods reproduce the oscillatory reference trajectory over the 500\-step prediction interval\. Figure[7](https://arxiv.org/html/2608.00400#S5.F7)reports the mean relative discreteℓ2\\ell\_\{2\}errors over the 100 test samples\.

\(a\)![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/wave_1d/error_RNN5_J5_long.png)\(b\)![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/wave_1d_modal/relative_error_M.png)

Figure 7:Average relativeℓ2\\ell\_\{2\}error over 100 test samples from nodal \(a\) and modal \(b\) FML\.
### 5\.3Mixed Nonlocal Diffusion Equation

Next, we consider the mixed nonlocal diffusion equation with periodic boundary conditions:

\(48\)\{∂tu​\(x,y,t\)=−c1​\(−Δ\)α2​u−c2​\(−Δ\)β2​u,for\(x,y\)∈\(0,2​π\)2,t\>0,u​\(x,y,0\)=u0​\(x,y\),for\(x,y\)∈\[0,2​π\]2,\\left\\\{\\begin\{aligned\} &\\partial\_\{t\}u\(x,y,t\)=\-c\_\{1\}\(\-\\Delta\)^\{\\frac\{\\alpha\}\{2\}\}u\-c\_\{2\}\(\-\\Delta\)^\{\\frac\{\\beta\}\{2\}\}u,\\quad\\mbox\{for\}\\ \\ \(x,y\)\\in\(0,2\\pi\)^\{2\},\\ \\ t\>0,\\\\ &u\(x,y,0\)=u\_\{0\}\(x,y\),\\quad\\mbox\{for\}\\ \\ \(x,y\)\\in\[0,2\\pi\]^\{2\},\\end\{aligned\}\\right\.where we chooseα=1\.5\\alpha=1\.5,β=0\.5\\beta=0\.5, andc1=c2=0\.05c\_\{1\}=c\_\{2\}=0\.05\. The numerical solutions are generated with the Fourier pseudospectral method in space and exact time integration withNx=Ny=80N\_\{x\}=N\_\{y\}=80, using a periodic grid with duplicated endpoints\.

We generateNtrain=106N\_\{\\rm train\}=10^\{6\}training sequences from random initial conditions of the form

\(49\)u0​\(x,y\)=∑n=0Ncx∑m=0Ncy\(A​\(n,m\)​cos⁡\(n​x\+m​y\)\+B​\(n,m\)​sin⁡\(n​x\+m​y\)\),u\_\{0\}\(x,y\)=\\sum\_\{n=0\}^\{N\_\{c\}^\{x\}\}\\sum\_\{m=0\}^\{N\_\{c\}^\{y\}\}\\big\(A\(n,m\)\\cos\(nx\+my\)\+B\(n,m\)\\sin\(nx\+my\)\\big\),whereNcx,Ncy∼U​\{0,1,2,3,4\}N\_\{c\}^\{x\},N\_\{c\}^\{y\}\\sim U\\\{0,1,2,3,4\\\}andA​\(n,m\),B​\(n,m\)∼U​\[−2−\(n\+m\),2−\(n\+m\)\]A\(n,m\),B\(n,m\)\\sim U\[\-2^\{\-\(n\+m\)\},2^\{\-\(n\+m\)\}\]\. For each initial condition, we select six consecutive snapshots, corresponding to five recurrent prediction steps:

\{𝐮\(i\)​\(ti\),𝐮\(i\)​\(ti\+Δ​t\),…,𝐮\(i\)​\(ti\+5​Δ​t\)\}i=1Ntrain,\\left\\\{\\mathbf\{u\}^\{\(i\)\}\(t\_\{i\}\),\\mathbf\{u\}^\{\(i\)\}\(t\_\{i\}\+\\Delta t\),\\ldots,\\mathbf\{u\}^\{\(i\)\}\(t\_\{i\}\+5\\Delta t\)\\right\\\}\_\{i=1\}^\{N\_\{\\rm train\}\},whereti∼U​\{0,Δ​t,…,15​Δ​t\}t\_\{i\}\\sim U\\\{0,\\Delta t,\\ldots,15\\Delta t\\\}andΔ​t=0\.05\\Delta t=0\.05\.

For nodal FML, we use a disassembly–assembly network withJ=5J=5,nd=na=2n\_\{d\}=n\_\{a\}=2, andnw=441n\_\{w\}=441\. For modal FML, we use a four\-block ResNet with three hidden layers per block and 50 neurons per layer\. Both models use a five\-step recurrent loss\.

For validation, we independently generate 100 test samples from the training distribution and evaluate five independently trained models for 1,000 time steps \(up toT=50T=50\)\. Figures[8](https://arxiv.org/html/2608.00400#S5.F8)and[9](https://arxiv.org/html/2608.00400#S5.F9)compare the reference and predicted solutions for one test sample\. Both methods reproduce the spatial decay pattern throughout the prediction interval\. Figure[10](https://arxiv.org/html/2608.00400#S5.F10)reports the mean relative discreteℓ2\\ell\_\{2\}errors over the 100 test samples\.

![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_2d/int_short_1.png)Figure 8:Reference and nodal FML predictions for the solution of \([48](https://arxiv.org/html/2608.00400#S5.E48)\)\.![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_2d_modal/id10_subs.png)Figure 9:Reference and modal FML predictions for the solution of \([48](https://arxiv.org/html/2608.00400#S5.E48)\)\.![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_2d/rel_error_RNN5.png)

![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_2d_modal/relative_l2_error.png)

Figure 10:Average relative discreteℓ2\\ell\_\{2\}error over 100 test samples from nodal \(left\) and modal \(right\) FML\.
### 5\.4Parameter recovery

We further assess the learned modal FML dynamics by recovering the fractional orders from the predicted trajectories\. To accomplish this, we recall the fractional Laplacian pseudo\-differential operator \([8](https://arxiv.org/html/2608.00400#S3.E8)\),

\(−Δ\)α/2​u​\(𝐱,t\)=ℱ−1​\[\|ξ\|α​ℱ​\[u\]\]\.\(\-\\Delta\)^\{\\alpha/2\}u\(\\mathbf\{x\},t\)=\\mathcal\{F\}^\{\-1\}\[\|\\xi\|^\{\\alpha\}\\mathcal\{F\}\[u\]\]\.
For the 1D diffusion problem, the modal solution satisfies

\(50\)u^​\(ξ,t\)=u^​\(ξ,0\)​e−\|ξ\|α​t,\\hat\{u\}\(\\xi,t\)=\\hat\{u\}\(\\xi,0\)e^\{\-\|\\xi\|^\{\\alpha\}t\},which provides a modewise estimateα​\(t,k\)\\alpha\(t,k\)\. On the periodic domain, the frequency variable is the integer Fourier modekk, soξ=k\\xi=k\. Figure[11](https://arxiv.org/html/2608.00400#S5.F11)shows the estimates from six nonzero modes of one predicted trajectory overt∈\[0,1\]t\\in\[0,1\]\. The estimates agree closely with the true valueα=1\.5\\alpha=1\.5at early times\. At later times, the modal amplitudes approach zero, making the estimate increasingly sensitive to numerical and prediction errors\.

![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_1d_modal/parafit_id30.png)Figure 11:Fractional orderα\\alpharecovered from a predicted trajectory of the 1D diffusion equation, plotted against the integer Fourier modekk\.For the 1D wave equation, the Fourier coefficientu^​\(ξ,t\)\\hat\{u\}\(\\xi,t\)satisfies the ODE

u^t​t​\(ξ,t\)\+c​\|ξ\|α​u^​\(ξ,t\)=0\.\\hat\{u\}\_\{tt\}\(\\xi,t\)\+c\|\\xi\|^\{\\alpha\}\\hat\{u\}\(\\xi,t\)=0\.We approximateu^t​t\\hat\{u\}\_\{tt\}from the predicted modal snapshots using the centered second\-order finite difference associated with the Leap\-Frog scheme and use the result to estimateα\\alpha\. Figure[12](https://arxiv.org/html/2608.00400#S5.F12)shows that the estimate remains close to the true value over a longer interval than in the diffusion example\. Unlike the diffusive modes, the wave modes do not undergo monotone amplitude decay, so their signal remains available for parameter recovery at later times\.

![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/wave_1d_modal/parafit_id20.png)Figure 12:Fractional orderα\\alpharecovered from a predicted trajectory of the 1D nonlocal wave equation\.For the 2D diffusion equation,ξ=\(k,l\)\\xi=\(k,l\)\. We recover\(α,β\)\(\\alpha,\\beta\)from the predicted modal coefficients using the MATLABfitfunction\. For eachk∈\{1,2,3,4\}k\\in\\\{1,2,3,4\\\}, we fit one parameter pair using the four nonzero modes withl∈\{1,2,3,4\}l\\in\\\{1,2,3,4\\\}\. Figure[13](https://arxiv.org/html/2608.00400#S5.F13)shows the four resulting estimates\. Becausec1=c2=0\.05c\_\{1\}=c\_\{2\}=0\.05, the modal amplitudes decay more slowly than in the 1D diffusion example, and useful estimates can be obtained over a longer time interval\.

![Refer to caption](https://arxiv.org/html/2608.00400v1/Figures/diffusion_2d_modal/parafit_id8_k4l4.png)Figure 13:Fractional ordersα\\alphaandβ\\betarecovered from a predicted trajectory of the 2D nonlocal diffusion equation\.

## 6Conclusion

We have extended flow map learning \(FML\) to unknown nonlocal PDE systems in both modal and nodal representations\. The proposed formulations learn the finite\-time evolution operator directly from solution data without explicit approximation of the underlying nonlocal operators\. Numerical experiments demonstrate accurate long\-time prediction for several representative fractional diffusion and wave equations\. These results suggest that FML provides a promising computational framework for data\-driven modeling of unknown nonlocal dynamics\. Future work will focus on more general classes of nonlocal operators and applications involving experimentally observed data\.

## References