Numerical methods for nonlocal and fractional models

Nonlocal diffusion models, including fractional models

• Nonlocal models represent distance-based interactions through integral equations and include fractional derivative models, while also enabling models beyond the fractional class.These models can represent physics differently from PDE models for the same applications.Source: Nonlocal diffusion models, including… · p. 3

1. General models for nonlocal diffusion

• Unlike a Green’s-function reformulation of the classical Poisson problem, the nonlocal model represents different physics from the local PDE model.The local reference problem describes steady-state diffusion with standard boundary conditions, whereas the nonlocal formulation uses interaction over a surrounding volume.Source: §1 · p. 4
• The model uses a nonlocal operator Lδ on Ω with Dirichlet, Neumann, Robin, or mixed volume constraints imposed through V and Nδ on ΩIδ.The operator is supported by interactions within distance δ and is a nonlocal analogue of ∇·(D∇u).Source: §1 · pp. 5–6
• The general nonlocal diffusion problem replaces boundary constraints on ∂Ω with volume constraints on the interaction domain ΩIδ.The interaction domain contains points outside Ω within horizon δ of some point in Ω, and δ may range from much smaller than diam Ω to infinity.Source: §1 · p. 5

1.1. A brief review of a nonlocal vector calculus

• A nonlocal vector calculus develops integral analogues of classical divergence, gradient, curl, and Green’s identities for nonlocal models.The framework was developed in Gunzburger and Lehoucq (2010), Du et al. (2012a), Du, Gunzburger, Lehoucq and Zhou (2013a), and subsequent work, with Du et al. (2013a) providing the main development reviewed here.Source: §1.1 · pp. 7, 9
• The horizon δ and antisymmetric kernel function αδ(x, y) are modeling choices that can produce operators with substantially different properties.The constitutive tensor Θδ(x, y) is another modeling choice governing media properties and solution regularity.Source: §1.1 · pp. 8, 10
• Composing nonlocal divergence and gradient operators defines Lδ as a nonlocal Poisson operator and naturally generalizes constant-coefficient fractional operators to variable coefficients.When Θδ is the identity tensor, Lδ is a nonlocal analogue of the Laplacian.Source: §1.1 · pp. 9–10
• Additional nonlocal operators yield vector identities, a nonlocal vector Laplacian, and Helmholtz-type decompositions, although some two-point operators do not recover classical counterparts as δ → 0.Related constructions support nonlocal Stokes equations and recover the classical Helmholtz decomposition through alternative one-point operators.Source: §1.1 · pp. 10–11

2. Weak formulations of nonlocal models

• The section introduces a weak formulation corresponding to the strong nonlocal problem and concludes with kernel examples and their associated energy spaces.The formulation is selected for inhomogeneous data and is later specialized to homogeneous data when convenient.Source: §2 · pp. 11–12, §2.1.2 · p. 15

2.1.1. Function spaces and norms

• The energy space V(Ω∪ΩIδ) is a kernel-dependent subspace of L2(Ω∪ΩIδ), while constrained and trace-like spaces support nonlocal volume constraints.Under suitable kernel conditions, a nonlocal Poincaré inequality makes the constrained seminorm an equivalent norm with constants independent of δ as δ →0.Source: §2.1.1 · pp. 12–13

2.1.2. Weak formulations and well-posedness

• The weak problem seeks u in V(Ω∪ΩIδ) satisfying Aδ(u,v)=⟨f,v⟩ for all constrained v while imposing u=g on ΩIδ.It follows from the nonlocal Green identity, and suitable kernels make Aδ symmetric, continuous, and coercive, yielding well-posedness by Lax–Milgram and an energy-norm estimate independent of δ.Source: §2.1.2 · pp. 13–14
• A nonlocal Neumann weak formulation uses the quotient space V(Ω∪ΩIδ) \ R for uniqueness and requires a compatibility condition on f and g for existence.The section notes alternative nonlocal Neumann formulations and fractional-Laplacian variants, but focuses primarily on the volume-constrained formulation.Source: §2.1.2 · p. 14

2.1.3. Relation to an energy minimization principle

• The weak formulation is the Euler–Lagrange equation of an energy minimization principle under suitable conditions on the kernel and data.This connects solving the nonlocal variational problem with minimizing its associated energy functional.Source: §2.1.3 · p. 15

2.2. Kernel choices and the corresponding energy spaces

• The kernel decomposition separates operator and constitutive roles: φδ determines first-order structure, while θδ represents constitutive behavior and together define γδ, Lδ, and Nδ.The construction assumes αδ is directed along y−x, with symmetric nonnegative constitutive factors; radial choices model homogeneous media, whereas inhomogeneity requires non-radial constitutive structure.Source: §2.2 · pp. 15–16

2.2.1. A list of kernel functions in common use

