First-Passage Time Statistics on Surfaces of General Shape: Surface PDE Solvers using Generalized Moving Least Squares (GMLS)
FFirst-Passage Time Statistics on Surfaces of General Shape:Surface PDE Solvers using Generalized Moving Least Squares (GMLS)
B. J. Gross , , P. Kuberry , P. J. Atzberger , [1] Department of Mathematics, University of California Santa Barbara (UCSB);[2] Department of Mechanical Engineering, University of California Santa Barbara (UCSB);[3] Sandia National Laboratories, Albuquerque, NM ∗∗ [email protected]; http://atzberger.org/ W e develop numerical methods for computing statistics of stochastic processes onsurfaces of general shape with drift-diffusion dynamics d X t = a ( X t ) dt + b ( X t ) d W t .We consider on a surface domain Ω the statistics u ( x ) = E x (cid:2)(cid:82) τ g ( X t ) dt (cid:3) + E x [ f ( X τ )] with the exit stopping time τ = inf t { t > | X t (cid:54)∈ Ω } . Using Dynkin’s formula, we com-pute statistics by developing high-order Generalized Moving Least Squares (GMLS)solvers for the associated surface PDE boundary-value problems. We focus particu-larly on the mean First Passage Times (FPTs) given by the special case f = 0 , g = 1 with u ( x ) = E x [ τ ] . We perform studies for a variety of shapes showing our methodsconverge with high-order accuracy both in capturing the geometry and the surfacePDE solutions. We then perform studies showing how FPTs are influenced by thesurface geometry, drift dynamics, and spatially dependent diffusivities. Introduction
First Passage Times (FPTs) [3, 12, 15] are of interest in many problems arising in fields includingbiology [28, 31, 34, 38, 40], physics [6, 8, 33, 42, 43, 49], engineering [20, 36, 46], finance [1, 22, 26,30], and machine learning [21, 51]. We consider here FPTs in the case of Ito stochastic processes X t with dynamics d X t = a ( X t ) dt + b ( X t ) d W t , [4, 12]. While stochastic numerical methods areavailable for approximating and sampling trajectories { X t } ≤ t ≤ T , [7], this can be computationallyexpensive given disparate dynamical time-scales and statistical sampling errors [4, 8, 34]. In practice,often of primary interest are a subset of statistics for characterizing underlying behaviors or formaking predictions of observations or measurements. We consider here the class of statistics for X t on an open domain Ω of the form u ( x ) = E x (cid:2)(cid:82) τ g ( X t ) dt (cid:3) + E x [ f ( X τ )] with the exit stoppingtime τ = inf t { t > | X t (cid:54)∈ Ω } . Through the choice of g , this provides information about the states x ∈ Ω realized by the stochastic trajectories. Through the choice of f , information is obtainedabout where stochastic trajectories encounter the boundary ∂ Ω. The Dynkin formula provides aconnection between these statistics and a collection of elliptic Partial Differential Equations (PDEs)with boundary-value problems of the form L u = − g, x ∈ Ω , u = f, x ∈ ∂ Ω, [12].For many problems, it is useful to consider stochastic processes X t in the manifold settingallowing for non-trivial topology or geometry [29, 33, 37, 43, 47]. We consider the case of stochasticprocesses X t within general surfaces, having potentially complicated shapes. We develop high-order methods for numerically computing statistics from the PDE boundary-value problems usingGeneralized Moving Least Squares (GMLS) approximations [17, 25]. Our approach provides meshless * Work supported by DOE Grant ASCR PhILMS DE-SC0019246 and NSF Grant DMS-1616353.** Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and EngineeringSolutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’sNational Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results andanalysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of theU.S. Department of Energy or the United States Government. SAND2021-XXXX J
Page 1 of 24 a r X i v : . [ m a t h . NA ] F e b ethods for solving PDEs on surfaces related to [19, 27, 45, 52, 53]. We focus particularly on themean first passage times given by the special case f = 0 , g = 1 with u ( x ) = E x [ τ ]. We performstudies for a variety of shapes to investigate how our methods converge. We show our methods canachieve high-order accuracy both in capturing the geometry and the solutions of the surface PDEs.We show how our methods can be utilized to investigate the behaviors of stochastic processes havingdrift-diffusion dynamics confined to surfaces. We perform studies for how FPTs are affected bysurface geometry, drift dynamics, and spatially dependent diffusivities. Our methods can also beused for computing other path-dependent statistics of stochastic processes confined to surfaces ofgeneral shape.
1. First-Passage Times and Path-Dependent Statistics on Surfaces
Figure 1.: (left) The mean First-Passage-Times (FPT) can be computed using our methods for manifoldsurfaces having a general shape with potentially complicated geometry or topology. Starting X t at any location X (0) = x on the surface, the FPT is obtained by computing the average arrival time at the boundary createdby slicing the surface at z = 0 . The manifold surface is represented by a collection of point samples. (right)The FPT is indicated by the color at location x on the surface. Our methods allow for computing FPTs forstochastic processes with general drift-diffusion dynamics d X t = a ( X t ) dt + b ( X t ) d W t . We show here resultsin the diffusive case with a ( x ) = 0 and b ( x ) = I . The mean First-Passage Time (FPT) is defined as u ( x ) = E X (0)= x [ τ ] , where τ = inf { t > | X t (cid:54)∈ Ω } . (1)This gives the time τ for a particle starting at location X (0) = x ∈ Ω to reach the boundary ofthe domain ∂ Ω. We denote by u ( x ) the average value of this passage time. In Figure 1, we showthe case of a manifold surface with boundary defined by a cut-plane at z = 0. The color at x indicates u ( x ) for the mean first passage time for a diffusion on the surface starting at x to reachthe boundary at z = 0, see Figure 1.We consider general drift-diffusion dynamics constrained to a surface M governed by theStochastic Differential Equations (SDEs), [12] d X t = a ( X t ) dt + b ( X t ) d W t , (2)subject to X t ∈ M . Page 2 of 24 he a term models drift in the dynamics and b the diffusive contributions. Throughout, we interpretour SDEs in terms of Ito stochastic processes [12].While in principle stochastic simulations can be performed to compute sample trajectories forMonte-Carlo estimates of statistics, in practice this can be expensive given disparate time-scales inthe dynamics and the slow decay of statistical errors [10, 34]. We develop alternative approachesusing Dynkin’s formula [12] which relates FPTs and other path-dependent statistics to PartialDifferential Equations (PDEs). We develop solvers for these associated boundary-value problems.For the stochastic process X t , we consider path-dependent statistics u of the form u ( x ) = E x (cid:20)(cid:90) τ g ( X t ) dt (cid:21) + E x [ f ( X τ )] . (3)The g ( x ) allows us to assign a value to individual trajectories based on the locations they traverseinside the domain x ∈ Ω and f ( x ) to assign a cost for the location where the process hits theboundary x ∈ ∂ Ω. The Dynkin formula for q ( x ) is given by [12] E x [ q ( X τ )] = q ( x ) + E x (cid:20)(cid:90) τ L q ( X s ) ds (cid:21) . (4)The L is the Infinitesimal Generator associated with the SDE in equation 2 which can be expressedas the differential operator L = a · ∂∂ x + 12 bb T : ∂ ∂ x . (5)Here, the · denotes the usual dot product with a gradient, and the : denotes the Hadamardelement-wise dot product (tensor contraction over indices). We take q ( x ) = f ( x ) , x ∈ ∂ Ω and L q ( x ) = − g ( x ) , x ∈ Ω. After rearranging terms, we have formally u ( x ) = q ( x ) , x ∈ ¯Ω. Using thisand equation 5, we can express the statistic u as the solution of a surface elliptic PDE of the form L u = − g, x ∈ Ω u = f, x ∈ ∂ Ω , (6)where L is given in equation 5. We require throughout f ∈ C ( ∂ Ω) and g ∈ C (Ω).In the special case with f = 0 and g = −
1, we obtain the mean first-passage time u ( x ) = E x [ τ ] , τ = inf { t | X t (cid:54) = Ω } . (7)This is the amount of time on average it takes for the drift-diffusion process starting at location x to hit the boundary of the domain ∂ Ω.To compute these and other path-dependent statistics for stochastic processes constrained tosurfaces, we develop solvers for the surface PDEs of equation 6. We show how these methods canbe used to investigate the roles played by geometry and the drift-diffusion dynamics modeled by a and b in first-passage times and the path-dependent statistics of equation 3. We consider general surface drift-diffusion dynamics associated with Ito stochastic processes of theform d X t = a dt + b d W t . (8) Page 3 of 24 he X t ∈ R n defines the location in terms of the ambient space. The a ( x ) ∈ R n and b ( x ) ∈ R n × m .In the embedding space we assume the drift a is in the tangent plane at each point a ( x ) ∈ T M x and that the range of the diffusion tensor b is in the tangent space range( b ) ⊂ T M x . As a result,the null-space of b T is the space of vectors orthogonal to the tangent space T M x . This has theconsequence that bb T has a null-space of vectors orthogonal to the tangent space and the range ofvectors spanned by b in the tangent space.In modeling systems, it is often convenient to specify the drift and diffusion using the embeddingspace representations of a and b . In numerical methods and when performing other practicalcalculations, it is often convenient to express the dynamics using local coordinate charts q = ( q , q )with embedding map x = σ ( q , q ). It will be convenient to have ways to convert between thesetypes of descriptions. The drift-diffusion dynamics in a local coordinate chart can be expressed asthe Ito stochastic process d q t = α dt + β d ˜ W t , (9)where α ( q ) ∈ R and β ( q ) ∈ R × ˜ m with ˜ m the number of noise sources.These descriptions are connected through the embedding map σ given by X = σ ( q ). Using theIto Lemma [12], we can connect the dynamics in the embedding space with those in the coordinatechart by d X t = (cid:20) α a + 12 Γ aij β i β j,T (cid:21) σ q a dt + σ q a β a d ˜ W t . (10)The Γ aij are the Christoffel symbols of the surface. Note that β a ∈ R is a row vector and d ˜ W t ∈ R ,so combining give a scalar for the coordinate component. In this notation we use the convention that ∂ q a = ∂ σ /∂q a = σ q a . We also use here the repeated index summation convention. This provides away to build up in terms of the components of the basis { ∂ q a } a the dynamics of X t in the embeddingspace.We obtain relations between descriptions using the expressions for the dynamics in the em-bedding space R and the dynamics expressed in the coordinates q . These can be summarized asfollows bb T = β a β b,T ∂ q a ∂ Tq b = β a β b,T σ q a σ Tq b , a c = g cb (cid:104) a , σ q b (cid:105) g , α c = a c −
12 Γ cab β a β b,T . (11)The β a β b,T is a scalar, since β c is a row vector. The (cid:104)· , ·(cid:105) g denotes the inner-product in theembedding space. The a c is useful so that we can represent the vector as a = a c σ q c . We stillsubscript with g even in the embedding space to highlight that this corresponds to the metricinner-product when using a coordinate chart. In practice, this is computed readily if we representthe vectors in the ambient space as the usual Euclidean inner-product between vectors. While thereare many possible choices of the noise terms β and b consistent with these equations, each areweakly equivalent since they have the same marginal probability densities. For convenience, wemake the specific choice b = σ q a β a . (12)This can be verified to satisfy the relations above. Page 4 of 24 he σ q a is a column vector and β a is a row vector so that b is a 2 × W t = ˜ W t with m = ˜ m = 2. With this correspondence, wealso can express β a in terms of b as σ Tq b b = σ Tq b σ q a β a = g ab β a . (13)This can be solved using the inverse metric tensor as β c = g cb σ Tq b b . (14)The σ Tq b is a row vector, b is a matrix, and this yields the row vector β c .We can substitute this in equation 11 to relate the coordinate chart drift and diffusion termsto the embedding space quantities as β c = g cb σ Tq b b , a c = g cb (cid:104) a , ∂ q b (cid:105) g , α c = a c −
12 Γ cab β a β b,T . (15)The Infinitesmal Generator L for the surface drift-diffusion process in equation 5 can be expressedin the local coordinate chart as L u = α c ∂u∂q c + 12 β c β d,T ∂ u∂q c ∂q d . (16)In practice, this is used to compute the action of L on u when evaluating terms in the surface PDEsin equation 6.We remark that a convenient feature of these calculations is that in the final expressionsno derivatives of the drift and diffusion coefficients were needed. This is particularly helpful oncomplicated surfaces for practical calculations. For modeling and simulation in practice, we usethe convention for the data structures that a is given as a vector field on the point cloud. We alsospecify that b is given in terms of components b [1], b [2], b [3], where b [ i ] is the i th column of thematrix. In this way, we can represent readily the tensors a , b and terms needed in the calculationsof L using equations 15 and 16 when solving equation 6.
2. GMLS Solvers for Surface Partial Differential Equations (PDEs)
To obtain first passage times and other path-dependent statistics, we need numerical methods forsolving the elliptic surface PDEs in equation 6. This requires approximation of the surface geometryand associated differential operators arising in equations 6, 15, and 16. This poses challenges giventhe high-order derivatives required of geometric quantities such as the local surface curvature andmetric tensor. We address this by developing meshless approaches based on Generalized MovingLeast Squares (GMLS) [17, 25]. We use this to obtain high-order approximations for geometricquantities associated with the local surface geometry facilitating the calculation of the neededderivatives. In meshless methods a collection of point samples is used to represent the manifoldgeometry and location of approximated surface field values, see Figure 2. We develop GMLSapproaches for estimating both local geometric quantities and differential operators to approximatethe action of the surface operators and to build collocation methods for the PDEs.
Page 5 of 24 .1. Generalized Moving Least Squares (GMLS) Approximation
The manifold is represented as a point cloud { x i } Ni =1 from which we need to approximate associ-ated differential operators and geometric quantities. Generalized Moving Least Squares (GMLS)approximates operators of the underlying surface fields by solving a collection of local least-squaresproblems [17, 25]. For a Banach space V with dual space V ∗ defined at each x i , an approximationis sought for a target operator τ x i [ u ], where τ x i ∈ V ∗ and u ( x ) ∈ V for x ∈ Ω ⊂ R d where Ω isa compact domain. In practice, we take V = V n to be a collection of polynomials up to degree n . We relate u ∈ V to a representative p ∗ ∈ V n , by considering a collection of linear functionalsΛ = { λ j } Nj =1 that serve to characterize u by Λ[ u ] = ( λ [ u ] , λ [ u ] , . . . , λ N [ u ]). We construct theapproximation p ∗ by solving the (cid:96) -optimization problem p ∗ = arg min p ∈ V n N (cid:88) j =1 ( λ j [ u ] − λ j [ p ]) ω ( λ j , τ x i ) . (17)The ω ( λ j , τ x i ) provides weights characterizing the importance of the particular sampling function λ j in estimating τ x j . For example, we can take λ j = δ ( x − x j ) with λ j [ u ] = u ( x j ) and we can take ω ( λ j , τ x i ) = ω ( x j − x i ) to be a function that decays as the distance increases between x j and x i .We can further take ω = ω ( (cid:107) x j − x i (cid:107) ) to be a radial function that decays to zero when r > (cid:15) . Inpractice, we use ω ( (cid:107) x j − x i (cid:107) ) = (1 − (cid:107) x j − x i (cid:107) /(cid:15) ) p + where ( · ) + = max( · , { φ i } d n i =1 for V n so that V n = span { φ , φ , . . . , φ d n } with d n = dim V n . Any function p ∈ V n now can be expressed as p ( x ) = Φ T a = d n (cid:88) i =1 a i φ i ( x ) . (18)Let τ [Φ] denote a vector where the components are the target operator acting on each of the basiselements. For the the target operator τ , we obtain the GMLS approximation ˜ τ by considering how τ acts on the optimal reconstruction p ∗ = Φ T a ∗ of u ,˜ τ [ u ] := τ [ p ∗ ] = τ [Φ] T a ∗ . (19)Conditions ensuring the existence of solutions to equation 17 depend primarily on the unisolvencyof Λ over V and distribution of the point samples { x i } . For theoretical results related to GMLSsee [17, 25]. GMLS has been primarily used to obtain approximations of constant coefficient lineardifferential operators [17]. Page 6 of 24 igure 2.:
Generalized Moving Least Squares (GMLS) on Surfaces. (left) For general surfaces, weapproximate the geometry and the action of operators τ [ u ] using variants of Moving Least Squares (MLS).To approximate a target function u ( x ) in the neighborhood of a given point x ∗ , Generalized Moving LeastSquares (GMLS) uses data samples { ( x j , u ( x j )) } mj =1 and probing functionals { λ k [ u ] } to select a local fit p ∗ = p ∗ ( x ; x ∗ ) ∈ V . For the function space V , we use throughout multivariate polynomials. (middle) A localapproximation to operators is obtained by solving a local least-squares problem with data weighted by thefunction w ( x ∗ , x j ) using equations 17 and 19. By using Principle Component Analysis (PCA) to find at x ∗ an approximate local tangent plane and using MLS in the Monge-Gauge, we can estimate surface geometricquantities such as the Gaussian curvature. (right) This approach is generalized further to estimate the actionof surface operators τ [ u ] , including the exterior derivative d , co-differential δ , Hodge star (cid:63) , and other surfaceoperators, see Appendix A. For our first-passage time problems, the target operators are technically non-linear given theirdependence on the geometry of the underlying manifold which also must be estimated from thepoint cloud. We handle this using a two stage approximation approach. In this first stage, weuse the point cloud to estimate two basis vectors ψ and ψ for the local tangent plane usingPrinciple Component Analysis (PCA) [2, 13]. We use these basis vectors to construct a localcoordinate chart ( q , q ). In this chart, we fit a function p ∗ ∈ V to the local point cloud to obtain aMonge-Gauge [14] representation of the surface ( q , q , p ∗ ( q , q )). We use GMLS to estimate targetgeometric quantities of interest, such as the local Gaussian curvature or high-order derivatives. In thesecond stage, we use the estimated geometric quantities to specify the target operators τ [ u ] for thefields u . We then use again GMLS to obtain an approximation of ˜ τ and to compute numerically thetarget operator values at x i . We have used related procedures for solving hydrodynamic equationson manifolds in [52]. We illustrate the approximation approach in Figure 2. To compute numerically the first passage times, we develop solvers for the elliptic PDE boundaryvalue problems given in equation 6. This is organized by representing the f and g function inputsand the a and b fields specifying the drift-diffusion dynamics. To avoid complications with localcoordinate charts, we numerically represent all input data globally using the ambient embedding Page 7 of 24 pace coordinates x ∈ R . The a ∈ R with only the tangential components playing a role inpractice. Similarly, the b tensor is represented by three vector fields b , b , and b . We use labelson the point cloud to determine which regions are to be considered interior to Ω and which arepart of the boundary of ∂ Ω. We only require f to be evaluated on ∂ Ω, while g must return reliablevalues for all x ∈ Ω.We use our GMLS methods in section 2.1 to estimate the surface geometric quantities and theaction of the operator in equation 6. This allows us to construct at each x i an equation for relating L u ( x i ) = − g ( x i ). This provides our collocation method for determining u ( x i ). Let [˜ u ] i = u ( x i ) and[˜ g ] i = g ( x i ). Collecting these equations together gives a sparse linear system A ˜ u = ˜ g . We solve thislarge sparse linear system using GMRES with algebraic multigrid (AMG) preconditioning usingTrilinos [18]. Our solvers have been implemented within a framework for GMLS problems usingthe Compadre library and PyCompadre [50]. The toolkit provides domain decomposed distributedvector representation of fields as well as global matrix assembly. The capability of the librarywere also extended for our surface geometry calculations by implementing symbolically generatedtarget operators. Our methods also made use of the iterative block solvers of (Belos [24]), blockpreconditioners of (Teko) and the AMG preconditioners of (MueLu [32, 35]) within the Trilinossoftware framework [18]. The framework facilitates developing a scalable implementation of ourmethods providing ways to use sparse data structures, parallelization, and hardware accelerations.
3. Convergence Results
We investigate the convergence of our GMLS solvers developed in section 2. The target operatorsthat arise in the surface PDE boundary-value problems involve a non-linear approximation. Thisarises from the coupling between the contributions to the error from the differential terms of thesurface operators and the GMLS estimations used for the surface geometry.
We perform studies using four different surface geometries: (i) ellipsoid, (ii) radial manifold I,(iii) radial manifold II, and (iv) torus. We label these as Manifold A–D, see Figure 3. We studyconvergence as the surface point sampling is refined. We take each refinement to have approximatelyfour times the number of points as the previous level. This aims to have the fill distance h halveunder each refinement.The manifolds can be described by the following implicit equations. Manifold A is an ellipsoiddefined by the equation x /a + y /b + z = s with a = 1 . , b = 1 . , s = 1. Manifold B is aradial manifold defined in spherical coordinates by ( θ, φ, r ( θ, φ ) where r ( θ, φ ) = 1 + r sin(3 φ ) cos( θ )with r = 0 .
1. Manifold C is a radial manifold defined in spherical coordinates by ( θ, φ, r ( θ, φ )where r ( θ, φ ) = 1 + r sin(7 φ ) cos( θ ) with r = 0 .
1. Manifold D is a torus defined by the equation( s − (cid:112) x + y ) + z = s with s = 0 . , s = 0 .
3. Each of the manifolds shown are representedby quasi-uniform point sets with approximately n = 10 samples. For quasi-uniform sampling weexpect the fill-distance h to scale as h ∼ / √ n . We report our results throughout using the notation¯ h − = √ n . Additional information on the number n used in samplings can be found in Appendix D.In our first passage time calculations, unless indicated otherwise, we generate the boundary ∂ Ω as the points at z = 0. In practice, we treat any points with | z | < × − as boundary points. Page 8 of 24 his serves to thicken the boundary region and at points x i near the boundary ∂ Ω helps achieveunisolvency in the local least-squares problems.
Figure 3.:
Point Sample Representations of the Surface Manifolds.
We study the accuracy of the our GMLS approximations of the surface operators by investigatingtheir action on test fields ˆ u ( x ). We generate u by using smooth functions ˆ u ( x, y, z ) parameterized inthe embedding space R . By the smoothness of the surface manifolds M , evaluation at x ∈ M givesa smooth surface field u . More formally, this corresponds to using the inclusion map ι : R (cid:44) −→ M toobtain u ( x ) = ι x ˆ u ( · ). We also use this inclusion map approach to generate test surface drift tensors a and diffusion tensors b for the stochastic process X t arising in equation 2 and 5.For our validation studies, we generate test fields u for the Manifolds A–D using u ( x ) = z ( x + y − x y ), which is an extension of a degree 5 Spherical Harmonic to R . For thedrift-diffusion tensors for Manifold A, we use a A = y ˜ σ + xz ˜ σ , b A = (cid:2) xz ˜ σ + y ˜ σ , xyz ˜ σ + cos( y ) ˜ σ , z ˜ σ + ( y + x ) ˜ σ (cid:3) . (20)The ˜ σ , ˜ σ provide a local orthogonal tangent basis at every point on the Manifold A. In practicein our numerical calculations, this does not need to be smooth in x and we construct these asconvenient using our local tangent plane approximation and Gram-Schmidt orthogonalization [23].For the Manifolds B and C, we use for the drift-diffusion tensors tangential projections of thevector fields in R given by a B,C = (cid:2) y, xz, x yz (cid:3) T , b B,C = (cid:2) ( xz, xyz, z ) , ( y , cos( y ) , y + x ) , ( x + y, e − z , e y ) (cid:3) T . (21)For Manifold D, we use drift-diffusion tensors a D = y σ u + xz σ v , b D = (cid:2) sin( xy ) σ u + y σ v , xyz σ u + cos( y ) σ v , e z σ u + ( y + x ) σ v (cid:3) . (22)The σ u , σ v provide a basis for the tangent space smooth in x given by the directional derivatives ofthe global parameterization σ ( u, v ) of the torus [14]. We characterize the accuracy of the GMLSapproximation ˜ L of the surface operators L using the (cid:96) -error (cid:15) op (cid:16) ˜ L (cid:17) = (cid:13)(cid:13)(cid:13) ˜ Lu − L u (cid:13)(cid:13)(cid:13) n, = 1 n n (cid:88) i =1 (cid:16) ˜ Lu ( x i ) − L u ( x i ) (cid:17) . (23) Page 9 of 24 n these studies, we evaluate to high precision the action of the operators L by symbolic calculationsusing SymPy [44]. In general, we emphasize that such calculations of expressions symbolically canbe prohibitive. Using this approach, we investigate the accuracy of the GMLS approximation of theoperator for each of the manifolds and with approximation spaces with multivariate polynomials ofdegree m ∈ { , , } . These results are reported in Tables 1– 3. We performed convergence studies for the Manifolds A–D shapes in Figure 3 with the test fields u and drift-diffusion tensors a , b discussed in Section 3.2. As the manifold resolution is refined,we study the convergence of the GMLS approximations of the surface operators using for spaces V multivariate polynomials of degree m ∈ { , , } . We estimate the convergence rates using thelog-log slope of the error (cid:15) op as the fill distance parameter h is varied between levels of refinement.We report these results in Tables 1– 3.We find in our empirical results that when using GMLS with polynomial spaces of degree m to evaluate the elliptic PDEs operator L of order k , we obtain convergence results of order m .We remark that there is no theory for our surface operators, since there is a non-linearity fromthe coupling between the surface geometry and the operator differential terms. Our results aresuggestive that our surface operator approximations achieve results exhibiting a trend consistentwith simpler differential operators approximated by GMLS [25]. This suggests our methods areresolving the surface geometric contributions with sufficient precision to minimize their contributionsto the error. Since the operator L is of differential order k = 2 we use throughout approximatingpolynomial spaces of at least degree m ≥ k ≥ Manifold A Manifold B Manifold C Manifold Dh (cid:96) -error Rate (cid:96) -error Rate (cid:96) -error Rate h (cid:96) -error Rate Table 1.:
GMLS Convergence Rates for L . For our GMLS methods using approximation spaces based onmultivariate polynomials of degree m = 2 , we find convergence rate of ∼ nd order. Manifold A Manifold B Manifold C Manifold Dh (cid:96) -error Rate (cid:96) -error Rate (cid:96) -error Rate h (cid:96) -error Rate Table 2.:
GMLS Convergence Rates for L . For our GMLS methods using approximation spaces based onmultivariate polynomials of degree m = 4 , we find convergence rate of ∼ th order. Page 10 of 24 anifold A Manifold B Manifold C Manifold Dh (cid:96) -error Rate (cid:96) -error Rate (cid:96) -error Rate h (cid:96) -error Rate Table 3.:
GMLS Convergence Rates for L . For our GMLS methods using approximation spaces based onmultivariate polynomials of degree m = 6 , we find convergence rate of ∼ th order.
4. Results for First Passage Times on Surfaces
For first passage times on surfaces, we investigate the roles played by the geometry, drift dynamics,and spatial-dependence of diffusivity. For our studies, we consider the general surface Langevindynamics d X t = − γ ∇ X U ( X ) dt + √ D I d W t (24)subject to X t ∈ M . The dynamics are expressed here in the embedding space X t ∈ R and constrained to be withinthe manifold surface M . We could also express these dynamics using local coordinate charts onthe surface as discussed in Section 1.1. The U is the potential energy, γ the friction coefficient, D = K B T /γ the diffusivity, and K B T the thermal energy with K B the Boltzmann constant and T the temperature [9]. In terms of the drift-diffusion tensors, a ( X t ) = γ − ∇ X U ( X ) , b ( X t ) = √ D I . Using our methods to capture drift in the surface dynamics, we study the role that can be playedby surface potentials in influencing stochastic trajectories and the first passage time. We considerthe case of a double-well potential U on a surface which from the conservative force F = −∇ X U in equation 24 introduces a drift into the dynamics influencing the mean first passage time. Weconsider geometries Ω described by the truncated torus, parameterized by X ( u, v ) = (cid:2) cos( u )( . . v )) , sin( u )( . . v )) , . v ) (cid:3) T , u ∈ (cid:20) π , π (cid:21) , v ∈ [0 , π ) . The double-well potential U is generated using the parameterization ( u, v ) as U ( x ) = k sin ( u ); x = X ( u, v ) . (25)We study how the mean first passage time changes as we vary the energy barrier k = ˜ k · k B T with ˜ k ∈ { . , . , . , . , . , . , . } . The case ˜ k = 0 . u = π and has energy barriers at u = π/ u = 3 π/ Page 11 of 24 oundaries at u = π/ u = 7 π/
4. We report the results of the first passage times starting atthe specific locations X i given by X : ( u, v ) = ( π, X : ( u, v ) = (cid:16) π, π (cid:17) X : ( u, v ) = ( π, π ) X : ( u, v ) = (cid:16) π , (cid:17) X : ( u, v ) = (cid:16) π , π (cid:17) X : ( u, v ) = (cid:16) π , π (cid:17) X : ( u, v ) = (cid:16) π , (cid:17) X : ( u, v ) = (cid:16) π , π (cid:17) X : ( u, v ) = (cid:16) π , π (cid:17) . We report our results in Figure 4.In our studies the starting locations X , X , X are at the minimum of the potential energy U .From these locations each trajectory must surmount at least one of the energy barriers to reach theboundary. This results in the largest first passage times. We find as the energy barrier is increased k (cid:38) k B T there is approximately exponential increase in the first passage times. This is in agreementwith theory for first passage times involving energy barriers based on asymptotic analysis usingKramer’s approximation [4]. Our numerical methods allow for computing the full solution u ( x ) ofequation 3, including first passage times, over a wide range of regimes k (cid:28) k B T , k ∼ k B T , and k (cid:29) k B T , and when starting at locations that are not only at the potential energy minimum, seeFigure 4. Figure 4.: (left) Double-well potential U on the surface of equation 25 influencing the drift of the stochasticdynamics in 24. (middle) The starting locations for first passage times with X (0) = X i . (right) First passagetimes from the selected starting points as the strength of the energy barrier k of the double-well potentialis varied. The locations X , X , X start at the minimum of the potential energy and exhibit the longestfirst passage times since they must surmount at least one of the energy barriers to reach the boundary. Asthe energy barrier is increased the first passage times increase, most rapidly when k (cid:38) k B T , and with anexponential trend in agreement with Kramer’s theory [4]. Page 12 of 24 .2. Role of Diffusion in Mean First Passage Times: Spatially-Dependent Diffusivities
We study the role of spatially-dependent heterogeneous diffusivities b = b ( X ) in first passage times.We consider the surface geometry Ω obtained from the truncated torus X ( u, v ) = (cid:2) . v ) , sin( u )(1 + 0 . v )) , − cos( u )(1 + 0 . v )) (cid:3) T , u ∈ U ( z = 0 , v ) , v ∈ [0 , π ) . This geometry is obtained by slicing the torus at z = 0. This defines the boundary ∂ Ω and forthe parameterization ( u, v ) a range U for u that depends on each value of v . We consider thespatially-dependent diffusivity b ( x ) = (cid:112) D ( x ) I = (cid:112) D (cid:18) − c exp (cid:18) −(cid:107) x − x c (cid:107) r (cid:19)(cid:19) I , (26) x c : ( u c , v c ) = (cid:18) π , π (cid:19) , r ≥ , ≤ c < . (27)This corresponds in equation 24 to D = D ( x ) , γ = γ ( x ) , and zero drift with U = 0 so a ( x ) = 0.The exponential is centered at x c with decay over length-scale r and influence amplitude c . For r sufficiently small the variations in diffusivity will occur primarily inside the domain Ω in a localizedregion (cid:107) x − x c (cid:107) (cid:46) r . The diffusivity elsewhere will be approximately constant ∼ D . We considerfor our initial starting locations X : ( u, v ) = ( π, X : ( u, v ) = (cid:16) π, π (cid:17) X : ( u, v ) = ( π, π ) X : ( u, v ) = (cid:18) π , (cid:19) X : ( u, v ) = (cid:18) π , π (cid:19) X : ( u, v ) = (cid:18) π , π (cid:19) X : ( u, v ) = (cid:18) π , (cid:19) X : ( u, v ) = (cid:18) π , π (cid:19) X : ( u, v ) = (cid:18) π , π (cid:19) . We study the role of spatially-dependent diffusivity in two cases: (i) when the depth c is variedwhile leaving the area scale r fixed, and (ii) when the area scale r is varied while leaving the depth c fixed. We use throughout the baseline parameters c = 0 . r = 0 .
5. We show the geometry andsolution for a typical spatially dependent diffusivity in Figure 5. We report our results in Figure 6.
Figure 5.:
Spatially Dependence Diffusivities. (left) The heterogeneous spatially dependent diffusivity b ( x ) inequation 26. (left-middle) The solution u ( x ) for the first passage time problem of equation 6. (right-middle),(right) The starting locations X (0) = X i . We report results of studies in Figure 6. We find varying the spatial extent r of the diffusivity plays a particularly strong role in thefirst passage times. The underlying mechanisms are different than the double-well potential. When Page 13 of 24 trajectory approaches x c the motion of the stochastic process slows down significantly as a resultof the smaller diffusivity. There is no barrier for X t to approach such regions but once within thisregion exhibits a type of temporary trapping behavior from the slow diffusion. We see from Figure 5,such a region influences the first passage time over a much larger range than the direct variations in b ( x ) given the high probability of a large fraction of trajectories encountering this trapping regioneven when a modest distance away. We see increases in r have a strong influence on increasingthe first passage times, see Figure 6. We also see in the limit that the localized diffusivity b ( x )approaches zero, even a relatively small probability of encountering this region can result in a largefirst passage time. We see this increasing influence as c approaches one, see Figure 6. Figure 6.:
Spatially Dependence Diffusivity and First Passage Times. (left) The first passage times for thespatially dependent diffusivity b ( x ) in equation 26 when varying the size r of the region of low diffusivity.(right) The first passage times when increasing the strength c of the diffusivity reduction in equation 26. We show how our methods can be used to investigate how the shape of the manifold surface influencesfirst passage times. In these studies, we consider purely diffusive dynamics. This corresponds to thecase with U = 0 with zero drift a ( x ) = 0 and a spatially homogeneous diffusivity b ( x ) = √ D I inequation 24. For the geometry we use a surface of revolution with an adjustable shape that forms aneck region near the bottom at z = 0, see Figure 7. This geometry is generated by considering firsta cylinder of radius r and height h = 0 .
05 which is capped at the top by a sphere. We then usea radial profile r ( z ) to connect the spherical cap with the bottom of the cylinder of radius r at z = 0. For this purpose, we use a bump function r ( z ) = (1 − . b ) + 0 . b exp (cid:16) − b b − b + . − z ) (cid:17) .For z ∈ [0 . , b + 0 . b so that this has arc-length π/
2, and the final radius r at z = 0.This serves to smoothly connect the unit hemisphere cap to the bottom of the cylinder with radius r . This ensures the geodesic distance from X (0) = X i to the boundary is always the same as wevary the shapes, see Figure 7. We report our results in Figure 8.We find as the neck region becomes smaller it acts as a hindrance for trajectories to reachthe boundary and first passage times increase. From the log-log plot we see that the first passage Page 14 of 24 ime appears over many radii to follow a power-law trend ¯ τ ∼ r α with values for α respectively for X (0) = X , X , X having α ≈ . , . , . X , X starting at points above the neck region, compared to the point X startingwithin the neck region. This indicates that as the radius shrinks the neck region increasingly acts asa hindrance for reaching the boundary. Interesting, we also see that starting at X also has FPTsthat significantly increase, since as r → X , X . We can characterize this mechanism by using a reaction-coordinate and notionof free energy (entropy contributions) arising from the constricting geometry. Figure 7.:
For the FPT studies when varying the geometry of the neck region, we show the select startinglocations X (0) = X i (red points). On the left we show the largest neck domain which does not present ageometric bottleneck for reaching z = 0 . On the right, we show the smallest neck region which presents asignificant geometric bottleneck potentially inhibiting diffusion in reaching z = 0 . The geometries used inthese studies were designed so that as the neck width is varied the geodesic distance from the starting pointsremains constant to reach the boundary at z = 0 . Page 15 of 24 igure 8.: (left) The first passage time as the neck radius r becomes small. We see as the neck becomesnarrower (towards the right) the geometry of the manifold restricts the diffusion and reaching the boundarybecomes more difficult. The r indicates the narrowest radius of the shape. (right) The shape acts to give ageometrically induced free energy Q ( z ) = − log (2 πr ( z )) for the dynamics when projecting to the z -axis. Wesee as the neck radius r narrows (towards the right) the geometry creates an increasingly large free-energybarrier and increases the first-passage times. Given the radial symmetry we can project the stochastic dynamics onto the z -component asan effective reaction-coordinate for reaching the boundary at z = 0. This provides a statisticalmechanics model [9] where the shape of the manifold gives an effective free energy having entropiccontributions proportional to Q ( z ) = − log(2 πr ( z )). This contributes to an effective mean-forcefor Z t yielding the drift a ( z ) = − ∂Q/∂z = − Q (cid:48) ( z ). This suggests the radially-averaged stochasticdynamics dZ t = − Q (cid:48) ( Z t ) dt + C dW t . From this perspective, the drift-diffusion process Z t must crossthe geometrically-induced energy barrier Q ( z ) to reach the boundary at z = 0. Similar to our resultson double-well potentials, this could be studied analytically with asymptotic Kramer’s theory [4].Our numerical methods allow for capturing directly these effects and we see as the neck narrowsthis significantly increases the first passage times. This correlates well with how the effective freeenergy barrier for Z t increases with the shape change, see Figure 8. Our results show how geometriceffects captured by our methods can provide mechanisms that influence significantly the behaviorsof stochastic processes and their first passage times. Conclusions
We have developed computational methods based on Generalized Moving Least Squares (GMLS) forinvestigating first passage times on surfaces of general shape and topology. We showed for a varietyof shapes that our methods converge with high-order accuracy both in capturing the geometryand the PDE solutions. We showed how our methods could be used to investigate how statisticsare influenced by the surface geometry, drift dynamics, and spatially dependent diffusivities. Thesolvers we have developed can be used more generally also for computing other path-dependentstatistics and solutions to elliptic boundary-value problems on surfaces of general shape.
Page 16 of 24 cknowledgements
The authors would like to acknowledge support from research grants DOE ASCR PhILMS DE-SC0019246 and NSF DMS-1616353.
AppendixA. Generalizing Vector Calculus to Manifolds using Exterior Calculus
Many of the operators used in vector calculus can be generalized naturally to the manifold setting.For this purpose, we utilize the operators of exterior calculus given by the Hodge star (cid:63) , exteriorderivative d , and vector to co-vector isomorphisms (cid:91), (cid:93) (definitions below) [5]. We denote vectorfields and tensors using the notation a = a i ...i k ∂ i · · · ∂ i k . We use ∂ i k to denote the basis vector ∂ i k = ∂ σ /∂ x ik and tensor product these together to represent vectors and tensors for the choice ofcoordinates x = ( x , x , . . . , x d ). We denote a differential k -form as α = (1 /k !) α i ,...,i k dx i ∧ · · · dx i k .The ∧ denotes the wedge-product of a tensor. We use the convention here with 1 /k ! to allowsummations over all permutations of the index values for i , . . . , i k . A more detailed discussion oftensor calculus on manifolds can be found in [5].We formulate the generalized operators in terms of the co-vectors (differential forms) f (cid:91) and F (cid:91) . We use that in the case of a scalar field we have quantitatively at each point f = f (cid:91) . Theisomorphisms (cid:91), (cid:93) mapping between the vector and co-vector spaces are given by a (cid:91) = (1 /k !) g i ,(cid:96) · · · g i k ,(cid:96) k a (cid:96) ...(cid:96) k dx i ∧ · · · dx i k (28) α (cid:93) = (1 /k !) g i ,(cid:96) · · · g i k ,(cid:96) k α (cid:96) ...(cid:96) k ∂ x i · · · ∂ x ik . (29)The exterior derivative d of a differential k -form α is defined in terms of the coordinates x as d α = 1 k ! ∂∂x j α i ,...,i k dx j ∧ dx i ∧ · · · dx i k . (30)The Hodge star (cid:63) is defined in terms of the coordinates x as (cid:63) α = (cid:112) | g | ( n − k )! k ! α i ,...,i k (cid:15) i ,...,i k ,j ,...,j n − k d x j ∧ · · · ∧ d x j n − k . (31)Note the indices have been raised here for the k -form with α i ,...,i k = g i (cid:96) · · · g i k (cid:96) α (cid:96) ,...,(cid:96) k . The (cid:15) (cid:96) ,...,(cid:96) n denotes the Levi-Civita tensor which gives the sign of the permutation of the indices (cid:96) , . . . , (cid:96) n and is otherwise zero [5].We can also generalize operators of vector calculus in the context of manifolds acting on scalarfields f and vector fields F asgrad M ( f ) = [ d f ] (cid:93) , div M ( F ) = − ( − (cid:63) d (cid:63) F (cid:91) ) = − δ F (cid:91) , curl M ( F ) = − (cid:63) d (cid:2) F (cid:91) (cid:3) , curl M ( f ) = [ − (cid:63) d f ] (cid:93) . (32)We define δ = ( − (cid:63) d (cid:63) ) which is referred to as the co-differential. To define d the exterior derivativeand (cid:63) the Hodge star, we consider the tangent bundle T M of the manifold and its dual co-tangent
Page 17 of 24 undle
T M ∗ . The tangent bundle defines the spaces for scalar fields, vector fields, and moregenerally rank m tensor fields over the manifold. The co-tangent bundle is the space of duals tothese fields. The co-tangent bundle can be viewed as the space of differential forms of order 0, 1,and m . More details on the associated differential geometry can be found in [5, 11, 14]. Additionalinformation on our approaches for using these operators can be found in [39, 41, 48]. B. Exterior Calculus Operators Expressed in Coordinates
We use as conventions that for differential 0-forms (scalar functions) f j = ∂ x j f , for differential1-forms (co-vector fields) α = α j d x j , and for vector fields v = v j ∂ j . In each case we have j ∈ { u, v } .The isomorphisms (cid:93) and (cid:91) between vectors and co-vectors can be expressed explicitly as v (cid:91) = ( v u σ u + v v σ v ) (cid:91) (33)= v u g uu d u + v u g uv d v + v v g vu d u + v v g vv d v = ( v u g uu + v v g vu ) d u + ( v u g uv + v v g vv ) d v ( α ) (cid:93) = ( α u d u + α v d v ) (cid:93) (34)= α u g uu σ u + α u g uv σ v + α v g vu σ u + α v g vv σ v = ( α u g uu + α v g vu ) σ u + ( α u g uv + α v g vv ) σ v . We use the conventions for the notation that for the embedding map σ we have σ u = ∂ u and σ v = ∂ v as in Appendix C. The exterior derivatives on these k -forms can be expressed as d f = ( ∂ u f ) d u + ( ∂ v f ) d v = f u d u + f v d v (35) d α = ( ∂ u α v − ∂ v α u ) d u ∧ d v. (36)The generalized curl of a 0-form and 1-form can be expressed in coordinates as − (cid:63) d f = (cid:112) | g | ( f u g uv + f v g vv ) d u − (cid:112) | g | ( f u g uu + f v g vu ) d v (37) − (cid:63) d α = ∂ v α u − ∂ u α v (cid:112) | g | . (38)A further discussions of these operators of differential geometry can be found in [5, 11, 14]. Additionaldetails on our approach can be found in [39, 41, 48]. C. Monge-Gauge Surface Parameterization
In the Monge-Gauge we parameterize locally a smooth surface in terms of the tangent planecoordinates u, v and the height of the surface above this point as the function h ( u, v ). This givesthe embedding map x ( u, v ) = σ ( u, v ) = ( u, v, h ( u, v )) . (39) Page 18 of 24 e can use the Monge-Gauge equation 39 to derive explicit expressions for geometric quantities.The derivatives of σ provide a basis ∂ u , ∂ v for the tangent space as ∂ u = σ u ( u, v ) = (1 , , h u ( u, v )) (40) ∂ v = σ v ( u, v ) = (0 , , h v ( u, v )) . (41)The first fundamental form I (metric tensor) and its inverse I − (inverse tensor) are given by I = (cid:20) g uu g uv g vu g vv (cid:21) = (cid:20) σ u · σ u σ u · σ v σ v · σ u σ v · σ v (cid:21) = (cid:20) h u ( u, v ) h u h v ( u, v ) h u ( u, v ) h v ( u, v ) 1 + h v ( u, v ) (cid:21) . (42)and I − = (cid:20) g uu g uv g vu g vv (cid:21) = 11 + h u + h v (cid:20) h v − h u h v − h u h v h u (cid:21) . (43)We use throughout the notation for the metric tensor g = I interchangeably. For notationalconvenience, we use the tensor notation for the metric tensor g ij and for its inverse g ij . Thesecorrespond to the first fundamental form and its inverse as g ij = [ I ] i,j , g ij = (cid:2) I − (cid:3) i,j . (44)For the metric tensor g , we also use the notation | g | = det( g ) and have that (cid:112) | g | = (cid:112) det( I ) = (cid:112) h u + h v = (cid:107) (cid:126)σ u ( u, v ) × (cid:126)σ v ( u, v ) (cid:107) . (45)The provides the local area element as dA u,v = (cid:112) | g | dudv . Lastly, we must compute the Christoffelsymbols of the second kind in our derivation of drift-diffusion on manifolds. In the general case thisis given by Γ kij = ∂ j σ i · σ k . (46)This can become complicated in general coordinates. However, in the case of the Monge Patch theabove simplifies to Γ kij = h ij h k (cid:112) | g | = h ij h k (cid:112) h u + h v . (47)For further discussions of these tensors and more generally the differential geometry of manifoldssee [5, 11, 14]. We use these expressions as the basis of our calculations of the action of our surfaceoperators. D. Sampling Resolution of the Manifolds
We provide a summary of the sampling resolution h used for each of the manifolds in Table 4. Werefer to h as the target fill distance . For each of the manifolds and target spacing h , we achieve anearly uniform collection of the points (quasi-uniform samplings) using DistMesh [16]. In practice,we have found this yields a point spacing with neighbor distances varying by about ≈ ±
30% relativeto the target distance h . We summarize for each of the manifolds how this relates to the numberof sample points n in Table 4. Robustness studies for related solvers have been carried out whenapplying perturbations to such point samples in [52]. Page 19 of 24 efinement Level A : h n B : h n C : h n D : h n1 .1 2350 .1 2306 .1 2002 .08 19122 .05 9566 .05 9206 .05 7998 .04 74783 .025 38486 .025 36854 .025 31898 .02 294944 .0125 154182 .0125 147634 .0125 127346 .01 118942 Table 4.:
Sampling Resolution for each of the Manifolds A–D. Relation between the target distance h andthe number of sample points n used for each of the manifolds. In each case, the neighbor distances betweenthe points sampled were within ≈ ± of the target distance h . References [1] Louis Bachelier. “Th´eorie de la sp´eculation”. fr. In:
Annales scientifiques de l’ ´Ecole NormaleSup´erieure
3e s´erie, 17 (1900), pp. 21–86. doi : . url : .[2] Karl Pearson F.R.S. “LIII. On lines and planes of closest fit to systems of points in space”.In: The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science doi : . eprint: https://doi.org/10.1080/14786440109462720 . url : https://doi.org/10.1080/14786440109462720 .[3] G. Klein and Max Born. “Mean first-passage times of Brownian motion and related prob-lems”. In: Proceedings of the Royal Society of London. Series A. Mathematical and PhysicalSciences doi :
10 . 1098 / rspa . 1952 . 0051 . eprint: https :/ / royalsocietypublishing . org / doi / pdf / 10 . 1098 / rspa . 1952 . 0051 . url : https ://royalsocietypublishing.org/doi/abs/10.1098/rspa.1952.0051 .[4] C. W. Gardiner. Handbook of stochastic methods . Series in Synergetics. Springer, 1985.[5] R. Abraham, J.E. Marsden, and T.S. Ratiu.
Manifolds, Tensor Analysis, and Applications . v.75. Springer New York, 1988. url : https://books.google.com/books?id=dWHet_zgyCAC .[6] Derek Y. C. Chan and Donald A. McQuarrie. “Mean first passage times of ions betweencharged surfaces”. In: J. Chem. Soc., Faraday Trans.
86 (21 1990), pp. 3585–3595. doi : . url : http://dx.doi.org/10.1039/FT9908603585 .[7] Kloeden.P.E. and E. Platen. Numerical solution of stochastic differential equations . Springer-Verlag, 1992.[8] Reinhard M¨uller, Peter Talkner, and Peter Reimann. “Rates and mean first passage times”.In:
Physica A: Statistical Mechanics and its Applications issn :0378-4371. doi : https://doi.org/10.1016/S0378-4371(97)00390-7 . url : .[9] L. E. Reichl. A Modern Course in Statistical Physics . John Wiley and Sons, 1998.[10] M Newman and G Barkema.
Monte carlo methods in statistical physics . Vol. 24. OxfordUniversity Press: New York, USA, 1999.[11] Micheal Spivak.
A Comprehensive Introduction to Differential Geometry . Ed. by 3. Vol. 1.Publish or Perish Inc., 1999.
Page 20 of 24
12] B. Oksendal.
Stochastic Differential Equations: An Introduction . Springer, 2000.[13] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.
Elements of Statistical Learning .Springer Series in Statistics. New York, NY, USA: Springer New York Inc., 2001.[14] A. Pressley.
Elementary Differential Geometry . Springer, 2001. url : https://books.google.com/books?id=UXPyquQaO6EC .[15] Sidney Redner. A Guide to First-Passage Processes . Cambridge University Press, 2001. doi : .[16] Per-Olof Persson and Gilbert Strang. “A Simple Mesh Generator in MATLAB”. In: SIAMReview
Scattered data approximation . Vol. 17. Cambridge university press, 2004.[18] Michael A. Heroux, Roscoe A. Bartlett, Vicki E. Howle, Robert J. Hoekstra, Jonathan J. Hu,Tamara G. Kolda, Richard B. Lehoucq, Kevin R. Long, Roger P. Pawlowski, Eric T. Phipps,Andrew G. Salinger, Heidi K. Thornquist, Ray S. Tuminaro, James M. Willenbring, AlanWilliams, and Kendall S. Stanley. “An Overview of the Trilinos Project”. In:
ACM Trans.Math. Softw. issn : 0098-3500. doi : . url : http://doi.acm.org/10.1145/1089014.1089021 .[19] John B. Greer, Andrea L. Bertozzi, and Guillermo Sapiro. “Fourth Order Partial DifferentialEquations on General Geometries”. In: J. Comput. Phys. issn :0021-9991. doi : . url : http://dx.doi.org/10.1016/j.jcp.2005.11.031 .[20] Russ Tedrake, Katie Byl, and JE Pratt. “Probabilistic stability in legged systems: Metastabilityand the mean first passage time (FPT) stability margin”. In: arXiv (2006).[21] J´erˆome Callut and Pierre Dupont. “Learning Partially Observable Markov Models from FirstPassage Times”. In: Proceedings of the 18th European Conference on Machine Learning . ECML’07. Warsaw, Poland: Springer-Verlag, 2007, 91–103. isbn : 9783540749578. doi : . url : https://doi.org/10.1007/978-3-540-74958-5_12 .[22] J. Hull. Options, Futures and Other Derivatives . Prentice Hall finance series. Pearson/PrenticeHall, 2009. isbn : 9780136015864. url : https://books.google.com/books?id=sEmQZoHoJCcC .[23] Richard L. Burden and Douglas Faires. Numerical Analysis . Brooks/Cole Cengage Learning,2010.[24] Eric Bavier, Mark Hoemmen, Sivasankaran Rajamanickam, and Heidi Thornquist. “Ame-sos2 and Belos: Direct and iterative solvers for large sparse linear systems”. In:
ScientificProgramming
IMA Journal of Numerical Analysis issn : 0272-4979. doi : .[26] ´Oscar Guti´errez. “American option valuation using first-passage densities”. In: QuantitativeFinance doi : . eprint: https://doi.org/10.1080/14697688.2013.794387 . url : https://doi.org/10.1080/14697688.2013.794387 . Page 21 of 24
27] Jian Liang and Hongkai Zhao. “Solving partial differential equations on point clouds”. In:
SIAM Journal on Scientific Computing
Proceedings of the National Academy of Sciences issn : 0027-8424. doi : . eprint: . url : .[29] E. Ben-Naim and P. L. Krapivsky. “First Passage in Conical Geometry and Ordering ofBrownian Particles”. In: First-Passage Phenomena and Their Applications . World Scientific,2014. Chap. Chapter 1, pp. 252–276. doi : . eprint: . url : .[30] R´emy Chicheportiche and Jean-Philippe Bouchaud. “Some Applications of First-PassageIdeas to Finance”. In: First-Passage Phenomena and Their Applications . World Scientific,2014. Chap. Chapter 1, pp. 447–476. doi : . eprint: . url : .[31] T. Chou and M. R. D’Orsogna. “First Passage Problems in Biology”. In: First-PassagePhenomena and Their Applications . World Scientific, 2014. Chap. Chapter 1, pp. 306–345. doi : . eprint: . url : .[32] Jonathan J. Hu, Andrey Prokopenko, Christopher M. Siefert, Raymond S. Tuminaro, andTobias A. Wiesner. MueLu multigrid framework . http://trilinos.org/packages/muelu .2014.[33] Remy Kusters and Cornelis Storm. “Impact of morphology on diffusive dynamics on curvedsurfaces”. In: Phys. Rev. E
89 (3 2014), p. 032723. url : https://link.aps.org/doi/10.1103/PhysRevE.89.032723 .[34] Ava J. Mauro, Jon Karl Sigurdsson, Justin Shrake, Paul J. Atzberger, and Samuel A. Isaacson.“A First-Passage Kinetic Monte Carlo method for reaction–drift–diffusion processes”. In: Journal of Computational Physics
259 (2014), pp. 536 –567. issn : 0021-9991. doi : https:/ / doi . org / 10 . 1016 / j . jcp . 2013 . 12 . 023 . url : .[35] Andrey Prokopenko, Jonathan J. Hu, Tobias A. Wiesner, Christopher M. Siefert, and RaymondS. Tuminaro. MueLu User’s Guide 1.0 . Tech. rep. SAND2014-18874. Sandia National Labs,2014.[36] Cenk Oguz Saglam and Katie Byl. “First Passage Value”. In: arXiv (2014). arXiv: . Page 22 of 24
37] O B´enichou, T Gu´erin, and R Voituriez. “Mean first-passage times in confined media: fromMarkovian to non-Markovian processes”. In:
Journal of Physics A: Mathematical and The-oretical doi : . url : https://doi.org/10.1088/1751-8113/48/16/163001 .[38] Nicholas F. Polizzi, Michael J. Therien, and David N. Beratan. “Mean First-Passage Times inBiology.” eng. In: Israel journal of chemistry
56 (9-10 2016), pp. 816–824.[39] Jon Karl Sigurdsson and Paul J. Atzberger. “Hydrodynamic coupling of particle inclusionsembedded in curved lipid bilayer membranes”. In:
Soft Matter issn : 1744-683X. url : http://dx.doi.org/10.1039/C6SM00194G .[40] Khem Raj Ghusinga, John J. Dennehy, and Abhyudai Singh. “First-passage time approachto controlling noise in the timing of intracellular events”. In: Proceedings of the NationalAcademy of Sciences (2017). issn : 0027-8424. doi :
10 . 1073 / pnas . 1609012114 . eprint: . url : .[41] B. Gross and P. J. Atzberger. “Spectral Numerical Exterior Calculus Methods for DifferentialEquations on Radial Manifolds”. In: Journal of Scientific Computing (2017). issn : 1573-7691. doi : . url : https://doi.org/10.1007/s10915-017-0617-2 .[42] C. Hohenegger, R. Durr, and D.M. Senter. “Mean first passage time in a thermally fluctuatingviscoelastic fluid”. In: Journal of Non-Newtonian Fluid Mechanics
242 (2017), pp. 48 –56. issn : 0377-0257. doi : https : / / doi . org / 10 . 1016 / j . jnnfm . 2017 . 03 . 001 . url : .[43] Alan E. Lindsay, Andrew J. Bernoff, and Michael J. Ward. “First Passage Statistics for theCapture of a Brownian Particle by a Structured Spherical Target with Multiple Surface Traps”.In: Multiscale Modeling & Simulation doi : .eprint: https : / / doi . org / 10 . 1137 / 16M1077659 . url : https : / / doi . org / 10 . 1137 /16M1077659 .[44] Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondˇrej ˇCert´ık, Sergey B. Kirpichev,Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, ThilinaRathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta,Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, ˇStˇep´anRouˇcka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and AnthonyScopatz. “SymPy: symbolic computing in Python”. In: PeerJ Computer Science issn : 2376-5992. doi : . url : https://doi.org/10.7717/peerj-cs.103 .[45] Ka Chun Cheung and Leevan Ling. “A kernel-based embedding method and convergenceanalysis for surfaces PDEs”. In: SIAM Journal on Scientific Computing . 2018, pp. 7063–7070. doi : . Page 23 of 24
47] Denis S. Grebenkov, Ralf Metzler, and Gleb Oshanin. “Strong defocusing of molecular reactiontimes results from an interplay of geometry and reaction control”. In:
CommunicationsChemistry issn : 2399-3669. url : https://doi.org/10.1038/s42004-018-0096-x .[48] B.J. Gross and P.J. Atzberger. “Hydrodynamic flows on curved surfaces: Spectral numericalmethods for radial manifold shapes”. In: Journal of Computational Physics
371 (2018),pp. 663 –689. issn : 0021-9991. doi : https://doi.org/10.1016/j.jcp.2018.06.013 . url : .[49] Adam Kells, Zsuzsanna ´E. Mih´alka, Alessia Annibale, and Edina Rosta. “Mean first passagetimes in variational coarse graining using Markov state models”. In: The Journal of ChemicalPhysics doi : . eprint: https://doi.org/10.1063/1.5083924 . url : https://doi.org/10.1063/1.5083924 .[50] Paul Kuberry, Peter Bosler, and Nathaniel Trask. Compadre Toolkit . Tech. rep. SandiaNational Laboratories, Jan. 2019. doi : . url : https://github.com/SNLComputation/compadre .[51] Junhong Xu, Kai Yin, and Lantao Liu. “Reachable Space Characterization of Markov DecisionProcesses with Time Variability”. In: (2019). arXiv: .[52] B.J. Gross, N. Trask, P. Kuberry, and P. J. Atzberger. “Meshfree methods on manifoldsfor hydrodynamic flows on curved surfaces: A Generalized Moving Least-Squares (GMLS)approach”. In: Journal of Computational Physics
409 (2020). url : https://doi.org/10.1016/j.jcp.2020.109340 .[53] Nathaniel Trask and Paul Kuberry. “Compatible meshfree discretization of surface PDEs”.In: Computational Particle Mechanics issn : 2196-4386. url : https://doi.org/10.1007/s40571-019-00251-2 ..