pith. sign in

arxiv: 2605.23252 · v1 · pith:7B54I4EKnew · submitted 2026-05-22 · 🧮 math.NA · cs.NA

A matrix-based spectral method for the numerical approximation of the fractional Laplacian and the fractional p-Laplacian of functions defined on mathbb R^n

Lo\"ic Constantin , Carlota M. Cuesta , Francisco de la Hoz This is my paper

Pith reviewed 2026-05-25 04:01 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords fractional Laplacianfractional p-Laplacianspectral methodBalakrishnan representationdifferentiation matricesnumerical approximation
0
0 comments X

The pith

Spectrally accurate differentiation matrices diagonalize the Laplacian so the fractional p-Laplacian reduces to an exact evaluation on the eigenvalue spectrum.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper introduces a matrix-based spectral method for approximating the fractional Laplacian and the fractional p-Laplacian on R^n. It starts from the Balakrishnan integral representation of these operators and discretizes the second-order Laplacian inside the integral using differentiation matrices that are spectrally accurate. These matrices admit a well-conditioned diagonalization, so the fractional operator can be applied directly to the eigenvalues and the integral collapses to an analytical spectral formula. The resulting scheme works in any dimension, requires no domain truncation and no multidimensional quadrature, and is demonstrated on the evolution equation whose solutions approach self-similar profiles.

Core claim

Discretizing the Laplacian with spectrally accurate differentiation matrices that can be diagonalized in a well-conditioned way lets the Balakrishnan integral for both the linear fractional Laplacian and the nonlinear fractional p-Laplacian be replaced by an exact analytical evaluation performed on the eigenvalue spectrum, producing a stable numerical scheme that extends immediately to arbitrary spatial dimensions without truncation or variational formulations.

What carries the argument

Well-conditioned diagonalization of spectrally accurate differentiation matrices for the standard Laplacian, which moves the fractional operator from the spatial grid to direct action on the eigenvalue spectrum.

If this is right

  • The scheme extends without change to any spatial dimension n.
  • The evolution equation partial_t u + (-Delta)_p^s u = 0 can be integrated stably in one and two dimensions.
  • Self-similar solutions are recovered as t tends to infinity.
  • No domain truncation or variational weak form is required.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The same diagonalization step could be reused for other nonlocal operators that possess a Balakrishnan-type integral representation.
  • The computational saving from replacing quadrature by a spectral sum grows rapidly with dimension.
  • The method supplies a direct way to compute steady states of related fractional equations by setting the time derivative to zero.

Load-bearing premise

The spectral differentiation matrices for the ordinary Laplacian remain accurate and well-conditioned when they are substituted into the Balakrishnan integral, including for the nonlinear term involving Phi_p.

What would settle it

A direct numerical comparison, in one or two dimensions, between the computed long-time solutions of the evolution equation and the known exact self-similar profiles would confirm or refute whether the scheme reproduces the correct asymptotics.

Figures

Figures reproduced from arXiv: 2605.23252 by Carlota M. Cuesta, Francisco de la Hoz, Lo\"ic Constantin.

