Source-linked AI summary
A Radial Basis Function (RBF)-Finite Difference (FD) Method for Diffusion and Reaction-Diffusion Equations on Surfaces
Varun Shankar, Grady B. Wright, Robert M. Kirby, Aaron L. Fogelson
TL;DR
Diffusion and reaction-diffusion equations on embedded surfaces require methods that handle scattered nodes without intrinsic-coordinate distortions or the cost of solving in the embedding space. The paper introduces a local RBF-FD method with shape-parameter optimization for stable surface operators, and demonstrates convergence and nonlinear-PDE applications on varied surfaces.
Problem
Existing surface methods can face intrinsic-coordinate singularities or distortions, while global RBF and embedded approaches can impose substantial computational costs.
Method
The paper uses local RBF-FD discretization in extrinsic coordinates with scattered surface nodes and stencil-dependent shape parameters selected through a global target condition number.
Results
The method provides a stable surface Laplace-Beltrami discretization without artificial-viscosity stabilization and is demonstrated through convergence studies and nonlinear PDE simulations on parametrically defined and point-cloud surfaces.
Takeaways & Limitations
Shape-parameter optimization can compensate for irregular surface sampling while combining the local sparsity of RBF-FD with stable surface PDE simulations.
Takeaways & Limitations
The method targets smooth boundaryless codimension-one manifolds, primarily surfaces embedded in R^3, and shape-parameter optimization adds root-finding cost.
Abstract
from arXiv · showhide
In this paper, we present a method based on Radial Basis Function (RBF)-generated Finite Differences (FD) for numerically solving diffusion and reaction-diffusion equations (PDEs) on closed surfaces embedded in $\mathbb{R}^d$. Our method uses a method-of-lines formulation, in which surface derivatives that appear in the PDEs are approximated locally using RBF interpolation. The method requires only scattered nodes representing the surface and normal vectors at those scattered nodes. All computations use only extrinsic coordinates, thereby avoiding coordinate distortions and singularities. We also present an optimization procedure that allows for the stabilization of the discrete differential operators generated by our RBF-FD method by selecting shape parameters for each stencil that correspond to a global target condition number. We show the convergence of our method on two surfaces for different stencil sizes, and present applications to nonlinear PDEs simulated both on implicit/parametric surfaces and more general surfaces represented by point clouds.
1 Introduction
The paper motivates a local RBF-FD approach for diffusion and reaction-diffusion equations on general surfaces, avoiding intrinsic-coordinate distortions and the high costs of dense global or embedded methods. It modifies an earlier surface RBF-FD formulation to address instability without hyperviscosity.
- Intrinsic methods can be efficient but may suffer from coordinate singularities or distortions, while embedded narrow-band methods incur costs that grow with embedding dimension.
- Global RBF surface discretizations achieve high convergence rates but require O(N^3) setup and O(N^2) matrix application for N surface nodes.
- RBF-FD uses local RBF interpolation to produce sparse differentiation matrices that handle irregular geometries and scattered node layouts.
- A direct extension of an earlier 1D-surface RBF-FD method to 2D surfaces was unstable and required hyperviscosity-based stabilization.
- The proposed formulation develops numerical and algorithmic strategies for general surfaces and appears to avoid hyperviscosity stabilization.
- The paper focuses on smooth boundaryless codimension-one manifolds, primarily surfaces embedded in R^3.
2 A review of RBF interpolation
RBF interpolation represents sampled functions with radial kernels and can be applied componentwise to vector data. The paper uses Euclidean distances on surface nodes and highlights the accuracy–conditioning trade-off controlled by the kernel shape parameter.
- RBF interpolation constructs a function from sampled values at scattered nodes using radial kernels based on Euclidean distances.
- The interpolation coefficients are obtained from a linear system whose matrix is invertible for distinct nodes under positive-definite or suitable conditionally positive-definite kernels.
- The method uses straight-line Euclidean distances in the embedding space for nodes lying on a lower-dimensional surface, avoiding explicit surface information during interpolation.
- Vector-valued functions are interpolated componentwise, allowing one factored interpolation matrix to support multiple components.
- The paper uses the positive-definite inverse multiquadric kernel for RBF-FD computations.
- Smaller shape parameters generally improve accuracy for smooth targets but make standard RBF-Direct interpolation ill-conditioned.
3 Surface Laplacian in Cartesian coordinates
The surface Laplacian is expressed using Cartesian operators and tangent-plane projection, avoiding singularities associated with intrinsic surface coordinates. The proposed discretization mirrors this formulation with local RBF interpolants.
- Cartesian formulation avoids singularities associated with intrinsic surface coordinates such as spherical-coordinate pole singularities.
- The projection operator maps vectors in R^3 onto the tangent plane using the surface unit normal.
- The surface gradient is formed by combining tangent-plane projection with the standard Cartesian gradient.
- The surface Laplacian is the surface divergence of the surface gradient and is written as GxGx + GyGy + GzGz.
- The proposed approximation follows the Cartesian Laplacian formulation while replacing global RBF interpolation with local RBF interpolants.
4 RBF-FD approximation to the surface Laplacian
The method approximates the surface Laplacian locally by repeatedly interpolating and differentiating RBFs on nearest-neighbor stencils, then assembles the resulting weights into a sparse differentiation matrix for method-of-lines simulations.
- Local stencils: Each surface node uses itself and its n−1 nearest Euclidean neighbors as a local stencil for approximating the surface Laplacian.The stencil is selected in the embedding space R^3.
- RBF-FD construction: RBF-FD weights are obtained by interpolating f on a stencil, computing surface-gradient components, interpolating them, and applying surface divergence.This repeated interpolation-and-differentiation procedure avoids analytically differentiating surface normal vectors.
- Sparse operator assembly: The center row of each local operator supplies the weights for its stencil center, and all center rows form the sparse matrix LX.The procedure is repeated over all N nodes.
- Computational cost: Building each stencil operator costs O(n^3), so with fixed n the total weight-construction cost scales as O(N) and is embarrassingly parallel.Including k-d-tree construction and lookup gives an overall building cost of O(N log N) + O(N).
- Implementation limitations: The approach uses Euclidean distances in the embedding space, requiring caution near thin surface features where cross-feature neighbors can corrupt the Laplacian approximation.The authors address this only by using sufficiently dense sampling around thin features.
- Method of lines: The sparse RBF-FD Laplacian enters a method-of-lines system of N coupled ODEs for diffusion and reaction-diffusion equations.Explicit time stepping can evaluate LX in O(N) operations.
5 Shape Parameter and Eigenvalue Stability
The method stabilizes RBF-FD differentiation matrices by selecting a shape parameter per stencil to target a common interpolation condition number. Tests show that this optimization places eigenvalues in the left half plane where a single mean shape parameter can produce a large positive eigenvalue.
- Stability requirement: Stability requires the eigenvalues of LX to lie in the ODE solver’s stability domain, generally at least in the left half plane.The basic RBF-FD procedure does not guarantee this property.
- Per-stencil optimization: Each stencil receives a shape parameter εk chosen to induce a target condition number κT for its RBF interpolation matrix.Because the condition number increases as ε decreases toward zero, εk is found from the zero of a scalar condition-number function.
- Irregular nodes: The optimization acts as a homogenization that compensates for irregular node distributions by equalizing stencil conditioning.This generalizes a sphere-based strategy in which quasi-uniform nodes allow one shape parameter to produce similar conditioning.
- Red Blood Cell test: For Red Blood Cell nodes, per-stencil optimization produced only left-half-plane eigenvalues, whereas the single mean ε produced one large positive eigenvalue.The comparison used N = 10000 mapped nodes and n = 31 stencils.
- Double-torus test: For the irregular double-torus nodes, optimized shape parameters again produced left-half-plane eigenvalues, while the single-ε alternative did not.The optimized test used N = 5041 nodes, n = 31, and κT = 10^11; the single parameter was ε = 2.47.
- Trade-offs: Shape optimization adds root-finding cost and can sacrifice local accuracy when a common target condition number is below a stencil’s naturally preferred value.The reported tests found no resulting impact on convergence.
6 Convergence studies
The convergence studies evaluate the RBF-FD method on sphere and torus diffusion problems while varying node count, stencil size, and condition-number strategies. Results show increasing convergence order with larger stencils, while fixed target condition numbers can cause saturation for high-order methods.
- Experimental setup: The studies vary stencil size n and total surface nodes N, using BDF4 with spatial errors intended to dominate temporal errors.Errors are measured in ℓ2 and ℓ∞ norms, with BICGSTAB solving the implicit systems.
- Increasing condition number: On the sphere, convergence ranges from orders two to three for n = 11 to slightly above order five for n = 31 in both norms.The sphere diffusion test uses icosahedral node sets and compares against exact solutions.
- Numerical conditioning: Allowing condition numbers to grow with N can create nearly singular or singular matrices in double precision, so the largest sphere case required quad-precision weights.The n = 31, N = 40962 weights used quad precision, although simulations using them ran in double precision.
- Increasing condition number: On the torus, convergence is order two for n = 11, close to order four for n = 17, and between orders five and six for n = 31 in both norms.The torus results are comparable to sphere diffusion and slightly better than forced diffusion on the sphere.
- Fixed condition number: For forced diffusion on the sphere, n = 11, 17, and 31 produce increasing convergence orders, with high-order cases eventually requiring larger target condition numbers.At κT = 10^14, n = 31 saturates for N > 2562, indicating that small shape parameters are needed for high-order RBF-FD methods.
- Fixed condition number: At κT = 10^20 on the sphere, n = 11 and n = 17 show no saturation, while n = 31 converges close to order six without saturation over the tested N values.The results support using smaller shape parameters to obtain higher convergence.
7 Application: Turing patterns
The RBF-FD method is applied to two-species Turing reaction-diffusion systems on parametrically defined, implicit, and point-cloud surfaces. Simulations produce qualitatively consistent spot and stripe patterns across several geometries.
- System and setup: The Turing system models coupled activator u and inhibitor v equations whose patterns depend on diffusion and coupling parameters.τ1 favors stripes, while τ2 favors spots; stripe patterns require longer to reach steady state.
- System and setup: The experiments use SBDF2 time stepping with Δt = 0.01 and RBF-FD discretization parameters reported for each surface.The diffusion terms are handled implicitly.
- Turing patterns on known surfaces: RBC, double-torus, and Bumpy Sphere simulations generate qualitatively similar spot and stripe patterns to previous results.The Bumpy Sphere experiment uses a sampled point set with N = 10000 vertices and interpolated unit normals.
- Turing patterns on point-cloud surfaces: Frog and Stanford Bunny models are preprocessed with Poisson surface reconstruction and Poisson disk sampling before simulation.The preprocessing smoothes approximate normals and produces point clouds with desired sampling density.
- Turing patterns on point-cloud surfaces: The Frog produces a Turing spot simulation and the Stanford Bunny a Turing stripe simulation, with results qualitatively similar to those on the earlier surfaces.Color maps encode concentration levels for visualizing the patterns.
8 Discussion
The paper presents a stable RBF-FD discretization for the Laplace-Beltrami operator on embedded surfaces, using optimized local shape parameters without artificial-viscosity stabilization. The discussion identifies thin features, evolving surfaces, and variable diffusion as important boundaries or extensions.
- Contributions: The method approximates the Laplace-Beltrami operator on codimension-one surfaces in R3 using scattered nodes without expanding into the embedding space.It is designed for surface PDE computations using only the surface point set and associated geometric data.
- Contributions: An optimization algorithm selects shape parameters for each interpolation patch and compensates for irregular surface sampling.The paper relates condition number, shape parameter, and convergence order through error and convergence studies.
- Limitations and extensions: The current method is limited to static objects; evolving-surface simulations would require dynamic k-d trees and likely parallelization for efficiency.The intended evolving-surface setting involves fluid-object interaction.
- Limitations and extensions: The method is not innately robust to thin surface features, motivating adaptive refinement or efficient surface-distance measurement.Robust sampling must sufficiently resolve both sides of a thin feature.
- Limitations and extensions: Spatially variable or anisotropic diffusion is identified as a future extension for applications such as gels and viscoelastic materials on surfaces.The authors describe this extension as conceptually straightforward but defer its realization to future work.