We present a general framework, treating Lipschitz domains in Riemannian manifolds, that provides conditions guaranteeing the existence of norming sets and generalized local polynomial reproductions—a powerful tool used in the analysis of various mesh-free methods and a mesh-free method in its own right. As a key application, we prove the existence of smooth local polynomial reproductions on compact subsets of algebraic manifolds in ℝ^n with Lipschitz boundary. These results are then applied to derive new findings on the existence, stability, regularity, locality, and approximation properties of shape functions for a coordinate-free moving least squares approximation method on algebraic manifolds, which operates directly on point clouds without requiring tangent plane approximations. There are two appendices: the first derives high order Markov inequalities for polynomials on algebraic manifolds and the second gives instructions for calculating the dimension of the space of degree m polynomials restricted to a real algebraic variety.

We present a high-order radial basis function-generated finite difference (RBF-FD) method for partial differential equations on moving manifolds ℳ(t)⊂ℝ^3 of co-dimension one. Our method builds on the tangent-plane formulation of surface RBF-FD and combines Lagrangian and Eulerian treatments: the manifold and material derivative are evolved in a Lagrangian fashion, while the remaining surface differential operators are reconstructed on the instantaneous point cloud and stabilized, when necessary, by a quasi-analytical hyperviscosity formulation. Lagrangian marker drift is handled by adaptive rearrangement using a compact global parametric model of the moving surface, which also supplies accurate normals and geometry-based quadrature. After rearrangement, we reconstruct the multistep history by backward semi-Lagrangian tracing and interpolation before resuming the Lagrangian time discretization. We exploit temporal coherence through three update strategies: defect correction for the local RBF-FD weights, a curvature-based update of the hyperviscosity coefficients between spectral recomputations, and a global defect-correction iteration that reuses an incomplete LU (ILU) factorization before preconditioned generalized minimal residual (GMRES) iterations. Finally, we enforce the prescribed global mass balance through a scalar projection based on the evolving surface quadrature. Numerical experiments demonstrate high-order convergence, conservation to roundoff in source-free problems, stable long-time integration, and substantial savings from the proposed update strategies.

We analyze rates of uniform convergence for a class of high-order semi-Lagrangian schemes for first-order, time-dependent partial differential equations on embedded submanifolds of Rd (including advection equations on surfaces) by extending the error analysis of Falcone and Ferretti [SIAM J. Numer. Anal., 35 (1998), pp. 909--940]. A central requirement in our analysis is a remapping operator that achieves both high approximation orders and strong stability, a combination that is challenging to obtain and of independent interest. For this task, we propose a novel mesh-free remapping operator based on \ell1 minimizing generalized polynomial reproduction, which uses only point values and requires no additional geometric information from the manifold (such as access to tangent spaces or curvature). Our framework also rigorously addresses the numerical solution of ordinary differential equations on manifolds via projection methods. We include numerical experiments that support the theoretical results and also suggest some new directions for future research.

Semi-implicit semi-Lagrangian (SISL) methods are commonly used for the shallow water equations (SWE) because they allow for larger time steps than those permitted by the Courant-Friedrichs-Lewy (CFL) stability condition in Eulerian schemes. In these methods, the semi-Lagrangian treatment of advection is typically performed using lower-order interpolation, such as tensor-product Lagrange interpolation with cubic or quintic polynomials. However, operational SISL schemes routinely employ spectrally accurate spatial discretizations, such as spherical harmonics or the double Fourier sphere (DFS) method, for computing horizontal derivatives of the prognostic variables. This creates a mismatch in numerical accuracy, making the use of low-order interpolation less clearly justified. In this work, we present the first numerical investigation of spectrally accurate interpolation in SISL schemes for the SWE. Our approach builds upon the recently developed DFS-based SWE model, incorporating a spectral interpolation scheme that is accelerated using the nonuniform fast Fourier transform (NUFFT) to maintain the same overall computational complexity as the original model. Using several standard SWE test cases, we evaluate the accuracy, conservation, and numerical diffusion of the new model, particularly over long integration times. Compared to an equivalent SISL model with low-order interpolation, the new model achieves higher accuracy, improved mass and energy conservation, and reduced numerical diffusion, demonstrating the potential benefits of incorporating spectrally accurate interpolation into SISL schemes.

