E F F I C I E N T A L G O R I T H M S F O R D I F F U S I O N - G E N E R A T E D M O T I O N B Y M E A N C U R V A T U R E By Steven J. Ruuth B M A T H , University of Waterloo, 1991 MSc, University of British Columbia, 1993 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F M A T H E M A T I C S A N D I N S T I T U T E O F A P P L I E D M A T H E M A T I C S We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A August 1996 © Steven J. Ruuth, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: Abst rac t This thesis considers the problem of simulating the motion of evolving surfaces with a normal velocity equal to mean curvature plus a constant. Such motions arise in a variety of applications. A general method for this purpose was proposed by Merriman, Bence and Osher, and consists of alternately diffusing and sharpening the front in a certain manner. This method (referred to as the MBO-method) naturally handles complicated topological changes with junctions in several dimensions. However, the usual finite dif-ference discretization of the method is often exceedingly slow when accurate results are sought, especially in three spatial dimensions. We propose a new, spectral discretization of the MBO-method which obtains greatly improved efficiency over the usual finite difference approach. These efficiency gains are obtained, in part, through the use of a quadrature-based refinement technique, by in-tegrating Fourier modes exactly, and by neglecting the contribution of rapidly decaying solution transients. The resulting method provides a practical tool, not available hitherto, for accurately treating the motion by mean curvature of complicated surfaces with junc-tions. Indeed, we present numerical studies which demonstrate that the new algorithm is often more than 1000 times faster than the usual finite difference discretization. New analytic and experimental results are also developed to explain important prop-erties of the MBO-method such as the order of the approximation error. Extrapolated algorithms, not possible when using the usual finite difference discretization, are proposed and demonstrated to achieve more accurate results. We apply our new, spectral method to simulate the motion of a number of three dimensional surfaces with junctions, and we visualize the results. We also propose and study a simple extension of our method to a nonlocal curvature model which is impractical to treat using the previously available finite difference discretization. Table of Contents Abstract ii List of Figures vi Acknowledgements ix 1 Introduction 1 1.1 Curvature-Dependent Motion 1 1.2 Methods for Curvature-Dependent Motion 4 1.3 Overview 8 2 Diffusion-Generated Motion by Mean Curvature Algorithm 10 • 2.1 The Two Phase Problem 10 2.2 Multiple Junctions 13 2.3 Selection of r 16 2.4 Finite Difference Discretizations of the MBO-Method 17 2.4.1 Selection of a Time-Stepping Method 17 2.4.2 Limitations of Finite Difference Discretizations 18 3 A New, Spectral Method 24 3.1 Discretization of the Heat Equation 24 3.2 Calculation of the Fourier Coefficients 26 3.3 Approximation of the Finest Subregions 30 3.3.1 Trivial Treatment of the Finest Subregions 30 iii 3.3.2 Piecewise Linear Approximation for Two-Phase Problems 31 3.3.3 Piecewise Linear Approximations for Junctions 39 3.4 Refinement Techniques 43 3.4.1 The Original Refinement Algorithm 43 3.4.2 A Method for a Gradual Refinement 47 3.5 Fast, Transform-Based Algorithms 49 3.5.1 Overview 52 3.5.2 The Unequally Spaced Fast Fourier Transform 53 3.6 Comparison to the Usual Finite Difference Discretization 57 4 Theoretical and Numerical Studies 60 4.1 Smooth Interfaces 60 4.1.1 Truncation Error Analysis 61 4.1.2 Extrapolation 63 4.1.3 Numerical Experiments 64 4.2 Nonsmooth Boundaries 66 4.3 Singularities in the Solution as Regions Disappear 69 4.4 Junctions in Two Dimensions 70 4.4.1 Error Analysis 71 4.4.2 Numerical Experiments 77 4.5 Summary 80 5 Numerical Experiments and Visualization 82 5.1 Three Dimensional, Two-Phase Problems 82 5.1.1 Visualization 82 5.1.2 Numerical Experiments 86 5.2 Junctions in Three Dimensions 87 i v 5.2.1 Visualization 87 5.2.2 Numerical Experiments 90 5.3 A Nonlocal Model 91 6 Conclusions 97 6.1 Summary 97 6.2 Future Research Directions 99 Bib l iography 101 A A n Es t imate for the N u m b e r of Basis Functions 107 v Lis t of Figures 1.1 A Region Evolving According to Mean Curvature Motion 2 .1.2 Evolution of Multiple Grains 3 1.3 A 3-state Potts Model 6 2.1 Initial Motion 12 2.2 Characteristic Set 12 2.3 After a Time, r 13 2.4 Sharpened Region 13 2.5 Initial Regions 15 2.6 Characteristic Sets 15 2.7 After a Time r 15 2.8 Sharpened Regions 15 2.9 Final Regions 16 2.10 Motion by Mean Curvature Results for Crank-Nicolson at Times, t . . . 19 2.11 A Smooth Shape with Widely Varying Local Curvatures 22 2.12 A Banded, Finite Difference Mesh 22 2.13 Sharpening a Shape 23 2.14 Sharpening a Perturbed Shape 23 3.1 Subdividing the Domain into its Coarsest Subregions 28 3.2 Dividing a Subregion 29 3.3 Subdividing a Square into Triangles 33 3.4 A Shape Represented by a Triangle 33 vi 3.5 A Shape Represented by Triangles 34 3.6 Errors Approximating Curved Segments 35 3.7 The Interpolation Step 36 3.8 Subdividing a Cube into Tetrahedrons 39 3.9 A Shape Approximated by a Tetrahedron 40 3.10 A Shape Approximated by Tetrahedrons 40 3.11 A Shape Represented as a Difference 41 3.12 Junction for which all Corner Phases Differ 42 3.13 Two Phases Represented at Corners 43 3.14 Refinement of a Smooth Region 45 3.15 Original Refinement can miss Slivers of the Region 45 3.16 A Neglected Section 46 3.17 A Problem with Sharp Corners 46 3.18 Refinement Methods for Corners 47 3.19 Gradual Refinement Captures the Entire Region 50 3.20 Refining a Cell 50 3.21 A Smooth Interface at Time, t 59 4.1 The Initial Interface 61 4.2 Extrapolated, Semi-Discrete Result 66 4.3 Initially Nonsmooth Interface at Times, t 67 4.4 The Initial Junction 72 4.5 A Smooth Three-Phase Problem 77 4.6 Evolution of a Junction Through a Singularity 79 5.1 From Q the Entire Curve is Visible 84 5.2 A Matrix Representation of the Surface 84 vii 5.3 Splitting a Shape into Easily Parameterized Portions 85 5.4 Piecewise Constant and Gouraud Shading of the Surface 85 5.5 Thin-Stemmed Barbell Moving by Mean Curvature Motion 88 5.6 Thick-Stemmed Barbell Moving by Mean Curvature Motion 89 5.7 Composition of a Junction 90 5.8 Three Phase Example Moving by Mean Curvature Motion 92 5.9 Four Phase Example Moving by Mean Curvature Motion 93 5.10 A Test Problem for the Nonlocal Curvature Algorithm 95 5.11 Nonlocal Model Which Preserves Area 96 A . l Contributing Rectangles 110 vm Acknowledgements First and foremost, many thanks to my supervisors Dr. Uri Ascher and Dr. Brian Wetton for their many helpful suggestions while working on this thesis. I would also like to thank Dr. Leslie Greengard for suggesting the use of the un-equally spaced fast Fourier transform and NSERC for financially supporting me during my graduate work. ix Chapter 1 Introduction 1.1 Curvature-Dependent Motion The topic of curve and surface evolution has recently generated tremendous interest in the mathematical sciences community and in several areas of application such as in grain growth and image enhancement. The computation of such motion and its rigorous justification have proved to be difficult tasks. This thesis considers the numerical approx-imation of curvature-dependent motion of surfaces with multiple phases. In particular, we study and develop fast methods for the case where the normal velocity, vn, of a surface is given by the sum of its principal curvatures (or mean curvature), K, We shall refer to this type of motion as motion by mean curvature. As illustrated for the two dimensional region of Figure 1.1, such a motion causes the most highly curved portions of a curve to smooth most rapidly. Indeed, any simple, closed curve in the plane shrinks to a small circle and disappears, regardless of the initial shape [31, 36]. For higher dimensional shapes, the most curved portions also move most quickly; however, topological merging and breaking is also possible [20]. Certain idealized models for grain growth fit into the framework of mean curvature motion. For example, when a liquid metal solidifies, crystallization begins in many locations with random orientations. As crystals grow, grains with different orientations meet, to form interfaces. To a good approximation, the surface energy of such a material . 1 Chapter 1. Introduction 2 Initial Motion Later Region (Solid) Figure 1.1: A Region Evolving According to Mean Curvature Motion is isotropic or independent of the orientation of the boundary of each grain [9]. By annealing the metal, it becomes warmed so that boundary atoms can change their phase. This produces an interface motion proportional to mean curvature motion [48, 51, 9]. In the interesting case where three or more grain orientations are present, junctions of moving surfaces can occur. This is illustrated in Figure 1.2 for the case of three grains. Although these important problems have been the subject of some study in two dimensions (e.g., [13, 14]), little is known about their numerical solution, especially in three dimensions. A n extension of the mean curvature model that we shall consider arises when the thermodynamic driving force of the interface motion depends on the volume phase change (i.e. bulk effects) as well as surface effects [71]. In this case, the normal velocity of the surface is given by its mean curvature plus a constant, We shall refer to this type of motion as affine velocity front motion. Curve evolution has also been studied extensively for image analysis and the enhance-ment of planar shape [43, 50, 42, 67, 41, 72]. Important application areas for this subject include noise suppression, image recognition and image interpretation. Mean curvature Chapter 1. Introduction 3 Figure 1.2: Evolution of Multiple Grains motion and affine motion are among the many curve evolution models that have been proposed for these applications. For example, affine velocity motion has been used to produce topological and other shape changes for characterizing shapes [42]. For other applications of these motions, see [42, 67] and references therein. Another important model that we will be considering fits into the framework of affine motion (1.2) and occurs when c(t) = - K a v ( t ) where K a v ( t ) is the average mean curvature over the surface at time t. Such motion preserves phase volumes (or areas in two dimensions) and occurs as a limit of a nonlocal model for binary alloys [57]. In the context of image enhancement, this nonlocal motion has also been suggested as a possible smoothing which preserves the area of shapes [64]. Theoretical aspects of mean curvature motion have also been studied extensively (see, e.g., [31, 36]). Some areas of computational interest include the approximation of bifurca-tion values [53] and the determination of self-similar solutions [21] under curvature flow. Computations of minimal surfaces have also been carried out by evolving surfaces accord-ing to mean curvature motion [20, 74]. Such surfaces find application in numerous areas [20] including soap film shapes, relativity theory, medical technology and architecture. Other applications for affine motion occur in certain flame propagation problems (e.g., Chapter 1. Introduction 4 grassfire flow) and automatic grid generation. See [54, 63, 66] for further details. 1.2 Methods for Curvature-Dependent Motion To study the phenomena outlined in the previous section, several numerical methods have been developed. Most of these can be divided into one of two groups: direct, or front tracking methods, where the motion of the interface is explicitly considered; and indirect methods, where the interface is given implicitly as the level set of some function. Several direct or front tracking methods have been proposed. In [14], a direct dis-cretization of the evolution equation for each interface is used to produce curvature-dependent motion for junctions in two dimensions. A number of other methods based on heuristic arguments have also been proposed for two dimensional junctions (see, e.g., [29, 30]). These methods are typically quite efficient for curves that never cross because they explicitly approximate the motion of the interface rather than a level set of some higher dimensional function. When line or planar segments interact, however, decisions must be made as to whether to insert or delete segments. Because very complicated topo-logical changes can occur in three dimensions, implementation of front tracking methods is often impractical in more than two dimensions. Another direct method, Brakke's method [9], minimizes surface energy to produce a variety of motions including motion by mean curvature. For three dimensional problems, however, this method requires user intervention when topological changes occur [9]. Also, there is no proof that the method actually approximates motion by mean curvature [71], although it is known that the Allen-Calm equation (see below) yields Brakke's motion in the limit for two-phase problems [40]. Phase field methods (see [16, 71] and references therein) give the interface as a level Chapter 1. Introduction 5 set of a reaction-diffusion equation such as the Allen-Cahn equation, ut = e A u — -f'(u). Here / ' is the derivative of a double well potential (e.g., f'(u) = u(u — l)(u — |)) and e is a small parameter. For the two-phase problem, it has been proven that certain phase field models produce mean curvature motion when e —> 0 [59, 58, 12, 26]. (Indeed, an important motivation for studying mean curvature motion is that it arises as a limit of certain phase field models.) By letting e be proportional to 7 ^ 7 , these methods provide a means of handling anisotropy (albeit in an ad hoc fashion) [71]. Junctions have also been treated by considering a vector-valued u [13, 56]. Unfortunately, phase field methods are often inherently too expensive for practical computation [48] because they represent the interface as an internal layer and thus require an extremely fine mesh (at least locally) to resolve this layer. Monte-Carlo methods for Q-state Potts models have also been used to simulate curvature-dependent motion with junctions (see, e.g., [33]). The Q-state Potts model assigns an integer state between 1 and Q to each point on a lattice to form regions (see Figure 1.3). The energy of each lattice site is taken to be a weighted average of the number of neighboring sites which are at a different state. To evolve regions, a random site and state are chosen. The selected site is set to the new state with a probability depending on temperature if the energy increases, and with probability 1 otherwise. It is unclear, however, what type of continuous motion this stochastic model approximates [71]. This method also introduces unwanted anisotropy into the motion due to the spatial mesh [71]. Furthermore, our numerical experiments indicate that these statistical methods are typically too slow to find accurate approximations to mean curvature motion. The Hamilton-Jacobi level set method of Osher and Sethian [54] is an appropriate Chapter 1. Introduction 6 1 1 1 1 1 1 1 1 1 1 1 1 i i i i 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 1 1[3 3 3 3 3 3 3 3 1 1 1 lT|3 3 3 3 3 3 3 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 1 1 2 2 2 2 2 213 3 3 3 ljT 2 2 2 2 2 2 2 213 3 2 2 2 2 2 2 2 2 2 2 2 2 Figure 1.3: A 3-state Potts Model choice for a wide variety of two-phase problems. This method naturally handles topo-logical merging and breaking in any number of spatial dimensions. To carry out the Osher-Sethian method for an initial surface T°, a continuous function is selected such that By solving r ° - {x e : $°(x) = o>. = F ( « ) | V $ | $(x,0) = for an arbitrary function of curvature, F{K), we obtain a surface T(t) = {x : $(x,t) = 0} which moves with a normal velocity F(K). It has been rigorously proven that this method converges for the case of mean curvature motion [27]. Furthermore, by assigning a func-tion, <fr;, to each region and evolving according to the method an extension to problems with junctions is possible [48]. This coupled Osher-Sethian approach allows for the Chapter 1. Introduction 7 specification of arbitrary velocities on each front. Unfortunately, this method can be excessively slow in three dimensions, in part because each <&; must be set to the signed distance to the interface at each step of the method. A very recent variational level set approach [74] has also been proposed for treating affine velocity motion. This method uses the level set formulation of Osher and Sethian [54] to arrive at an algorithm for a theoretical variational problem posed in [56]. This approach is reported to treat problems exhibiting topological changes with junctions in three dimensions [74]. The usefulness of the method is currently limited, however, by a time-stepping restriction similar to that associated with an explicit treatment of the heat equation. Such a restriction can produce an inefficient method for mean curvature flows because very small time steps must be used whenever a fine spatial mesh is applied. A method based on the model of diffusion-dependent motion of level sets has recently been proposed by Merriman, Bence and Osher [47, 48]. We shall refer to this method as the MBO-method (cf. [25]) although the name DGCDM-algorithrn has also been used [48]. The specifics of this method are elaborated upon in subsequent chapters, so will not be repeated here. This method naturally handles complicated topological changes with junctions in several dimensions. It also produces accurate approximations to mean curvature flow more efficiently than phase field models, the recent variational approach or Monte-Carlo methods. Thus, this method is often the only practical choice for three dimensional problems involving junctions. Even so, the method is often exceedingly slow for three dimensional and affine velocity motions. Similarly to all other methods for multiple-phase problems, no convergence results are known for the MBO-method. However, [25, 4] do give rigorous convergence proofs for two-phase mean curvature motion and [45] gives some further asymptotic results. In this thesis, we describe a fast, new algorithm which improves upon the speed of the usual discretization of the MBO-method, often by a factor of a thousand or more. We also Chapter 1. Introduction 8 develop analytic and experimental results to explain important properties of the method, such as the order of the approximation error and its enhancement by extrapolation. Finally, the improved utility of the new algorithm is demonstrated by approximating a number of two and three dimensional motions including a nonlocal curvature model. The resulting method provides a practical tool, not available hitherto, for accurately approximating the motion by mean curvature of complicated surfaces with junctions. 1.3 Overview An outline of the rest of the thesis follows. In Chapter 2, algorithms describing the MBO-method for two phase and multiple phase problems are given. This is followed by a discussion on how to select the step size, r , of the method. For the case of the finite difference discretizations originally proposed [47], the selection of an appropriate time-stepping scheme is discussed and several limitations of the method are identified. In Chapter 3, a new, spectral method for the realization of the MBO-method is pro-posed and described in detail. A spatial discretization is given and an efficient quadrature for calculating the corresponding Fourier coefficients is provided. This quadrature ob-tains accurate approximations to the front using a piecewise linear approximation to the surface and a gradual refinement technique. Unequally spaced transform methods for the rapid evaluation of the Fourier coefficients are also applied. This chapter concludes with a comparison of the proposed method and the usual finite difference approach. In particular, numerical experiments are presented to illustrate the efficiency gains which arise from our method. In Chapter 4, the order of accuracy of the MBO-method is studied for a variety of two dimensional shapes which move according to mean curvature motion. At first, local Chapter 1. Introduction 9 error expansions and numerical studies are carried out to demonstrate that the method produces a first order error in r for smooth, two-phase problems. Next, numerical stud-ies for nonsmooth corners and for singularities in the solution are discussed. The error here appears to be C(rlog(r)) . For three-phase problems, an expansion is derived which suggests that junction angles are stable and remain within 0(y/r) of the correct config-uration. Numerical studies are also given which confirm that an 0(y/r) error arises for junctions. Throughout this chapter, extrapolated algorithms are proposed and demon-strated to achieve more accurate results. While the discussion in Chapter 4 is independent of the discretization method used to implement the MBO-method, the spectral method of Chapter 3 was found crucial to realize the claimed errors at an affordable cost. In Chapter 5, the new method for discretizing the MBO-method is applied to a number of three dimensional surfaces. In particular, two-phase barbell-shaped regions are evolved according to mean curvature motion and the results are visualized. Mul-tiple phase motions are also approximated for examples involving three and four-phase junctions in three dimensions. A simple extension of the MBO-method to a nonlocal curvature model is also proposed and studied. Overall, the results are rather encouraging. Conclusions and suggestions of directions for future research are presented in Chapter 6. Chapter 2 Diffusion-Generated Motion by Mean Curvature Algorithm An algorithm for following interfaces propagating according to mean curvature mo-tion (1.1) was introduced in a paper by Merriman, Bence and Osher [47, 48]. This algorithm is one of the only methods of treating both topological changes and junc-tions. This chapter describes the method for the two phase and multiple phase problems and discusses some serious limitations of the usual finite difference discretization of the method. Subsequent chapters discuss a new method which greatly improves upon the usual finite difference approach. In this and later chapters, we will make use of a simplified version of the Von Neumann-Mullins parabolic law [51]. Specifically, we will use the fact that the area, A , enclosed by a simple curve which moves by mean curvature motion obeys % = -*<• <2"3' 2.1 The Two Phase Problem Suppose we wish to follow an interface moving by mean curvature motion (see, e.g., Fig-ure 2.1). To carry out this motion over a domain, T>, the MBO-method uses a diffusion-generated motion: MBO-Method (Two Regions) B E G I N (1) Set U equal to the characteristic function for the initial region. 10 Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 11 1 if x belongs to the initial region i.e., set U(x,0) = 1 0 otherwise. R E P E A T for all steps, j , from 1 to the final step: B E G I N (2) Apply diffusion to U for some time, T . Ut = AU, i.e., find U(x,jr) using | ^ = 0 on dV on starting from U(x,(j — l ) r ) . (3) "Sharpen" the diffused region by setting 1 i f < 7 ( x , j r ) > ! 0 otherwise. U(X,JT) = < E N D E N D For any time t, the level set {x : U(x,t) — |} gives the location of the interface. In step (2) of the method, zero flux conditions dU — =0ondV (2.4) on are selected. These boundary conditions cause the curve to meet the boundary at right angles, as is appropriate for the case of grain growth [14]. It has also been argued that these boundary conditions (2.4) are the most natural for image processing since they do not impose any value at the boundary [2]. Although our experiments will concentrate on the important zero flux case, other types of boundary conditions are sometimes consid-ered. In particular, Dirichlet conditions have been used for computing minimal surfaces [20]. Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 12 Figure 2.1: Initial Motion Figure 2.2: Characteristic Set To illustrate this algorithm for the problem of Figure 2.1, a grey-scaled representation of U after each of the steps (1) to (3) is given below: After step (1) U = 1 for the black region of Figure 2.2 and U = 0 elsewhere, (2) U ranges between 0 and 1 as represented by the greyscale image1 of Figure 2.3, (3) U = 1 for the black region of Figure 2.4 and U = 0 elsewhere. An extension to the case of affine velocity front motion (1.2) is also possible. To obtain such a motion, we track the level set instead of the usual level set of | [45]. 1 T h e " r ing ing" i n F igure 2.3 and others is an artifact which arises from pr in t ing . These features are not present in the or ig ina l grayscale images. Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 13 Figure 2.3: After a Time, r Figure 2.4: Sharpened Region 2 .2 M u l t i p l e J u n c t i o n s An extension of the previous method to multiple phases has also been made [48]. To treat such problems, the characteristic function for each region is diffused for a time r , as is outlined below for the case of 120-120-120 degree junction angles2. MBO-Method (Multiple (r) Regions) B E G I N (1) For i = 1,..., r Set Ui(x,0) equal to the characteristic function for the ith. region. R E P E A T for all steps, j, from 1 to the final step: B E G I N 2See [45, 48] for an extension to nonsymmetric junction angles. Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 14 (2) For i — 1,..., r, starting from U(x, (j — 1)T), Apply diffusion to Ui for some time slice, r . i.e., find Ui(x,jr) using < ^ = AU- ' ^ = o on ax>. k on (3) "Sharpen" the diffused regions by setting the largest Ui equal to 1 and the others equal to 0 for each point on the domain: For i = 1,...,r Q . T T ^ . , j lit $ e {ye T>:Ui(y,JT)>Uj(yjT),j = l, Set Ui[x,jT) — s. [ 0 otherwise E N D E N D For any time t, the interfaces are given by M {x : Ui(x,t) = max{Uj(x,t)}}. (2.5) i=l,...,r To illustrate this algorithm for the problem of Figure 2.5, a grey-scaled representation of Ui after each of the steps (1) to (3) is given below: After step (1) Ui = 1 for the black regions of Figure 2.6 and Ui = 0 elsewhere, (2) Ui ranges between 0 and 1 as represented by the greyscale images of Figure 2.7, (3) Ui = 1 for the black regions of Figure 2.8 and Ui = 0 elsewhere. Reconstruction of the interfaces by (2.5) then gives the final regions displayed in Fig-ure 2.9. Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm Figure 2.8: Sharpened Regions Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 16 Figure 2.9: Final Regions 2.3 Selection of r To accurately resolve the motion of features of the interface, it is also important to select r appropriately. In particular, r should be small enough so that diffusive information does not travel a distance comparable to the local radius of curvature, ^ (see [48]). By Equation (2.3), it is straightforward to show that a shrinking circular interface of curvature, disappears at time t = This gives a restriction on T , V ^ < - - (2-6) Diffusion must also proceed long enough so that the motion of the interface over each step can be resolved by the spatial discretization. For the case of a finite difference discretization, the level set U = | must move at least one grid point, otherwise the interface remains stationary. This produces the restriction that (speed of motion of the interface) x r > grid spacing Letting K be the curvature and h the grid spacing, we arrive at a second restriction for the finite difference approach, ACT » h. (2.7) As we shall see, the restriction (2.7) does not appear for the new, spectral method that we propose in Chapter 3. Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 17 2.4 F in i t e Difference Discret izat ions of the M B O - M e t h o d We now discuss several aspects of the usual finite difference approach for simulating the MBO-method. In particular, this section considers how to select a time-stepping scheme and outlines several limitations of finite difference discretizations of the MBO-method. 2.4.1 Selection of a Time-Stepping M e t h o d Explicit time-stepping for the diffusive step of the MBO-method is very expensive because of the usual stability time step restriction arising from the heat equation. This stability time step restriction can be overcome by using implicit time-stepping schemes. When applied to the semi-discrete heat equation, some implicit methods, such as backward Euler, give a strong decay of high frequency error modes. Other methods, such as Crank-Nicolson, and ADI methods [69] give a very weak decay of these error modes [3, 69]. After the sharpening step of the MBO-method, the solution is discontinuous. Because high fre-quency modes make a very important contribution to such a result, it is essential to use a time-stepping scheme which produces appropriate damping (i.e., strong damping) of these components. To illustrate this idea, we shall consider the motion by mean curvature of a shrinking circular interface with initial radius | . By Equation (2.3), the exact solution of this U = AhU Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 18 problem at time t is a circle of radius \J{\)2 — It. Using a finite difference simulation of the MBO-method, this circular interface was evolved to time t = 0.02. The spatial discretization was carried out using the standard 5-point Laplacian with grid spacing h = 1/256. For the time discretization, both backward Euler and Crank-Nicolson were considered using a time step, A i = 0.0005. The value of r was chosen to be a multiple of A i such that easy to show that this choice satisfies both restrictions (2.6) and (2.7). For general problems, however, an appropriate value of r is difficult or impossible to determine using Equations (2.6) and (2.7). In Figure 2.10, numerical results are given for Crank-Nicolson at various times, t. These results demonstrate that high frequency error modes can linger on to produce disastrous results when a weakly damping time-stepping scheme such a Crank-Nicolson is used. Application of the strongly damping backward Euler, however, gives the cor-rect phase areas to within 3%. This simple time-stepping technique is often adequate for diffusion-generated motion by mean curvature because other, larger sources of error dominate (see Chapter 4). 2.4.2 L imi t a t ions of F in i t e Difference Discret izat ions We have seen from the previous section that there are restrictions on our choice of r and h for finite difference discretizations of the MBO-method. To produce an accurate result, r should be small enough (see Equation (2.6)) so that diffusive information does not travel a distance on the order of the local radius of curvature. Once a sufficiently small T is selected, the mesh spacing, h, must be chosen small enough (see Equation (2.7)) so Because the solution to this problem at time i is a circle of radius J{\)2 — 2i, it i is Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 19 Figure 2.10: Motion by Mean Curvature Results for Crank-Nicolson at Times, t Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 20 that the level set moves at least one grid point, otherwise the sharpening step leaves the front stationary. Since the MBO-method does not explicitly evaluate curvature, we would prefer to select r and h without using Equations (2.6) and (2.7). One simple method for adap-tively choosing r for a given mesh, is to set r to a multiple of the time taken for the level set to move one grid point. Such a choice ensures that the front propagates and produces r-values which are roughly inversely proportional to the maximum speed of propagation of the front. Unfortunately, numerical tests indicate that an appropriate choice of the multiple is not always possible because it depends on factors such as the local curvature of the interface and the local mesh spacing. Thus, this adaptive method is often inappropriate for finding accurate results. Satisfying the restriction (2.7) can be computationally impractical even for smooth, two dimensional problems. Consider, for example, the motion by mean curvature of the boundary of the spiral region given in Figure 2.11. (The curvature-dependent motion of similar shapes has been considered in biological models [73].) Since the local curvature of the boundary of this problem varies tremendously, it is impractical to satisfy (2.7) everywhere using a uniform mesh. To achieve a more efficient finite difference algorithm, one might consider discretizing the MBO-method using a local mesh refinement at the level of the P D E . However, car-rying out local mesh refinement is rather involved for level set methods when curvature terms arise (see [49]). An alternative approach is to place a narrow band of grid points around the front (cf. [1]). Even this optimized, finite difference approach can lead to a prohibitive number of operations per step when an accurate solution is sought. For example, consider the motion by mean curvature of a smooth curve. For such a curve, each step of the MBO-method produces an C ( r 2 ) error in the position of the front (see Section 4.1.1). To preserve the overall accuracy of the method, grid points must be Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 21 at most a distance 0(r2) apart since each step produces an error which is comparable to the mesh spacing. Noting that the front travels a distance 0(T) per step of the method (see Section 4.1.1), it is clear that a minimum of O grid points are needed to safely band a curve (see, e.g., Figure 2.12). Thus, a minimum of O operations per step are required to preserve the overall accuracy of the method, which is often prohibitively expensive when accurate results are sought. A further limitation of the finite difference approach is that the error is not regular. Specifically, very small differences in the position of the level set 1/2 before sharpening can produce jumps in the front location after sharpening. This type of error is unde-sirable because it makes the construction of higher order accurate, extrapolated results impractical. Figures 2.13 and 2.14 illustrate how a small change in the position of the level set 1/2 can lead to a jump in the front location after sharpening. (We shall see in the next two chapters that our proposed method essentially eliminates the spatial error to allow for higher order accurate extrapolations in r.) To avoid the limitations outlined in this section, we introduce a new, spectral method for realizing the MBO-method in the next chapter. Figure 2.12: A Banded, Finite Difference Mesh Chapter 2. Diffusion-Generated Motion by Mean Curvature Algorithm 23 t + + * +i f +• M t +• +1 +- H + + h + + + + + + + + + + + + 4-u * I Sharpening 9 • • * • • • • • o — o o o o 6 o o 6 o o 6 o o 6 ..©.——©. -6 Figure 2.13: Sharpening a Shape 7 4 *• t * + 4* — - A H + + + + + + + + + + + + Sharpening o o o o o 6 o o o 6 o o o 6 o o 6 o o -b Figure 2.14: Sharpening a Perturbed Shape Chapter 3 A New, Spectral Method As we shall see later in this chapter, accurate computation of solutions using the usual finite difference discretization of the MBO-method can be expensive, even for simple two dimensional problems. Since we are mainly interested in three dimensional problems or problems involving more than two phases, a faster method is desired. This chapter describes a new, spectral method for realizing the MBO-method which is typically much faster than the usual finite difference approach. For notational simplicity, the algorithm focuses on the two dimensional case over the domain, T> = [0,1] x [0,1]. Certain extensions to three spatial dimensions and more phases are also discussed. 3.1 Discretization of the Heat Equation As we have seen in Section 2.1, carrying out diffusion-generated motion by mean curvature requires us to solve the heat equation ut = Au, (3.8) du „ — = 0 on dV On repeatedly over time intervals of (possibly variable) length r , starting from the charac-teristic function of the region to be followed. Over any of these time intervals, u may be 24 Chapter 3. A New, Spectral Method 25 approximated by the Fourier cosine tensor product, n -1 U(x,y,t) = ^2 dij(t) cos(nix) cos(7rj'y) (3-9) i,j=o because zero flux boundary conditions are considered. Substitution of (3.9) into (3.8), then gives n -1 U(x,y,t) = cij exp(-ir2\i2 + j2][t-tstart])cos(irix)cos(TTJy) (3.10) for t s t a r t < t < tstart + r, where Cij = dij{tstart) and tstart is the time when the current interval starts. One might expect that a Fourier spectral approximation for u would be unwise because u is initially discontinuous at interfaces. We are only interested in the solution after a time r , however. After a sufficiently large time T , high frequency modes have dissipated. Since the problem is linear, different modes do not interact (they are eigenfunctions) and thus there is never a need to approximate high frequency modes (not even near tstart-, when high frequency modes make an important contribution to the solution). For this reason, an accurate approximation to (3.8) at time r can be obtained using far fewer basis functions than might otherwise might be expected. Indeed, Appendix 1 shows that for smooth problems in two dimensions, selecting the number of Fourier modes in each direction to satisfy V TT T for a curve of length L produces an error in the position of the front which is at most e + O(CT). In practice, however, our implementations simply select an n satisfying » > f f (3-11) V TT T and verify the corresponding results by repeating the calculation with a different n (cf. [8])-Chapter 3. A New, Spectral Method 26 3.2 C a l c u l a t i o n of the F o u r i e r Coefficients The values of the Fourier coefficients, c,-j, of equation (3.10) must still be determined at the beginning of each time step (i.e., immediately following the sharpening of the previous step). In fact, we carry out sharpening (at t = tstart) as part of the calculation for the Fourier coefficients. These coefficients are found using an adaptive quadrature method rather than a pseudospectral method. Begin by defining to be the approximation of the phase we are following. By multiplying equation (3.10) at time t = tstart by cos(irix) cos(7rjy) and integrating over the domain we obtain - i ri n-1 Rt = {{x,y) : U(x,y,t) > -} which simplifies via the usual orthogonality conditions to give c0j Jo lo1 U(x, y> tstart) dx dy, 2 Jo Jo U{x,y, tstart) cos(7rza;) dx dy 2 lo lo U(x,y, tstart) cos(njy) dx dy for i ^ 0, for j ^ 0, 4 / 0 1 / o 1 ^ ( a ; 5 y 5 ^ ^ r O c o s ( 7 r ^ ) c o s ( 7 r i y ) dx dy for i,j ^ 0. Immediately after sharpening, U{x,y,t) 0 otherwise which implies that (3.12) Chapter 3. A New, Spectral Method 27 where 1 if i = j = 0 Hi = \ 4 if % ± 0 and j ^ 0 (3.13) 2 otherwise Thus, simple functions must be integrated over a complicated, non-rectangular region, Rt. This may be accomplished by recursively subdividing the domain (cf. [62, 61]), as we illustrate for the region, R, given in Figure 3.1a. We begin by evaluating U at the corners of a number of equally-sized subregions, so as to capture the large-scale features of the shape. Typically, n x n subregions are selected because the corresponding [/-values can be evaluated in just (9(n2log(n)) operations using a fast Fourier transform (see, e.g., [17]). If the phase at all four corners of any subregion corresponds to white, then we assume that the subregion does not intersect with R and hence no contribution to the Fourier coefficients is made. This case corresponds to the subregions of Figure 3.1b which have at least one dashed edge. If all four corners of a subregion, Q, correspond to grey, however, we assume that Q C R and add a contribution to each of the Fourier coefficients, C ; J , for 0 < i,j < n — 1. This case corresponds to the subregions of Figure 3.1b which have at least one thin, solid edge. Finally, if two phases occur, further subdivisions are carried out. We demonstrate this subdivision procedure for the subregion, Q, of Figure 3.1b. Because Q is a mixed region, we divide it into quadrants, as shown in Figure 3.2b. Since the phase color at all corner points of quadrant Q\ is white, we assume that this quadrant does not intersect with 7^ and hence does not contribute to the Fourier coefficients. For each of the remaining quadrants, Q\, Q\ and Q\, two phases occur, so Q Chapter 3. A New, Spectral Method 28 —I Fig. 3.1a. Initial Region, R Fig. 3.1b. Coarsest Subdivisions Figure 3.1: Subdividing the Domain into its Coarsest Subregions further subdivision is required. See Figure 3.2c. Focusing on the refinement of the subregion, Q\, we find that the phase of the upper right hand corner of Q\ is different than that of the other corners. Thus, Q\ is also subdivided. Corner points of the remaining subregions are grey, so we assume Q\ C R for k = 2,3,4 and add contributions to each of the Fourier coefficients, c,j, for 0 < i,j < n — 1. Recursive subdivisions of the domain continue (see, e.g., Figure 3.2d) until regions containing multiple phases can be safely approximated by some simple numerical tech-nique. The next section discusses methods for approximating the regions at the finest grid subdivisions. Chapter 3. A New, Spectral Method (0.375, 0.875] (0.5, 0.S75] (0.375, 0.875) (0.375, 0.75) (0.5, 0.75) Fig. 3.2a. Initial Subregion (0.5, 0.875} cr Q 1 Q 4 (0.375, 0.75) (0.5, 0.75) Fig. 3.2b. One Subdivision (0.375, 0.875) (0.5, 0.875) Q 3 2 Q 3 Q 1 2 Q' (0.375, 0.75) (0.5, 0.75) (0.375, 0.875) (0.5, 0.875) (0.375, 0.75) Fig. 3.2c. Two Subdivisions Fig. 3.2d. Four Subdivisions Figure 3.2: Dividing a Subregion Chapter 3. A New, Spectral Method 30 3.3 A p p r o x i m a t i o n of the F i n e s t S u b r e g i o n s In the previous section, a method was introduced for recursively dividing the domain into rectangles. At some point, however, we must stop subdividing and treat the finest cells. This section discusses how to approximate the contributions to the Fourier coefficients at the finest grid subdivisions. 3.3.1 T r i v i a l T r e a t m e n t of the F i n e s t S u b r e g i o n s A n easy method for treating the finest grid subdivisions is to add one half the contribution that would occur for the whole subregion to each of the Fourier coefficients. In other words, a contribution Often, this approach is unsatisfactory when accurate results in three dimensions are sought. If the finest subregions are of width h, then an O(h) error in the position of next chapter), leading to O ( T 2 i - 2 J cells at the finest level of subdivision in d dimensions. Such a mesh is often impractical to treat, especially when three dimensional results are sought. To achieve a more accurate approximation of regions and hence much faster results, we next consider piecewise linear approximations of the interface. is added to c tj,0 < i,j < n — 1, for each of the finest subregions, Qk. the front occurs. For smooth curves, we seek an 0(T2) approximation of the front (see Chapter 3. A New, Spectral Method 31 3.3.2 Piecewise Linear Approx ima t ion for Two-Phase Problems In the last subsection, we saw that a trivial treatment of the integrals (3.12) at the finest grid subdivision produced an 0{h) error in the position of the front, where h is the width of the finest grid subdivision. To produce an O (j^- + / i 2 ) approximation of the interface, a simplicial decomposition of the region, R, with a piecewise linear approximation to the boundary can be used. We now describe such a method for two-phase problems in two and three dimensions. Two Dimens iona l Problems There are three main steps for approximating the integrals (3.12) over the finest grid subdivisions for two-phase problems in two dimensions. These are detailed below. Step 1. Divide the Square Cell into Two Triangles. We begin by breaking the square subdomain into two triangles and consider each separately (see, e.g., Figure 3.3). In two dimensions, this subdivision step is optional because it only slightly simplifies the implementation of Step 2. Step 2. Approximate Regions Using Triangles. We next approximate the desired phase with a number of triangular subregions. Details for this approximation method are now given for each of the four possible cases. (Five cases arise if Step 1 is omitted.) Case 0. If none of the corners of the triangle belong to R, then we assume that R and the triangular subdomain do not overlap. No contribution to the Fourier coefficients is made. Case 1. If one corner is in R, then linear interpolation is used to determine a trian-gular approximation to the subregion. For example, consider approximating Chapter 3. A New, Spectral Method 32 the gray phase in AABC of Figure 3.4. Letting U(A), U{B) and U(C) be the function values at points A, B and C, we approximate where the curve crosses the edges of the triangle by points a and b, a = B + \ ~ U{B) (C-B), (3.14) b = A + U(C) - U{B) (C-A). U(C) - U(A) This gives a triangular approximation, AbaC, to the desired subregion which can be used to approximate the contributions to the Fourier coefficients (see Step 3 below). Case 2 . If two corners are in R, then we linearly interpolate to estimate the loca-tion of the interface and break the subregion into two triangles. For example, the shaded region of APQR in Figure 3.5 is approximated by APpR and APrp, respectively. Case 3. If three corners are in R, then we assume that the entire subdomain be-longs to R, and we approximate the integrals (3.12) over the entire subdomain. We seek an estimate of the error produced by this step for a smooth curve1. One source of error occurs when smooth curves are approximated by line segments. By Figure 3.6, this approximation produces an 0{h2) error in the position of the front, since the curvature is independent of h. We also produce errors by replacing the actual front position with the interpola-tion (3.14). To determine this error, we consider the one dimensional analogue of curvature motion given by Figure 3.7. (We expect a similar result to hold in 1 F o r affine veloci ty mot ions , nonsmooth corners may arise from singulari t ies i n the so lu t ion (see, e.g., [42]). E a c h corner can produce an 0{h?) error i n the phase areas. However, because these corners are r ap id ly smoothed away, they typ ica l ly do not affect the overall order of accuracy of the method . Chapter 3. A New, Spectral Method Figure 3.4: A Shape Represented by a Triangle Chapter 3. A New, Spectral Method 34 Figure 3.5: A Shape Represented by Triangles two dimensions.) We seek an estimate of how well the interpolated value, xm, approximates the front position, xm; i.e., we seek an estimate of \xm — xm\. Treating r as a constant, an elementary Taylor series expansion gives us that U(xm) = U(xi) + U'(x1)(xm - X l ) + ^U"(Cl)(xm - Xif where £i £ ( x i , x m ) . Replacing U'(xi) by a one-sided difference, we find U(xm) = U(Xl)+U{f>^~^l)(xm-x0-^//(6)(^2-x1)(xm-x1)+^//(e1)(xm-x1)2 {x2 — xi) 2 2 where £2 £ (xi,x2). Isolating xm and subtracting the interpolated value, xm, gives \xm xm\ W"^2)(x2 - xi){xm - xi) + \U"{ii){xm - Xi) U{x2) - U(xi) < max £e [x 1,3:2] u"(0 h3 U(x2) - U{xi) where h = x2 — x\. Since U(x) = \ + |erf-(^"/F*) ^ 1S e a s y t o show that (3.15) U M - U M = ^ e x p ^ - x ^ ) ] ^ ) ^ ^ ) U"(x) = -2 ^ 3 T 2 Chapter 3. A New, Spectral Method 35 approximately JL = 0(h ) \ 1 /K=radius of curvature Figure 3.6: Errors Approximating Curved Segments Applying these results to Equation (3.15) yields 1 (h2\ ,, \ , \ x m - x m \ < - — max (£ - xm) + h.o.t. Thus, \xm xm = o 'h3 Taking into account both of the contributions to the error, we find that this tri-angular approximation of regions produces an O {^r + h2^ error in the position of the front. S t e p 3. Integrate over each Triangular Subregion. We are now left with the task of adding a contribution 7 Y / / « - ( « * ) C O B C ^ C M (3.16) Ujk to each Fourier coefficient, c,-j, for each triangular subregion, Tk. Expanding the integrand about some point, (£, y), in Tk yields hjk = JtJArea(Tk) cos(iTix) COS(TTjy) Chapter 3. A New, Spectral Method 36 U = l u=o x m x 2 U ( x 2 ) 1/2 Fig. 3.7a. Sharpened Intervals Fig. 3.7b. After a Time, r Figure 3.7: The Interpolation Step Tk + 0(\i2+j2]hA) where Area(Tfc) is the area of triangle Tk. Because the integral term cancels when we expand about the centroid of the triangle2, we choose to approximate each contribution (3.16) by lijk ~ 7jj Area(Tfc) cos^nix) cos(7rjy) (3.17) where (x,y) is the centroid of Tk- This approximation is preferred over the direct evaluation of the integrals (3.16) because it is much faster (it only requires two trigonometric evaluations) and it produces errors which are typically small relative to those arising in Step 2. 2 T h e centroid of a tr iangle w i t h corners A, B and C is A+B+c. Chapter 3. A New, Spectral Method 37 Three Dimens iona l Problems There are also three steps for approximating the contributions to the Fourier coefficients over the finest grid subdivisions for two-phase problems in three dimensions. These are outlined below. Step 1. Divide the Cube into Six Tetrahedrons. We begin by breaking cube-shaped subdomains into six tetrahedrons and consider each separately. For example, the cube in Figure 3.8 would be split into the tetra-hedrons, TABFD,TAEFD,TDHEF,TCBFD,TCGFD and TDHGF- In three dimensions, this subdivision step is highly recommended because it significantly simplifies the implementation of Step 2. Step 2 . Approximate Regions Using Tetrahedrons. We next approximate the desired phase with a number of tetrahedrons. For smooth surfaces3, this step produces an O {^r + h2) error in the position of the front where h is the width of the finest grid subdivision. A n outline of this approximation method is now given for each of the five possible cases. (Nine cases arise if Step 1 is omitted.) Case 0. If none of the corners of the tetrahedron belong to R, no contribution to the Fourier coefficients is made. Case 1. If one corner is in R, then the region is approximated by a tetrahedron. For example, the shaded region in tetrahedron TAB CD of Figur.e 3.9 would be approximated by TabCD. 3Nonsmooth corners may arise from singularities in the solution (see, e.g., [35]). Each corner can produce an 0(h3) error in the phase volumes. However, because these corners are rapidly smoothed away, they typically do not affect the overall order of accuracy of the method. Chapter 3. A New, Spectral Method 38 Case 2 . If two corners are in R, then we linearly interpolate to estimate the lo-cation of the interface and break the subregion into three tetrahedrons. For example, the shaded region of TABCD in Figure 3.10 is approximated by the tetrahedrons Tseg^T^jgh and TshOg-, respectively. Case 3. If three corners are in R, then we represent the shape as the difference of shapes which are treated using Cases 1 and 4. See Figure 3.11 for an example. Case 4. If four corners are in R, then we assume that the entire subdomain belongs to R, and we approximate the integrals (3.18) over the entire subdomain. Step 3. Integrate over each Tetrahedron. For each tetrahedron, T(, a contribution /«« = 7 i * / / / o s ^ x ) c e * , , , ) c o * * , ) dV (3.18) T, W here 1 if i = j = k = 0 2 if exactly two of i,j or k equals 0 4 if exactly one of i,j or k equals 0 8 otherwise must be added to each Fourier coefficient, c ^ . Similar to the two dimensional case, we expand the integrand about the centroid of the tetrahedron4, to arrive at lijke ~ 7ijfcVolume(T^) cos(7rzx) cos(7rjy) cos(irkz) where (x,y,z) is the centroid of Ti and Volume(T^) is its volume. This approxi-mation is preferred over the direct evaluation of the integrals (3.18) because it is 4 T h e centroid of a tetrahedron w i t h corners A, B, C and D is A + B + c + D Chapter 3. A New, Spectral Method 39 tetrahedrons Figure 3.8: Subdividing a Cube into Tetrahedrons much faster (it only requires three trigonometric evaluations) and it produces errors which are typically small relative to those arising in Step 2. 3.3.3 Piecewise Linear Approx imat ions for Junct ions The previous subsections do not explain how to approximate contributions to the Fourier coefficients when junctions occur. Several methods for carrying out such approximations are possible. This subsection describes three of these methods for two-dimensional prob-lems; extensions to three dimensions are straightforward. M e t h o d 1. Trivial Treatment of the Junction. If all corners of the triangular subregion correspond to a different phase (e.g., A ABC of Figure 3.12 or AACD or Figure 3.13) then § the contribution (3.12) that would occur for the entire subregion is added to each of the three sets of Fourier coefficients. Other regions are assumed to contain only two phases and hence are treated according to the discussion of the previous subsection. Chapter 3. A New, Spectral Method 40 Figure 3.10: A Shape Approximated by Tetrahedrons Chapter 3. A New, Spectral Method 41 Figure 3.11: A Shape Represented as a Difference Although the method is somewhat crude, it only produces an 0(h2) error in the phase areas for cells which possess three phases. Thus, the overall order of accu-racy of the method is not degraded since the number of junctions present in two dimensions is independent of h. M e t h o d 2 . Triangular Approximation of Regions. We may also use several triangles to approximate regions that contain junctions. If different phases occur at each corner of a triangular region, then that region can be broken into three triangles which have at most two phases each. These triangular subregions are treated according to the discussion of the previous subsection. For example, AABC or Figure 3.12 can be broken into AbaC, AAah and AABa where a and b are found by linear interpolation, a UX{C)-U2{C) B, b = U2(B)-U2(C) + U1(C)-Ul(B) U2(A) - U3(A) C + U2(B)-U2(C) + U1(C)-U1(B) U3{C)-U2{C) A, U2(A) - U2{C) + U3{C) - U3(A) U2(A) - U2{C) + U3(C) - U3{A) and Ui(X), 1 < i < 3, is the value of Ui at the point X. Figure 3.12: Junction for which all Corner Phases Differ The other possible junction occurs when exactly two different phases are present at the corners of the subdomain (e.g., AABC of Figure 3.13). In this case, we break the region into three triangles, as we would for the two-phase case (see, e.g., Figure 3.13). The triangles that arise contain two phases, and hence should be treated according to the discussion of the previous subsection. M e t h o d 3. Subdividing Regions Containing Multiple Phases. The final approach that we mention is simply to subdivide any region containing more than two phases. After a few iterations, the smallest subregions that arise can be treated by Method 1 to give an accurate approximation of each integral. Our implementations have used the first method (i.e., the trivial treatment), since this straightforward approach typically produces errors which are very similar to those for the other, more involved methods. For example, in the three-phase problem displayed in Figure 4.6, cells containing multiple phases contribute less than 2% of the total error arising from the spatial discretization. Chapter 3. A New, Spectral Method 43 Figure 3.13: Two Phases Represented at Corners 3.4 Refinement Techniques In Section 3.2, a recursive algorithm for subdividing the domain was introduced. We now carry out a more detailed study of the method and introduce a gradual refinement which overcomes certain limitations of the original algorithm. For illustrative purposes, all examples set the width of the coarsest grid to be H = |. Similar results arise for the usual choice of H = -. n 3.4.1 The Or ig ina l Refinement A l g o r i t h m The original refinement algorithm of Section 3.2 is effective for a variety of problems. Application of this method to the two-phase shape of Figure 3.14a, for example, produces a refinement which captures the entire interface at the level of the finest grid subdivision. See Figure 3.14b. For certain smooth regions, however, small slivers can be missed when the algorithm is applied. Consider, for example, a translation by ( 7 5 ^ , -^j of the shape found in Figure 3.14a. Applying the subdivision algorithm to the translated shape gives the Chapter 3. A New, Spectral Method 44 mesh displayed in Figure 3.15a. A magnification of the leftmost part of the shape (see Figure 3.15b) indicates that a small, thin region is missed by the algorithm. We would like at least a crude estimate of the error arising from these slivers. Suppose that the local curvature of the boundary for such a region is K. By Figure 3.16, this leads to an 0(KH3) error in the corresponding phase area. The number of coarse-grid cells that neglect these slivers will be'problem-dependent. For a convex shape, however, only the cells at the extreme top, bottom, left and rightmost parts of the curve need to be considered. As the shape shrinks under mean curvature motion, we expect that the proportion of time that slivers of width 0(H2) arise in any of these four cells will be 0(H). Assuming H = ^ and n = O(-y^) (see inequality (3.11)), the total neglected area over O(^) steps of the method will be 0(T), on average. The original refinement algorithm also produces errors when applied to nonsmooth shapes. Consider, for example, the region displayed in Figure 3.17. Such a shape may arise for an affine velocity motion when a topological breaking occurs. (Note that such topological changes do not occur for mean curvature motion in two dimensions [31, 36].) Applying the original subdivision algorithm to the shape gives the mesh displayed in Figure 3.18a. Clearly, an 0(H2) error in the phase area is produced at the cell containing the sharp corners. This corresponds to an 0(T) error when H = ^ and n is chosen according to (3.11). For the numerical experiments that we have considered, this flaw in the refinement technique produces errors which are no larger than the 0(T) errors which arise from the MBO-method for smooth, two phase problems and are smaller than the (9(A/F) which arise for junctions (see next chapter). Nonetheless, we introduce a more accurate refinement in the next subsection to give us a greater confidence in our results, especially for complicated, non-convex initial regions. A more accurate refinement is also of value when higher order, extrapolated methods are considered (see next chapter). Chapter 3. A New, Spectral Method Fig. 3.15a. Refined Domain Fig. 3.15b. Zoom-in Near a Neglected Section Figure 3.15: Original Refinement can miss Slivers of the Region Chapter 3. A New, Spectral Method 46 H • H T 1/K = radius of curvature y ( i / K ) 2 - ( H / 2 ) 2 approximately K H 2 / 8 Figure 3.16: A Neglected Section Figure 3.17: A Problem with Sharp Corners Chapter 3. A New, Spectral Method 47 Fig. 3.18a. Original Refinement Fig. 3.18b. Gradual Refinement Figure 3.18: Refinement Methods for Corners 3.4.2 A M e t h o d for a Gradua l Refinement As we have seen in the previous subsection, the original subdivision algorithm can miss small pieces of both smooth and nonsmooth shapes. We seek a refinement which captures the entire interface at the level of the finest grid subdivision, even for nonsmooth shapes. To achieve this objective, a gradual refinement was implemented. This method pro-ceeds according to the original subdivision algorithm of Section 3.2, with the following additional consideration: Whenever any cell is refined, check the subdivision level of the neighboring cells. Subdivide neighbors which are two or more levels of refinement coarser. This method accurately represents the narrow, sliver-shaped regions that were missed using the original refinement. By using a fine subdivision in a small neighborhood of the interface, this method even captures the rapid variations' in the front that arise from corners. See Figures 3.19 and 3.18b for examples. Certainly, this gradual refinement produces more cells than the original approach. The order of the number of cells is unchanged, however. To see this, note that cells of Chapter 3. A New, Spectral Method 48 width he = 2~eH, where 0 < I < log 2 (^j form a band at most two cells wide on each side of the interface. The length of each band can be bounded by a constant, K, independent of h (e.g., bands for a convex region are shorter than the perimeter of the domain). Letting n-h be the number of cells of width h, we observe that Total number of cells = n-h + «2/i + • • • + ns. + rijj, 4K . 4K . . AK ~~h 8K < -h+-2-h+-- + T + n2' < — + n2. h Thus, 0(j^ + n 2) cells are required, which matches in order the result for the original refinement. Implementation of this gradual refinement is somewhat more involved than the orig-inal approach because cell neighbors must be found. Many data structures appropriate for this task have been considered [62, 61]. Our implementation defines the grid as a list of vertices (cf. [28]), each of which is described by a data structure: s t r u c t u r e ver tex B E G I N x, y; % Coordinates of the vertex, u; % Function value at this vertex. n, s, e, w; % Pointers to the vertices north (above), south (below), east (right) and west (left) of this vertex, h; % If there is no cell northeast of this vertex, set to (-1). Otherwise, set to the width of the cell northeast of this vertex. E N D Chapter 3. A New, Spectral Method 49 Access of cells and their neighbors is carried out indirectly by traversing their vertices. For example, the refinement of region Q in Figure 3.20a can be carried out as follows: 1. Determine the cell width, h, from the data structure defining vertex, u 4 . Starting at V4, traverse the perimeter of Q by going north, east, south and finally west a distance h. 2. From the traversal in Step 1, it is clear that no vertex should be added to edge (v4,vj). Thus, four vertices, u g , U i o , U n and vi2, are added to the grid as shown in Figure 3.20b. Updates to the data structures for all boldfaced vertices in Fig-ure 3.20b must also be carried out. 3. Finally, we check if neighboring cells need to be refined. Since there are vertices immediately north of vj and v$, the cell north of Q is at most one level coarser than the refined regions. Thus, refinement north of Q is unnecessary. Similarly, no refinement occurs east or west of Q. There is no vertex immediately south of v5, however, so a refinement of the region south of Q is needed. This is carried out by applying steps 1 to 3 to the region northeast of vertex v\. 3.5 F a s t , T r a n s f o r m - B a s e d A l g o r i t h m s The refinements of Sections 3.2 and 3.4 lead to a large number of function evaluations, n-l U(x,y)=- c i j - exp ( -7 r 2 [ i 2 + (i')2]7")cos(7rja:)cos(7rj'y). (3.19) j,j'=o Because these evaluations occur on an unequally spaced grid, a fast Fourier transform cannot be used. Direct evaluation of Equation (3.19) at Nq points, however, is often prohibitively expensive because 0(n2Nq) operations are required. Similarly, evaluation Chapter 3. A New, Spectral Method 0.64 0.62 0.6 0.58 0.56 0.54 0.52 Fig. 3.19a. Refined Domain Fig. 3.19b. Zoom-in near a Difficult Section Figure 3.19: Gradual Refinement Captures the Entire Region V7 Q V6 V3 V . V 2 V7 V 12 V8 V6 V 11 V3 V . V V 9 ! V 2 Fig. 3.20a. Initial Domain Fig. 3.20b. After Subdividing Q Figure 3.20: Refining a Cell Chapter 3. A New, Spectral Method 51 of the Fourier coefficients using Equations (3.12) and (3.17) leads to a Fourier sum of the form Np-1 CJJ' = E ^ cos(njxe) COS(TTj'ye) (3.20) 1=0 where 0 < < n — 1 and (x£,ye) are unequally spaced. Once again, a fast Fourier transform cannot be used, and direct evaluation leads to 0(n2Np) operations. In this section, we consider recent methods for the fast evaluation of (3.19) and (3.20). Specifically, we discuss an unequally spaced fast Fourier transform method [5] that we have applied to the new, spectral method. This transform leads to an algorithm that typically requires only . o(Ilog 2(r)) operations per r-step for the basic MBO-method. For consistency with [5], our discussion will assume that Equations (3.19) and (3.20) have been re-arranged as exponential sums. Specifically, we express Equation (3.19) as n-l U(x,y) = hi' exp(-iirjx)exp(-inj'y) (3.21) j,j'=—n+l where / i i ' = 7 — W I exp ( -7r 2 [ j 2 + (j') 2]r) 133' and 7jj/ is given by Equation (3.13). Equation (3.20) is re-written as en, = ^Real (5 j y ) + ( - l )^Real(40 (3.22) where Np-1 Sjj' = E deexp(-iirjx£)exv(-inj'yi), (3.23) 1=0 JVP-1 Sji' = deexp(-iirj[l - xt])exp{-i7rj'ye), (3.24) 1=0 and 0 < j,j' < n — 1. (Note that it is easy to show that Equation (3.22) is symmetrical with respect to x and y by factoring e~t7TJ = (—l)-7 from Equation (3.24).) Chapter 3. A New, Spectral Method 52 3.5.1 Overview Several methods for the fast evaluation of Equations (3.21) and (3.22) have been de-veloped. Lagrange interpolation has been used to replace function values at arbitrary points by several function values on an equally spaced grid [55, 10]. Taylor expansions have also been used to correct for deviations from an equally spaced grid [70]. Although these methods produce a significant speed-up over a direct approximation of the sum-mations, they are inherently inefficient because they must oversample or compute on a much finer than n x n mesh to accurately evaluate all n 2 Fourier modes. Furthermore, no studies have been undertaken to determine how the accuracy and speed of these methods depends on the oversampling factor and the choice of interpolation [5]. Recent methods do not need a significant amount of oversampling and have been proven to converge quickly. In [5], for example, a method based on multiresolution analysis was developed and implemented that evaluates Equation (3.21) at Nq points in 0(Vglog2Q) +n 2 log(n)) (3.25) operations, and evaluates all n2 Fourier coefficients of (3.22) in o ( /V p l og 2Q )+n 2 l og (n ) ) (3.26) operations, where e is the precision of the computation. Essentially the same bounds are obtained in [24] for an interpolation using Gaussian Bells and in [68] for an algorithm which uses Lagrange interpolation and Green's theorem. In practice, numerical experiments [24, 68, 5] suggest that Beylkin's algorithm [5] is the fastest of the three recent methods. For this reason, we next consider Beylkin's approach for evaluating Equations (3.21) and (3.22). Chapter 3. A New, Spectral Method 53 3.5.2 The Unequally Spaced Fast Fourier Transform In this subsection, we outline Beylkin's algorithm [5] for evaluating Equations (3.21) and (3.22). We have implemented this algorithm in two dimensions; extensions to three dimensions are possible [5]. The version of the algorithm given in [5] requires that grid points must not occur within a distance O of the boundary of the domain [0,1] x [0,1]. In the following implementation, we have scaled and translated points to overcome this limitation. A more elegant and faster version can be derived using periodic extensions of the Fourier coefficients given in Equations (3.21) and (3.22) [5]. Fast Evaluation of Fourier Coefficients Three steps are required for the evaluation of Fourier coefficients using Beylkin's algo-rithm. These are described below for Equation (3.23). Equation (3.24) can be treated by replacing xi in step 1 by (1 — xg). 1. The algorithm begins by projecting an integral representation of the sum onto the span of a number of translated central B-splines. Specifically, we evaluate for 0 < k, k' < 2n — 1. The m^-order central B-spline, ft(m\ can be stably evaluated using the recursion fkk1 = ) p(m\x) | ( m + l) + x — x 1 m m where 1 if < 0 otherwise. Chapter 3. A New, Spectral Method 54 See [7] for an algorithm that evaluates each of the components of / in only 0(m2) operations. 2. Using a fast Fourier transform, we next evaluate 2 n - l 2 n - l Fjj' = E E fkk>e k=0 k'=0 -2irikj/(2n) -2mk'j'/ (In) for — n < j,j' < n. 3. A final scaling gives an approximation to the desired coefficients Sj+s. j'+s. Srf —, Fiji 2 2 ^)(j/(2n))a^)(j'/(2n)) where t=m e=-m a n d - f < j , / < f - l . Summing over steps 1 to 3, we see that a total of 0(m2Np + n2 log(n)) operations are taken. A proof of the operation count (3.26) can be derived using the relationship between m and the precision of the calculation, e [5]. Our implementations set m = 23, since this choice has been experimentally shown to give errors which are similar to those arising from roundoff in double precision calculations (see [5]). Fast Eva lua t ion of Fourier Summations Beylkin also gives a three step algorithm for the fast evaluation of Fourier sums [5]. This method is described below for the evaluation of Equation (3.21). Chapter 3. A New, Spectral Method 55 1. The algorithm begins by modifying the Fourier coefficients according to fp Wkk> f k k ' — A—A— bkb'k where { fkk' for — n + l<k,k'<n — 1, 0 otherwise (m-l)/2 fc=-(m-l)/2 and - 2 n <k,k' <2n- 1. 2. Using a fast Fourier transform and the result from step 1, 2n- l 2n - l k=-2n k'=-2n is evaluated for —2rc < j,j' < 2n — 1. 3. Having completed the pre-processing steps 1 and 2, an arbitrary number of function evaluations may be carried out using U(x,y)= F33,^m\2nx-])f3^\2ny-3'). (3.27) j,j'=-2n Using the algorithm discussed in [7], each of these evaluations requires only 0(m2) operations. Summing over all three steps, we find that O(m2N0 + n2 log(n)) operations are taken to produce Nq function evaluations. This corresponds to the work estimate given in (3.25). Chapter 3. A New, Spectral Method 56 A p p l i c a t i o n to the N e w , Spectral M e t h o d We seek an estimate of how much work arises from an iteration of the new, spectral method when these transform methods are applied: 1. We begin by carrying out the first two steps of Beylkin's method for the fast evaluation of Fourier-sums. This step produces a total of 0(n2 logn) operations. 2. The domain is refined according to Sections 3.2 and 3.4. Since each evaluation of the function, U, can be carried out in 0(m2) operations using Equation (3.27), refinement produces a total operation count of (3.25). 3. We next collect the contributions to each Fourier coefficient into a sum (3.20) in 0(NP) operations. The number of operations to evaluate this sum according to Beylkin's method is given by (3.26). Summing over all three steps, we find that O (Np log2(e) + Ng log2(e) + n2 log(n)) operations are taken where e is the precision of the calculation. Using the fact that O (j^j refined cells arise (see Section 3.4.2), it is clear that Nq = O (j^j and Np — O {j^j. The remaining 0{n2) coarse grid cells may be treated with a fast Fourier transform in 0(n 2 log(n)) operations. Applying these relationships, along with h <C ^, we see that a total of C ( ( £ ) log a(/i) + n 2log(n)) operations arise at each iteration of the spectral discretization of the MBO-method. As we shall see in the next chapter, the basic MBO-method produces an 0{r2) error in the position of a smooth curve at each step of the method. To avoid degrading this accuracy Chapter 3. A New, Spectral Method 57 (see Section 3.3.2), we select h = 0(T) to arrive at o ( i l o g 2 ( r ) ) operations per step. For the case of junctions, we may apply the same considerations to determine that operations are required per step to avoid degrading the overall accuracy of the method. 3.6 Comparison to the Usual Finite Difference Discretization There are several reasons why the spectral method described in this chapter is preferred over the usual finite difference approach. These reasons are outlined below. 1. As has been discussed in Section 3.1, only low frequency modes need to be approx-imated provided r is not taken very small. A large amount of computational work is saved by only treating these low frequency modes. 2. The new, spectral method does not require any time-stepping between tstart and tstart + T. This eliminates a possible source of error and produces large savings in computational work. 3. Local refinement is much simpler to implement for the new, spectral approach because it is done in the context of a quadrature, rather than a discretization of a differential equation. 4. By using a spectral method, the error arising from discretizing the heat equation can be nearly eliminated. This is an attractive feature, because it makes extrapolation in T practical (see next chapter), which in turn allows for larger T-steps. When Chapter 3. A New, Spectral Method 58 larger r-steps are taken, even fewer basis functions are required to solve the heat equation to a given accuracy. 5. The original finite difference algorithm must satisfy (2.7) globally, or part of the front may erroneously remain stationary. By recursively refining near the interface, the new, spectral approach can essentially eliminate this restriction. 6. The new, spectral method also gives an O {^y + h2^j approximation of the location of the front, which is greatly superior to the first order approximation arising for finite differences. As we saw in the previous section, this improved accuracy, in part, explains why operations are needed per step for the basic method. This compares very favorably to the idealized finite difference result for smooth curves, O(y), which was derived in Section 2.4.2. These are indeed formidable advantages for the new, spectral method over the usual finite difference approach. Even when the fast transform-based methods of the previous section are not used, large computational savings are typically observed. To illustrate the performance improvement, consider the motion by mean curvature of the kidney-shaped region displayed in Figure 3.21. Using the new, spectral method 5 and a finite difference approach6, we compare the area lost over a time t — 0.0125 with the exact answer, 0.0125 x IT = 0.0785398 (see [51]). From Table 1, we see that the new, spectral method is adequate for finding solutions to within a 1% error. As we shall see in the next chapter, even more accurate results are practical using the transform methods 5 A direct evaluat ion of the Fourier summat ions was carried out. 6 T h e difference a lgo r i thm uses an adaptive r -s tepping method (see Section 2.4.2) on a un i fo rm mesh. A m u l t i g r i d technique [3] was used to solve the i m p l i c i t equations which arose f rom a backward Eu le r t ime-stepping scheme. Chapter 3. A New, Spectral Method 59 Figure 3.21: A Smooth Interface at Time, t described in Section 3.5. The finite difference approach, however, is impractical when accurate solutions are sought (see Table 2). Numerical tests for the problems described in the next chapter also found that the new, spectral method often requires less than 0.1% of the computational time of the usual finite difference approach. For this reason, the numerical studies in the following two chapters are carried out using the new, spectral method. T h Error Time7 0.003125 2 - 9 4% 0.4 s 0.00078125 2 - i i 1% 8 s At A x Error Time7 1 x 10- 4 2 x 10~5 l 128 1 512 4% 3% 85 s 10341 s Table 1. New, Spectral Method Table 2. Finite Difference Discretization 7A11 t imings were carr ied out on an H P 7 3 5 / 1 0 0 works ta t ion . Chapter 4 Theoret ical and Numer i ca l Studies In this chapter, the order of accuracy of the MBO-method is studied for a variety of two dimensional initial shapes which move according to mean curvature motion. In particular, theoretical and numerical studies are presented to show that the method produces a first order error in r for smooth, two-phase problems. For three-phase problems, an expansion is derived to help explain why junction angles remain within 0(-^/r) of the correct configuration and numerical experiments are presented to show that the basic method demonstrates 0(y/r) errors in actual computations. Numerical studies for initial corners and for singularities in the solution are also presented. Throughout this chapter, extrapolation in r is used to produce higher order methods. Our efforts fall short of a global convergence analysis. What we present instead is a combination of heuristics: Local error analysis coupled with numerical experiments which suggest that the derived local error orders may also hold globally. 4.1 Smooth Interfaces In this section, we consider the application of the MBO-method to smooth curves. Specif-ically, we determine the local order of accuracy of the method and we propose an extrap-olation to provide, a higher order method. 60 Chapter 4. Theoretical and Numerical Studies 61 g(x) Figure 4.1: The Initial Interface 4.1.1 Truncat ion E r r o r Analys is We begin this chapter with a demonstration that the MBO-method is locally first order accurate in r for smooth curves. Specifically, we show that the level set F(x, Z,T) = | for 3F_ dt = A F , F(x,z,0) = < 1 if z > g[x) 0 otherwise where g(x) is the initial interface, gives an 0(T2) approximation to the position of the front after a time, r. Begin by considering a smooth (4 times differentiable) interface which initially passes through the origin, tangent to the x-axis1 as in Figure 4.1. By [45], the value of F(x, z,t) along the z-axis is given by 1 1 tz y1 1 t00 x2 r F(0,z,t) = - =/ e-*<dy + — K 1 2 Jhrt Jo y 4nt J-oo Jo 0 0 _ s i f9^ (*-y)2 e 4t I e 4t dy dx. Expanding two of the exponentials around 0, and replacing g(x) by the first few terms 1For consistency with [45], the x-z plane is considered rather than the usual x-y plane. Chapter 4. Theoretical and Numerical Studies 62 of the Taylor series gives, F(0,z,t) = L= f (l-^-)dy + V ' ' ; 2 y/Wth \ At) y^ f°° f / 6 4' J-oo JO 1 f°° / • | 3 < 2 » ( 0 ) : C 2 + i S ( 3 ) ( 0 ) 3 ; 3 + i 9 ( 4 ) ( 0 ) a g 4t 47Tt J-oo ( * - y ) 2 4* <fy + h.o.t. 2 JAfrt V I2t 1 + 1 r°° - — / e 4t Ant J-oo y + 12t y= b( 2) ( 0 ) s 3 + ± P < 3 ) (0)rr 3 + ± f l<«> (0)x4 (fa + h.o.t., y=o - - -±= ( * _ £ - ) + 2 12*/ i f 1 /-00 _s£ e it -co y + -3z 2y + 3zy 2 - y 3 12t dx + h.o 3/=0 Noting that terms which are odd powers of x cancel in the integral, and writing for y ( n ) (0) we find mz,t) = l -/ ¥ t \ 12* , / T 1 y 0 0 -;— / e 4 t 47rt J-oo ^ y ( 2 ) ^ 2 + ^9™x* + \z2g™x2 + 3z (y^x2)2 - (\g^x2y 12t dx + h.o.t Applying integration by parts and the well-known identity, e-»dt = f-J-oo V a yields I z zc F(0,Z,t) = --^= + —-y= + 2 2Vnt Uiy/nt (4.28) y ( 4 ) , * (^ ( 2 )) 21 4\ + m 8ty^F + h.o.t. To find where the interface intersects the z-axis, set F(0,z,t) = \ to obt am: 4! 16* 21 t2--(gW)3t2 + h.o.t. Chapter 4. Theoretical and Numerical Studies 63 Substitution of the first order approximation of z into the right hand side gives a higher order estimate which simplifies to *(*) = {9{2})312 + 9{2)t - \ (g{2))3t2 + 7 > * 2 + \ (s(2)) V - \ (g(2))3t2 + 0(t% Thus, Z(T) = 5r ( 2 )r + T2 + 0(T3). Using a known result (see [54]) for the exact position, z*. zt — TT77 . \ 2 it is straightforward to show that z*(T) = g ^ r + [2 g{4) - ( 9 { 2 ) y T2 + 0(T3). Thus, the MBO-method is locally first order in r for smooth curves. (4.29) (4.30) 4.1.2 E x t r a p o l a t i o n From Equations (4.29) and (4.30) we expect that extrapolation in r can be used to produce higher order accurate results since the error, \ (g{2))3 r2 + 0 ( r 3 ) varies smoothly with r and depends only on the local properties of the curve. Assuming that the sharpened Fourier coefficients at time t are given by {cij(t)}, an extrapolated method for smooth curves may be obtained as follows: Chapter 4. Theoretical and Numerical Studies 6 4 Extrapolation I 1. Carry out one step of the MBO-method with a step size 2r to obtain the sharpened coefficients, c2j(t + 2 T ) . 2. Carry out two steps of the method with a step size r to obtain sharpened coeffi-cients, c-j(t + 2T ) . 3. Set ClJ(t + 2r) = 2cT.(i + 2r) - c2J(t + 2r) to obtain the extrapolated coefficients. • We hope to obtain second order accurate results in r using this extrapolation. To un-derstand why, suppose Mir) produces an 0(T) approximation to some unknown quantity M*. Assuming the error for the approximation of Mir) to M* can be expressed as M* = M(T) + KT + 0(T2) where K is independent of r, we find that the extrapolation 2M(r) - M(2r) = M* + G(T2) is second order accurate. See [32, 34, 37] for treatments of extrapolation for initial value problems. The extrapolation process may be repeated, of course, as in Romberg's integration to obtain methods which we expect are even higher order accurate. See, for example, [15]. 4.1.3 Numerical Experiments From the previous two sections, we expect that the basic and extrapolated methods will give globally first and second order convergence rates, respectively. To experimentally Chapter 4. Theoretical and Numerical Studies 65 verify these convergence estimates, we consider the motion by mean curvature of a shrink-ing ellipse with principle axes 0.8 and 0.6. (Similar results are obtained for the shapes of Figure 3.21 (cf. [60]) and Figure 3.14a.) Using the basic method and Extrapolation I, the area lost over a time t = 0.025 was compared to the exact answer2 for several r . The results from a number of experiments are given in Table 3, below. Basic Method Extrapolation I T Error Conv. Rate3 Error Conv. Rate 0.00125 -8.31e-04 1.0 -4.05e-05 -0.000625 -4.11e-04 1.0 7.71e-06 2.39 0.0003125 -2.05e-04 1.0 2.02e-06 1.93 0.00015625 -1.02e-04 1.0 4.97e-07 2.02 Table 3. Extrapolated Results for Two Phases These results support the conclusion that the MBO-method is first order in r and suggest that extrapolation can be used in conjunction with the new, spectral method to produce higher order results. The extrapolated errors in Table 3 change sign (and the convergence rate fluctuates) because the first step of the method uses the characteristic set of the initial shape in the diffusive step, whereas later steps use extrapolated results which can take on values between —1 and 2 (see, e.g., Figure 4.2). Clearer estimates of the convergence rate are obtained by computing the errors that arise away from t = 0. For example, Table 4 gives the error in the area lost from time t — 0.02 to t = 0.025 for several values of r . 2 T h e exact result is easily derived using Equa t ion (2.3). 3 I f the error for a step of size r is ET, then we estimate the convergence rate as l o g 2 Chapter 4. Theoretical and Numerical Studies 66 Figure 4.2: Extrapolated, Semi-Discrete Result T Error Conv. Rate 0.00125 9.66e-06 -0.000625 2.28e-06 2.08 0.0003125 5.53e-07 2.04 0.00015625 1.34e-07 2.04 Table 4. Extrapolated Method for Two Phases These results support the conjecture that the extrapolated method is second order in r. 4.2 Nonsmoo th Boundaries In the previous section, we demonstrated that the MBO-method produces a first order er-ror in r for smooth curves. We now consider the application of the method to nonsmooth initial shapes. In order to determine the form of the error that arises for nonsmooth cor-ners, we consider a fixed r. To produce an optimized code, however, variable sized steps in r should be considered. Consider, for example, the motion by mean curvature of the initially nonsmooth interface given in Figure 4.3. Using the basic method and Extrapolation I, the area lost Chapter 4. Theoretical and Numerical Studies 67 0.2 0,3 0.4 0.5 0.6 0.7 0.8 Figure 4.3: Initially Nonsmooth Interface at Times, t over a time t — 0.0125 was compared to the exact answer2 for several r . The results over a number of experiments are given in Table 5, below. Basic Method Extrapolation I T Error Conv. Rate Error Conv. Rate 0.0015625 6.40e-03 0.71 -1.98e-03 1.16 0.00078125 3.75e-03 0.77 -8.94e-04 1.15 0.000390625 2.15e-03 0.80 -4.23e-04 1.08 0.0001953125 1.21e-03 0.83 -2.06e-04 1.04 Table 5. Basic and Extrapolated Results for a Nonsmooth Shape These results suggest that the extrapolation is first order in r . To obtain a clearer understanding of the form of these errors, let c(r) be the error which arises for the basic method using a step size r. Assuming that the extrapolated result is indeed first order, i.e. 2e(r) - e(2r) = cr + O(T) for some constant c, we write Chapter 4. Theoretical and Numerical Studies 68 ^Hr)]-Y = - c + o(l), Integrating with respect to r , we find that e(r) = -c r log( r ) + o(rlog(r)), = 0(r log(r)) . Thus, we suspect that the basic method is order r log(r) for the case of nonsmooth initial corners4. To derive a higher order method for nonsmooth curves, we assume that the sharpened Fourier coefficients at time t are given by {cij} and carry out Extrapolation I twice to eliminate the leading order error term. This repeated extrapolation produces the follow-ing method: Extrapolation II 1. Carry out one step of the MBO-method with a step size AT to obtain the sharpened coefficients, cfj(t + AT). 2. Carry out two steps of the MBO-method with a step size 2r to obtain the sharpened coefficients, c2J(t + AT). 3. Carry out four steps of the method with a step size r to obtain sharpened coeffi-cients, cjj(t + AT). A. Set + AT) = AcTl3{t + AT) - Ac}j{t + AT) + c%{t + AT) 4 Note that if the order is indeed 0(rlog(r)) then the quantity measured as convergence rate gives 1 + l o g 2 ( l + (log 2 r ) - 1 )—*• 1 as r -> 0 which is consistent with the results in Table 5. Chapter 4. Theoretical and Numerical Studies 69 to obtain the extrapolated coefficients. • Applying this extrapolated method to the nonsmooth problem yields the results given in Table 6. T Error | Conv. Rate 0.0015625 4.91e-04 -0.00078125 ' 2.60e-04 0.91 0.000390625 7.46e-05 1.80 0.0001953125 1.89e-05 1.98 Table 6. Results for a Nonsmooth Shape Using Repeated Extrapolation These results show that repeated extrapolation can be very effective, even when nons-mooth initial corners are present. Indeed, Table 6 suggests that an approximately second order method arises from the use of Extrapolation II. Topological mergings and breakings in three dimensional or affine flows can also pro-duce nonsmooth corners. For this reason, we expect no better than an 0(Tlog(r)) error to arise in these situations. To obtain a deeper understanding of the form of these errors, further studies are needed since an additional dependency on nonlocal properties of the interface occurs. 4.3 Singularities in the Solution as Regions Disappear To conclude our study of the MBO-method for two-phase problems, we consider the method's treatment of regions as they shrink to points and disappear under mean curva-ture motion. Similar to the case of nonsmooth corners, this motion causes condition (2.6) to be violated. For this reason, we suspect that an C(rlog(r)) error may arise whenever the basic method is applied. Chapter 4. Theoretical and Numerical Studies 70 To test if an 0(r log(r)) error occurs in practice, we consider the motion by mean curvature of a shrinking ellipse with principle axes 0.2 and 0.25. Using the basic method, the disappearance time TT was compared to the exact answer, TE = 0.025, for several r . A n extrapolation in the disappearance time, 4TT — 4T2r + T 4 t , was also computed. The results for a number of experiments are given in Table 7, below. T Conv. Rate (\rpr — 4T2T + T 4 T ) — TE 0.005 -3.15e-03 - 1.49e-03 0.0025 -1.90e-03 0.73 -3.16e-04 0.00125 -1.12e-03 0.76 -1.85e-05 0.000625 -6.42e-04 0.80 8.34e-06 Table 7. Basic and Extrapolated Errors for a Shrinking Ellipse These results show a slow tendency upwards to first order for the basic method and indicate that extrapolation can be very effective as phase regions disappear. (Unfor-tunately, we were unable to determine a clear estimate of the convergence rate for the extrapolation.) Similar to the case of nonsmooth corners, these results are suggestive of an order r log(r) error for the basic method. 4 . 4 Junct ions in Two Dimensions In this section, we consider the application of the MBO-method to triple junctions in two dimensions. Specifically, we present numerical experiments to estimate the order of accuracy of the basic method and propose an extrapolated method for improved accuracy. We also derive asymptotic results which explain the stability of junction angles and suggest a source of the 0(^/r) error which arises in numerical experiments. Chapter 4. Theoretical and Numerical Studies 71 4.4.1 E r r o r Analys i s Under mean curvature motion, two dimensional triple junctions form a stable 2f-^3z-2f angle configuration (see, e.g., [13]). In this subsection, we show that each step of the MBO-method produces an 0(y/r) error in the junction angles which is rapidly dissipated during subsequent steps. We now derive an expansion for the angles of a two dimensional triple junction after one step of the MBO-method. We shall assume that each angle approximates | 7 r , initially. The In i t ia l Junc t ion We begin by orienting a polar coordinate system so that the junction angle furthest from | 7r is centered about 9 = 0. We denote the initial interfaces by £o,£i and £ 2 and the initial regions by Ro, Ri and R2 as in Figure 4.4. To represent the small deviations from the 2 ^ - 2 ^ - 2 ^ junction configuration we define e = l£0£i -3 3 3 2TT ce L£x£2 3 ' 2TT where L£{£j is the angle between and £j. Since l£0£i + L£\£2 + l£2£0 = 2ir it is easy to see that 2ir . 27r _ | || | = k i l l + c i . 6 Hence, the constant c satisfies — 1 < c < 0. In order to carry out our expansions, we want an expression for each interface, L£dx 2TT ~ T > y L£Q£\ 2TT i£2£o -2TT ~ T > T 4 = { ( r A ( r ) ) : r > 0 } Chapter 4. Theoretical and Numerical Studies 72 Figure 4.4: The Initial Junction for some function, 0e,(r)- Using the above definitions it is straightforward to show that M O = -\n-\t+\KQr + 0{r% M O = ^ + i £ + I « i r + (9 (r 2 ) , 0l2 (r) = TT + (c + t + i / c 2 r + Q(r2) where is the curvature of line at the origin. Chapter 4. Theoretical and Numerical Studies 73 A p p r o x i m a t i o n of Ui We now want to estimate U = (Uo, Ui, U2) after a time, r . At t = 0, we know that 1 if (r, 9) e Rl Ui(r,e,0) = 0 otherwise, for 0 < i < 2. The Green's function representation of Uo(r,9,r) gives \ /-co U0(r,9,r) = 4lTT JO 6'AR) ( [RcosU) -rcoa(9)Y + [RsmU) - r s m ( 0 ) r \ „ . . . „ exp — — — -J— I R dip dR, AT J-exp (-£) Texp f - f ) /"•'"'exp (^*^>) * # ^ 47TT \ AT J Jo \ 4T J Je(o (R) \ 2 r J Replacing the exponential in the inner integral by its series and integrating term by term yields, 1 / r2\ 0 0 T 1 / r R \ n ( R2\ retAR) AlTT 1 47TT exp , 4 r ) 1 7? exp *i2 n! V2-4r (4.31) 1 (rR + 4\2T +E(r,9,r 0tl(R)-9eo(R)+1- sin (2 (0tl (fl) - *)) _ I sin (2 (9£o(R) - 9)) dR where En(r,9,r) n—3 1 ( exp - , ATTT \ AT ) We seek an estimate of En(r,9,r). Noting that 8e0(R) [6 " ' cosn(<j>-9) d<j> < 2?r, Je. Chapter 4. Theoretical and Numerical Studies 74 it is clear that < < Assuming ^ < 1, we find r V / r2\ T e x p " 4 7 ' E{r,0,T) = O\[-= ' T , Similar to the above result, the second and higher order terms of 9^ (R) make an O yr + ^ - + (^77) contribution to Uo{r,9,r) in Equation (4.31). Thus, (7„M,r) 47T ri2 + 27 ( R2\ (2 1, Rexp I - — I I - 7 r + e + - ( « i - «o)# sin ( - 7 T H—e H — K i i ? — 9\ — sin ( ——TT — —e + -K0R — 6 1 3 2 2 / V 3 2 2 1 /7 \R + 4 V2r 2 • 1 . /2 -7T H— sin -7T — 29 3 2 V3 1 / 2 - sin — 7 T - 20 2 V 3 dR M 1 E R " 2 f r We now expand the integrand (using Maple [18]) and apply integration by parts with the well-known identity r°° ( R2^ / exp Jo \ a , dR = -Va7r to obtain £ / o M,r) 1 1 I, . Fr ry/3cos(9) o + 7Te + 7(»i - « ° ) \ / - + A 3 2f 4 V 7T 4-V/T7T Chapter 4. Theoretical and Numerical Studies 75 1 \/3 e + — ( « i - K0)r cos(8) + —— (/ci + Ko)^ sin(0) -f- ,—r cos(0) 47r 47r o v / 7 r r + ^ ( - W r ' + 0 ( r + ^ + ( ^ ) " ) . Re-writing in Cartesian coordinates yields U0(x,y,r) = ~ + ~e+ \(KX - K0)J^ + y^=x (4.32) o ZTT 4 V 7T 4yT7T 1 . . \/3 . x e "v/3 / n 2 A + T ~ ( K l - + -7— (« i + «o)y + Q ,—x + —— [x -y ) Air 4TT S A / T T T 1O7TT V ' +0 (T + e- (x2 + y2) + (x2/r + y2/r) f) . Expansions for Ui and C/2 may be obtained via rotations of (4.32) to give 1 1 1 [Y 3y — \^3x 1 -v/3^1 Ui{x,y,r) = - + —ce + -{K2 - « i ) y ~ + g ,— — (2/e2 + KI)Z y 0 Z7T 4 V 7T Q\/T7r 47T 47T 3 + 4c y/3 3 \/3 / 2 2 \ ==ea; 7=eV — ~ T ~ — x y + 777;— (y — 2; ) 16A/7TT 16y 7TT 1 6 7 T T 3 2 7 T T V ' +0 (T + -r (x2 + y 2 ) + (X2/T + y2/r) f) . U2{x,y,r) = 1 - U0(x,y,T) - Ui(x,y,r). Angle Expansions We now seek expansions for the angle configuration of the junction after a time r . Begin by letting lT,l\ and £ 2 be the MBO-approximations to the branches of the triple junction after a time T . T O approximate the angle between IQ and £[, we require the location of the triple junction at time r. This can be found by solving the system U0(xT,yT,T) = U1(xT,yT,r), U0{xT, y T , r ) = <7 2 (x T ,y T ,T) for (xT, y T) to give xT = - 2 y ^ e - -J=(KI ~ K o ) r + O(er) + O ( r f ) , Chapter 4. Theoretical and Numerical Studies 76 yT = -l(2c+l)^e+l(K0-2K2-rK1)T-rO{eT) + O(rl) . 6 V 7T 6 v ' Our next task is to find the slope of £Q at the triple junction. Since £Q is given by {(x,y) : r / 0 (x,y,r) = U2(x,y,r)}, the slope of £rQ is given implicitly by Dx [U0(x, y, r) - f/2(x, y, r)] = 0. Solving for the slope at the triple junction, (x T ,y T ) , yields 2 { 2\ Pr~ m0 = -V3 + - (c - 1 + (c + 2 ) - ) e + 2 (2K0 + «I + K2) \ / ~ + C?(r) + C (e2) . Similarly, the slope of £\ at the triple junction is given by mi = V3 + % (c + 2 + (c - 1)-) t + 2(K0 + 2/ci + « 2 ) \ / ^ + O(r) + 0 (e2) . Thus, the value of the desired angle is given by L£rG£\ = 7r — arctan ( — — ) , VI + m 0 m i / = h + {l-^^l^-Ko) + ° { T ) + ° (e2) • ( 4 - 3 3 ) Similarly, L£\£\ = 2-n+(^-^ce+1-(K2-K1)]f^+0(T) + 0(e2), L£V\ = \^-{^-\){1 + C)C + \ ( K 0 - K 2 ) ^ + O{T) + O{C2). Thus, each step of the MBO-method produces an 0(y/r) error in the junction angles which is rapidly dissipated during subsequent steps. Summing up such contributions5 over many r steps, we expect to obtain a rapidly converging geometric sum which gives rise to an 0(\/T) error in total. This is an interesting result because it gives an explana-tion for the stability of junction angles and suggests a source of the 0(y/r) error which arises in numerical experiments (see next section). 5 T h i s s u m m a t i o n step is non-rigorous because it assumes, among other things, that KQ, K\ and K2 are bounded independent of T. Chapter 4. Theoretical and Numerical Studies 77 1 I y y o l 0 1 x x t = 0.0 t = 0.04 Figure 4.5: A Smooth Three-Phase Problem 4.4.2 N u m e r i c a l Exper iments We have just seen that the MBO-method is expected to produce an 0(y/r) error in the angle configuration of junctions. For this reason, we suspect that the basic method may produce an 0(y/r) error when junctions are present. To test if this is indeed the case, we consider the motion by mean curvature of the three-phase problem given in Figure 4.5. Using the basic method, the area lost for the central region, A T , over a time t = 0.04 was compared to the exact answer6, 0.0837758, for several r . Because an 0(y/r) error seems plausible from our asymptotic results, an extrapolation in the area, (y/2AT — A2T^J , was also computed to eliminate the conjectured leading order error term. The results for a number of experiments are given in Table 8, below. 6 A p p l y i n g the V o n N e u m a n n - M u l l i n s parabol ic law [51], gives us that the area of the central phase obeys Chapter 4. Theoretical and Numerical Studies 78 T Error in AT Conv. Rate Error in ^7 [y/2AT - A 2 T ) Conv. Rate 0.00625 2.60e-03 0.67 -1.08e-03 1.33 0.003125 1.70e-03 0.61 -4.63e-04 1.22 0.0015625 1.14e-03 0.57 -2.12e-04 1.13 0.00078125 7.79e-04 0.55 -1.00e-04 1.08 0.000390625 5.37e-04 0.54 -4.86e-05 1.05 Table 8. Basic MBO-method for Three Phases These results support the conjecture that the MBO-method is 0(y/r) for the case of junctions and suggest that extrapolation can be used in conjunction with the new, spec-tral method to produce higher order estimates of quantities of interest such as phase areas. Accurate, extrapolated estimates of the disappearance time of the smallest phase for a more complicated three-phase problem (see Figure 4.6) have also been determined in this manner [60]. Using the basic method, the disappearance time, TT, was compared to an estimate of the exact answer7 for several r . A n extrapolation in the disappearance time, (V2Tr — T 2 t ) , was also computed. The results for a number of experiments are reported in Table 9, below. T Error in TT Conv. Rate Error in (v / 2T T - T 2 T ) Conv. Rate 0.04 0.16435 - 0.07782 -0.02 0.12103 . 0.44 0.01645 2.24 0.01 0.08802 0.46 0.00833 1.38 0.005 0.06354 0.47 0.00444 0.91 Table 9. Results for the Disappearance Time of a Phase Region 7 T h i s result, T = 0.33051, was obtained using B r i a n Wet ton ' s front t racking code; see [14]. Chapter 4. Theoretical and Numerical Studies 79 t = 0.0 * = 0.6 < = 1.2 * = 1.3 Figure 4.6: Evolution of a Junction Through a Singularity These results are suggestive of an 0 ( A / T ) or 0{yjr log(r)) error for the basic method. We also find that a marked improvement in the error occurs when extrapolation is used. An extrapolated method for approximating the entire interface is also sought. Assum-ing that the sharpened Fourier coefficients at time t are given by {cij(t)}, an extrapolated method for junctions may be obtained as follows: Extrapolation III 1. Carry out one step of the MBO-method with a step size 2r to obtain the sharpened coefficients, c2J(t + 2r). Chapter 4. Theoretical and Numerical Studies 80 2. Carry out two steps of the method with a step size r to obtain sharpened coeffi-cients, cjAt + 2r). 3. Set Cij(t + 2r) = V2cUt + 2r)-c2J(t + 2T) y/2-l to obtain the extrapolated coefficients. • Using this extrapolated method, approximations to the three-phase problem given in Figure 4.5 were computed. A comparison of the exact area change for the central region with the computed values yields the results given in Table 9. T Error Conv. Rate 0.00625 -1.10e-03 1.25 0.003125 -5.07e-04 1.11 0.0015625 -2.56e-04 0.98 0.00078125 -1.38e-04 0.89 0.000390625 -7.76e-05 0.83 Table 10. Extrapolated Method for Three Phases Based on these results, it is not clear if the extrapolated algorithm possesses a higher order of accuracy than the basic method. Nonetheless, the extrapolated method may be of practical use since it typically results in a marked improvement in the error (cf. Tables 8 and 10). 4.5 S u m m a r y Based on a number of two dimensional studies, we have determined the dominant local error term of the MBO-method for certain classes of problems and have suggested ex-trapolations for improved accuracy. We have observed corresponding behaviour of the Chapter 4. Theoretical and Numerical Studies 81 error in selected actual computations. These results are summarized in Table 11 below. Type of Interface Conjectured Error Order Extrapolation Choice Two Phases: Smooth T Extrapolation I Nonsmooth T l o g ( r ) Extrapolation II Singularities rlog(r) Extrapolation II Junctions Extrapolation III Table 11. Summary of Theoretical and Numerical Studies in 2-D For the remainder of this thesis, we focus on the approximation of three dimensional and nonlocal models of curvature-dependent motion. Chapter 5 Numer i ca l Exper iments and Visua l i za t ion In this chapter, we apply the MBO-method to a number of three dimensional problems. In particular, we evolve surfaces with junctions according to mean curvature motion and visualize the results. A n extension of the MBO-method to a nonlocal model that preserves phase areas is also proposed and studied. To obtain results at an affordable cost, all problems are solved using the new, spectral method rather than with a finite difference discretization. 5.1 Three Dimens iona l , Two-Phase Problems We begin this chapter by considering the motion by mean curvature of two-phase prob-lems. Results developed here will also be used in the next section for the visualization of three dimensional junctions. 5 .1 .1 V i sua l i za t ion To obtain a clearer understanding of the motion by mean curvature of surfaces, we seek a method for visualizing our results. This section describes a simple approach for generating movies of evolving shapes using Matlab [46]. See [6, 44] for discussions on more advanced visualization techniques. Our principal task is to construct each frame of the movie. A relatively straightfor-ward approach for generating surfaces with diffuse, ambient and specular lighting effects is to use Matlab's surfl(X,Y, Z) command [46]. This command produces a surface by 82 Chapter 5. Numerical Experiments and Visualization 83 interpolating between the points given by the matrices X, Y, and Z. Although it is pos-sible to display small portions of a surface, we prefer to use larger segments. This avoids the shading variations that occasionally arise between segments when the full surface is displayed using a number of smaller fragments. Certain smooth shapes are easy to represent as an m x m matrix. Suppose, for exam-ple, there is a point Q in the region which has a direct line of sight to the entire surface (e.g., point Q in Figure 5.1, but not point Q). Then, the following algorithm can be used to represent the shape: Visua l i za t ion I Select a Q satisfying the "line of sight" property listed above. (Currently Q is user selected.) For 0 < i , j < m — 1, Set Pij = (Xij,Y{j, Zij) equal to the intersection of the surface with the line £, where £ passes through Q and has an azimuth or horizontal rotation of and a vertical elevation of - ^ V — IT. This approach has been used to represent a variety of smooth surfaces (e.g., Figure 5.2). More complicated shapes have also been considered by dividing the shape into subregions and treating each separately (e.g., Figure 5.3). Having constructed an appropriate m x m matrix, we need only call surflQ and specify a shading model to display the surface. Matlab provides for piecewise constant and Gouraud (piecewise bilinear) shading. For the surfaces we have considered, superior results arose from Gouraud shading (see, e.g., Figure 5.4). Chapter 5. Numerical Experiments and Visualization 84 Figure 5.1: From Q the Entire Curve is Visible Figure 5.2: A Matrix Representation of the Surface Chapter 5. Numerical Experiments and Visualization 85 Figure 5.3: Splitting a Shape into Easily Parameterized Portions Fig. 5.4a. Piecewise Constant Shading Fig. 5.4b. Gouraud Shading Figure 5.4: Piecewise Constant and Gouraud Shading of the Surface Chapter 5. Numerical Experiments and Visualization 86 5 . 1 . 2 Numerical Experiments We now report on several experiments for two-phase mean curvature motion. Throughout this section, a piecewise linear fit to the interface is used (see Section 3.3.2) with a finest cell width of h = rib-Collapsing Sphere From the previous chapter, we saw that the basic MBO-method was first order in r for smooth, two dimensional problems. To test if this result also holds in three dimensions, we consider the motion by mean curvature of a collapsing sphere with initial radius 0.4. Using the new, spectral discretization of the MBO-method the volume lost over a time t = 0.03 was compared to the exact answer, 0.167T, for several r . The results for a number of experiments are given in Table 11, below. T Error Conv. Rate 0.01 0.0600 -0.005 0.0258 1.22 0.0025 0.0119 1.12 0.00125 0.00582 1.03 Table 11. The Basic MBO-Method for the Shrinking Sphere These results suggest that the basic MBO-method is first order in r for smooth surfaces without junctions. Thin-Stemmed Barbell Examples involving topological changes are also naturally handled by the method. For example, Figure 5.5 displays the motion of a thin-stemmed barbell using a step size, Chapter 5. Numerical Experiments and Visualization 87 r = 0.0004. From these images, it is clear that the center handle pinches off to form two pieces. As expected from [39], these convex shapes become nearly spherical as they disappear. Thick-Stemmed Barbell A wider stem can produce a qualitatively different motion. For example, Figure 5.6 displays the motion of a thick-stemmed barbell using a step size, r = 0.0005. From these images, we see that no topological changes arise, and that the shape eventually becomes ellipsoidal and more spherical as it disappears. 5.2 Junctions in Three Dimensions In this section, we consider the motion by mean curvature of surfaces with junctions. Specifically, we report on some numerical experiments and give a visualization method for certain multiple phase problems. 5 .2 .1 Visualization The visualization of multiple phases can also be carried out in a straightforward manner for certain problems. Suppose, for example, that a clear phase surrounds several opaque regions. Then the following algorithm can be used to display the visible surfaces: Visualization II For each opaque phase, 0 < j< m, Draw the surface Uj(x,y, z) = Uc\ear(x, y, z) before sharpening using Matlab's surfl() command and Visualization I. Chapter 5. Numerical Experiments and Visualization 88 Chapter 5. Numerical Experiments and Visualization 89 t = 0.0150 t = 0.0200 Figure 5.6: Thick-Stemmed Barbell Moving by Mean Curvature Motion Chapter 5. Numerical Experiments and Visualization 90 Figure 5.7: Composition of a Junction This method draws each of the visible surfaces between the opaque and clear regions. It also constructs several artificial surfaces which are not displayed by surfl() since they are hidden. See, for example, Figure 5.7. It is an interesting fact that these artificial surfaces are often nonsmooth since the corresponding values of Uj and c7ciear are zero to within the tolerance of the line search algorithm used in Visualization I. Visualization II has been applied to a variety of surfaces with junctions. Indeed, each of the movie frames displayed in the next section was constructed using this approach. 5.2.2 N u m e r i c a l Exper iments We now report on experiments for the motion by mean curvature of surfaces with junc-tions. Throughout this section, a trivial treatment of the finest subregions was used (see Section 3.3.1) with a finest cell width of h = yfg • Three Phase Example From Section 5.1.2, we saw that the MBO-method naturally treats two-phase problems in three dimensions. Multiple phase problems have also been evolved using the method. Chapter 5. Numerical Experiments and Visualization 91 For example, Figure 5.8 displays the motion by mean curvature of a cylindrical, three phase shape using a step size T = 0.0005. From these images, we see that the central M u l t i p l e P h a s e E x a m p l e The previous example considered the motion by mean curvature of three phases meeting at 120-120-120 degree angles. The evolution of four-phase junctions may also be studied using our new, spectral discretization of the MBO-method. For example, Figure 5.9 displays the motion of a spherical four-phase shape using a step size, r = 0.0004. From these images, we see that the four-phase junction is stable under mean curvature motion, as is expected from experimental studies of recrystallized metal [38]. 5 . 3 A N o n l o c a l M o d e l The bulk of this thesis has been concerned with the motion by mean curvature of inter-faces. We now consider an extension of this motion to a nonlocal model. Consider a collection of disjoint interfaces, which separate a number of regions. If we evolve these interfaces using a normal velocity where n(x,t) and Ka(t) are the mean curvature and average mean curvature1 respec-tively, then we obtain a nonlocal volume-preserving flow [57]. This motion is of interest 1Here, the average mean curvature is given by blue phase pinches off to form two spherical regions, which eventually disappear. (5.34) where |r,j is the perimeter of the interface of r,-. Chapter 5. Numerical Experiments and Visualization 92 t = 0.0220 t = 0.0260 Figure 5.8: Three Phase Example Moving by Mean Curvature Motion Chapter 5. Numerical Experiments and Visualization 93 t = 0.0144 t = 0.0180 Figure 5.9: Four Phase Example Moving by Mean Curvature Motion Chapter 5. Numerical Experiments and Visualization 94 because it arises as a nonlocal model for binary alloys [57] and because it gives a possible smoothing for generating multiscale representations of planar curves [64]. We now develop a simple modification of the MBO-method to treat this nonlocal model in two dimensions. An extension to the three dimensional case is straightforward. We begin by noting that an approximation to the desired motion (5.34) is obtained at each step if we track the level set \ " \Ka{t){^ ( 5 " 3 5 ) rather than the usual level set of one-half [45]. Thus, phase areas remain approximately constant provided we follow the contour (5.35). Since U is continuous at any time t before sharpening, the area enclosed by a level set c, A(c,t) = Area({x : U(x,t) < c}) is a strictly increasing and hence invertible function of c. Thus, the desired contour (5.35) can be approximated by the level set that preserves phase areas. Applying this result to the MBO-method gives a simple algorithm for computing solutions to the nonlocal model: Non loca l Curva ture A l g o r i t h m To obtain an approximation to the nonlocal curvature model (5.34), we may carry out the MBO-method for two phases given in Section 2.1 using the following as a replacement for step (3): (3a) Find the level set that preserves phase areas; i.e., determine the value c satisfying A(c,t) = A ( i o ) . (5.36) Solving (5.36) for c may be accomplished by a variety of line search algorithms. For example, [11] gives an efficient and reliable approach based on a combination Chapter 5. Numerical Experiments and Visualization 95 0 0.2 0.4 0.6 0.8 1 Figure 5.10: A Test Problem for the Nonlocal Curvature Algorithm of bisection, secant and inverse quadratic interpolation methods. Note that this step is relatively simple to implement since the area of the phase we are following is given by the Fourier coefficient coo(i). (3b) Carry out step (3) of the MBO-method using the contour c rather than the usual choice of one-half. To test how well this method approximates the nonlocal model, we consider the motion of two circles with initial radii 0.2 and 0.15 (see Figure 5.10). Using the nonlocal curvature algorithm, the area of the smaller circle after a time t = 0.02 was compared to the exact answer, 0.044508, which was found by integration of (5.34). The results for a number of experiments are given in Table 12, below. r Error Conv. Rate 0,004 -0.0105 -0.002 -0.00389 1.43 0.001 -0.00182 1.10 0.0005 -0.00104 0.80 Table 12. Results for the Nonlocal Curvature Algorithm Chapter 5. Numerical Experiments and Visualization 96 0 0.2 04 0.6 OB 1 0 02 04 06 0.8 1 t = 0.00875 t = 0.01250 Figure 5.11: Nonlocal Model Which Preserves Area These results indicate that the error decreases with r. Unfortunately, our tests did not obtain a clear convergence rate. Examples involving topological changes can be naturally handled by the method. For example, Figure 5.11 displays the motion of three regions using a step size of r = 0.00025 and a finest cell width of h = From these images, it is clear that the large elliptical regions grow at the expense of the smaller circle. Before long, the two elliptical regions merge and the circle disappears. As we expect, this final shape slowly smoothes to become more circular. C h a p t e r 6 C o n c l u s i o n s 6.1 S u m m a r y In this thesis, we have considered a recent method [47, 48] for the motion by mean curvature of curves and surfaces. This method (MBO) naturally handles complicated topological changes with junctions. Because it produces approximations to these prob-lems more efficiently than other methods, it is often the only practical choice for the accurate treatment of three dimensional surfaces involving junctions. The usual finite difference discretization of the MBO-method has a number of limita-tions. For example, a very fine mesh must be used (see restriction (2.7)) to prevent the front from becoming stationary. Satisfying such a restriction globally is often computa-tionally impractical when accurate results are needed. Furthermore, each step produces an error in the position of the front which is comparable in size to the mesh spacing. To preserve the overall accuracy of the method for smooth curves, this leads to a work estimate of O ) operations per step for an optimized finite difference approach. (Here, T <§C 1 represents the time-length of each step of the MBO-method.) The usual finite dif-ference approach also produces irregular errors, which precludes the use of extrapolation in T for obtaining a higher order accuracy. To overcome these, and other, limitations we have proposed a new, spectral discretiza-tion of the MBO-method which we have found is often more than 1000 times faster than the usual finite difference approach. This method uses a Fourier cosine tensor product 97 Chapter 6. Conclusions 98 to discretize the heat equation which arises at each r-step of the MBO-method. The corresponding Fourier coefficients are determined using a quadrature with a piecewise linear approximation to the surface. This approach leads to several advantages over the usual finite difference discretization especially when used in combination with the recent, unequally spaced fast Fourier transform method given in [5]. Firstly, local refinement is much easier to carry out within a quadrature, rather than within a discretization of a differential equation (cf. [49]). Indeed, our proposed method recursively refines near the front to essentially eliminate the restriction (2.7), which plagues finite difference discretizations of the method. Furthermore, the piecewise linear approximation of the front which we use is more accurate than the crude approximation which arises for the finite difference approach. Indeed, we show that our new, spectral method requires only (9( log 2(r)) operations per step to preserve the overall accuracy of the MBO-method for smooth curves. This compares very favourably to the O operation count that arises for an idealized finite difference approach. In practice, how-Our proposed method also has the advantage that it allows for higher order extrapola-tions in r since it essentially eliminates irregular spatial errors. By essentially eliminating these spatial errors we also produce an algorithm which is nearly Euclidean invariant and hence more attractive for certain image enhancement applications (e.g., [42]). Further improvements in efficiency are obtained for the new method by neglecting high order Fourier modes which correspond to rapidly decaying solution transients. Gains also arise by carrying out the time integration exactly, rather than by a time-stepping method. New analytic and experimental results are also given to explain important properties of the MBO-method, such as the approximation error. In particular, an asymptotic expansion for the position of the front is derived to show that the method is first order Chapter 6. Conclusions 99 in r for smooth curves. For nonsmooth corners and singularities, numerical experiments are carried out to demonstrate that an order T log( r ) error arises. Higher order results for two-phase problems are obtained using extrapolation in the Fourier coefficients An asymptotic expansion for the angles of a two dimensional triple junction is also derived. This expansion indicates that each step of the MBO-method produces an 0(y/r) error in the junction angles which rapidly dissipates during subsequent steps. This result is of interest because it explains the stability of junction angles and suggests that an O(yfr) error arises for problems with junctions. Numerical experiments confirm that an 0(\/T) error occurs, and indicate that extrapolation can be used in conjunction with the proposed method to produce a marked improvement in the error. Finally, the improved utility of the new method is demonstrated by approximating and visualizing the motion by mean curvature of a number of three dimensional surfaces. Specifically, we give examples of two-phase, barbell-shaped regions which can undergo topological breaking depending on the width of the initial stem. Multiple phase motions are also approximated for examples involving three and four-phase junctions. We con-clude our results with a simple extension of the MBO-method to a nonlocal curvature model that preserves phase volumes. 6.2 Future Research Direct ions There are many directions for future work in diffusion-generated motion by mean curva-ture. A detailed theoretical investigation of the method would be desirable. In particular, a convergence proof for the case of junctions would be of great interest. Unfortunately, the mathematical tools which have been applied to produce rigorous convergence proofs for two-phase problems with topological changes cannot be readily extended to junctions Chapter 6. Conclusions 100 (see, e.g., [25]). However, because no method has been proven to converge for junctions we feel that even a formal convergence proof away from singularities would be an important achievement. Application of the MBO-method to different curvature-dependent motions would also be of interest. Currently, arbitrary junction angles can be treated by the method [45]. Affine velocity motions for two-phase problems can also be carried out [45]. Extending the method to affine velocity problems with nonsymmetric junctions (cf. [74]) would be useful for modelling a combination of surface and bulk effects in idealized grain growth applications [71]. Extensions to anisotropic mean curvature flows [71], complicated do-main geometries [20] and mixed boundary conditions [20] would also be of interest. Experimental studies of mean curvature motion for surfaces with junctions might also be carried out using our realization method for the MBO-method. For example, self-similar solutions (cf. [21]) and singularities (cf. [22]) arising from mean curvature flow might be studied. Further studies of nonsmooth corners might also be undertaken by comparing the results of the MBO-method with known similarity solutions [65]. Certainly, effective visualizations of two-phase surfaces have been carried out [21]. Methods for the improved visualization of junctions should also be considered. We have found the effective display of these surfaces to be rather challenging since features of interest (e.g., moving junctions) can become obscured by nearby regions. Finally, novel applications of junctions could also be investigated. In particular, the motion by mean curvature of junctions may be of value in certain image enhancement applications (cf. [19, 52, 23]). Bibl iography [1] D. Adalsteinsson and J .A. Sethian. A fast level set method for propagating inter-faces. Journal of Computational Physics, 118(2):269-277, 1995. [2] L. Alvarez, P. Lions, and J . Morel. Image selective smoothing and edge detection by nonlinear diffusion: II. SIAM J. Numerical Analysis, 29(3):845-866, 1995. [3] U . M . Ascher, S.J. Ruuth, and B.T.R. Wetton. Implicit-explicit methods for time-dependent PDE's . SIAM J. Numerical Analysis, 32(3):797-823, 1995. [4] G. Barles and C. Georgelin. A simple proof of convergence for an approximation scheme for computing motions by mean curvature. SIAM Journal of Numerical Analysis, 32(2):484-500, 1995. [5] G. Beylkin. On the fast Fourier transform of functions with singularities. Applied and Computational Harmonic Analysis, 2:363-381, 1995. [6] J . Bloomenthal. Polygonization of implicit surfaces. Computer Aided Geometric Design, 5(4):341-355, November 1988. [7] C. De Boor. A Practical Guide to Splines. Springer-Verlag, 1978. [8] J.P. Boyd. Chebyshev & Fourier Spectral Methods. Springer-Verlag, 1989. [9] K . A . Brakke. The surface evolver. Experimental Mathematics, 1(2):141—165, 1992. [10] A . Brandt. Multilevel computations of integral transforms and particle interactions with oscillatory kernels. Computer Physics Communications, 65:24-38, 1991. [11] R. Brent. Algorithms for Minimization Without Derivatives. Prentice-Hall, 1973. [12] L. Bronsard and Kohn. Motion by mean curvature as the singular limit of Ginzburg-Landau dynamics. Journal of Differential Equations, 90(2):211-237, 1991. [13] L. Bronsard and F. Reitich. On three-phase boundary motion and the singular limit of a vector-valued Ginzburg-Landau equation. Arch. Rat. Mech., 124:355-379, 1993. [14] L. Bronsard and B.T.R. Wetton. A numerical method for tracking curve networks moving with curvature motion. Journal of Computational Physics, 120(l):66-87, 1993. 101 Bibliography 102 R.L . Burden and J.D. Faires. Numerical Analysis. PWS-Kent Pub. Co., 1993. G. Caginalp. The dynamics of a conserved phase field system: Stefan-like, Hele-Shaw, and Cahn-Hilliard models as asymptotic limits. IMA J. applied Math., 44(l):77-94, 1990. C. Canuto, M . Y . Hussaini, A . Quarteroni, and T .A . Zang. Spectral Methods in Fluid Dynamics. Springer-Verlag, 1987. B . W . Char, K . O . Geddes, G.H. Gonnet, B .L . Leong, M . B . Monagan, and S.M. Watt. Maple V Language Reference Manual. Springer-Verlag, 1991. Y . Chen. Segmentation and association among lines and junctions for a line image. Pattern Recognition, 27(9):1135-1157, 1994. D. L. Chopp. Computing minimal surfaces via level set curvature flow. Journal of Computational Physics, 106(1):77—91, 1993. D.L. Chopp. Computation of self-similar solutions for mean curvature flow. Exper-imental Mathematics, 3(1):1—15, 1994. D.L. Chopp and J .A. Sethian. Flow under mean curvature: singularity formation, minimal surfaces and geodesies. Experimental Mathematics, 2(4):235-255, 1992. M . Demi. Contour tracking by enhancing corners and junctions. Computer Vision and Image Understanding, 63(1):118—134, 1996. A. Dutt and V . Rokhlin. Fast Fourier transform for nonequispaced data. SIAM Journal on Scientific and Statistical Computing, 14(6):1368—1393, 1993. L .C . Evans. Convergence of an algorithm for mean curvature motion. Indiana University Mathematics Journal, 42:553-557, 1993. L .C . Evans, H . M . Soner, and P.E. Souganidis. Phase transitions and generalized motion by mean curvature. Communications on Pure and Applied Mathematics, 45(9):1097-1123, 1992. L .C . Evans and J. Spruck. Motion of level sets by mean curvature, I. J. Differential Geometry, 23:69-96, 1986. R. E. Ewing. Adaptive mesh refinements in large-scale fluid flow simulation. In International Conference on Accuracy Estimates and Adaptive Refinements in Finite Element Computations, pages 299-314, Lisbon, Portugal, 1984. V . E . Fradkov, M . E . Glicksman, M . Palmera, J . Nordberg, and K . Rajan. Topological rearrangements during 2D normal grain growth. Physica D, 66:50-60, 1993. Bibliography 103 H . Frost, C. Thompson, C. Howe, and J . Whang. A two-dimensional computer sim-ulation of capillarity-driven grain growth: preliminary results. Scripta Metallurgica, 22:65-70, 1988. M . Gage and R.S. Hamilton. The heat equation shrinking convex plane curves. J. Differential Geometry, 23:69-96, 1986. C. W. Gear. Numerical initial-value problems in ordinary differential equations. Prentice-Hall, 1971. J . Glazier, M . Anderson, and G. Grest. Coarsening in the two-dimensional soap froth and the large-Q Potts model: a detailed comparison. Philosophical Magazine B, 62:615-645, 1990. W . B . Gragg. On extrapolation algorithms for ordinary initial-value problems. SIAM Journal on Numerical Analysis, 2:384-403, 1965. M . A . Grayson. A short note on the evolution of a surfaces by its mean curvature. Duke Mathematical Journal, 58(3):555-558, 1989. M . A . Grayson. Shortening embedded curves. Annals of Mathematics, 129:71-111, 1989. E. Hairer, S.P. Norsett, and G. Wanner. Solving Ordinary Differential Equations I. Springer-Verlag, 1987. D. Harker and E.R. Parker. Grain shape and grain growth. In Transactions of the A. S.M, pages 156-161, Cleveland, 1945. G. Huisken. Flow by mean curvature of convex surfaces in spheres. J. Differential Geometry, 20:237-266, 1984. T. Ilmanen. Convergence of the Allen-Cahn equation to Brakke's motion by mean curvature. Journal of Differential Geometry, 38(2):417—461, 1993. M . Kass, A . Witkin, and D. Terzopoulos. Snakes: Active contour models. Interna-tional Journal of Computer Vision, 1(4):321—331, 1987. B. B. Kimia, A .R . Tannenbaum, and S.W. Zucker. Shapes, shocks, and deforma-tions I: The components of two-dimensional shape and the reaction-diffusion space. International Journal of Computer Vision, 15(3):189—224, 1995. J.J. Koenderink and A . J . van Doom. Dynamic shape. Biological Cybernetics, 53:383-396, 1986. Bibliography 104 [44] W.E . Lorensen and H.E. Cline. Marching cubes: A high resolution 3D surface construction algorithm. In M . C . Stone, editor, Computer Graphics (SIGGRAPH '87 Proceedings), volume 21(4), pages 163-169, Anaheim, California, July 1987. [45] P. Mascarenhas. Diffusion generated by mean curvature. C A M Report 92-33, Uni-versity of California, Los Angeles, 1992. [46] The MathWorks. MATLAB. High-Performance Numeric Computation and Visual-ization Software. Springer-Verlag, 1992. [47] B . Merriman, J . Bence, and S. Osher. Diffusion generated motion by mean curva-ture. In J .E. Taylor, editor, Computational Crystal Growers Workshop, pages 73-83. American Mathematical Society, Providence, Rhode Island, 1992. [48] B . Merriman, J . Bence, and S. Osher. Motion of multiple junctions: a level set approach. Journal of Computational Physics, 112(2):334—363, 1994. [49] B . Milne. Adaptive Level Set Methods Interfaces. PhD thesis, University of Califor-nia, Berkeley, 1995. [50] F. Mokhtarian and A . K . Mackworth. A theory of multiscale, curvature-based shape representation for planar curves. IEEE Transactions on Pattern Analysis and Ma-chine Intelligence, 14(8):789-805, 1992. [51] W.W. Mullins. Two-dimensional motion of idealized grain boundaries. Journal of Applied Physics, 27(8):900-904, 1956. [52] J .A. Noble. Finding half boundaries and junctions in images. Images and Vision Computing, 10(4):219-232, 1992. [53] R . H . Nochetto, M . Paolini, and C. Verdi. A dynamic mesh algorithm for curvature dependent evolving interfaces. Journal of Computational Physics, 123(2):296—310, 1996. [54] S. Osher and J .A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79:12-49, 1988. [55] W . H . Press and G.B. Rybicki. Fast algorithm for spectral analysis of unevenly sampled data. Astrophysical Journal, 338:227-280, 1989. [56] F. Reitich and H . M . Soner. Three-phase boundary motions under constant veloci-ties. I: The vanishing surface tension limit. Technical Report 94-NA-015, Centre for Nonlinear Analysis, Carnegie Mellon University, 1994. Bibliography 105 J. Rubinstein and P. Sternberg. Nonlocal reaction-diffusion equations and nucle-ation. IMA J. Appl. Math., 48:248-264, 1992. J . Rubinstein, P. Sternberg, and J.B. Keller. Fast reaction, slow diffusion and curve shortening. SIAM J. Appl. Math., 49(1):116-133, 1989. J . Rubinstein, P. Sternberg, and J .B. Keller. Reaction-diffusion processes and evo-lution to harmonic maps. SIAM J. Appl. Math., 49(6):1722-1733, 1989. S.J. Ruuth. An algorithm for generating motion by mean curvature. In Proc. 12th International Conference on Analysis and Optimization of Systems Images, Wavelets and PDE's, Paris, France, 1996. H . Samet. Applications of spatial data structures: computer graphics, image pro-cessing, and GIS. Addison-Wesley, 1990. H . Samet. The design and analysis of spatial data structures. Addison-Wesley, 1990. G. Sapiro and A. Tannenbaum. Affine invariant scale-space. International Journal of Computer Vision, ll(l):25-44, 1993. G. Sapiro and A. Tannenbaum. Area and length preserving geometric invariant scale-space. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17(l):67-72, 1995. D. W. Schwendeman. A front dynamics approach to flow by mean curvature. Techni-cal report, Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, New York, 12180-3590, 1996. J .A. Sethian. Curvature flow and entropy conditions applied to grid generation. Journal of Computational Physics, 115:440-454, 1994. J .A. Sethian. Theory, algorithms, and applications of level set methods for propa-gating interfaces. Acta Numerica, 5:309-395, 1996. E. Sorets. Fast Fourier transforms of piecewise constant functions. Journal of Computational Physics, 116:369-379, 1995. J.C. Strikwerda. Finite Difference Schemes and Partial Differential Equations. Wadsworth L Brooks/Cole, 1989. T.D. Sullivan. A technique for convolving unequally spaced samples using fast Fourier transforms. Technical report, Sandia Report, 1990. J .E. Taylor, J .W. Cahn, and C A . Handwerker. I-Geometric models of crystal growth. Acta Metall. Meter., 40(7): 1443-1474, 1992. Bibliography 106 [72] D. Terzopoulos, A . Witkin, and M . Kass. Constraints on deformable models: Re-covering 3D shape and nonrigid motion. Artificial Intelligence, 36(1):91—123, 1988. [73] J.J. Tyson and J.P. Keener. Singular perturbation theory of traveling waves in excitable media. Physica D, 32:327-361, 1988. [74] H . Zhao, T. Chan, B. Merriman, and S. Osher. A variational level set approach to multiphase motion. C A M Report 95-36, University of California, Los Angeles, 1995. A p p e n d i x A A n Es t imate for the N u m b e r of Basis Functions We now estimate how many basis functions are required to accurately represent the position of a smooth, two dimensional closed curve (e.g., Figure 1.1 or Figure 2.1) after a time, r. Specifically, we compare the position of the level set | for n-l U?(x,y) = Y °ij exp ( -7r 2 [ i 2 + j 2 ] r ) cos(7u:r) COS(TTjy) i,j=o and oo U™(x,y) = Y Cij exp(-n2\i2 + j2]r) cos(Trix) cos(njy) to find a bound, nc, such that for all n > nc the approximate interface locations for U™ and differ by at most e + C(er). There are three steps in the derivation: 1. Lemma A . l demonstrates that if 1 7 7 ° ° _ jp\ < ! Z\/TTT then the approximate interface locations for U™ and differ by at most e + C(er). 2. Lemma A.2 derives an estimate of the size of the Fourier coefficients, c,-j. 3. Finally, Theorem A . l estimates the sum of the neglected Fourier terms of [/" and shows that a safe choice of n is given by ln(: n > nc = \r Vl96"L; 107 Appendix A. An Estimate for the Number of Basis Functions 108 where L is the length of the curve at the beginning of the current step. L e m m a A . l If the initial interface is smooth (four times continuously differentiable) and \ur- ^ n i o o < 2 ^ ' then the \ level sets for and [/" are at most a distance e + C(er) apart. Proof. Begin by orienting the shape according to Figure 4.1 for some arbitrary point on the initial interface. Following the notation of Section 4.1.1, we see from Equa-tion (4.28) that after a time r the level set ( | + 2 J ^ r ) is a distance (2), Z(T) = g [ 2 ) T - e + i„(4) _ £r„(2h3 (gW)3 r2 + 0(er) + 0(T3) _T 3 above the initial interface position 1. Comparing to the result for the level set | , g{2)r + Z(T) T2 + 0(T3), and noting that all terms not involving a factor of e are identical, we see that \Z(T)-Z(T)\ < e + 0(er). • L e m m a A . 2 If i ~> n, then the Fourier coefficients satisfy I Cij I ^ n where L is the length of the curve at the beginning of the current step. 1 W e have assumed e < T. T h i s w i l l be true for any reasonable choice of e since we expect to take steps i n to ta l . Appendix A. An Estimate for the Number of Basis Functions 109 Proof. From Equation (3.12), we know that h,-1 < 4 Rt < 4 / / c o s ( ^ ) dA Rt Over any rectangle, R C Rt, of width -{ ^ cosiixix) dA R (A.37) since J ^ ° + T cos(irix)dx = 0. Dividing the domain into rectangles of width | , we see that only subregions in contact with the interface contribute to the bound (A.37). See, for example, Figure A . l . Indeed, letting the height of each rectangle tend to zero, it is easy to see that only contributions within a horizontal distance | from the interface contribute. Using this fact, and integrating only over half of each rectangle so the sign of cos(7rza;) does not change, I Cij I <-AL AL i ~ n where L is the length of the interface at the beginning of the current step. • We now give the main result of this section. Theorem A . l Suppose that U™ approximates and that the initial interface is smooth. If n > nc = 9 ' TT T then the interface approximations for U™ and are within a distance e + 0(er) of one another. Appendix A. An Estimate for the Number of Basis Functions 110 Proof. We begin by estimating oo \U™ - U?\oo = | Yl Cij exp(-n2\i2 + j2]r) cos(irix) cos^jy)]^, i,j = o max(i,j) > n oo < J2 |c,-i|exp(-7r2[i2+j2]r), i,j = o max(z, j) > n oo oo Using Lemma A . 2 , this simplifies to give us O T OO OO Wr-UlU < - £ £ e x p ( _ V [ z 2 + j 2 ] r ) , " i =n j=0 = — ^ e x p ( - ^ 2 i 2 r ) ^ e x p ( - 7 r 2 j 2 r ) , n i=n j=0 < — (exp(—7r 2n 2r) + J exp(—n2nsr) dsj ( 1 + j exp(—n2s2r) ds = — {1 + (1 + --4=) exp(-^Vr) Now, assuming that irn2T > 1 and r < 1 (this is consistent with the final result for any reasonable error tolerance), it is easy to show that |C /r- t/ T n |oo < ^ e x p ( - A V ) . Appendix A. An Estimate for the Number of Basis Functions 111 By Lemma A . l , we want | ^ T ° ° - ^ l o o < 6 2 ^ ' so we seek an n such that r- exp(—7r n T) < N , Solving for ra, we arrive at a bound, w > „ = ,/iMifa)i; •
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Efficient algorithms for diffusion-generated motion...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Efficient algorithms for diffusion-generated motion by mean curvature Ruuth, Steven J. 1996
pdf
Page Metadata
Item Metadata
Title | Efficient algorithms for diffusion-generated motion by mean curvature |
Creator |
Ruuth, Steven J. |
Date Issued | 1996 |
Description | This thesis considers the problem of simulating the motion of evolving surfaces with a normal velocity equal to mean curvature plus a constant. Such motions arise in a variety of applications. A general method for this purpose was proposed by Merriman, Bence and Osher, and consists of alternately diffusing and sharpening the front in a certain manner. This method (referred to as the MBO-method) naturally handles complicated topological changes with junctions in several dimensions. However, the usual finite dif-ference discretization of the method is often exceedingly slow when accurate results are sought, especially in three spatial dimensions. We propose a new, spectral discretization of the MBO-method which obtains greatly improved efficiency over the usual finite difference approach. These efficiency gains are obtained, in part, through the use of a quadrature-based refinement technique, by in-tegrating Fourier modes exactly, and by neglecting the contribution of rapidly decaying solution transients. The resulting method provides a practical tool, not available hitherto, for accurately treating the motion by mean curvature of complicated surfaces with junc-tions. Indeed, we present numerical studies which demonstrate that the new algorithm is often more than 1000 times faster than the usual finite difference discretization. New analytic and experimental results are also developed to explain important prop-erties of the MBO-method such as the order of the approximation error. Extrapolated algorithms, not possible when using the usual finite difference discretization, are proposed and demonstrated to achieve more accurate results. We apply our new, spectral method to simulate the motion of a number of three dimensional surfaces with junctions, and we visualize the results. We also propose and study a simple extension of our method to a nonlocal curvature model which is impractical to treat using the previously available finite difference discretization. |
Extent | 8198422 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079751 |
URI | http://hdl.handle.net/2429/6142 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1996-148254.pdf [ 7.82MB ]
- Metadata
- JSON: 831-1.0079751.json
- JSON-LD: 831-1.0079751-ld.json
- RDF/XML (Pretty): 831-1.0079751-rdf.xml
- RDF/JSON: 831-1.0079751-rdf.json
- Turtle: 831-1.0079751-turtle.txt
- N-Triples: 831-1.0079751-rdf-ntriples.txt
- Original Record: 831-1.0079751-source.json
- Full Text
- 831-1.0079751-fulltext.txt
- Citation
- 831-1.0079751.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079751/manifest