Grady B. Wright

Grady B. Wright, Bengt Fornberg

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 3D spherical shell.

Edward J. Fuselier, Grady B. Wright

In this paper we present a high-order kernel method for numerically solving diffusion and reaction-diffusion partial differential equations (PDEs) on smooth, closed surfaces embedded in $\mathbb{R}^d$. For two-dimensional surfaces embedded in $\mathbb{R}^3$, these types of problems have received growing interest in biology, chemistry, and computer graphics to model such things as diffusion of chemicals on biological cells or membranes, pattern formations in biology, nonlinear chemical oscillators in excitable media, and texture mappings. Our kernel method is based on radial basis functions (RBFs) and uses a semi-discrete approach (or the method-of-lines) in which the surface derivative operators that appear in the PDEs are approximated using collocation. The method only requires nodes at "scattered" locations on the surface and the corresponding normal vectors to the surface. Additionally, it does not rely on any surface-based metrics and avoids any intrinsic coordinate systems, and thus does not suffer from any coordinate distortions or singularities. We provide error estimates for the kernel-based approximate surface derivative operators and numerically study the accuracy and stability of the method. Applications to different non-linear systems of PDEs that arise in biology and chemistry are also presented.

Alex Townsend, Heather Wilber, Grady B. Wright

A collection of algorithms is described for numerically computing with smooth functions defined on the unit sphere. Functions are approximated to essentially machine precision by using a structure-preserving iterative variant of Gaussian elimination together with the double Fourier sphere method. We show that this procedure allows for stable differentiation, reduces the oversampling of functions near the poles, and converges for certain analytic functions. Operations such as function evaluation, differentiation, and integration are particularly efficient and can be computed by essentially one-dimensional algorithms. A highlight is an optimal complexity direct solver for Poisson's equation on the sphere using a spectral method. Without parallelization, we solve Poisson's equation with $100$ million degrees of freedom in one minute on a standard laptop. Numerical results are presented throughout. In a companion paper (part II) we extend the ideas presented here to computing with functions on the disk.

Edward J. Fuselier, Grady B. Wright

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.

Michael Chiwere, Daniel Fortunato, Grady B. Wright

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.

Grady B. Wright, Mohsin Javed, Hadrien Montanelli et al.

Algorithms and underlying mathematics are presented for numerical computation with periodic functions via approximations to machine precision by trigonometric polynomials, including the solution of linear and nonlinear periodic ordinary differential equations. Differences from the nonperiodic Chebyshev case are highlighted.