• Common kernel classes include translation-invariant, critical, peridynamic, and fractional functions, each associated with characteristic energy spaces and singularity behavior.The fractional class has the fractional-Laplacian singularity and yields 2s derivatives of smoothing under f∈H−s, whereas integrable kernels may preserve data discontinuities.Source: §2.2.1 · pp. 16–18, Table 2.1 · p. 17
• Kernel scaling constants are modeling choices that can enforce local δ→0 or whole-space δ→∞ limits, or instead be selected from the modeled physics.The listed classes also distinguish smooth, nonsmooth, singular-integrable, critical, peridynamic, and fractional behavior; square-integrable kernels may instead have L2 as their energy space.Source: §2.2.1 · pp. 18–19

3. Fractional diffusion models

• Fractional diffusion uses the fractional Laplacian −(−∆)^s, s ∈(0, 1), as a nonlocal model associated with long-range Lévy-flight interactions and superdiffusion.The operator approaches the identity as s →0+ and the PDE Laplacian −∆ as s →1−.Source: §3 · pp. 19–20

3.1. Integral fractional Laplacian models for diffusion on bounded domains

• On bounded domains, the integral fractional Poisson problem imposes volume constraints outside Ω, is variationally well posed, and converges strongly to the local diffusion solution as s →1−.The bilinear form is H^s_c(Ω)-coercive and continuous, yielding well-posedness through the Lax–Milgram theorem.Source: §3.1 · pp. 20, 22

3.1.1. Truncated interaction domains

• Truncating the infinite interaction domain produces a finite-horizon model whose solution approximates the full-space solution, while potentially better representing practical large but finite horizons.The truncation error is controlled by an estimate depending on norms of f and g.Source: §3.1.1 · p. 23

3.2. Spectral fractional Laplacian models for diffusion on bounded domains

• The spectral fractional Laplacian instead fractionalizes the bounded-domain PDE Laplacian through its homogeneous Dirichlet eigenpairs, recovering the identity at s →0 and the integer-order Laplacian at s →1.Unlike the integral definition, the spectral and integral operators differ on bounded domains, although they coincide on ℝ^d.Source: §3.2 · pp. 23–25

3.3.1. Other integral representations of the fractional Laplacian

• Alternative formulations support computation through Dunford–Taylor and heat-kernel integral representations, while extension methods replace nonlocality with a local weighted problem in one additional dimension.The extension approach introduces a singular weight and an additional spatial dimension.Source: §3.3.1 · p. 25, §3.3.2 · pp. 25–26

3.3.3. Regularity of solutions

• Fractional Poisson solutions lack the classical integer-order smoothing property and exhibit formulation-dependent regularity, with boundary behavior characterized separately for integral and spectral operators.For the integral problem, regularity is established under smooth-boundary and Sobolev-data assumptions; the spectral problem satisfies a corresponding r+2s regularity result.Source: §3.3.3 · p. 26

3.3.4. Analytic solutions in one and two dimensions

• Closed-form solutions on the unit ball provide verification cases for numerical methods and reveal boundary terms behaving like δ(x)^s, where δ(x) is the distance to the boundary.Figure 3.1 displays example analytic solutions for one- and two-dimensional settings.Source: §3.3.4 · p. 27, Figure 3.1 · p. 28

3.4. Inhomogeneous constitutive functions in fractional models

• Inhomogeneous fractional media are modeled by replacing the standard kernel with one containing a constitutive tensor Θ_s(x, y), which must be nonradial to represent spatially varying properties.Radial Θ_s yields a homogeneous medium, while scalar radial Θ_s corresponds to a constant-coefficient fractional analogue.Source: §3.4 · pp. 28–29

4. Introductory remarks

• Nonlocal discretizations involve horizon δ and resolution parameters h or N, with continuous and discrete limits considered as these parameters vary.Relevant limits include δ →0, δ →∞, h or N refinement, and combinations of these limits.Source: §4 · pp. 29–30
• For some algorithms, taking δ and h limits in different orders or simultaneously can produce different, potentially incorrect, limiting solutions.Consistency therefore requires examining how error-estimate constants depend on δ, h, or N.Source: §4 · p. 30

4.1. Asymptotically compatible schemes

• Asymptotically compatible schemes converge to nonlocal models at fixed δ and to corresponding local discrete schemes as δ vanishes, independently of limit order.This property applies with fixed numerical resolution and with increasing resolution for continuum models; Figure 4.1 illustrates the relationships among uδ, u0, uδ,h, and u0,h.Source: §4.1 · p. 30, Figure 4.1 · p. 31
• When δ is proportional to h, some peridynamic discretizations converge to incorrect limits, including piecewise constant conforming finite elements and simple Riemann-sum quadrature.AC alternatives include conforming Galerkin, discontinuous Galerkin, collocation, and quadrature-based strong-form approximations (Chen and Gunzburger 2011, Tian and Du 2014, Xu, Gunzburger, Burkardt and Du 2016b).Source: §4.1 · pp. 30–31

5. Finite element methods for nonlocal models

