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$

Abstract: Given a function $u$ defined on $\mathbb R<sup>n$, its fractional $p$-Laplacian is given by $$(-Δ)<em>p^su(\vec x)=C_1(n,s,p)\int</em>{\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<sup>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 $$(-Δ)_p^su(\vec x)=C_4(n,s,p)\int_0^{\inftyΔ(t-Δ)^{{-1}\left[Φ_p(u(\vec}} x)-u(\cdot))\right](\vec x)\frac{dt}{t<sup>{1-sp/2}},$$ where $C_4(n,s,p)$ is another normalization constant, and $Φ_p(t)=|t|<sup>{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<sup>n$. Our approach builds on the Balakrishnan representation, where we discretize the 2nd-order derivatives in $Δ$ 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}+(-Δ)<sup>s_pu=0,$in one and two spatial dimensions, being able to capture the self-similar solutions that arise as $t\to\infty$.

• The paper introduces a matrix-based pseudospectral method that transforms unbounded domains into finite ones for efficient computation of fractional Laplacians.
• It exploits Chebyshev mappings and diagonalizable differentiation matrices to achieve spectral accuracy with errors as low as 10⁻¹¹.
• The method extends to the nonlinear fractional p-Laplacian, enabling robust simulations of evolution equations and self-similar asymptotics.

A Matrix-Based Spectral Method for the Numerical Approximation of the Fractional Laplacian and Fractional $p$-Laplacian on $\mathbb{R}^n$

---

Introduction and Motivation

The fractional Laplacian and fractional $p$-Laplacian are nonlocal operators with significant roles in mathematical analysis, PDEs, and many applied disciplines including physics and machine learning. While the fractional Laplacian ($p=2$) features in linear diffusion and stochastic processes, the fractional $p$-Laplacian ($p\neq2$) governs nonlinear, nonlocal phenomena and arises in modeling anomalous diffusion and nonlinear nonlocal flows.

Traditional numerical schemes for these operators encounter substantial obstacles: domain truncation, intricate quadrature rules, poor scalability in higher dimensions, and limited spectral accuracy. The present work addresses these limitations by developing a matrix-based pseudospectral scheme leveraging the Balakrishnan (and related Bochner) integral representations, enabling efficient spectral computation of both fractional Laplacian and $p$-Laplacian in arbitrary dimension on unbounded domains.

---

Spectral Differentiation Matrix Construction

The method begins with a rational Chebyshev mapping $x = L \cot \xi$ transforming $\mathbb{R}$ to a finite domain $(0, \pi)$, with nodes corresponding to mapped Chebyshev points. The function $\mathbb{R}^n$0 is composed as $\mathbb{R}^n$1, and spectral differentiation matrices $\mathbb{R}^n$2 and $\mathbb{R}^n$3 are constructed via the classical formulas for periodic grids, with careful even extension to ensure compatibility with MATLAB FFT-based spectral algorithms. Second derivatives $\mathbb{R}^n$4 are recoverable from these matrices by explicit, stable transformations.

Substantial attention is given to optimal matrix construction: only half-rows are needed due to symmetries; Toeplitz/circulant structure is exploited for computational efficiency; and careful asymptotic manipulation ensures numerical stability, even for large node numbers (e.g., $\mathbb{R}^n$5). Empirical tests on analytic functions demonstrate discrete $\mathbb{R}^n$6 errors as low as $\mathbb{R}^n$7 for derivative approximation, validating spectral precision and stability for extreme $\mathbb{R}^n$8.

---

Diagonalization and Operator Representation

The differentiation matrices are shown to be stably diagonalizable, a critical advance enabling evaluation of fractional powers. Empirically, the condition number growth is sublinear in $\mathbb{R}^n$9, and the eigenstructure is well-behaved, permitting robust spectral analysis for operator powers—unlike standard Chebyshev matrices, which are often ill-conditioned.

For multidimensional Laplacians, tensor products of the $p$0 matrices are formed with analogous diagonalization. The eigenvalue spectrum permits reduction of the Balakrishnan (or Bochner) integral representation for the fractional operator to analytical evaluation at the spectral level, avoiding quadrature altogether.

The paper proves (for any $p$1) that the spectral fractional Laplacian action reduces to a matrix product involving Hadamard powers of $p$2 applied to the eigenvector projections of $p$3, yielding:

$p$4

This formulation generalizes efficiently to arbitrary dimension, is free of domain truncation, and sidesteps variational approaches.

---

Extension to the Fractional $p$5-Laplacian

For $p$6, the fractional $p$7-Laplacian is nonlinear and introduces the function $p$8, with additional subtleties for $p$9. The matrix-based spectral approximation works by forming, for each point, the difference array, applying $p=2$0 pointwise, and then using the spectral fractional Laplacian machinery (with order $p=2$1). The appropriate scaling constants are derived analytically, linking the Balakrishnan and Bochner formulations and ensuring proper normalization.

The $p=2$2-dimensional fractional $p=2$3-Laplacian is then computed as a sequence of spectral projections, Hadamard powers, and matrix multiplications applied to the nonlinear transformed difference data. MATLAB implementations exploit vectorization and symmetry for both loop-based and loopless methods, balancing memory and speed.

---

Numerical Experiments: Convergence, Accuracy, and Scaling

Extensive benchmarks are presented for both operators. For radial functions with known analytic fractional Laplacians (involving hypergeometric functions), the spectral errors in multi-dimensional settings rapidly decay to $p=2$4 or better for modest $p=2$5, confirming high accuracy. The accuracy improves with increased $p=2$6 and $p=2$7, and is particularly high for the half-Laplacian ($p=2$8).

For $p=2$9, the method exhibits coherent numerical behavior even in the absence of analytic solutions. Singularities and structural transitions (number of peaks) are observed as $p$0 passes critical values ($p$1), with two distinct regimes: single-peaked for $p$2, dual-peaked for $p$3. This is demonstrated both in 1D and 2D, agreeing with theoretical predictions and previous qualitative studies.

---

Evolution Equations: Asymptotic Self-Similarity

The method is applied to time-dependent fractional $p$4-Laplacian evolution equations. The Runge-Kutta integration is combined with spectral spatial discretization, and mass conservation and asymptotic self-similarity are numerically confirmed. The simulated profiles in self-similar variables converge tightly to the predicted analytic scaling forms for large $p$5, corroborating theoretical results on the fundamental solution and asymptotic regime for $p$6 above critical values.

Mass is spectrally computed with high precision; accurate self-similarity is achieved for both $p$7 and $p$8 across multiple parameter regimes. Challenging cases requiring large $p$9 and $p\neq2$0 are handled efficiently, and no domain truncation is necessary. This capability to handle fast diffusion ($p\neq2$1), normal diffusion, and slow diffusion ($p\neq2$2) in a single framework is a notable computational advantage.

---

Practical and Theoretical Implications

The matrix-based spectral scheme offers:

• High spectral accuracy for derivatives and fractional powers in unbounded domains;
• Computational efficiency and scalability in dimensions up to at least $p\neq2$3 without domain truncation or costly quadrature;
• Stable diagonalization enabling fractional power computation even for very large $p\neq2$4;
• Direct extension to nonlinear operators, supporting fractional $p\neq2$5-Laplacian computations;
• Robust simulation of evolution equations, elucidating self-similar asymptotics and mass conservation.

These features facilitate advanced numerical studies in nonlocal nonlinear PDEs, probabilistic models, and applied domains requiring precise computation in unbounded geometry. From a theoretical standpoint, the approach provides computational validation for qualitative and asymptotic analysis, and opens doors to broader generalizations (possibly variable order, mixed operators, or non-Euclidean geometries).

---

Conclusion

This work presents an efficient, stable, and spectrally accurate matrix-based pseudospectral method for computing the fractional Laplacian and fractional $p\neq2$6-Laplacian in arbitrary dimensions on unbounded domains. By leveraging analytical diagonalization, the Balakrishnan and Bochner integral representations are reduced to algebraic spectral manipulation, sidestepping domain truncation and expensive quadrature. Numerical results confirm high precision, robust scaling, and compatibility with nonlinear and time-dependent settings, marking an advance in computational nonlocal analysis and PDE discretization (2605.23252).