Abstract
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 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.
Similar content being viewed by others
Notes
The condition \(\tau > 3/2\) ensures that functions within the kernel’s native space are continuous on \(\mathbb{R }^3\).
See the discussion preceeding Proposition 2 in “Appendix A” for details.
See Theorem 2 for details.
With this integral operator comes the pseudodifferential operators \(T^{-r}, \,r>0\). A function \(f\) is in the native space if and only if \(T^{-1/2}f \in L_2(\mathbb{M })\) [35, Proposition 4.9], so we expect functions such that \(T^{-1}f\in L_2(\mathbb{M })\) to be twice as smooth.
References
3DXM Consortium: 3D-XplorMath. Website. http://3D-XplorMath.org. Accessed 27 July 2011
Adalsteinsson, D., Sethian, J.A.: Transport and diffusion of material quantities on propagating interfaces via level set methods. J. Comput. Phys. 185(1), 271–288 (2003). doi:10.1016/S0021-9991(02)00057-8
AIM@SHAPE: “Bumpy sphere” from AIM@SHAPE shape repository. Website. http://shapes.aimatshape.net. (2011)
Ascher, U.M., Ruuth, S.J., Wetton, B.T.R.: Implicit-explicit methods for time-dependent partial differential equations. SIAM J. Numer. Anal. 32, 797–823 (1995)
Barkley, D.: A model for fast computer simulation of waves in excitable media. Phys. D 49, 61–70 (1991)
Barkley, D.: Barkley model. Scholarpedia 3(10), 1877 (2008)
Barrio, R., Varea, C., Aragón, J., Maini, P.: A two-dimensional numerical study of spatial pattern formation in interacting Turing systems. Bull. Math. Biol. 61, 483–505 (1999)
Bergdorf, M., Sbalzarini, I., Koumoutsakos, P.: A Lagrangian particle method for reaction-diffusion systems on deforming surfaces. J. Math. Biol. 61, 649–663 (2010). doi:10.1007/s00285-009-0315-2
Bertalmío, M., Cheng, L., Osher, S., Sapiro, G.: Variational problems and partial differential equations on implicit surfaces. J. Comput. Phys. 174, 759–780 (2001)
Bollig, E.F., Flyer, N., Erlebacher, G.: Solution to PDEs using radial basis function finite-differences (RBF-FD) on multiple GPUs. J. Comput. Phys. 231(21), 7133–7151 (2012)
Calhoun, D.A., Helzel, C.: A finite volume method for solving parabolic equations on logically Cartesian curved surface meshes. SIAM J. Sci. Comput. 31(6), 4066–4099 (2009). doi:10.1137/08073322X
Chaplain, M.A., Ganesh, M., Graham, I.G.: Spatio-temporal pattern formation on spherical surfaces: numerical simulation and application to solid tumour growth. J. Math. Biol. 42, 387–423 (2001). doi:10.1007/s002850000067
Davydov, V.A., Manz, N., Steinbock, O., Zykov, V.S., Muller, S.C.: Excitation fronts on a periodically modulated curved surface. Phys. Rev. Lett. 85, 868–871 (2000). doi:10.1103/PhysRevLett.85.868
Deckelnick, K., Dziuk, G., Elliott, C.M., Heine, C.: An $h$-narrow band finite-element method for elliptic equations on implicit surfaces. IMA J. Numer. Anal. 30(2), 351–376 (2010). doi:10.1093/imanum/drn049
Dziuk, G.: Finite elements for the Beltrami operator on arbitrary surfaces. In: Hildebrandt, S., Leis, R. (eds.) Partial Differential Equations and Calculus of Variations, Lecture Notes in Mathematics, vol. 1357. Springer, Berlin (1988)
Dziuk, G., Elliot, C.: Finite elements on evolving surfaces. IMA J. Numer. Anal. 27, 262–292 (2007)
Dziuk, G., Elliot, C.M.: Surface finite elements for parabolic equations. J. Comput. Math. 25, 385–407 (2007)
Epstein, I.R., Pojman, J.A.: An Introduction to Nonlinear Chemical Dynamics. Oxford University Press, Oxford (1998)
Evans, E., Fung, Y.C.: Improved measurements of the erythrocyte geometry. Microvasc. Res. 4(4), 335–347 (1972). doi:10.1016/0026-2862(72)90069-6
Fasshauer, G.E.: Meshless Approximation Methods with MATLAB. Interdisciplinary Mathematical Sciences, vol. 6. World Scientific Publishers, Singapore (2007)
Fasshauer, G.E., McCourt, M.J.: Stable evaluation of Gaussian radial basis function interpolants. SIAM J. Sci. Comput. 34(2), A737–A762 (2012). doi:10.1137/110824784
FitzHugh, R.: Fitzhugh-nagumo simplified cardiac action potential model. Biophys. J. 1, 445–466 (1961)
Flyer, N., Lehto, E.: Rotational transport on a sphere: Local node refinement with radial basis functions. J. Comput. Phys. 229, 1954–1969 (2010)
Flyer, N., Lehto, E., Blaise, S., Wright, G.B., St-Cyr, A.: A guide to RBF-generated finite differences for nonlinear transport: shallow water simulations on a sphere. J. Comput. Phys. 231, 4078–4095 (2012)
Flyer, N., Wright, G.B.: Transport schemes on a sphere using radial basis functions. J. Comput. Phys. 226, 1059–1084 (2007)
Flyer, N., Wright, G.B.: A radial basis function method for the shallow water equations on a sphere. Proc. R. Soc. A 465, 1949–1976 (2009)
Fornberg, B., Larsson, E., Flyer, N.: Stable computations with Gaussian radial functions. SIAM J. Sci. Comput. 33, 869–892 (2011)
Fornberg, B., Lehto, E.: Stabilization of RBF-generated finite difference methods for convective PDEs. J. Comput. Phys. 230, 2270–2285 (2011)
Fornberg, B., Piret, C.: A stable algorithm for flat radial basis functions on a sphere. SIAM J. Sci. Comput. 200, 178–192 (2007)
Fornberg, B., Piret, C.: On choosing a radial basis function and a shape parameter when solving a convective PDE on a sphere. J. Comput. Phys. 227, 2758–2780 (2008)
Fornberg, B., Wright, G.: Stable computation of multiquadric interpolants for all values of the shape parameter. Comput. Math. Appl. 48, 853–867 (2004)
Fornberg, B., Wright, G., Larsson, E.: Some observations regarding interpolants in the limit of flat radial basis functions. Comput. Math. Appl. 47, 37–55 (2004)
Fuselier, E., Narcowich, F., Ward, J., Wright, G.: Error and stability estimates for divergence-free RBF interpolants on the sphere. Math. Comput. 78, 2157–2186 (2009)
Fuselier, E.J., Wright, G.B.: Stability and error estimates for vector field interpolation and decomposition on the sphere with rbfs. SIAM J. Numer. Anal. 47(5), 3213–3239 (2009). doi:10.1137/080730901
Fuselier, E.J., Wright, G.B.: Scattered data interpolation on embedded submanifolds with restricted positive definite Kernels: Sobolev error estimates. SIAM J. Numer. Anal. 50(3), 1753–1776 (2012). doi:10.1137/110821846
Fuselier, E.J., Wright, G.B.: Order preserving approximations of derivatives with periodic radial basis functions (2013). Submitted.
Gia, Q.T.L.: Approximation of parabolic PDEs on spheres using spherical basis functions. Adv. Comput. Math. 22, 377–397 (2005)
Gomatam, J., Amdjadi, F.: Reaction-diffusion equations on a sphere: meandering of spiral waves. Phys. Rev. E 56, 3913–3919 (1997)
Greer, J.B.: An improvement of a recent Eulerian method for solving PDEs on general geometries. J. Sci. Comput. 29, 321–352 (2006)
Grindrod, P., Lewis, M.A., Murray, J.D.: A geometrical approach to wave-type solutions of excitable reaction-diffusion systems. Proc. R. Soc. A 433(1887), 151–164 (1991). doi:10.1098/rspa.1991.0040
Hardin, D.P., Saff, E.B.: Discretizing manifolds via minimum energy points. Notices Am. Math. Soc. 51, 1186–1194 (2004)
Kansa, E.J.: Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—ii: solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput. Math. Appl. 19, 147–161 (1990)
Keener, J., Sneyd, J.: Mathematical Physiology. Springer, New York (1998)
Kondo, S., Miura, T.: Reaction-diffusion model as a framework for understanding biological pattern formation. Science 329(5999), 1616–1620 (2010). doi:10.1126/science.1179047
Landsberg, C., Voigt, A.: A multigrid finite element method for reaction-diffusion systems on surfaces. Comput. Visual Sci. 13, 177–185 (2010)
Larsson, E., Lehto, E., Heryudono, A., Fornberg, B.: Stable computation of differentiation matrices and scattered node stencils based on Gaussian radial basis functions. Report No. 2012–020, Department of Information Technology, Uppsala University, Uppsala, 2012
Liu, D., Xu, G., Zhang, Q.: A discrete scheme of Laplace–Beltrami operator and its convergence over quadrilateral meshes. Comput. Math. Appl. 55, 1081–1093 (2008)
Liu, Y., Liu, W.K.: Rheology of red blood cell aggregation by computer simulation. J. Comput. Phys. 220, 139–154 (2006)
MacDonald, C.B., Ruuth, S.J.: Level set equations on surfaces via the closest point method. J. Sci. Comput. 35, 219–240 (2008)
MacDonald, C.B., Ruuth, S.J.: The implicit closest point method for the numerical solution of partial differential equations on surfaces. SIAM J. Sci. Comput. 31, 4330–4350 (2009)
Maini, P.K., Othmer, H.G. (eds.): Mathematical Models for Biological Pattern Formation. Springer, New York (2000)
Manz, N., Davydov, V.A., Müller, S.C., Bär, M.: Dependence of the spiral rotation frequency on the surface curvature of reaction-diffusion systems. Phys. Lett. A 316(5), 311–316 (2003). doi:10.1016/S0375-9601(03)01148-4
Murray, J.: Mathematical Biology II: Spatial Models and Biomedical Applications. Springer, New York (2003)
Narcowich, F.J., Ward, J.D., Wendland, H.: Sobolev bounds on functions with scattered zeros, with applications to radial basis function surface fitting. Math. Comput. 74(250), 743–763 (2005). (electronic)
Narcowich, F.J., Ward, J.D., Wright, G.B.: Divergence-free RBFs on surfaces. J. Fourier Anal. Appl. 13(6), 643–663 (2007)
Palais, R., Palais, B., Karcher, H.: Point clouds: distributing points unifomly on a surface. Technical report (2011)
Piret, C.: The orthogonal gradients method: a radial basis functions method for solving partial differential equations on arbitrary surfaces. J. Comput. Phys. 231(20), 4662–4675 (2012)
Ruuth, S.J., Merriman, B.: A simple embedding method for solving partial differential equations on surfaces. J. Comput. Phys. 227, 1943–1961 (2008)
Sbalzarini, I.F., Hayer, A., Helenius, A., Koumoutsakos, P.: Simulations of (an)isotropic diffusion on curved biological surfaces. Biophys. J. 90, 878–885 (2006)
Schwartz, P., Adalsteinsson, D., Colella, P., Arkin, A.P., Onsum, M.: Numerical computation of diffusion on a surface. Proc. Nat. Acad. Sci. 102(32), 11151–111156 (2005)
Shankar, V., Wright, G.B., Fogelson, A.L., Kirby, R.M.: A study of different modeling choices for simulating platelets within the immersed boundary method. Appl. Numer. Math. 63, 58–77 (2013)
Shelley, M.J., Baker, G.R.: Order-preserving approximations to successive derivatives of periodic functions by iterated splines. SIAM J. Numer. Anal. 25(6), 1442–1452 (1988). doi:10.1137/0725084
Stam, J.: Flows on surfaces of arbitrary topology. In: SIGGRAPH ’03: ACM SIGGRAPH 2003 Papers, pp. 724–731. ACM, New York, (2003). doi:10.1145/1201775.882338
Turing, A.: The chemical basis of morphogenesis. Philos. Trans. R. Soc. B 237, 37–52 (1952)
Turk, G.: Generating textures on arbitrary surfaces using reaction-diffusion. In: SIGGRAPH ’91: Proceedings of the 18st Annual Conference on Computer Graphics and Interactive Techniques, pp. 289–298. ACM, New York, (1991)
Tyson, J.J., Keener, J.P.: Singular perturbation theory of traveling waves in excitable media (a review). Physica D 32(3), 327–361 (1988). doi:10.1016/0167-2789(88)90062-0
Varea, C., Aragon, J., Barrio, R.: Turing patterns on a sphere. Phys. Rev. E 60, 4588–4592 (1999)
Wendland, H.: Scattered data approximation, Cambridge Monographs on Applied and Computational Mathematics, vol. 17. Cambridge University Press, Cambridge (2005)
Womersley, R.S., Sloan, I.H.: Interpolation and cubature on the sphere. Website. http://web.maths.unsw.edu.au/rsw/Sphere/. Accessed 17 Jan 2013
Wright, G.B., Flyer, N., Yuen, D.: A hybrid radial basis function - pseudospectral method for thermal convection in a 3D spherical shell. Geochem. Geophys. Geosyst. 11, Q07,003 (2010)
Xu, G.: Discrete Laplace–Beltrami operators and their convergence. Comput. Aided Geom. Des. 21, 767–784 (2004)
Yagisita, H., Mimura, M., Yamada, M.: Spiral wave behaviors in an excitable reaction-diffusion system on a sphere. Physica D 124(1–3), 126–136 (1998). doi:10.1016/S0167-2789(98)00182-1
Acknowledgments
Research support for G. B. Wright was provided, in part, by grants DMS-0934581, DMS-0540779, and DMS-1160379 from the National Science Foundation.
Author information
Authors and Affiliations
Corresponding author
Appendices
Appendix A: Convergence Results
In this section we present the convergence results for our discrete differential operators. The results presented here are for the special case when \(\mathbb{M }\subset \mathbb{R }^3\) is a two dimensional surface, although similar results hold in higher dimensions. The arguments will depend heavily on error estimates for kernel interpolants on smooth, embedded submanifolds of \(\mathbb{R }^d\) given recently in [35]. We will frequently reference results from that paper throughout this section. Sobolev error estimates for the surface gradient and divergence operators will immediately follow the results in [35]. However, the estimates for the discrete surface Laplacian will require more work.
Proposition 1
Let \(1\le q \le \infty , \,\phi \) be a positive definite kernel satisfying (25) with \(\tau > 3/2 + 1\), and define \(s = \tau - 1/2\). Then there is a constant \(h_{\mathbb{M }}\) depending only on \(\mathbb{M }\) such that if a finite node set \(X\subset \mathbb{M }\) satisfies \(h\le h_\mathbb{M }\), for all \(f\in \mathcal N _\phi (\mathbb{M })\) and \(\mathbf{f }\in (\mathcal N _\phi (\mathbb{M }))^3\) we have
where \((x)_+ = x\) if \(x\ge 0\) and is zero otherwise.
Proof
We will focus on the discrete gradient operator. We have
where \(W^1_q(\mathbb{M })\) is the \(L_q\) Sobolev space of order 1 (see Sect. 2, [35]). The assumptions on \(\tau \) allow us to apply Theorem 4.6 in [35], and the results follow. The error bounds in Theorem 4.6 in [35] carry over to the vector-valued case, thus bounds for the divergence operator error are obtained in a similar manner. \(\square \)
When \(\phi \) satisfies (26), i.e. \(\phi \) has finite smoothness, the estimates come in two types: the first concerns targets that may be too rough to be within the native space, and the second applies to functions with additional smoothness. This additional smoothness is measured with the inverse of the integral operatorFootnote 4
We denote a vector version of this operator by \(\mathbf{T }\), which simply applies \(T\) to each component of a \(3\)-dimensional vector field. Similar to the proof above, the result below follows from Theorem 4.12 and Corollary 4.10 in [35].
Proposition 2
Let \(1\le q \le \infty \), and let \(\phi \) be a positive definite kernel satisfying (26) with \(\tau > 3/2 + 1\), and define \(s = \tau - 1/2\). Then there is a constant \(h_{\mathbb{M }}\) depending only on \(\mathbb{M }\) such that if a finite node set \(X\subset \mathbb{M }\) satisfies \(h\le h_\mathbb{M }\), we have the following:
-
1.
Rough target functions. Let \(\beta \) be such that \(s \ge \beta > 2\). Then for all \(f \in H^\beta (\mathbb{M })\) and \(\mathbf{f }\in \mathbf{H }^\beta (\mathbb{M })\) we have
$$\begin{aligned} \Vert G_\mathbb{M }f-\nabla _{\mathbb{M }} f\Vert _{\mathbf{L }_q(\mathbb{M })}&\le C h^{\beta -1-2(1/2-1/q)_+}\rho ^{s -\beta }\Vert f\Vert _{H^{\beta }(\mathbb{M })},\\ \Vert D_{\mathbb{M }}\mathbf{f }-\nabla _{\mathbb{M }}\cdot \mathbf{f }\Vert _{L_q(\mathbb{M })}&\le C h^{\beta -1-2(1/2-1/q)_+}\rho ^{s - \beta }\Vert \mathbf{f }\Vert _{\mathbf{H }^{\beta }(\mathbb{M })}. \end{aligned}$$ -
2.
Smooth target functions. For all \(f\in L_2(\mathbb{M })\) such that \(T^{-1}f\in L_2(\mathbb{M })\) and \(\mathbf{f }\in \mathbf{L }_2(\mathbb{M })\) such that \(T^{-1}\mathbf{f }\in \mathbf{L }_2(\mathbb{M })\), we have
$$\begin{aligned} \Vert G_\mathbb{M }f-\nabla _{\mathbb{M }} f\Vert _{\mathbf{L }_q(\mathbb{M })}&\le C h^{2s-1-2(1/2 - 1/q)_+}\Vert T^{-1}f\Vert _{L_2(\mathbb{M })},\\ \Vert D_{\mathbb{M }}\mathbf{f }-\nabla _{\mathbb{M }}\cdot \mathbf{f }\Vert _{\mathbf{L }_q(\mathbb{M })}&\le C h^{2s-1-2(1/2 - 1/q)_+}\Vert \mathbf{T }^{-1}\mathbf{f }\Vert _{\mathbf{L }_2(\mathbb{M })}. \end{aligned}$$
Obtaining convergence rates for the discrete Laplace operator is not as straightforward, and will require extra tools. The first is the “zeros lemma” from Narcowich, Ward and Wendland [54]. We will use the manifold version stated in [35, Lemma 4.5], with parameters adapted to our situation.
Proposition 3
Let \(1\le q\le \infty , \,t\in \mathbb{R }\) with \(t > 1\). Let \(\mu \in \mathbb{N }\) satisfy \(0\le \mu \le \lceil t - 2(1/2 - 1/q)_{+}\rceil -1\). Also, let \(X\subset \mathbb{M }\) be a discrete set with mesh norm \(h_{X,\mathbb{M }}\). Then there is a constant depending only on \(\mathbb{M }\) such that if \(h_{X,\mathbb{M }}\le C_{\mathbb{M }}\) and if \(u\in H^{t}(\mathbb{M })\) satisfies \(u|_{X}=0\), then
Proposition 3 also holds with full Sobolev norms, and a similar result holds for vector fields.
Since the operator \(L_{\mathbb{M }}\) is a mixture of differential and interpolation operators, it will be beneficial to bound the norm of the interpolation operator and its associated error in several different Sobolev spaces. In particular, we will use the following.
Lemma 1
Let \(\phi \) satisfy (26) with \(\tau > 3/2\), and define \(s = \tau - 1/2\). Let \(\beta \) be such that \(s \ge \beta > 1\). Then there is a constant \(h_{\mathbb{M }}\) depending only on \(\mathbb{M }\) such that if a finite node set \(X\subset \mathbb{M }\) satisfies \(h\le h_\mathbb{M }\) then for all \(f\in H^{\beta }(\mathbb{M })\) we have
Proof
The last estimate in the proof of Theorem 4.12 in [35] is
Applying a triangle inequality and observing that \(1\le \rho \), we get
\(\square \)
Now we are poised to present the following theorem.
Theorem 1
Let \(1\le q \le \infty , \,\phi \) satisfy (26) with \(\tau > 3/2 + 2\), and define \(s = \tau - 1/2\). Let \(\beta \) be such that \(s \ge \beta > 3\). Then there is a constant \(h_{\mathbb{M }}\) depending only on \(\mathbb{M }\) such that if a finite node set \(X\subset \mathbb{M }\) satisfies \(h\le h_\mathbb{M }\), then for all \(f\in H^\beta (\mathbb{M })\) we have the estimate
Proof
First, we have
For the rightmost term, the assumption \(\beta > 3\) allows us to use [35, Theorem 4.12] to get:
For the other term in (33), we have
Note that the last estimate is the interpolation error for the target function \(\mathbf{g }:= \nabla _\mathbb{M }I_{\phi }f\).
We claim that \(\mathbf{g }\) is in \(\mathbf{H }^{t}(\mathbb{M })\) for all \(t < 2s - 2\). Since \(\phi \) satisfies (26), \(I_{\phi }f \in H^{\nu }(\mathbb{R }^3)\) for all \(\nu < 2\tau - 3/2\). By the trace theorem restricting to \(\mathbb{M }\) puts \(I_{\phi }f\) in \(H^{\nu - 1/2}(\mathbb{M })\). Thus \(\mathbf{g }\in \mathbf{H }^{\nu - 3/2}(\mathbb{M })\) for all \(\nu < 2\tau - 3/2\), which is equivalent to \(\mathbf{g }\in \mathbf{H }^{t}(\mathbb{M })\) for all \(t < 2s - 2\). In particular, \(\mathbf{g }\in \mathbf{H }^{\beta - 1}(\mathbb{M })\). Also, since \(\beta > 2\) then we have
and that \(\beta -1 > 1\). This allows us to estimate (34) with Proposition 3 (with parameter \(\mu = 1\) and target smoothness \(t = \beta -1\)) to get
Further, since \(\beta - 1>1\), we can also apply Lemma 1, and with two applications of it we get
This completes the proof. \(\square \)
Continuing with our error analysis, we now shift our attention to target functions that are very smooth. First we will need a lemma.
Lemma 2
Let \(\phi \) satisfy (26) with \(\tau > 3/2\), and let \(s = \tau - 1/2\). Then for all \(f\in L_2(\mathbb{M })\) such that \(T^{-1}f \in L_2(\mathbb{M })\) we have the following estimate
Also, for all \(\mathbf{f }\in \mathbf{L }_2(\mathbb{M })\) such that \(\mathbf{T }^{-1}\mathbf{f }\in \mathbf{L }_2(\mathbb{M })\) we have
Proof
The first estimate is established in the proof of Corollary 4.10 in [35]. By working component-wise the second estimate follows from the first. \(\square \)
Theorem 2
Let \(1\le q \le \infty , \,\phi \) satisfy (26) with \(\tau > 3/2 + 2\), and define \(s = \tau - 1/2\). Let \(f\in L_2(\mathbb{M })\) be such that \(T^{-1}f \in L_2(\mathbb{M })\) and \(\mathbf{T }^{-1}\nabla _{\mathbb{M }}f \in \mathbf{L }_2(\mathbb{M })\). Then there is a constant \(h_{\mathbb{M }}\) depending only on \(\mathbb{M }\) such that if a finite node set \(X\subset \mathbb{M }\) satisfies \(h\le h_\mathbb{M }\), we have the estimate
Proof
We begin as in the proof of the previous theorem. We have
To estimate the rightmost term we can use the error estimates in [35, Corollary 4.10] to get:
For the other term in (36), we proceed as in the proof of Theorem 1 to get
where \(\mathbf{g }= \nabla _\mathbb{M }I_{\phi }f\), and as before we know that this function is in \(\mathbf{H }^{t}(\mathbb{M })\) for all \(t < 2s - 2\). In particular, \(\mathbf{g }\in \mathbf{H }^{s-1}(\mathbb{M })\) and we can estimate (37) with Proposition 3 (with \(t=s-1\) and \(\mu =1\)) to get
This is where the proof detours from that of the previous theorem. To estimate \(\Vert I_{\Phi }\mathbf{g }- \mathbf{g }\Vert _{\mathbf{H }^{s-1}(\mathbb{M })}\), note that it is bounded it by following quantity:
We will bound each term individually. First we concentrate on \(I\). An application of Lemma 1, which we may apply since \(s - 1 > 1\), gives us
Thus a bound for \(III\) will result in a bound for \(I\). To bound \(III\), we can apply Lemma 2 to get
To bound \(II\), we again employ Lemma 2:
This completes the proof. \(\square \)
Appendix B: Surfaces and Node Sets from the Numerical Experiments
1.1 B.1 Unit Sphere
This manifold is, of course, described implicitly by
The node sets we use for discretizing the unit sphere are the minimum energy (ME) node sets of Womersley and Sloan [69]. These node sets are approximately uniformly distributed over the surface of the sphere and have the nice property that the mesh norm \(h\) and the separation radius \(q\) decrease uniformly like the inverse of the square root of the number of nodes \(N\), i.e. \(h,q \sim \frac{1}{\sqrt{N}}\). Additionally, the nodes in these sets are not oriented along any vertices or lines, which emphasizes the ability of our method to handle arbitrary node layouts. They have been used quite successfully in many other RBF applications, e.g. [25, 26, 28, 33, 34, 55, 70].
1.2 B.2 Red Blood Cell
This manifold is a mathematical model for human red blood cells in static equilibrium conditions. The model was first proposed in [19] and has been used in many subsequent studies (e.g. [48]). The model can be described parametrically as follows:
where \(-\pi /2 \le \theta \le \pi /2, \,-\pi \le \lambda < \pi , \,r_0 = 3.91/3.39, \,c_0 = 0.81/3.39, \,c_2 = 7.83/3.39\), and \(c_4=-4.39/3.39\). The node sets we used for discretizing this manifold were obtained using a radial projection of the ME points for the surface of the sphere described above.
1.3 B.3 Bumpy Sphere
This manifold is constructed from the “bumpy sphere” surface [3] using the following procedure. First, \(N=5256\) points on the bumpy sphere surface were obtained from http://shapes.aimatshape.net. This surface is homeomorphic to the unit sphere, so spherical coordinates for each of the \(N=5256\) points on the bumpy sphere were obtained. Upon normalizing the original nodes by the maximum distance of any of the nodes from the origin, a parametric model for the surface was constructed using the RBF geometric modeling technique from [61]. The original (normalized) \(N=5256\) points on the bumpy sphere are used as the computational nodes in all the numerical experiments and the parametric model is used for computing the normal vectors to the surface.
1.4 B.4 Torus
This manifold is given by the implicit equation:
The node sets used for discretizing this manifold were obtained by arranging the nodes so that their Reisz energy (with a power of 2) is near minimal as described in Hardin and Saff’s seminal article [41]. As with the ME sphere nodes, these ME torus nodes sets provide a near uniform discretization of the manifold with an approximately uniform decrease in the mesh norm \(h\) and the separation radius \(q\). The node sets used in the experiments were kindly given to us by Drs. Douglas Hardin and Edward Saff and Ms. Ayla Gafni from Vanderbilt University.
1.5 B.5 Dupin’s Cyclide
This manifold is given by the implicit equation:
where \(a=2, \,b=1.9, \,d=1\), and \(c^2 = a^2 - b^2\). The node set for discretizing this manifold was obtained from the algorithm of Palais, Palais, and Karcher (PPK) [56], which will be included in a future version of the amazing mathematical visualization software package 3-D-XplorMath [1]. This algorithm generates “uniformly random” point clouds for surfaces and works for both implicit and parametrically defined surfaces. To obtain the node set displayed in Fig. 1e we used the PPK algorithm to generate a point cloud consisting of 9562 points. We then thinned this point set by removing points that were too close together in order to increase the separation radius. The final node set ended up at \(N=4948\). Both choices for the starting and ending number of nodes were chosen somewhat arbitrarily.
1.6 B.6 Bretzel2
This manifold is given by the implicit equation:
Unlike the the other surfaces, this one does not have an explicit parameterization. The node set for this domain was also obtained from the PPK algorithm. In this case we started with a point cloud of 10278 points and thinned it to end up with a node set of \(N=5041\) nodes.
Appendix C: Example Matlab Code
The code below simulates the Turing system (30) on the red blood cell surface using the parameters for spots in Table 1. The results should be similar to those shown in the top right plot of Fig. 6. The code below uses SBDF2 as the numerical time integration method to simplify the presentation. For the code to work, the \(N=4096\) ME node set for the sphere must be downloaded from [69], or from the second authors’ webpage http://math.boisestate.edu/~wright/research/.
Rights and permissions
About this article
Cite this article
Fuselier, E.J., Wright, G.B. A High-Order Kernel Method for Diffusion and Reaction-Diffusion Equations on Surfaces. J Sci Comput 56, 535–565 (2013). https://doi.org/10.1007/s10915-013-9688-x
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10915-013-9688-x
Keywords
- Radial basis functions
- Mesh-free
- Manifold
- Collocation
- Method-of-lines
- Pattern formation
- Turing patterns
- Spiral waves