• Conforming finite element discretizations construct piecewise-polynomial spaces over Ω∪ΩIδ, impose volume constraints, and yield well-posed Galerkin problems by continuity and coercivity.Singular interaction kernels require specialized quadrature for stiffness-matrix evaluation; for some integrable radial kernels, DG functions are already conforming because nonlocality couples neighboring elements.Source: §5 · pp. 31–33

5.1. Asymptotically compatible conforming finite element methods for nonlocal diffusion

• Asymptotically compatible conforming Galerkin methods with continuous piecewise linear functions recover the correct local limit as both δ and h decrease, without requiring h<δ.For piecewise constants, the condition h=o(δ) is required, whereas fixed bandwidth can produce an incorrect local limit.Source: §5.1 · pp. 34–35

5.2. Nonconforming and DG FEMs for nonlocal models with sufficiently singular kernels

• For sufficiently singular kernels, Tian and Du (2015) use kernel cutoffs to make discontinuous Galerkin spaces conforming, with convergence to the original nonlocal solution as the cutoff is removed and h decreases.With piecewise constants, convergence requires h=o(1/n).Source: §5.2 · pp. 35–37

5.3. Adaptive mesh refinement for nonlocal models

• Adaptive finite element methods provide reliable and efficient residual estimators, while anisotropic refinement around solution jumps achieves h^2 convergence with degrees of freedom comparable to a regular local grid.The strategy uses elongated elements of thickness O(h^4) across discontinuities and length O(h) along them, alongside abrupt grid transitions and localized discontinuous bases.Source: §5.3 · pp. 37–39

5.4. Approximations of fractional models as limit of nonlocal models with a finite range of interactions

• Conforming Galerkin approximations of finite-horizon fractional-type nonlocal models converge to fractional equations as δ→∞, while the same framework also supports the local limit as δ→0.The fractional-limit convergence holds as h→0 and δ→∞ without requiring h to depend on δ.Source: §5.4 · pp. 39–41

6. Finite element methods for the integral fractional Laplacian

• The finite element space is X_h for s < 1/2 and the interior-node subspace vanishing on ∂Ω for s ≥ 1/2, with dimension determined accordingly.The discrete solution satisfies A_s(u_h, v) = ⟨f, v⟩ for every v ∈ V_h.Source: §6 · p. 41
• For u ∈ H2(Ω), globally quasi-uniform meshes yield the expected convergence rate h2−s = O(N_h^(s−2)/d), while in two dimensions refined meshes achieve at most O(N_h−1/2+ε).Boundary refinement restores the optimal rate in d = 1 but not in higher dimensions because element shape-regularity becomes limiting.Source: §6 · p. 42
• Because solutions generally have boundary singularities, general geometries require adaptively refined meshes, supported by a posteriori error estimators for the integral fractional Laplacian.Appropriately graded meshes may suffice for simple geometries.Source: §6 · p. 42

6.1. Quadrature rules

• The fractional Poisson discretization produces a dense system with singular-integral matrix entries, so element-pair and element–external-edge contributions must be decomposed for quadrature.Non-disjoint element pairs require boundary-element-style treatments because their singular integrals are not directly amenable to numerical quadrature.Source: §6.1 · p. 42
• Non-uniform-order Gauss-type rules and transformed weighted quadrature handle fractional-Laplacian singularities, with O(log N_h^2d) points per element pair making quadrature consistency error subordinate to discretization error.The transformed approach maps rules on [0, 1]2d to element pairs and incorporates singularity through the quadrature weight.Source: §6.1 · p. 43

7. Finite element methods for the spectral fractional Laplacian

• The extension problem is posed on Ω × [0, ∞), with the spectral fractional Poisson solution recovered as the trace u = trace_Ω U.Its weak formulation uses a weighted solution space on the extension domain.Source: §7 · pp. 43–44

7.1. Truncation in the extended direction

• Truncating the extended cylinder at z_trun yields an exponentially small truncation error, after which graded-mesh finite elements approximate the problem.With M ∼ h^-1, γ > 3/(1 − β), and z_trun ∼ |log N|, the resulting method is optimal in regularity but sub-optimal in complexity; sparse grids and hp-finite elements address this limitation, with hp methods achieving exponential convergence even for low domain regularity and incompatible forcing.Source: §7.1 · pp. 44–45

7.2. Spectral method in the extended direction

• A hybrid alternative uses finite elements in the spatial direction and a spectral method in the extension direction, approximating eigenvalues through asymptotic laws, coarse finite elements, and decimation.Only |log h|^p eigenvalue approximations are required for some p > 0, while tensor structure reduces the linear solve to integer-order reaction–diffusion problems handled by conjugate-gradient or multigrid solvers.Source: §7.2 · pp. 45–46