This article treats kernel approximation and interpolation on embedded manifolds of ℝ^Nusing restrictions of positive and conditionally positive definite kernels. The main challenge is to develop an approximation theory that treats error measured in highly regular smoothness spaces relative to the kernel. This means that the order of smoothness is higher than that of the kernel's associated native space (in the positive definite case, the reproducing kernel Hilbert space generated by the kernel). This prevents the use of standard techniques for controlling error in this setting, especially RKHS space arguments like orthogonality of the interpolation projector, or bounds using the power function. We generalize an approximation scheme introduced by DeVore and Ron which treats target functions that are in the range of the kernel's integral operator. In the case of embedded manifolds, this generalization is now feasible due to recently developed local polynomial reproductions for certain submanifolds of ℝ^N. Furthermore, we give sufficient conditions on kernel and manifold which allow the range of the integral operator to be precisely identified: in particular, guaranteeing that the range is a Sobolev space. Finally, we provide new kernel-based Bernstein inequalities for embedded manifolds which lead to estimates for interpolation in Sobolev spaces compactly contained in the native space.

Spherical and polar geometries arise in many important areas of computational science, including weather and climate forecasting, optics, and astrophysics. In these applications, tensor-product grids are often used to represent unknowns. However, interpolation schemes that exploit the tensor-product structure can introduce artificial boundaries at the poles in spherical coordinates and at the origin in polar coordinates, leading to numerical challenges, especially for high-order methods. In this paper, we present new bivariate trigonometric barycentric interpolation formulas for spheres and bivariate trigonometric/polynomial barycentric formulas for disks, designed to overcome these issues. These formulas are also efficient, as they only rely on a set of (precomputed) weights that depend on the grid structure and not the data itself. The formulas are based on the double Fourier sphere method, which transforms the sphere into a doubly periodic domain and the disk into a domain without an artificial boundary at the origin. For standard tensor-product grids, the proposed formulas exhibit exponential convergence when approximating smooth functions. We provide numerical results to demonstrate these convergence rates and showcase an application of the spherical barycentric formulas in a semi-Lagrangian advection scheme for solving the tracer transport equation on the sphere.

We investigate the spectrum of differentiation matrices for certain operators on the sphere that are generated from collocation at a set of scattered points $X$ with positive definite and conditionally positive definite kernels. We focus on the cases where these matrices are constructed from collocation using all the points in $X$ and from local subsets of points (or stencils) in $X$. The former case are called global methods (e.g., the Kansa or radial basis function (RBF) pseudospectral method), while the latter are referred to as local methods (e.g., the RBF finite difference (RBF-FD) method). Both techniques are used extensively for numerically solving certain partial differential equations on spheres, as well as other domains. For time-dependent PDEs like the diffusion equation, the spectrum of the differentiation matrices and their stability under perturbations are central to understanding the temporal stability of the underlying numerical schemes. In the global case, we present a perturbation estimate for differentiation matrices which discretize operators that commute with the Laplace-Beltrami operator. In doing so, we demonstrate that if such an operator has negative (non-positive) spectrum, then the differentiation matrix does, too. For conditionally positive definite kernels this is particularly challenging since the differentiation matrices are not necessarily diagonalizable. This perturbation theory is then used to obtain bounds on the spectra of the local RBF-FD differentiation matrices based on the conditionally positive definite surface spline kernels. Numerical results are presented to confirm the theoretical estimates.

Recovering pressure fields from image velocimetry measurements has two general strategies: (i) directly integrating the pressure gradients from the momentum equation and (ii) solving or enforcing the pressure Poisson equation (divergence of the pressure gradients). In this work, we analyze the error propagation of the former strategy and provide some practical insights. For example, we establish the error scaling laws for the pressure gradient integration (PGI) and the pressure Poisson equation. We explain why applying the Helmholtz–Hodge decomposition (HHD) could significantly reduce the error propagation for the PGI. We also propose to use a novel HHD-based pressure field reconstruction strategy that offers the following advantages or features: (i) effective processing of noisy scattered or structured image velocimetry data on a complex domain; (ii) using radial basis functions (RBFs) with divergence/curl-free kernels to provide divergence-free correction to the velocity fields for incompressible flows and curl-free correction for pressure gradients; and (iii) enforcing divergence/curl-free constraints without using Lagrangian multipliers. Complete elimination of divergence-free bias in measured pressure gradient and curl-free bias in the measured velocity field results in superior accuracy. Synthetic velocimetry data based on exact solutions and high-fidelity simulations are used to validate the analysis as well as demonstrate the flexibility and effectiveness of the RBF-HHD solver.