Figure 1
Figure 1. Figure 1: Errors in semilogarithmic scale corresponding to (34), in on [PITH_FULL_IMAGE:figures/full_fig_p020_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Errors in semilogarithmic scale corresponding to (43) with [PITH_FULL_IMAGE:figures/full_fig_p021_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Errors in semilogarithmic scale corresponding to (43) with [PITH_FULL_IMAGE:figures/full_fig_p021_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: (−∆)s p e −x 2 1, for s = 0.8, N = 2000, L = 20. Left: p ∈ {1.01, 1.02, . . ., 2.49}. Right: p ∈ {1.01, 1.02, . . ., 9.99}. x -5 -4 -3 -2 -1 0 1 2 3 4 5 (−∆) s p e −x 2 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 s = 0.8, N = 2000, L = 20 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 x -5 -4 -3 -2 -1 0 1 2 3 4 5 (−∆) s p e −x 2 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 s = 0.8, N = 2000, L = 20 2.05 2.10 2.15 2.20 2.… view at source ↗
Figure 5
Figure 5. Figure 5: (−∆)s p e −x 2 1, for s = 0.8, N = 2000, L = 20. Left: p ∈ {1.60, 1.65, . . ., 1.95}. Right: p ∈ {2.05, 2.10, . . ., 2.45} [PITH_FULL_IMAGE:figures/full_fig_p022_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: (−∆)s p e −x 2 1−x 2 2 , for s = 0.8, N = 250, L = 4. Left: p = 1.1. Right: p = 2.4. 22 [PITH_FULL_IMAGE:figures/full_fig_p022_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Error in the spectral approximation of the mass [PITH_FULL_IMAGE:figures/full_fig_p025_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Solutions in the self-similar variables v(x, t) = u(x, t)M(t) −spβt α and r(x, t) = M(t) (2−p)β t −βx, for n = 1 and p ∈ (pc, p1). Left: s = 0.2, p = 1.74, with N = 1001, L = 2 × 104 , ∆t = 10−2 , t ∈ {18.9, 18.91, . . ., 19}. Right: s = 0.8, p = 1.4, with N = 501, L = 10, ∆t = 10−5 , t ∈ {4.8, 4.81, . . ., 4.9}. and 11 correspond to the case p ∈ (pc, p1); Figures 9 and 12 to the case p ∈ (p1, 2); and Figu… view at source ↗
Figure 9
Figure 9. Figure 9: Solutions in the self-similar variables v(x, t) = u(x, t)M(t) −spβt α and r(x, t) = M(t) (2−p)β t −βx, for n = 1 and p ∈ (p1, 2). Left: s = 0.2, p = 1.9, with N = 5001, L = 500, ∆t = 0.1, t ∈ {19, 19.1, . . ., 20}. Right: s = 0.8, p = 1.8, with N = 501, L = 10, ∆t = 10−3 , t ∈ {1.8, 1.81, . . ., 1.9}. r(x, t) -5 -4 -3 -2 -1 0 1 2 3 4 5 v(x, t) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 n = 1, s = 0.2, p = 2.5, N = 2001… view at source ↗
Figure 10
Figure 10. Figure 10: Solutions in the self-similar variables v(x, t) = u(x, t)M(t) −spβt α and r(x, t) = M(t) (2−p)β t −βx, for n = 1 and p ∈ (2, 2/s). Left: s = 0.2, p = 2.5, with N = 2001, L = 10, ∆t = 1, t ∈ {180, 181, . . ., 190}. Right: s = 0.8, p = 2.2, with N = 1501, L = 10, ∆t = 5 × 10−4 , t ∈ {18, 18.1, . . ., 19}. -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 r(x, t) ×10-13 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 v(x, t) ×1023n = 2, s = 0.… view at source ↗
Figure 11
Figure 11. Figure 11: Solutions in the self-similar variables v(x, t) = u(x, 0, t)M(t) −spβt α and r(x, t) = M(t) (2−p)β t −βx, for n = 2 and p ∈ (pc, p1). Left: s = 0.2, p = 1.84, with N = 151, L = 104 , ∆t = 10−2 , t ∈ {38.9, 38.91, . . ., 39}. Right: s = 0.8, p = 1.6, with N = 81, L = 15, ∆t = 2 × 10−4 , t ∈ {4.8, 4.81, . . ., 4.9}. 26 [PITH_FULL_IMAGE:figures/full_fig_p026_11.png] view at source ↗
Figure 12
Figure 12. Figure 12: Solutions in the self-similar variables v(x, t) = u(x, 0, t)M(t) −spβt α and r(x, t) = M(t) (2−p)β t −βx, for n = 2 and p ∈ (p1, 2). Left: s = 0.2, p = 1.95, with N = 201, L = 75, ∆t = 0.1, t ∈ {18, 18.1, . . ., 19}. Right: s = 0.8, p = 1.7, with N = 101, L = 10, ∆t = 0.001, t ∈ {4.8, 4.81, . . ., 4.9}. r(x, t) -5 -4 -3 -2 -1 0 1 2 3 4 5 v(x, t) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 n = 2, s = 0.2, p = 2.5, N = 2… view at source ↗
Figure 13
Figure 13. Figure 13: Solutions in the self-similar variables v(x, t) = u(x, 0, t)M(t) −spβt α and r(x, t) = M(t) (2−p)β t −βx, for n = 2 and p ∈ (2, 2/s). Left: s = 0.2, p = 2.5, with N = 201, L = 3, ∆t = 1, t ∈ {80, 81, . . ., 90}. Right: s = 0.8, p = 2.2, with N = 101, L = 1, ∆t = 5×10−4 , t ∈ {1.8, 1.81, . . ., 1.9}. 27 [PITH_FULL_IMAGE:figures/full_fig_p027_13.png] view at source ↗
read the original abstract

Given a function $u$ defined on $\mathbb R^n$, its fractional $p$-Laplacian is given by $$(-\Delta)_p^su(\vec x)=C_1(n,s,p)\int_{\mathbb R^n}\frac{|u(\vec x)-u(\vec y)|^{p-2}(u(\vec x)-u(\vec y))}{\|\vec x-\vec y\|_2^{n+sp}}d\vec y,\quad\vec x\in\mathbb R^n,$$where the integral is understood in the principal value sense, $p\in(1,\infty)$, $s\in(0,1)$, and $C_1(n,s,p)$ is a normalization constant. A formally equivalent nonlinear Balakrishnan formulation is given by $$(-\Delta)_p^su(\vec x)=C_4(n,s,p)\int_0^\infty\Delta(t-\Delta)^{-1}\left[\Phi_p(u(\vec x)-u(\cdot))\right](\vec x)\frac{dt}{t^{1-sp/2}},$$ where $C_4(n,s,p)$ is another normalization constant, and $\Phi_p(t)=|t|^{p-2}t$. In this paper, we present a matrix-based spectral method to approximate numerically the fractional Laplacian (i.e., the linear case, where $p = 2$) and the fractional $p$-Laplacian for functions defined on $\mathbb R^n$. Our approach builds on the Balakrishnan representation, where we discretize the 2nd-order derivatives in $\Delta$ using spectrally accurate differentiation matrices. A key advantage is that these matrices can be diagonalized in a well-conditioned manner, enabling a stable and robust numerical scheme that naturally extends to arbitrary spatial dimensions $n$. In particular, this diagonalization allows the fractional operator to act directly on the eigenvalue spectrum, effectively reducing the Balakrishnan integral to an analytical evaluation at the spectral level and thereby avoiding costly multidimensional quadrature. The resulting method also avoids domain truncation and variational formulations, making it both computationally efficient and conceptually straightforward. As a practical application, we simulate the evolution of $$\frac{\partial u}{\partial t}+(-\Delta)^s_pu=0,$$in one and two spatial dimensions, being able to capture the self-similar solutions that arise as $t\to\infty$.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 1 minor

Summary. The paper proposes a matrix-based spectral method for numerically approximating the fractional Laplacian (p=2) and fractional p-Laplacian on R^n. It discretizes the Laplacian via spectrally accurate differentiation matrices within the Balakrishnan representation, claims that diagonalization reduces the integral representation to an analytical spectral-level evaluation (avoiding multidimensional quadrature and domain truncation), and applies the scheme to simulate the evolution equation ∂u/∂t + (-Δ)_p^s u = 0 in 1D and 2D, capturing self-similar solutions as t→∞.

Significance. A stable, dimension-independent scheme that avoids truncation and quadrature would be useful for fractional PDEs if the claimed reduction and accuracy hold. The extension from linear to nonlinear p-Laplacian and the reproduction of self-similar behavior are potentially valuable, but the absence of quantitative validation limits the assessed impact.

major comments (2)
  1. [Abstract] Abstract (method description paragraph): the claim that diagonalization 'effectively reduc[es] the Balakrishnan integral to an analytical evaluation at the spectral level' does not hold for the nonlinear fractional p-Laplacian (p≠2). The term Φ_p(u(x)−u(·)) produces a general vector; each resolvent (t−Δ)^{-1} must still be applied to this vector in the eigenbasis and the outer integral over t remains a numerical quadrature whose error is not analyzed.
  2. [Numerical results] Numerical experiments (inferred from abstract claim of capturing self-similar solutions): no error tables, convergence rates, or comparisons against known exact solutions or alternative discretizations are provided to support the assertions of stability, robustness, and accuracy for either the linear or nonlinear operator.
minor comments (1)
  1. [Abstract] The explicit forms or literature references for the normalization constants C_1(n,s,p) and C_4(n,s,p) are not supplied in the abstract or method outline.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the constructive comments, which help clarify the presentation of our method. We address each major comment below.

read point-by-point responses

  1. Referee: [Abstract] Abstract (method description paragraph): the claim that diagonalization 'effectively reduc[es] the Balakrishnan integral to an analytical evaluation at the spectral level' does not hold for the nonlinear fractional p-Laplacian (p≠2). The term Φ_p(u(x)−u(·)) produces a general vector; each resolvent (t−Δ)^{-1} must still be applied to this vector in the eigenbasis and the outer integral over t remains a numerical quadrature whose error is not analyzed.

    Authors: We agree that the abstract statement is imprecise for p≠2. For the linear case (p=2), the operator commutes with the eigenbasis, allowing the Balakrishnan integral to be evaluated analytically on the eigenvalues after diagonalization. For the nonlinear case, Φ_p is applied pointwise in physical space to produce a vector; this vector is then transformed to the eigenbasis where the resolvents act by diagonal scaling, after which the result is transformed back. The outer integral over t nevertheless requires numerical quadrature. We will revise the abstract and method description to distinguish the two cases explicitly, remove the claim of full analytical reduction for the nonlinear operator, and add a brief discussion of the quadrature procedure and its error control. revision: yes

  2. Referee: [Numerical results] Numerical experiments (inferred from abstract claim of capturing self-similar solutions): no error tables, convergence rates, or comparisons against known exact solutions or alternative discretizations are provided to support the assertions of stability, robustness, and accuracy for either the linear or nonlinear operator.

    Authors: The manuscript prioritizes the derivation of the matrix-based spectral discretization and the qualitative demonstration that self-similar profiles emerge in the evolution equation. We acknowledge that quantitative evidence (error tables, observed convergence rates, and benchmark comparisons) is needed to substantiate claims of accuracy and stability. In the revised manuscript we will add a dedicated numerical validation section containing grid-refinement studies, L^∞ and L^2 error tables for the linear operator against known Fourier multipliers, and, where feasible, comparisons with existing discretizations for both the linear and nonlinear cases. revision: yes

Circularity Check

0 steps flagged

No circularity: standard Balakrishnan representation and spectral matrices applied without self-referential reduction

full rationale

The paper applies the known Balakrishnan integral representation (both linear and nonlinear forms) to discretize via differentiation matrices that are diagonalized using standard linear algebra. The claimed reduction to spectral-level evaluation holds exactly for the linear p=2 case by eigenmode decoupling; for p≠2 the method still performs matrix-vector operations inside a numerical t-integral but does not redefine any quantity in terms of itself or rename a fit as a prediction. No uniqueness theorems, ansatzes smuggled via self-citation, or self-definitional steps appear in the given text. The scheme is a direct numerical construction whose accuracy is benchmarked externally rather than forced by its own inputs.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The approach relies on the equivalence of the integral definition and the Balakrishnan formula (a standard result in the field) together with the spectral accuracy of differentiation matrices for the ordinary Laplacian; no new free parameters or invented entities are introduced.

axioms (2)
  • domain assumption The Balakrishnan representation is formally equivalent to the principal-value integral definition of the fractional p-Laplacian.
    Invoked in the abstract as the starting point for the numerical scheme.
  • domain assumption Spectrally accurate differentiation matrices for the standard Laplacian can be stably diagonalized and used inside the parameter integral.
    Central to reducing the integral to an analytical spectral evaluation.

pith-pipeline@v0.9.0 · 6000 in / 1510 out tokens · 26338 ms · 2026-05-25T04:01:21.773095+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

18 extracted references · 18 canonical work pages · 1 internal anchor

  1. [1]

    J. P. Boyd. The Optimization of Convergence for Chebyshev Poly nomial Methods in an Unbounded Domain. Journal of Computational Physics , 45(1):43–79, 1982

    work page 1982

  2. [2]

    J. P. Boyd. Chebyshev and Fourier Spectral Methods . Dover Publications, (2nd revised edition) edition, 2001

    work page 2001

  3. [3]

    Cayama, C

    J. Cayama, C. M. Cuesta, and F. de la Hoz. Numerical approximat ion of the fractional Laplacian on R using orthogonal families. Applied Numerical Mathematics , 158:164–193, 2020

    work page 2020

  4. [4]

    Cayama, C

    J. Cayama, C. M. Cuesta, and F. de la Hoz. A pseudospectral me thod for the one-dimensional fractional Laplacian on R. Applied Mathematics and Computation , 389:125577, 2021

    work page 2021

  5. [5]

    Cayama, C

    J. Cayama, C. M. Cuesta, F. de la Hoz, and C. J. Garc ´ ıa-Cerver a. A fast convolution method for the fractional Laplacian in R. Revista de la Real Academia de Ciencias Exactas, F ´ ısicas y Naturales. Serie A. Matem´ aticas, 119(80):1–45, 2025

    work page 2025

  6. [6]

    C. M. Cuesta, F. de la Hoz, and I. Girona. Numerical computation of the half Laplacian by means of a fast convolution algorithm. Electronic Transactions on Numerical Analysis , 60:59–98, 2024

    work page 2024

  7. [7]

    A Non-Recursive, Dimension-Independent Schur-Decomposition Algorithm for $N$-Dimensional Sylvester Tensor Equations

    Carlota M. Cuesta and Francisco de la Hoz. A non-recursive Schu r-Decomposition Algorithm for N -Dimensional Matrix Equations. arXiv:2412.15840, 2024

    work page internal anchor Pith review Pith/arXiv arXiv 2024

  8. [8]

    Jakobsen

    F´ elix del Teso, Jø rgen Endal, and Espen R. Jakobsen. Robust n umerical methods for nonlocal (and local) equations of porous medium type. Part II: Schemes and expe riments. SIAM J. Numer. Anal. , 56(6):3611–3647, 2018

    work page 2018

  9. [9]

    Jakobsen

    Felix del Teso, Jørgen Endal, and Espen R. Jakobsen. Robust Nu merical Methods for Nonlocal (and Local) Equations of Porous Medium Type. Part I: Theory. SIAM Journal on Numerical Analysis , 57(5):2266–2299, 2019

    work page 2019

  10. [10]

    Higher-order a symptotic expansions and finite difference schemes for the fractional p-Laplacian

    F´ elix del Teso, Mar ´ ıa Medina, and Pablo Ochoa. Higher-order a symptotic expansions and finite difference schemes for the fractional p-Laplacian. Math. Ann. , 390(1):157–203, 2024

    work page 2024

  11. [11]

    Three Representations of the Fractional p-Laplacian: Semigroup, Extension and Balakrishnan Formulas

    F´ elix del Teso, David G´ omez-Castro, and Juan Luis V´ azquez. Three Representations of the Fractional p-Laplacian: Semigroup, Extension and Balakrishnan Formulas. Fractional Calculus and Applied Analysis, 24(4):966–1002, 2021

    work page 2021

  12. [12]

    Kwa´ snicki

    M. Kwa´ snicki. Ten equivalent definitions of the Fractional Lapla ce Operator. Fractional Calculus and Applied Analysis , 20(1):7–51, 2017

    work page 2017

  13. [13]

    Fast Fourier-like Mapped Chebyshev Spectral-Galerkin Methods for PDEs with Integral Fra ctional Laplacian in Unbounded Domains

    Changtao Sheng, Jie Shen, Tao Tang, Li-Lian Wang, and Huifang Yuan. Fast Fourier-like Mapped Chebyshev Spectral-Galerkin Methods for PDEs with Integral Fra ctional Laplacian in Unbounded Domains. SIAM Journal on Numerical Analysis , 58(5):2435–2464, 2020

    work page 2020

  14. [14]

    MATLAB, Version R2025b, 2025

    The MathWorks Inc. MATLAB, Version R2025b, 2025

    work page 2025

  15. [15]

    Spectral Methods in MATLAB

    Lloyd Nicholas Trefethen. Spectral Methods in MATLAB . SIAM, 2000

    work page 2000

  16. [16]

    The evolution fractional p-Laplacian equation in RN : Fundamental solution and asymptotic behaviour

    Juan Luis V´ azquez. The evolution fractional p-Laplacian equation in RN : Fundamental solution and asymptotic behaviour. Nonlinear Analysis, 199:112034, 2020

    work page 2020

  17. [17]

    The fractional p-Laplacian evolution equation in RN in the sublinear case

    Juan Luis V´ azquez. The fractional p-Laplacian evolution equation in RN in the sublinear case. Calculus of Variations and Partial Differential Equations , 60(140):1–59, 2021

    work page 2021

  18. [18]

    Growing solutions of the fractional p-Laplacian equation in the Fast Diffusion range

    Juan Luis V´ azquez. Growing solutions of the fractional p-Laplacian equation in the Fast Diffusion range. Nonlinear Analysis, 214:112575, 2022. 29

    work page 2022