• For periodic one-dimensional diffusion, Fourier modes are eigenfunctions of the nonlocal operator, and their eigenvalues converge to the Laplacian eigenvalues n^2 as δ tends to zero.The method uses projection onto the finite Fourier space X_N of zero-mean periodic functions.Source: §8 · p. 47
• The Fourier spectral-Galerkin scheme is asymptotically compatible, with local-solution convergence at least quadratically in 1/N and uniform error estimates for sufficiently small δ and large N.The estimates do not restrict the relative sizes of δ and N, and their constants are independent of δ and N.Source: §8 · pp. 47–48
• Because nonlocal eigenvalues generally lack analytic expressions and involve highly oscillatory integrals, a hybrid algorithm uses series truncation for small nδ and ODE integration for large nδ.For large nδ, high-order Runge–Kutta methods evaluate the auxiliary function K(nδ).Source: §8 · pp. 48–49
• On the unit sphere, spherical harmonics diagonalize the nonlocal Laplace–Beltrami operator, while modified Clenshaw–Curtis quadrature and asymptotics reduce eigenvalue-computation costs.The cost is O(l^2) per eigenvalue generally and O(l log l) for sufficiently large l, followed by fast spherical harmonic transforms.Source: §8 · p. 49
• Spectral approximations extend to nonlocal Allen–Cahn phase-field equations, retaining spectral accuracy for smooth minimizers and convergence even when minimizers are discontinuous.Such discontinuities can arise when interaction kernels lack sufficiently strong singularities, despite conventional phase-field models having smooth profiles.Source: §8 · pp. 49–50
• Spectral methods can capture solution discontinuities but exhibit spurious Gibbs phenomena, while extensions beyond periodic boundary conditions remain open because suitable eigenvalue estimates are unavailable.The analysis currently provides uniform convergence rates in δ and N primarily for periodic settings.Source: §8 · p. 50

9. Spectral-Galerkin methods for fractional diffusion

• Fractional diffusion solutions decay algebraically on R^d and develop weak boundary singularities on bounded domains, making conventional discretizations difficult and motivating specialized spectral bases.On bounded domains, such bases may achieve higher accuracy than finite difference and finite element methods.Source: §9 · p. 51

9.1. Spectral-Galerkin methods in unbounded domains

• In unbounded domains, avoiding truncation with orthogonal functions addresses the large errors caused by slow power-law decay and avoids unresolved artificial-boundary constructions.The discussion restricts the presentation to one dimension, with higher dimensions treated by tensor products.Source: §9.1 · pp. 51–52

9.1.1. Approximation by Hermite polynomials

• Hermite–Galerkin discretization yields the energy-norm estimate ∥u −uN∥Hs(R) ≤CNs−m|u| b Hm(R) + CN1/6−ℓ/2|f| b Hℓ(R), but not exponential convergence.The method has a diagonal mass matrix and efficiently computable stiffness matrix via the Fourier transform.Source: §9.1.1 · p. 53

9.1.2. Approximation by modified Gegenbauer polynomials

• Modified Gegenbauer bases with nonlinear singular mappings are better suited than classical Hermite or Laguerre bases for algebraically decaying functions and support corresponding approximation estimates.Their fractional-Laplacian action is less compact than for Hermite polynomials but can be approximated by finite hypergeometric-function series.Source: §9.1.2 · pp. 53, 55

9.1.3. Approximation by modified Chebyshev polynomials

• Mapped Chebyshev spectral-Galerkin methods diagonalize the reformulated system, use FFT and sinc quadrature, extend by tensor products, and achieve O(N(log N)^d) complexity.The relevant error estimate is similar to the mapped Gegenbauer estimate with λ = 0.Source: §9.1.3 · pp. 55–56

9.2. Spectral-Galerkin methods in bounded domains

• For bounded domains, Jacobi-weighted Petrov–Galerkin bases encode boundary behavior and yield exponential convergence when the source term f is smooth.The analysis extends to s ∈(0, 1/2] and multi-interval domains, while Jacobi and spherical-harmonic bases provide fast, accurate methods on unit balls in two and three dimensions.Source: §9.2 · pp. 57–59
• Boundary singularities and insufficient regularity limit bounded-domain spectral rates for fractional diffusion–reaction models, and optimal L2(R^d) estimates remain poorly understood.For unbounded problems, existing spectral methods also fail to converge exponentially because solutions generally decay algebraically at infinity.Source: §9.1.3 · pp. 56–57, §9.2 · p. 59
• One-dimensional fractional Poisson problems without lower-order terms dominate spectral-method studies because relation (9.14) makes Jacobi polynomials suitable for solution approximation.Developing rigorous techniques for convection and reaction terms is therefore an important direction.Source: §9.2 · p. 59
• Higher-dimensional fractional diffusion requires spectral methods that address irregular domains and unclear regularity theory, even on convex polygons.These difficulties make the numerical analysis non-trivial and motivate further study in bounded higher-dimensional domains.Source: §9.2 · p. 59
• Regularity in suitable spaces, appropriate basis functions, and optimal convergence rates remain largely open theoretical problems.Source: §9.2 · p. 60