Approximating differential operators defined on two-dimensional surfaces is an important problem that arises in many areas of science and engineering. Over the past ten years, localized meshfree methods based on generalized moving least squares (GMLS) and radial basis function finite differences (RBF-FD) have been shown to be effective for this task as they can give high orders of accuracy at low computational cost, and they can be applied to surfaces defined only by point clouds. However, there have yet to be any studies that perform a direct comparison of these methods for approximating surface differential operators (SDOs). The first purpose of this work is to fill that gap. For this comparison, we focus on an RBF-FD method based on polyharmonic spline kernels and polynomials (PHS+Poly) since they are most closely related to the GMLS method. Additionally, we use a relatively new technique for approximating SDOs with RBF-FD called the tangent plane method since it is simpler than previous techniques and natural to use with PHS+Poly RBF-FD. The second purpose of this work is to relate the tangent plane formulation of SDOs to the local coordinate formulation used in GMLS and to show that they are equivalent when the tangent space to the surface is known exactly. The final purpose is to use ideas from the GMLS SDO formulation to derive a new RBF-FD method for approximating the tangent space for a point cloud surface when it is unknown. For the numerical comparisons of the methods, we examine their convergence rates for approximating the surface gradient, divergence, and Laplacian as the point clouds are refined for various parameter choices. We also compare their efficiency in terms of accuracy per computational cost, both when including and excluding setup costs.

We develop a new meshfree geometric multilevel (MGM) method for solving linear systems that arise from discretizing elliptic PDEs on surfaces represented by point clouds. The method uses a Poisson disk sampling-type technique for coarsening the point clouds and new meshfree restriction/interpolation operators based on polyharmonic splines for transferring information between the coarsened point clouds. These are then combined with standard smoothing and operator coarsening methods in a V-cycle iteration. MGM is applicable to discretizations of elliptic PDEs based on various localized meshfree methods, including RBF finite differences (RBF-FD) and generalized finite differences (GFD). We test MGM both as a standalone solver and preconditioner for Krylov subspace methods on several test problems using RBF-FD and GFD, and numerically analyze convergence rates, efficiency, and scaling with increasing point cloud sizes. We also perform a side-by-side comparison to algebraic multigrid (AMG) methods for solving the same systems. Finally, we further demonstrate the effectiveness of MGM by applying it to three challenging applications on complicated surfaces: pattern formation, surface harmonics, and geodesic distance.

Surface reconstruction from a set of scattered points, or a point cloud, has many applications ranging from computer graphics to remote sensing. We present a new method for this task that produces an implicit surface (zero-level set) approximation for an oriented point cloud using only information about (approximate) normals to the surface. The technique exploits the fundamental result from vector calculus that the normals to an implicit surface are curl-free. By using a curl-free radial basis function (RBF) interpolation of the normals, we can extract a potential for the vector field whose zero-level surface approximates the point cloud. We use curl-free RBFs based on polyharmonic splines for this task, since they are free of any shape or support parameters. Furthermore, to make this technique efficient and able to better represent local sharp features, we combine it with a partition of unity (PU) method. The result is the curl-free partition of unity (CFPU) method. We show how CFPU can be adapted to enforce exact interpolation of a point cloud and can be regularized to handle noise in both the normal vectors and the point positions. Numerical results are presented that demonstrate how the method converges for a known surface as the sampling density increases, how regularization handles noisy data, and how the method performs on various problems found in the literature.

