Carlota M. Cuesta

Loïc Constantin, Carlota M. Cuesta, Francisco de la Hoz

Given a function $u$ defined on $\mathbb R^n$, its fractional $p$-Laplacian is given by $$(-Δ)_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 $$(-Δ)_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^{1-sp/2}},$$ where $C_4(n,s,p)$ is another normalization constant, and $Φ_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 $Δ$ 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}+(-Δ)^s_pu=0,$$in one and two spatial dimensions, being able to capture the self-similar solutions that arise as $t\to\infty$.

Carlota M. Cuesta, Francisco de la Hoz

In this paper we present a non-recursive direct solver, based on the Bartels-Stewart algorithm, for $N$-dimensional Sylvester tensor equations. The method relies only on Schur decompositions of the coefficient matrices and reduces the computation to a single sequential sweep over tensor entries, making it entirely independent of the dimension $N$. Its main advantages are simplicity, a dimension-independent formulation, and the ability to solve very high-dimensional problems limited only by available memory, which is used efficiently. We successfully solve cases up to $N=29$ on a standard laptop with $32$ GB RAM. Compared with the recursive blocked method of Chen and Kressner (state of the art), both approaches achieve identical accuracy. The recursive method is faster for large coefficient matrices, whereas our solver is competitive or superior when matrices are small, especially for large $N$, where recursive methods cannot effectively exploit BLAS-3 kernels. It also uses memory more efficiently: for near-capacity problems (e.g., matrices of order $19$ with $N=7$, where solution and right-hand side occupy $\approx 28$ GB), the Chen-Kressner method exceeds available memory, while ours succeeds. The method is also significantly simpler to implement, fully independent of $N$, and correctly handles singleton dimensions. We detail the algorithm and derive accurate cost estimates. For reproducibility, we provide pseudocode and complete MATLAB implementations. As an application, we compute solutions of linear $N$-dimensional ODE systems with constant coefficients at arbitrary times, and thus of evolutionary PDEs after spatial discretization, including highly accurate solutions of an advection-diffusion equation on $\mathbb{R}^N$.