Quadrature-based finite differences

• AC quadrature finite differences depend on uniform truncation errors independent of δ; Tian and Du (2013) established the property, with extensions by Du et al. (2019c).The multidimensional construction uses W(z) = |z|2/|z|1, while the one-dimensional case uses W(z) = |z|.Source: §10 · p. 60
• The Cartesian-grid scheme is second-order accurate and AC, with O(h2) truncation error independent of δ and solution approximation rate O(h2 + δ2).Quadratic exactness is the key mechanism, and the rate follows once numerical stability is established.Source: §10 · p. 61
• Nonnegative stencil coefficients give the Cartesian-grid method a discrete maximum principle and numerical stability, whereas higher-order interpolants can lose this property.Higher-order stability remains open except for specialized kernels, including work by Leng et al. (2019).Source: §10 · pp. 61–62

Extensions

• AC strong-form methods have also been developed for coupled local–nonlocal, time-dependent, infinite-domain, nonlocal-in-time, nonlinear, and conservation-law models.For nonlocal conservation laws, AC discretization demonstrates convergence with or without singular solutions.Source: §10 · p. 62

RK collocation

• RK collocation reproduces multilinear polynomials, preserves quadratic exactness, and satisfies the same uniform consistency result, but all RK schemes fail the discrete maximum principle.Leng et al. (2019) establish stability for the linear scheme by comparing Fourier symbols with a stable Galerkin approximation; higher-order stability remains unresolved.Source: §10 · pp. 63–64

11. Numerical methods for the strong form of fractional diffusion

• Strong-form fractional diffusion is discretized through quadrature-based finite differences, Monte Carlo methods, and radial basis function methods.These approaches contrast with the weak-form discretizations studied earlier.Source: §11 · p. 64

11.1. Quadrature rule-based finite difference methods

• Quadrature-based finite differences split the singular fractional-Laplacian integral into singular and regular parts, enabling specialized treatment alongside standard quadrature.This widely used approach relies on the simplicity of directly discretizing the integral operator and has variants for fractional diffusion models.Source: §11.1 · pp. 64–65
• The resulting finite-difference schemes admit maximum-principle-based stability and error estimates, while higher-dimensional extensions improve convergence through adapted vanishing-viscosity approximations.Huang and Oberman (2014) established the maximum-principle argument for k = 1 and 2, and Del Teso, Endal and Jakobsen (2018, 2019) extended and improved the approach.Source: §11.1 · p. 66
• Uniform-grid discretizations can be translation-invariant and FFT-accelerated, but graded meshes for boundary singularities produce dense unstructured matrices requiring O(N2) storage and O(N3) computation.Hierarchical matrices are proposed as an alternative for unstructured meshes with linear storage and computational complexity (Xu and Darve 2018a).Source: §11.1 · pp. 66–67

11.2. Monte Carlo method by Feynman–Kac formula

• The fractional Feynman–Kac representation replaces Brownian motion with an isotropic α-stable Lévy process, whose paths exit domains by jumps rather than boundary hitting.For (−∆)^s with α = 2s, the representation uses the first exit time from the domain and the associated stable-process distribution.Source: §11.2 · p. 68
• Direct Feynman–Kac Monte Carlo is inefficient because evaluating u(x) requires many stable-process paths, motivating walk-on-spheres methods that exploit distributional symmetry without simulating complete paths.Walk-on-spheres iteratively samples exits from the largest inscribed spheres until leaving a convex domain.Source: §11.2 · pp. 68–69

11.3. Radial basis function methods

• Radial basis functions offer a mesh-free route that may require only minor changes to local-PDE algorithms, especially in high dimensions, but fractional-diffusion studies remain scarce.Rosenfeld, Rosenfeld and Dixon (2019) used compactly supported Wendland functions in a mesh-free pseudo-spectral method.Source: §11.3 · pp. 69–70
• RBF interpolation estimates require smoothness and parameter conditions; experiments found no measurable cutoff-function effect, while optimal convergence without cutoffs and broader collocation stability questions remain open.Open issues include suitable RBF selection, stability, optimal convergence rates, adaptivity, and fast accurate fractional-Laplacian evaluation in high dimensions.Source: §11.3 · p. 70

12. Conditioning and fast solvers

• Nonlocal stiffness-matrix conditioning depends generally on both the mesh size h and horizon δ, so nonlocal systems are not automatically better conditioned than local diffusion systems.For fractional-type translation-invariant kernels, the condition number is typically O(h^-2sδ^2s−2).Source: §12 · p. 71

12.1. Fast algorithm for kernels with non-smooth truncation

• Tian and Engquist (2019) decompose kernels into components that enable FMM or hierarchical-matrix compression away from singularities and non-smooth support boundaries.The interaction boundary is finely resolved with small boxes, while smooth regions use larger boxes.Source: §12.1 · pp. 71–72, Figure 12.1 · p. 72, Figure 12.2 · p. 72
• The kernel-decomposition algorithm achieves O(N log N) complexity in one dimension and O(N^2−1/d) complexity in spatial dimension d.The higher-dimensional cost reflects the increasing number of small boxes needed to resolve co-dimension 1 interaction boundaries.Source: §12.1 · p. 72