Divergence-free (div-free) and curl-free vector fields are pervasive in many areas of science and engineering, from fluid dynamics to electromagnetism. A common problem that arises in applications is that of constructing smooth approximants to these vector fields and/or their potentials based only on discrete samples. Additionally, it is often necessary that the vector approximants preserve the div-free or curl-free properties of the field to maintain certain physical constraints. Div/curl-free radial basis functions (RBFs) are a particularly good choice for this application as they are meshfree and analytically satisfy the div-free or curl-free property. However, this method can be computationally expensive due to its global nature. In this paper, we develop a technique for bypassing this issue that combines div/curl-free RBFs in a partition of unity framework, where one solves for local approximants over subsets of the global samples and then blends them together to form a div-free or curl-free global approximant. The method is applicable to div/curl-free vector fields in $\R^2$ and tangential fields on two-dimensional surfaces, such as the sphere, and the curl-free method can be generalized to vector fields in $\R^d$. The method also produces an approximant for the scalar potential of the underlying sampled field. We present error estimates and demonstrate the effectiveness of the method on several test problems.

We present a high-order radial basis function finite difference (RBF-FD) framework for the solution of advection-diffusion equations on time-varying domains. Our framework is based on a generalization of the recently developed Overlapped RBF-FD method that utilizes a novel automatic procedure for computing RBF-FD weights on stencils in variable-sized regions around stencil centers. This procedure eliminates the overlap parameter δ, thereby enabling tuning-free assembly of RBF-FD differentiation matrices on moving domains. In addition, our framework utilizes a simple and efficient procedure for updating differentiation matrices on moving domains tiled by node sets of time-varying cardinality. Finally, advection-diffusion in time-varying domains is handled through a combination of rapid node set modification, a new high-order semi-Lagrangian method that utilizes the new tuning-free overlapped RBF-FD method, and a high-order time-integration method. The resulting framework has no tuning parameters and has O(N logN) time complexity. We demonstrate high-orders of convergence for advection-diffusion equations on time-varying 2D and 3D domains for both small and large Peclet numbers. We also present timings that verify our complexity estimates. Finally, we utilize our method to solve a coupled 3D problem motivated by models of platelet aggregation and coagulation, once again demonstrating high-order convergence rates on a moving domain.

The direct method used for calculating smooth radial basis function (RBF) interpolants in the flat limit becomes numerically unstable. The RBF-QR algorithm bypasses this ill-conditioning using a clever change of basis technique. We extend this method for computing interpolants involving matrix-valued kernels, specifically surface divergence-free RBFs on the sphere, in the flat limit. Results illustrating the effectiveness of this algorithm are presented for a divergence-free vector field on the sphere from samples at scattered points.

The Hierarchical Equal Area isoLatitude Pixelation (HEALPix) scheme is used extensively in astrophysics for data collection and analysis on the sphere. The scheme was originally designed for studying the Cosmic Microwave Background (CMB) radiation, which represents the first light to travel during the early stages of the universe's development and gives the strongest evidence for the Big Bang theory to date. Refined analysis of the CMB angular power spectrum can lead to revolutionary developments in understanding the nature of dark matter and dark energy. In this paper, we present a new method for performing spherical harmonic analysis for HEALPix data, which is a central component to computing and analyzing the angular power spectrum of the massive CMB data sets. The method uses a novel combination of a non-uniform fast Fourier transform, the double Fourier sphere method, and Slevinsky's fast spherical harmonic transform [38]. For a HEALPix grid with N pixels (points), the computational complexity of the method is O(Nlog2⁡N), with an initial set-up cost of O(N3/2log⁡N). This compares favorably with O(N3/2) runtime complexity of the current methods available in the HEALPix software when multiple maps need to be analyzed at the same time. Using numerical experiments, we demonstrate that the new method also appears to provide better accuracy over the entire angular power spectrum of synthetic data when compared to the current methods, with a convergence rate at least two times higher.

We present a new hyperviscosity formulation for stabilizing radial basis function-finite difference (RBF-FD) discretizations of advection-diffusion-reaction equations on manifolds $\mathbb{M} \subset \mathbb{R}^3$ of codimension 1. Our technique involves automatic addition of artificial hyperviscosity to damp out spurious modes in the differentiation matrices corresponding to surface gradients, in the process overcoming a technical limitation of a recently developed Euclidean formulation. Like the Euclidean formulation, the manifold formulation relies on von Neumann stability analysis performed on auxiliary differential operators that mimic the spurious solution growth induced by RBF-FD differentiation matrices. We demonstrate high-order convergence rates on problems involving surface advection and surface advection-diffusion. Finally, we demonstrate the applicability of our formulation to advection-diffusion-reaction equations on manifolds described purely as point clouds. Our surface discretizations use the recently developed RBF-least orthogonal interpolation method and, with the addition of hyperviscosity, are now empirically high-order accurate, stable, and free of stagnation errors.