12.2. Conditioning and solvers for finite element discretizations of the integral fractional Laplacian model

• For fractional finite element systems, diagonal-preconditioned conjugate gradients converge in a fixed number of steps for small s when the time step is sufficiently small relative to the mesh size.For larger fractional orders or steady-state problems, iteration counts depend on problem size; multigrid restores a uniform iteration bound.Source: §12.2 · pp. 73–74
• Dense fractional stiffness matrices require efficient matrix–vector products; fast transforms and matrix compression achieve quasi-optimal O(Nh log Nh) solve complexity.Fast transforms are limited to uniform meshes, whereas matrix compression is more difficult to implement.Source: §12.2 · p. 74

PART THREE Selected extensions

• The article frames these conditioning and solver results as part of broader extensions to other model settings and additional approximation approaches.The surrounding discussion also introduces applications of the considered nonlocal and fractional models.Source: PART THREE Selected extensions · p. 74

13.1. Indefinite and non-self-adjoint problems

• Indefinite and non-self-adjoint nonlocal problems require inf-sup conditions rather than Lax–Milgram’s symmetric coercivity framework, with Petrov–Galerkin and weakly coercive formulations covering distinct trial and test spaces.The continuous problem is well-posed when the bilinear form is bounded, satisfies an inf-sup condition, and has the stated nondegeneracy property; discrete inf-sup stability must generally be verified for the chosen finite element spaces.Source: §13 · p. 74, §13.1 · p. 75
• Conforming finite element discretizations are well-posed when their discrete spaces satisfy the required approximation and discrete inf-sup conditions, with parameter-independent constants needed for asymptotic compatibility.For fixed δ, parameter independence is unnecessary; it is imposed for the asymptotic compatibility analysis of Tian and Du (2020).Source: §13.1 · pp. 75–76

13.2. Asymptotically compatible schemes

• A scheme is asymptotically compatible when its discrete spaces satisfy the prior conformity and stability requirements and become asymptotically dense in the local PDE energy space as δ → 0 and h → 0.This extends the earlier asymptotic compatibility framework for symmetric, coercive bilinear forms to the discrete problems involving both horizon and grid-size parameters.Source: §13.2 · p. 76

13.3. Operator-splitting and mixed formulations

• Operator-splitting and mixed formulations broaden nonlocal modeling by recasting nonlocal diffusion and Stokes-type systems into first-order or mixed forms, with operator-splitting finite elements recovering local LDG discretizations as δ → 0.Mixed finite element methods require discretizing two-variable quantities in both x and y, whereas operator-splitting methods have been developed for radial, integrable kernels; nonlocal Stokes models also admit spectral and finite difference approximations under periodic boundaries.Source: §13.3 · pp. 76–77

14. Nonlocal convection–diffusion problems

• Fully nonlocal convection–diffusion models replace both diffusion and convection, producing generally non-symmetric kernels that can also represent non-symmetric diffusion.The general formulation treats the two nonlocal terms as separate phenomena and extends the kernel class considered by Du, Huang and Lehoucq (2014), based on D’Elia et al. (2017).Source: §14 · p. 78

14.1. Non-symmetric kernels and nonlocal convection–diffusion operators

• The nonlocal operator permits distinct diffusion and convection horizons and kernels, with tensor-valued diffusion and vector-valued convection parameters that generally yield non-radial, non-translationally invariant kernels.The framework includes specialized one-dimensional and conservative formulations as special cases, while existing asymptotically compatible analyses use specialized kernels and different techniques.Source: §14.1 · pp. 78–80

14.1.1. Steady-state nonlocal convection–diffusion problems

• Steady-state problems are well-posed under nonlocal Poincaré, coercivity, and data assumptions, with solution operators smoothing diffusion problems more strongly than convection problems.The stated analysis is effective for diffusion-dominated problems; convection-dominated cases may require the formulation of Tian and Du (2020).Source: §14.1.1 · pp. 81–82

15. Time-dependent nonlocal problems

• Time-dependent nonlocal problems use steady-state spatial discretizations and local-PDE temporal schemes, with weak-formulation, well-posedness, and error analyses following analogous local pathways.Examples include backward Euler, Crank–Nicolson, leap-frog, and other explicit methods.Source: §15 · p. 82

15.1. Time-dependent nonlocal diffusion

• For radial integrable kernels, nonlocal diffusion preserves diffusive decay properties but solutions may be no smoother than the data, unlike parabolic PDE solutions.Fractional heat equations and nonlocal-in-time variants are also considered.Source: §15.1 · p. 83

15.2. Time-dependent convection–diffusion problems

• Time-dependent nonlocal convection–diffusion formulations are well-posed under coercivity and continuity, requiring only bounded ∥D_c,δμ∥∞ under standard variational arguments.This condition is weaker than the assumptions initially suggested for the nonlocal convection velocity.Source: §15.2 · p. 84

15.3. Nonlocal wave equations

• Nonlocal wave equations conserve energy and exhibit dispersive wave velocities, which converge quadratically in δ to the nondispersive local-wave relation as δ →0.For local waves, ω/k = ±c is independent of frequency and wavenumber; nonlocal velocity is a nonlinear function of k.Source: §15.3 · pp. 84–85

16. Inverse problems

• Inverse problems use sparse or noisy state/output data to identify uncertain sources, volume constraints, coefficients, kernels, horizons, or fractional exponents through regularized control and optimization.The objective combines a state-control cost with regularization, although admissible state or control constraints can sometimes ensure well-posedness without regularization.Source: §16 · pp. 85–86

16.1.1. Distributed optimal control in nonlocal diffusion for a matching functional

• Distributed nonlocal optimal control admits unique optimal state-control pairs and converges to the local problem as δ →0; discontinuous targets are matched better because irregular nonlocal states are admissible.For large horizons, the optimal state can perfectly match a discontinuous target while using a lower-amplitude, smaller-L2-norm control.Source: §16.1.1 · pp. 87–88

16.1.2. Coefficient identification

• Coefficient-identification methods prove existence and discretization convergence, recover smooth and discontinuous diffusion coefficients, and neural-network approaches recover kernel parameters such as δ and variable fractional order s(x).Box constraints can provide well-posedness without regularization, while numerical tests demonstrate recovery for square-integrable and truncated fractional kernels.Source: §16.1.2 · pp. 89–91

16.2.1. Distributed optimal control in fractional diffusion for a matching functional

• Fractional optimal-control studies establish existence, optimality conditions, regularity, stability, and a priori error estimates for elliptic and parabolic problems using finite-element and backward-Euler discretizations.Related formulations also control volume-constraint data g for steady-state and time-dependent fractional problems.Source: §16.2.1 · pp. 91–92
• Antil et al. (2019b) establish well-posedness and regularity for time-dependent nonlocal Dirichlet and Robin optimal control, propose a discretization, and illustrate it numerically in two dimensions.Their main contribution is showing that nonlocal control can use information outside the domain, not only boundary information.Source: §16.2.1 · p. 92
• Antil et al. (2020) consider control problems for spectral fractional semilinear operators and integral fractional Laplacians, with spectral existence and regularity results extendable to the integral setting after small modifications.Source: §16.2.1 · p. 92
• Pang, Lu and Karniadakis (2019b) use physics-informed neural networks to identify the fractional power of the integral fractional Laplacian, as a special case of Pang et al. (2019a).Other cited kernel-identification works address the fractional power only for the spectral fractional-Laplacian definition.Source: §16.2.1 · p. 93
• Control problems involving the spectral fractional Laplacian have also been analyzed by Antil, Otarola and Salgado (2016, 2018) and Ot´arola and Salgado (2018).Source: §16.2.1 · p. 93

Problem formulation

• The nonlocal obstacle problem imposes −Lu ≥ f, u ≥ ψ, complementarity, and homogeneous volume constraints, modeling applications including elastic membranes, contact mechanics, and American put options.The obstacle function ψ specifies the lower bound, and the formulation is the nonlocal analogue of the local PDE obstacle problem.Source: §17 · p. 93

Analysis and related work

• Nonlocal obstacle problems require less obstacle smoothness for well-posedness than local problems, while their solution behavior can differ substantially from the local setting.The section reviews analyses covering regularity, well-posedness, and finite element approximation, including fractional and truncated kernels.Source: §17 · pp. 93–94

Mixed formulation

• The mixed formulation seeks u ∈ V and λ ∈ M satisfying A(u,v) − B(λ,v) = ⟨f,v⟩ and B(η − λ,u − ψ) ≥ 0 for admissible v and η.Here V is the energy space, Vd its dual, A is the nonlocal bilinear form, and M is the closed convex dual cone of nonnegative functionals.Source: §17 · p. 94

Regularity

• Under f ∈ L2(Ω) and (−Lδψ − f)+|Ω ∈ L2(Ω), improved fractional-obstacle regularity holds for every s ∈ (0,1), with β = min{s, 1/2 − ε}.The regularity results are obtained under the stated data and obstacle assumptions, with ε > 0.Source: §17 · p. 95

Finite element approximation

• Mixed finite elements use continuous piecewise linear trial functions and discontinuous linear multiplier functions; discrete inf-sup stability yields a well-posed discrete problem.The multiplier space is constructed with locally bi-orthogonal basis functions, and the stability is established for V_h^d × V_h.Source: §17 · pp. 95–96

Numerical behavior and truncation