The Cosmic Microwave Background Radiation (CMBR) represents the first light to travel during the early stages of the universe's development. This sphere of relic radiation gives the strongest evidence for the Big Bang theory to date, and refined analysis of its angular power spectrum can lead to revolutionary developments in understanding the nature of dark matter and dark energy. Satellites collect CMBR data over a sphere using a Hierarchical Equal Area isoLatitude Pixelation (HEALPix) grid. While this grid gives a quasiuniform discretization of a sphere, it is not well suited for doing fast \emph{and} accurate spherical harmonic analysis -- a central component to computing and analyzing the angular power spectrum of the massive CMBR data sets. In this paper, we present a new method that overcomes these issues through a novel combination of a non-uniform fast Fourier transform, the double Fourier sphere method, and Slevinsky's fast spherical harmonic transform (Slevinsky, 2017). The method has a quasi-optimal computational complexity of $\mathcal{O}(N\log^2 N)$ with an initial set-up cost of $\mathcal{O}(N^{3/2}\log N)$, where $N$ represents the number of points in the HEALPix grid. Additionally, we provide the first analysis of the method used in the current HEALPix software for computing the spherical harmonic coefficients. Numerical results illustrating the effectiveness of the new technique over the current method are also included.

We present three new semi-Lagrangian methods based on radial basis function (RBF) interpolation for numerically simulating transport on a sphere. The methods are mesh-free and are formulated entirely in Cartesian coordinates, thus avoiding any irregular clustering of nodes at artificial boundaries on the sphere and naturally bypassing any apparent artificial singularities associated with surface-based coordinate systems. For problems involving tracer transport in a given velocity field, the semi-Lagrangian framework allows these new methods to avoid the use of any stabilization terms (such as hyperviscosity) during time-integration, thus reducing the number of parameters that have to be tuned. The three new methods are based on interpolation using 1) global RBFs, 2) local RBF stencils, and 3) RBF partition of unity. For the latter two of these methods, we find that it is crucial to include some low degree spherical harmonics in the interpolants. Standard test cases consisting of solid body rotation and deformational flow are used to compare and contrast the methods in terms of their accuracy, efficiency, conservation properties, and dissipation/dispersion errors. For global RBFs, spectral spatial convergence is observed for smooth solutions on quasi-uniform nodes, while high-order accuracy is observed for the local RBF stencil and partition of unity approaches.

A radial basis function (RBF) method based on matrix-valued kernels is presented and analyzed for computing two types of vector decompositions on bounded domains: one where the normal component of the divergence-free part of the field is specified on the boundary, and one where the tangential component of the curl-free part of the field specified. These two decompositions can then be combined to obtain a full Helmholtz-Hodge decomposition of the field, i.e. the sum of divergence-free, curl-free, and harmonic fields. All decompositions are computed from samples of the field at (possibly scattered) nodes over the domain, and all boundary conditions are imposed on the vector fields, not their potentials, distinguishing this technique from many current methods. Sobolev-type error estimates for the various decompositions are provided and demonstrated with numerical examples.

One commonly finds in applications of smooth radial basis functions (RBFs) that scaling the kernels so they are 'flat' leads to smaller discretization errors. However, the direct numerical approach for computing with flat RBFs (RBF-Direct) is severely ill-conditioned. We present an algorithm for bypassing this ill-conditioning that is based on a new method for rational approximation (RA) of vector-valued analytic functions with the property that all components of the vector share the same singularities. This new algorithm (RBF-RA) is more accurate, robust, and easier to implement than the Contour-Padé method, which is similarly based on vector-valued rational approximation. In contrast to the stable RBF-QR and RBF-GA algorithms, which are based on finding a better conditioned base in the same RBF-space, the new algorithm can be used with any type of smooth radial kernel, and it is also applicable to a wider range of tasks (including calculating Hermite type implicit RBF-FD stencils). We present a series of numerical experiments demonstrating the effectiveness of this new method for computing RBF interpolants in the flat regime. We also demonstrate the flexibility of the method by using it to compute implicit RBF-FD formulas in the flat regime and then using these for solving Poisson's equation in a 3-D spherical shell.