• Numerical solutions show distinct local and nonlocal primal fields and multipliers: the nonlocal multiplier lies in L2(Ω), whereas the local multiplier consists of Dirac delta functions; truncated-kernel solutions converge to the untruncated solution.The same fractional-Laplacian results, including the truncated-kernel analysis, are proved for the truncated fractional Laplacian kernel.Source: Figure 17.1 · p. 95, §17 · p. 96

18. Reduced-order modelling

• Reduced-order modelling replaces high-dimensional nonlocal discretizations with low-dimensional surrogates, enabling many parameter studies at substantially lower online cost.This is especially important for nonlocal systems because their discrete matrices are less sparse or even full, creating a double burden of dimensionality and density.Source: §18 · p. 96
• The ROM framework targets parametrized nonlocal diffusion and constructs a basis with Nrom ≪ Nh whose Galerkin solutions remain accurate across selected parameter domains.Parameters enter the constitutive function and kernel, while the high-dimensional model is supplied by a finite-element or other variational discretization.Source: §18 · pp. 96–97
• GRB and POD approaches provide ROMs with error analyses and numerical illustrations for nonlocal diffusion models, while ROMs can be parametrized by horizon δ and fractional exponent s.Source: §18 · pp. 97–98
• A six-dimensional POD basis closely approximates the discontinuous Galerkin high-dimensional solution, while a one-dimensional basis performs poorly for the singular-kernel test.The illustration uses a time-dependent one-dimensional problem with a non-integrable 1/|y −x| kernel and a manufactured jump-discontinuous solution.Source: §18 · pp. 98–99, Figure 18.1 · p. 98

19. A turbulent flow application

• The fractional-viscosity turbulence model studies Richardson-type superdiffusive separation by adding σ(−∆)^s u to the standard Navier–Stokes viscous term.The model uses s ∈(0, 1), a bounded domain, and a divergence-free velocity field; Richardson’s classical prediction concerns mean-square separation proportional to t^3.Source: §19 · pp. 99–100
• For s = 1/3, the modified Navier–Stokes model deviates from Kolmogorov’s k^-5/3 energy-spectrum scaling, whereas other s ∈(0, 1) values remain consistent with it.Gunzburger et al. (2018) add a fractional Laplacian viscosity term and obtain k^-5/3+β decay for s = 1/3.Source: §19 · pp. 99–100
• An implicit backward-Euler IMEX scheme is unconditionally stable and first-order accurate, while two-stage splitting reduces coding effort despite requiring dense fractional-Laplacian solves.The fractional term must remain implicit because lagging it causes serious stability problems; Taylor–Hood finite elements and truncated fractional operators enable fully discrete schemes.Source: §19 · pp. 100–101

20. Peridynamics models for solid mechanics

• Peridynamics provides a nonlocal vector-valued displacement model for solid mechanics that can nucleate and propagate defects such as fractures.The linear state-based model of Mengesha and Du (2014) is a nonlocal analogue of the Navier equations and is defined through energy minimization.Source: §20 · p. 101
• The peridynamic weak problem is well-posed because its bilinear form is coercive and continuous on the associated energy space.Mengesha and Du (2014) establish these properties and apply the Lax–Milgram theorem.Source: §20 · p. 103
• As δ →0 with an appropriately chosen kernel, weak solutions of the nonlocal state-based peridynamic model converge to the classical Navier–Cauchy elasticity solution.This establishes the model’s consistency with local linear elasticity.Source: §20 · p. 103
• Conforming Galerkin approximations containing continuous piecewise linear functions recover the correct local limit as both δ and h diminish, even when δ decreases faster than h.For discontinuous piecewise constant conforming elements, the correct limit is recovered when h = o(δ), whereas fixed bandwidth can produce an incorrect limit.Source: §20 · pp. 104–105

21. Image denoising

• Nonlocal diffusion filters use neighbourhood intensity patterns to preserve image structures and textures more effectively than methods based only on infinitesimal neighbourhoods.The discussion adopts the additive noise model f = u + η and connects convex nonlocal regularization to nonlocal diffusion, including integral fractional problems.Source: §21 · pp. 105–106

21. Image deblurring via minimization of a nonlocal functional

• Nonlocal deblurring is formulated by minimizing a fidelity term plus a kernel-dependent nonlocal energy, whose Euler–Lagrange equation is a nonlocal reaction–diffusion equation.Well-posedness follows from the properties of the kernel γδ.Source: §21.1 · p. 106
• Nonlocal reconstruction shows superior ability to capture high-contrast features compared with the classical total variation method in the illustrated deblurring example.Figure 21.1 compares clean, noisy, nonlocal, and local reconstructions.Source: §21.1 · p. 106, Figure 21.1 · p. 107

21.1.1. Optimization of the denoising parameters

• Bilevel optimization identifies the denoising weight function and has a proved well-posed formulation, a second-order solution algorithm, and numerical tests on benchmark images.The framework constrains w in H1(Ω) with 0 ≤ w(x) ≤ w for finite w.Source: §21.1.1 · pp. 108–109