M O N O T O N E S C H E M E S F O R D E G E N E R A T E A N D N O N - D E G E N E R A T E P A R A B O L I C P R O B L E M S by MIRNA LIMIC M.B.A. Rochester Institute of Technology, 2001 B.Sc, Economics and Business, 1999 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS F O R T H E D E G R E E O F M a s t e r of Science in T H E F A C U L T Y OF G R A D U A T E STUDIES (Computer Science) T H E U N I V E R S I T Y OF BRITISH C O L U M B I A December 2006 © Mirna Limic. 2006 Abstract This thesis analyses several monotone schemes for degenerate second order parabolic PDEs over two-dimensional domains which is in most cases equal (0, l ) 2 . For such problems the degeneracy means that in addition to the axis spanned by the standard basis vectors, there exists another pair of orthogonal cordinate axes such that the spatial difference operator is of the second order in one of the directions, and of the first order in the remaining directions. The direction of one axis along which the spatial difference operator is of the second order we call the direction of diffusion. The thesis considers only constant coefficient PDE's and therefore the degeneracy can be easily determined by solving the eigenvalue problem for diffusion tensor. We analyse the impact of the direction of diffusion on the convergence of a scheme. Previous work on a second order elliptic problem has shown that central differences taken in the direction of diffusion produce a convergent scheme, wherease disalignment between the two result in non-convergent schemes. One of our aims was to check the validity of this finding and possibly improve the construction of schemes. As a result we present a novel approach to building monotone schemes from the diffusion tensor by taking spatial step sizes of different length in order to align the second order central differences with the direction of diffusion. We give a step by step algorithm and supply all of the findings. Our findings are based on M A T L A B m file, and C language implementations. ii Contents Abstract ii Contents iii List of Tables . . . v List of Figures ' vi Acknowledgments viii Dedication ix 1 Introduction 1 1.1 Introduction 1 1.2 Motivation 2 1.3 Outline 4 2 Theoretical background 6 2.1 Degenerate equations 6 2.2 Central difference approximations 9 2.2.1 Central differences denned 10 2.2.2 Differences of the Toolbox of Level Set Methods 1.1 11 2.2.3 Checking derivative approximations 12 2.2.4 Difference approximations and matrices of positive type 12 2.3 Courant-Friedrichs-Levy (CFL) condition for monotone schemes 14 2.4 Consistency of monotone schemes 16 3 Analysing Problems 2. and 3 17 3.1 Degenerate heat equation 17 3.1.1 Goal and problem setting 17 3.1.2 Deriving analytical solutions 18 3.1.3 Two-dimensional eigenvalue problems 21 3.1.4 Solutions on squares 22 3.1.5 Solution with the vortex velocity field 22 3.1.6 Findings 24 3.2 General degenerate parabolic equation 29 3.2.1 Goal and problem setting 29 iii 3.2.2 Monotonicity of the schemes 31 3.2.3 Rotation of coordinate axes for General degenerate parabolic equation 32 3.2.4 Findings 33 4 Building monotone schemes • . . . 46 4.1 Central difference approximations 46 4.2 Computational procedure 48 4.3 Problem definition 49 4.4 Findings ' 50 4.5 Analysing Problem 4 59 4.5.1 Problem Definition 59 4.5.2 Finding's 59 5 Software implementation 63 5.1 Implementation details of m files 63 5.2 m file v.s. C language implementations .' 65 5.3 Implementation of C files 66 6 Future work 67 Bibliography 69 Appendix 71 iv Lis t of Tables Table 4.1 Example of the IVP sin(67r(x - (y - v2t))) for final time 1.0, and v = (0,0.5) '. 60 Table A . l A list of m ( M A T L A B ) files implemented in the Toolbox of Level Set Methods 1.1 73 v List of Figures Figure 3.1 Analytical solution of Problem A computed for the final time 0.5 , boundary condition 0, velocity v = (0,1), grid size [100 x 100], and initial condition sm(2mj) s'm(nx)/\/2 Va- G domain D, and Vy G [0,1/2] 25 Figure 3.2 Numerical solution of Problem A computed for the final time 0.5 boundary condition 0, velocity v = (0,1), grid size [100 x 100], and initial condition s'm{2wy) sin(7ra)/\/2 V.T G domain D, and Vy G [0,1/2] 26 Figure 3.3 Analytical solution of Problem B computed for the final time 0.5 , boundary condition 0, vortex velocity field, grid size [100 x 100], and initial condition sin(7ra) sm{mj) Vz,y G domain D 27 Figure 3.4 Numerical solution of Problem B computed for the final time 0.5 boundary condition 0, vortex velocity field, grid size [100 x 100], and initial condition sin(7r.r) sin(7ra) Vx,y G domain D 28 Figure 3.5 Numerical solution for the final time 0.2 for the discretization (3.11), boundary condition sin(67r(.x — y)), grid size [100 x 100], and initial condition i) V.-/:.;//. . . . '. 38 Figure 3.6 Analytical and numerical solutions for the final time 1.0 for the dis-cretization (3.11), boundary condition sin(67r(.x — y)), grid size [100 x 100], and initial condition zero 39 Figure 3.7 Final time 1.2 for the discretization (3.12), boundary condition sin(67r(.x-— y)}, grid size [100 x 100], and initial condition 0 Vx,y 40 Figure 3.8 Final time 1.2 for the discretization (3.13), boundary condition sin(67r(a;— y)), grid size [100 x 100], and initial condition zero for all x,y 41 Figure 3.9 Final time 0.2 for the discretization (3.11), periodic boundary condi-tions, grid size [100 x 100], and initial condition $in(6n(x — y)).. 42 Figure 3.10 Final time 0.2 for the discretization (3.13), periodic boundary condi-tions, grid size [100 x 100], and initial condition sin(67r(a; — y)) 43 Figure 3.11 Numerical solution for the final time 1.0 for the discretization (3.12), boundary condition sin(67r(a; + '</)), grid size [100 x 100], and initial condition 0 Wx,y G D. . • 44 Figure 3.12 Numerical and analytical solutions for the final time 1.2 for the dis-cretization (3.11), boundary condition sin(67r(a.' + y)), grid size [100 x 100], and initial condition 0 Vx,y G D 45 Figure 4.1 The numerical neighbourhoods used in the calculation of the terms of the system matrix A in (4.2) 48 vi Figure 4.2 Numerical solution for the final time 0.1 using unequal grid step sizes discretization, boundary condition sin(67r(a; — 2y)), grid size [100 x 201], and initial condition 0 Vx,y G D 52 Figure 4.3 Analytical solution for the final time 0.1 using unequal grid step sizes discretization, and grid size [100 x 201] 53 Figure 4.4 Numerical solution for the final time 0.25 using equal grid step sizes discretization, boundary condition sin(67r(a; — 2y)), grid size [100 x 100], and initial condition 0 Vrr.y 6 £) 54 Figure 4.5 Numerical and analytical solutions for the final time 0.25 using unequal grid step sizes discretization, boundary condition sin(67r(3x — 2y)), grid size [100 x 68], and initial condition 0 Vx,y € D 55 Figure 4.6 Numerical solution for the final time 0.25 using equally spaced grid step sizes, boundary condition sin(67r(3a; — 2j/)), grid size [100 x 100], and initial condition 0 Vx, y S D 56 Figure 4.7 Numerical solution for the final time 0.05, with boundary condition —x (L~x) on the top and bottom ends of the domain, y——— on the left and an a22 right ends of the domain, and grid size [100 x 100] 57 Figure 4.8 Analytical solution for the final time 0.05, with boundary condition —x^~x^ on the top and bottom ends of the domain, y^—^- on the left and «11 1 J "22 right ends of the domain, and grid size [100 x 100]. . . , 58 Figure 4.9 Numerical solution for the final time 1.0, with boundary condition sin(67r(.x — (y — V2t))), where v-x = 0.5, grid size [100 x 100], and initial condition sin(67r(a: — y)) 61 Figure 4.10 Numerical solution for the final time 0.5, with boundary condition sin(67r(a; — (y — i^*))), wherein = 1, grid size [100 x 100], and initial condition sin(67r(:i; — y)) 62 vii Acknowledgments I would like to thank my supervisior for his guidance. In addition, I would like to express many thanks to Professor Anthony Peirce of the University of British Columbia Mathematics Department for the suggestions he gave as the second reader of my thesis. Finally, I thank everyone else who helped me with their knowledge and support in this thesis research. viii my family whose support enabled me to make this thesis. This is our success. ix Chapter 1 Introduction 1.1 Introduction We focus our study on monotone schemes for a class of degenerate linear parabolic partial dif-ferential equations. For the sake of comparison, monotone schemes for some non-degenerate problems are also analysed, such as the heat equation. Our goals are to (1) investigate the efficiency of monotone schemes in generating a "good1' numerical ap-proximation, and (2) confirm or disprove through numerical calculations and mathematical tools the claim that monotonicity has impact on the convergence of a scheme. We consider two-dimensional problems on bounded domains. We limit ourselves to a square two-dimensional domain, which is in most cases (0, l ) 2 , and sometimes (—1, l ) 2 . For such domains the degeneracy means that there exists another pair of orthogonal coordinate axes such that the spatial differential operator is of the second order in one of the directions, and of the first order in the remaining directions. Therefore, for degenerate equations on the defined domains we use terminology in the direction of diffusion, meaning, in the direction of one axis along which the diffusion operator is of the second order. The monotone schemes are interesting because of the following simple fact. If stable and consistent, they are convergent, as proved by Barles and Souganidis [5]. Our findings confirm the work of Oberman [2] which states that central difference approximations to second order terms taken in the direction of diffusion result in a convergent monotone scheme for the corresponding elliptic problem'. To that result our next goal was to create a method for building monotone schemes for degenerate and non-degenerate parabolic equations where the central difference approximations are taken in the direction of diffusion. The idea rests on taking discretization steps of different size along axes in order to "align" the derivative approximations with the diffusion direction. Spatial discretization of linear PDEs gives a system of linear ODEs ^u(t) + Au(t) = 0 , 1 where the matrix A is called the system matrix. For monotone schemes the system matrix has a particular structure. It has positive values along the diagonal and negative or zero values on the off-diagonal places. If we decide to solve the above system of ODEs using the explicit Euler method this special property of A allows us to compute the Courant-Friedrichs-Levy bound on the time step, r. In some cases (like the case of the two-dimensional heat equation) this bound can be found analytically, while in general it must be computed numerically, ff we decide to solve the system of ODEs using the implicit Euler, or some other implicit method, the described structure of A enables very easy implementation. It is important to note that there is a 1-to-l relationship between monotone schemes and matrices with above specified structure. Therefore, in order to build a monotone scheme for a linear P D E , one is in effect creating a system matrix of ODE's which has positive values along the diagonal and non-positive or zero values on the off-diagonal positions. A class of linear PDEs to which we limit our research are given in Section 1.2. The schemes for these equations are built from degenerate schemes of the individual terms in the respective equations. In [2] the author studied transformations preserving monotonicity. For instance, the sum of monotone functions is a monotone function. Therefore, we could also say that we built our monotone schemes by summing the monotone schemes of each individual term in a respective equation. For example, we sum a monotone scheme for a convective part of an equation with a monotone scheme for the difussive part of that equation thus obtaining a monotone scheme for the equation. The monotone sliemes we use are explained in Section 2.2.1. So far we have presented some benefits of using a degenerate, i.e. monotone differ-ence scheme. However clue to the monotonicity requirement the accuracy of a monotone difference scheme is at most of the first order [2]. This fact was first explained by Barles and Souganidis [5]. .Higher than first order schemes such as E N O , and W E N O as given in Osher, Fedkiw [15] result in matrices that do not have the special property forementioned, and are therefore not monotone. Since we study monotone schemes we restrict ourselves to first order integration and use Forward Euler in deriving our solutions. The software implementations are present in two programming languages. One is the M A T L A B m language, and the other is C programming language. The files written in the in language are designed to work with a software package, the Toolbox of Level Set Methods, version 1.1, [1]. The Toolbox of Level Set Methods 1.1 is a set of routines implementing level set methods which solve time-dependent Hamilton-Jacobi partial differential equations. The C implementations are designed as stand-alone programs. The difference and reasons for two different implementations are given in Section 5. 1.2 Motivation The analysis in this thesis is motivated by the research presented in [2], where degener-ate elliptic and parabolic equations of the second order are examined for the purpose of constructing their numerical or approximate solutions. The emphasis is put on numerical methods which are called monotone schemes. For their definition refer to Section 2.1. The analysis of monotone schemes in this thesis focusses on three initial value problems on the 2 square domain, D = (0,1) x (0,1). We define the problems for a function / as follows Problem 1. Heat equation Problem 2. Degenerate heat equation with advection ft - °~t fxx ~ °l fyy = 0 /,. + vV/ ~ a2 ./;,,, = 0 Problem 3. ,i. General degenerate f\ — ^ a,; <), f — 0 parabolic equation » j ' = 1 Problem 4. ,i General degenerate parabolic /t + v V / — ^ a»j did-jf = 0 equation with advection *>'=! We implemented an example of Heat equation to observe behaviour of numerical solution in time when the scheme used is monotone. However, since Heat equation is discussed in many differential equation books, for example Boyce, DiPrima [9], we skip its theoretical analysis. To execute the implemented example see Appendix. We supply initial and boundary conditions to the above problems in later chapters, yet we do not defer the definition of the parameters present in the above equations. The positive numbers a2, a2, a2 are usually called diffusion constants. Continuous functions v have two components vx,vy and define velocity fields on D — (0, l ) 2 . Continuous functions afj define a matrix a — {a%iSn which is called diffusion tensor. At this point we are obliged to explain our matrix notation. A matrix a = {<H'-j}n is a matrix with indices going from (1,1) to (2,2). The matrix a must be symmetric at each point of D and must be either positive definite or positive semidefinite. In the former case it has two positive eigenvalues and Problem 3. is non-degenerate. In the latter case the matrix a has one zero eigenvalue and one positive eigenvalue and the corresponding Problem 3. is degenerate. Actually the degeneracy does not need to be realized at each point of D. It is sufficient that the degeneracy happens at a single point of D = [0,1]2. Problem 3. can be rewritten in another form, which is utilized in Game theory, Problem 3. ft - tn\ce[LDR] = 0 where L, R are 2 x 2 matrices and D is Hessian of / . The matrices L, and R. cannot be arbitrary since trace[LZ).R] must coincide with Ylij=i aijdi&j.f-In order to introduce monotone schemes we start our analysis with Heat equation. Next, Degenerate heat equation is considered as a degenerate version of Heat equation. Monotone schemes which are developed for Heat equation are applied to Degenerate heat equation. An extension of investigation onto General degenerate heat equation is necessary 3 in order to draw conclusion on the efficiency of monotone schemes based on the results of our research. In addition, in the research of General degenerate heat equation we confirmed that the direction in which the central differences tire taken has an impact on the convergence of a scheme. Chapter 4 presents a new method, created by M. Limic, for building monotone schemes using spatial step sizes of different length. This method builds monotone schemes from a diffusion tensor of the degenrate elliptic differential operator. We conclude our analysis of monotone schemes with General degenerate parabolic equation with advection taking into account the importance of the direction of central differences. In the analysis we use the method The best way to demonstrate the efficiency of monotone schemes would be to com-pare the solutions / of the considered problem and its numerical approximation fapp. Here we are presented with an obvious difficulty. Solutions are not known and this fact is the basic reason which forces us to look for numerical or approximate solutions. Apart from the-oretical results telling us about the convergence of numerical solutions towards solutions / , the efficiency can be studied by using examples for which the unique solutions / are known. We say that solution / is known in a closed form. In this report we analyze numerical solutions using examples for which solutions / are known in closed forms. 1 .3 Outline We have organized the material into eight chapters. Here we give a brief outline of the issues discussed in each chapter. (i) In chapter 1 we give the motivation and goals of this thesis. We also introduce the topic of monotone schemes, and note some important facts regarding the system matrices associated to monotone schemes and how schemes are built. (ii) Chapter 2 contains the theory of the degenerate elliptic and parabolic equations and corresponding numerical schemes. Following the theoretical background is a descrip-tion of the central differences we had to construct to derive monotone schemes for problems considered. Finally we give our derivation of the.CFL condition for mono-tone schemes and its comparison to the C F L condition implemented in the Toolbox of Level Set Methods 1.1. (iii) Chapter 3 is devoted to a detailed study of monotone schemes by analysing the problems of Section 1.2. We present the derivation of analytical solution to three examples of Degenerate heat equation, as well as solutions to the examples of General degenerate parabolic equation. We consider each problem in turn giving our findings as the concluding discussion. While Degenerate heat equation serves exclusively for the analysis of the efficiency of montone schemes, General degenerate parabolic equation is additionally used for an analysis of the importance of the direction in which the central differences are taken on the convergence of a scheme. (iv) Chapter 4 presents a method we developed for computing monotone schemes for the case of degenerate and non-degenerate diffusion tensors. This is achieved by taking numerical step sizes of different length along two coordinate axes. We present the 4 algorithm for deriving the step size in y-axis direction given the number of nodes in the .u-axis direction. Our findings are based on the implementation of several examples of General degenerate parabolic equation to this setting. (v) Chapter 5 is concerned with the implementation details of m files only. Since the C programming language implementations were not required at any time during this thesis research, their implementation will not be discussed here. They were done as an aid to help overcome the disadvantages of slowness and memory overconsumption of the m language implementations. If however, the reader is interested in obtaining the implementation details of C files as well at the C implementations themselves, they can contact me at mirnalim@cs.ubc.ca. (vi) Based on the work completed, in Chapter 6 we propose some topics that can be taken as further research. (vii) Appendix serves as a small manual for compiling and executing the programs. It contains a listing of the programs that implement the problems considered in the course of the thesis research. As an end note to this chapter it should be said that this thesis is written to be a learning tool just as much as a presentation of the published, work and our research conducted in the area of degenerate monotone equations and their schemes. 5 Chapter 2 Theoretical background The theory we look at and develop in this thesis rests mainly on the published work in the field of degenerate elliptic and parabolic equations. Methods of the eigefunction theory are used for the purpose of deriving analytical solutions to PDEs. 2.1 Degenerate equations In the work of Barles and Souganidis [5] it is shown that monotone, stable, and consistent schemes are also convergent. The only drawback in building such schemes is the fact that they are of the first order of convergence. If the coefficients of a P D E are smooth enough, schemes with higher order of convergence can be successfuly used. However, at points where the coefficients are not smooth one prefers to use monotone schemes due to the stability they have. Monotone schemes are not only used for linear P D E problems, but also for non-linear PDEs. In effect, they are usually defined for non-linear PDEs. Spatial discretization of an initial value problem (IVP) for a P D E gives a system of ODEs of the form . ju^t) = F>(t,u(t)), i = l ,2 , . . . , iV , (2.1) where u(t, x) are the solutions of the IVP for a PDE. These solutions are at grid nodes .r f approximated by functions ui(t). These functions define a column vector u(t) of approximate solutions. The functions F'1 result from spatial discretization and can depend explicitly on t, as well as on the approximate solutions u;. The constructed system (2.1) of ODEs must be solved using discretization in the time variable. This means that we consider discrete times to = 0,ti,... ,t,M = T and approximate solutions «;(£fc) by numbers u\k\ Thus the columns u(tfc) are approximated by colums u1-^. In the case of the Forward Euler method applied to (2.1) we obtain u<*+1> = u f> + r f ( t t , u ( t > ) := ff'(r,tfc,u(fc)). (2.2) Based on the definition of monotonicity as given in [14] we give the following defi-nition. 6 D E F I N I T I O N 2.1.1 (Monotone scheme) We say that a scheme is monotone if the functions Hl(r, tk, u'fc') do not decrease when their variables increase. It is important to note that a scheme is not monotone for all values of time increment T. This fact implies a necessity to analyse the monotonicity of a scheme with respect to the choice of the time step. The aim here is to find the biggest time step for which the scheme is monotone. This issue is discussed in Section 2.3. For example, take the parabolic equation ut — uxx = 0. After carrying out the central difference spatial discretization of uxx we get the O D E system d where the functions are defined by dtu(t) - F[u(t)] = 0, Fl{u) -U; — U. i-l dx ui ~ui+i Then, in the present case, the system of difference equations (2.2) has the form ~ V dx*) 1 + dx* \U'+l + Ui-1)' defining a scheme. Looking at the definition of monotonicity 2.1.1 we see that this scheme is monotone for all 0 < T < dx2/2. We computed this result by noting that we must have 2r „ dx2 in order for the considered scheme to be monotone. Next we give the definition of degenerate ellipticity for which we note that it includes the definition of monotonicity. For each node i the symbol N(i) denotes the neighbourhood of node i. D E F I N I T I O N 2.1.2 We say that a scheme Ui — Uj ,.i = l N is degenerate elliptic if each F% is a non-decreasing function in each of its arguments. Given the above definition of degenerate ellipticity it follows that every degenerate scheme is also monotone. Next we would like to state an important theorem taken from [2] and class notes of [6] for our discussion on stability and monotonicity. 7 T H E O R E M 2.1 A scheme is monotone and non-expansive in the l°° norm if and only if it is degenerate elliptic. T H E O R E M 2.2 The accuracy of a monotone finite difference scheme is at most first order. The three problems defined in Section 1.2 must be supplied with boundary and initial conditions. So defined we call them initial boundary value problems for linear PDEs. In the remaining.part of this section we describe a procedure for constructing mono-tone schemes. First we have to define conditions on the coefficients ay, U;, and c of a general linear elliptic differential operator of the second order d ' d A(x) = - __ a:jXxU),;i, + Y^ViWdi + c(x). (2.3) i,j = l i=l D E F I N I T I O N 2.1.3 The differential operator is called degenerate elliptic if the following three conditions are valid (a) The functions ay = aji; vt (i, j = 1, 2, . . ., d) and c satisfy the following conditions: \vi\i lcl fi M, c < 0, aij,v-i,c uniformly continuous on Rd. (2.4) (b) The differential operator (2.3) is elliptic on R d, i.e there exists a positive number M such that d M\z\\> Y _ a . y ( x ) ^ ^ > 0 , x e M d , (2.5) where Zi £ Cd are complex numbers and \z\2 the corresponding l2-norm. In other words, diffusion tensor is positive semidefinite. (c) There is a, point x £ Rd such that j=i aij{x)zi^j = 0. We say that the diffusion tensor is defined by functions a.y, convection is defined by velocities v-t and the function c is called the killing rate. For the problems considered in this thesis the killing rate is zero. Examples of this thesis have either homogeneous or non-homogeneous Dirichlet boundary conditions or periodic boundary conditions. Homogeneous Dirichlet boundary condition is specified by zero value of solution, non-homogeneous by non-zero value of solu-tion. For an IVP defined by dtu(t,x) - A(x)u(t,x) = 0, u(0,x) = u0(x), • (2.6) u(t,x)dD = 0, we use those spatial discretizations for which we get monotone schemes. In Section 2.2.1 some of the possibilities are described in detail. A discretization of the differential operator A(x) results in a system matrix A with entries A y , where indices i,j refer to grid nodes. Discretizations leading to monotone schemes are easily characterized using the notion of matrices of positive type [4], which we define next. 8 D E F I N I T I O N 2.1.4 A matrix A with entries Aij is said to be of positive type if it has positive'diagonal entries, non-positive off.-diagonal entries and positive or zero row sums, An > 0, Ai:j < 0 for i ± j, Aj = ]T A{j > 0. ' i s / A matrix of positive type with zero row sums, Aj = 0 Vj £ / , is called conservative. If we wish to construct a monotone scheme for the IVP (2.6) we need to follow certain principles RULES OF DISCRETIZATION Spatial discretization of the differential operator A(x) must be carried out using the following two rules in order to ensure the monot'onicity of the resulting scheme 1) The convection term Yli=ivi(x-)^iu's discret'ized using the upwinding scheme. 2) The diffusion term —J2ij=iaij(x)^i^ju must be discretized in a way to make the resulting system matrix of positive type. If we abide to these rules, we build a system of ODEs of the form jtu{t) + Au(t) = 0, where the matrix A is of positive type. This matrix, as mentioned before, is called the system matrix. A system matrix of positive type results in a monotone scheme, and a monotone scheme always has a system matrix of positive type. Therefore we could laicaly say that system matrices of positive type and monotone schemes are in a 1—to—1 relationship. Let us define a diagonal matrix D with entries Du = An and a matrix B with entries Bij = —Ay i ^ j. Then A. = D - B. If we use Forward Euler method, the constructed O D E is converted into the following system of difference equations u<fc+1> = (/ - TD)U^ + T B U O This scheme is monotone for sufficiently small r. The computation of T is given in Section 2.3 where we use the same decomposition of the matrix A onto D and B in order to compute the expression for calculating the upper bound on T, i.e. the C F L condition. 2 . 2 Central difference approximations We would like to introduce a way of defining central differences that is different from the one available in the Toolbox of Level Set Methods 1.1 distribution. Central differences as defined below are the ones used by Oberman [2], and differ from the ones in the Toolbox of Level Set Methods 1.1 in their calculation of the second order mixed partial derivative. 9 2.2.1 C e n t r a l differences defined The orthogonal coordinate system in Rd is determined by unit vectors e; of the canonical basis. Points x = h Yl'i=x h^u h is the grid step, h £ Z, define a numerical grid G/i on R''. The partial differential operators of the first and second orders can be approximated in the usual way cV(x) -» a.,/(x) = ^ [/(x + eih) - / (x - eth)}, (2.7) • 9?/(x) - dufix) = i [ / ( x + e l / l ) - 2 / ( x ) + / ( x - e ^ ) ] , (2.8) where we use the symbol 3 to indicate the approximation to a first or second order derivative, and symbol —> to show that the operator to the left of the arrow is approximated by the discrete operator to the right of the arrow. Following our rules of discretization we discretize the operators Vidi by using the upwinding procedure rather than by using the central difference operator 3;/(x). Thus the first partial derivative in the a>axis direction would be approximated by one of the following two possibilities ^ [/(x + eih) - /(x)] , i [/(x) - / (x - e</0] • which are called the forward and backward finite difference operators respectively. The second order partial differential operator df is approximated by two'subsequent applications of the first order operators. For example, if the first step of the approximation is the forward finite difference operator i [/(x + eth) - /(x)], then the application of the backward finite difference operator gives us I [ l ( / ( x + e I / l ) - / ( x ) ) - i ( / ( x ) - / (x - e i / 0 ) ] = £ [/(x + eth) - 2 /(x) + / (x - e,/t)] • We would have arrived at the same answer had we applied the directions of the approxima-tion in the reversed order, i.e. first backward, and then forward. The mixed partial differential operator can be approximated by one of the following four possibilities OiOifix) -» a. i j/(x) = J _ f f(x±eih±ejh)-f{x±eih)-f(x±ejh) + f(x), (2.9) h? \ -f(x ± e i h - ejh) + f(x ± eth) + / (x =p ejh) - /(x). The above operators are derived by subsequent applications of forward finite difference operator and backward finite difference operator in the directions ei and e 2 . We now show how one of the four above discretizations for a mixed partial differential operator can be obtained. We take first the backward finite difference operator in the first variable, i. I[ / (x)- / (x- e i /0)] . 10 and then apply the backward finite difference operator in j i [ / ( x ) - / ( x - e ^ ) ) ] . The resulting expression is i [ i ( / ( x ) - / (x - eji)) - i ( / ( x ejh) -r / ( x - eih - ejh))] • = [/(x) - / (x - e.;/),) - / (x - ejh) + / (x - eji - e j / i ) ] . This is easily recognisable as one of the four possibilities given in (2.9). With a similar-procedure one can obtain the remaining three possibilities given in (2.9). In the problems considered in this thesis we have used upwinding for the approxi-mation of the first order derivatives. When approximating second order non-mixed terms we used operator 3a. For mixed partial derivative approximations we used half-sums of the operators defined in the expression (2.9). For example, a half sum of the first two operators of (2.9) is 2^5 [/(x + e-ih + ejh) - / (x + e . ; / i ) - / (x + ejh) + /(x) + / (x - eth - ejh) - / ( x - eth) - / (x - ejh) + / ( x ) j . Which half-sums are taken is explained in Section 2.2.4. The choice must give us a system matrix of positive type. 2.2.2 Differences of the T o o l b o x of L e v e l Set M e t h o d s 1.1 The difference approximations of second derivatives implemented in the Toolbox of Level Set Methods 1.1 differ from our central differences of Section 2.2.1 in their computation of the mixed partial derivative approximation. In place of upwinding, the Toolbox of Level Set Methods 1.1 implements - ^ ( / ( x + e f / i ) - - / ( x - e i / i ) ) as the approximation of the first order term in direction e;, which we recognize as our first order approximation operator B;/(x). Let the first order partial derivative approximation in the variable i be as above. Then applying to it the approximation in the variable j yields Wi [sTI + e i h + e J h ) ~ /(X ~ e i , } ' + eJh)) ~ 2ir(/(X + e i k ~ e J k ) ~ /(X - e i k - eJh))] = -^j | / (x + e . ; / i + ejh) - f(x - eji + ejh) - / (x + eth - ejh) + / (x - e-Ji - e J / i ) j . This difference approximation of the mixed partial derivative is not monotone. Therefore, schemes using these differences in general are not monotone. This is discussed further in Section 3.2.2. The non-monotonicity resulting from using the above non-monotone mixed partial derivative approximation can give a non-convergent scheme. This is discussed further for the case of General degenerate parabolic equation. 11 2.2.3 C h e c k i n g der iva t ive approx ima t ions It is well known that the second order difference operator 3 a is the correct approximation There exists a simple check which demonstrates that the four discretizations of d\d2 (2.9) given in Section 2.2.1 are correct and ensure the convergence of numerical solutions to the unique solution of considered IVPs for PDEs. Here we show the calculation for one of possible four approximation of (2.9). A check for the other three approximations can be done using exactly the same procedure. Our check is based on the the necessary and sufficient condition for convergence [8] below. For any polynomial of the form P(xi,X2) = t^ l-Tl + C02X2 + OJ3XiX2 expressions d\d2P(x\,X2) and 3I2P ( :EI, .T2) must have the same value equal to LO3. It is obvious that d\d2P(x\,X2) is LO3, so let us check that 3i2-P(^i, 2:2) has value UJ3 for the example of Section 2.2.1 where two applications of the backward mixed partial difference operator were applied to derive the operator 3y, i ^ j . First we apply the backward differential operator in the first variable l ( p ( x ) - P ( x - e 1 / l ) ) = . = i [UJ\XI + u)2X2 + U3X1X2 - (u>i(xi - h) + W2X2 + OJ3(XI - h)x2^j = i (u>1h + (AJ3X2I1J and then the backward difference operator in the second variable. This gives i \ji(oJih + u!3x2h).- ±(uih + u3(x2 - /i)/i)j A similar check can be done for any of the other-three possibilities from (2.9). Therefore, we conclude that the operator 3y is a good approximation of the operator dLdj for i ^ j. 2.2.4 Difference approx ima t ions and matr ices of pos i t ive t ype Using our discretizations, the quadratic operator a\i(d\)2 + a.22(<92)2 is approximated by the differences an 3 n + a-22 3 22 • Values of ay, i ^ j, dictate the choice of approximation for a^didj. If 0,12 > 0, then ai20ic*2 should be approximated by the half sum of the first two possibilities of (2.9), and otherwise by the half sum of the second two possibilities of (2.9). Thus when ai2 > 0 we use the following approximation of the mixed partial derivatives 2 ^ [/(x + eth + e,jh) - / (x + eth) - / (x + e,-/i)4-/ (x - eih - ejh) - / (x - eth) - / (x - ejh) + 2/(x)j, 12 while when a\2 < 0 we use 2^7 [ - / ( x + e{h - ejh) + / ( x + eji) + / ( x - ejh) - / ( x - e ; / i + ejh) + / ( x - eth) + / ( x + e ^ / i ) - 2 / ( x ) j . An approximation based on a different choice can result in a non-convergent method. For further discussion on this subject see Section 3.2. Discussed approximations of first and second order terms give us a system of ODEs 4u(t) = -Au(t), u(0) = u 0 at where A is the n x n system matrix, and uo is the initial condition. If we were to write out the terms of A, we would see that it has positive diagonal entries and non-positive off-diagonal entries. In other words matrix A is of positive type. For the two dimensional case where A = —^aijdidj using centered differences from (2.9) results in entries of the system matrix corresponding to the grid node (k, I) of the following form Akiki. - 2h-2[au + a 2 2 - K2I] Aklk±u = -h-2[an - \al2\] ^ ^ Akiki±i = -hr2[a22 - \a12W Akik±u±i = — h~2\ai2\ (= Akik±u^i). In the case that values in square brackets are positive the resulting matrix with entries Atj ki is of positive type, i.e. it has positive diagonal entries .and non-positive off-diagonal entries. Otherwise, schemes are not monotone. Motzkyn arid Wasov [4] suggest for non-monotone schemes to perform rotations. There is also another possibility. In order to obtain monotone schemes one can use stencils having different length in x\ and x2 directions. This idea is thoroughly pursued in this thesis. Let us apply the described approximation procedure for the general case of the differential operator A to Degenerate heat equation. For this equation the operator is 2 2 A = - ^ oij di dj + '"i Qi • (2-12) •y=i «=i In the present example the corresponding diffusion tensor a = {a^}^ is " 1 0 " 0 0 The off-diagonal entries in the system matrix corresponding to the grid node (/c, I) have the form /l _ /,-2„ u-i j ~vk+u f° r Vk+U Aki.k±u - ~h a.u - h < 1 c 2 I vk-ll fOT Vk-U < 0 A _ L - l / -Vll+l i m V t l + \ >0 Aklkl±l — —n, ' - v l ^ iorv2 >0 " ( 2 1 3 ) Jki-x 1U1 uki-i 13 where k, I, G N identify nodes on the grid k units in direction e\h, and I units in the direction e2h, and velocity at a node {k.l) is v^i = {vli,vlt). The above matrix A is of positive type if the velocities are chosen in such, a way to ensure that Akt k±u and Aki ki±i are non-negative. The diagonal entries of the system matrix must be equal to or greater than the negative sum of the off-diagonal entries in order to ensure monotonicity. In this thesis we have chosen to set the diagonal entries equal to the negative sum of the off-diagonal entries. Let us look at another example. On the elliptic differential operator of General degenerate parabolic equation 2 A = - "ij<h'h (2-14) ij=i one can apply the described approximation procedure when the corresponding diffusion tensor-has the form " 1 1 1 1 Observe that det(a) = 0, which means that the problem is degenerate. We can say that a discretization of General degenerate parabolic equation relies on the appropriate choice of the discretization of the mixed partial derivatives. What "a good choice" is, depends on the diffusion tensor. 2.3 Courant-Friedrichs-Levy (CFL) condition for mono-tone schemes We now present our derivation of the C F L condition for monotone schemes. Let us remind ourselves of the structure of a matrix A of positive type. It must have positive diagonal entries, negative or zero off-diagonal entries and the sum of entries in each row must be either zero or positive. Our differential operator (2.3) must be discretized by a matrix A, where A must be of positive type. The system matrix of the initial value problem for ODEs ^u(t) + Au(t) = 0 can be rewritten in the following form. Let A = D - B, where D is a diagonal matrix with entries du = an > 0, and B is a non-negative matrix with entries fey = —ay > 0, bu = 0. Set p = max; du and rewrite the matrix A as • A = PI - Q, 14 where the matrix Q has entries = bij for i ^ j and qa = p — du. Then the above initial value problem for ODEs has the form ju(t) - pu(t). + Qu(t) = 0. Using Forward Euler method we get a system for discretizations u^fc^ of u(ifc) U(fc+D = ( !_ p r ) u ( fc ) + T Q u ( * ) = o. (2.15) Then the C F L condition is determined by the maximal time-interval, TCFL, TCFL = P"1-For matrices A of positive type so determined TQFL is theoretical, and consequently optimal. The Toolbox of Level Set Methods 1.1 calculates the timestep by considering terms of a partial differential equation separately, meanwhile making sure not to incorrectly set the time step to the minimum of the advection time step and the diffusion time step, as shown in the following example. Let A = Ai + A2 where both matrices Ar are of positive type, Ar = prI — Qr, and p1 = p2. Then a bound of the maximal time interval r for the system (2.15) from the corresponding time intervals rr = l/pr is obtained as follows / 1 1 \ 1 . T < mm — , — I = —. Vpi p2J p However, the maximal time interval TCFL for A is 1 TCFL < 7T-2p So the actual maximal time step is twice smaller than the time step computed as a minimum of the Ti and T 2 . In other words the scheme is unstable for r = ^, since r > TCFL-Let us now correctly state the above computation of the maximal time step. L E M M A . 2.1 Let A = A\ + A2 where both matrices Ar are of positive type, Ar = prI — Qr. Then a bound of the maximal time interval TCF'L for the system (2.15) can be obtained from the corresponding time-intervals TR = l/pr as follows ( 1 1 TCFL > 1 W r 2 A proof of this result follows simply from p <pi +p2- Therefore, taking r = + ^•)7 1 makes the method stable. We have to point out that this is not the largest possible "stable" time step. It follows that the optimal value of TCFL should be calculated numerically from the diagonal entries du rather than using this lemma. Let us show an application of this lemma. In the case of level set equation in one dimension the matrix entries da are h~l\v{xi)\ so that the associated T\ must be h T\ = max v{xi)\' 15 This is the theoretical value of T \ . This T\ depends on the coefficients of differential operator and discretization procedure (upwinding in the considered case). In the case of level set equation in two dimensions we have theoretical values T\ , T2 and a bound on TCFL from the lemma. However, the theoretical TCFL is generaly larger and can be obtained directly from the entries da. Let us point out again that r depends on the coefficients of the elliptic differential operator A defined by (2.3) and the discretization procedure used for the approximation of the first and second order terms. For example, assume that we use discretization that results in the matrix entries (2.11) and .additionally that the system matrix is of positive type. Then TQFL = 1/p implies _ IS 2 maxfci [an(xki) + a22{xki) - |ai2(x f c/)|] ' where xki are grid nodes. The obtained expression coincides with the well-known expression TCFL = / l 2 (4a 2 )" 1 for A = - a 2 £ 2 = 1 5.2. 2.4 Consistency of monotone schemes We have shown in Section 2.3 how to ensure stability of a monotone scheme. Here we discuss the consistency of the monotone schemes investigated in this thesis starting with a definition from Strikwerda, [8]. D E F I N I T I O N 2.4.1 Given a-partial differential equation Pu = f and a finite difference scheme, P,u,dxV = f, we say the finite difference scheme is consistent with the partial differential equation if for any smooth function 4>(t, x\, x2) P4> — Pdt,<ix<t> —* 0 o.s dt —> 0, dx = (dxi, da;2) —> 0, the convergence being pointwise convergence at each grid point. In our case P = dt + J2ivi^i ~ J2i3 " o ' V ^ and Pdt.dx = Forward Euler difference operator + Upwinding difference operator — Second order difference operator. i In the class notes [6] it is shown that explicit Euler difference method, and Upwinding method are consistent. Therefore, we have to demonstrate that our Second order difference operator is consistent with Oijdidj. We denote them by Pdt,dx and P, respectively. We already know that both operators give the same result when applied to second order polynomials T2(.Ti, 2:2). Hence, if 4>(x1,x2) = T2(xi, x2) + R3(dxx, d,x2)., .Ti = hk 4- dxx, x2 = hi + dx2, at a grid node (hk, hi) where T2 is a second order Taylor polynomial, then P^-PdtM = PR3.- Pdt,dxR3-The remainder R3 is the series of monomials of arguments dx\ and dx2 of the order 3 and higher. Second order partial derivatives of R3 behave like 5 = \dxi\ + \dx2\ for small values of 5, proving the consistency. 16 Chapter 3 Analysing Problems 2. and 3 In this Chapter we focus on Degenerate heat equation and General degenerate parabolic equation using the tools developed in Chapter 2. 3.1 Degenerate heat equation We write Degenerate heat equation where v is a velocity field. This parabolic equation is interesting since it has a convection and a diffusion term, where the diffusion is present only in the direction of the x-axis. We study the equation for three separate velocity fields for which analytical solutions can be found. 3.1.1 G o a l and p r o b l e m set t ing Our goal is the implementation of monotone schemes for the three problems defined below. We consider these three cases in order to investigate the efficiency of monotone schemes for an equation involving constant and vortex velocity fields. We look at two examples of the constant velocity field, since in one of the examples the velocity is aligned with the coordinate axis, while in the other example it is not. We begin with the definition of our two examples for Degenerate heat equation which we will call Problem A, and Problem B. Problem A. Consider domain D = (0, l ) 2 , with initial condition '(/.,. (x) + v v u - uxx = 0 (3.1) {x,y) = g(y)f2sm(-Kx), where sin(27T7/) for y £ [0,1/2] 0 for y £ [ 0 , 1 / 2 ] . 17 The velocity is constant v = (0,1), and the solution to (3.1) is u(t,x,y) = exp(—-n2t)g{y — i^i) \/2 sin(7r:r). The function u(t,x,y) is actually defined on a strip (0,1) x ( — 0 0 , 0 0 ) C K 2 , and has zero values at x = 0 as well as at x = 1. In addition, the function u is different from zero only for 0 < y — V2't < 5 . If we consider the problem for t S [0, then we can claim that solution has zero values on the whole boundary. Problem B. Consider again domain D = (0, l ) 2 . where velocity is field given as in [7], (dyO, —dx8), by the so called stream function 0 = 8(x,y) = - sin2(7rrt) sm2{iry). n The initial condition is defined by oj(x.y) = sin(7r.r) sin(7rjy), and solution to (3.1) is given by u(t,x,y) = - exp(-i7T 2)w(.T, y). We have the boundary condition 0 everywhere on the boundary. The analytical solutions to the three problems are constructed by using methods of eigenfunctioiis [IS], [16]. 3.1.2 D e r i v i n g ana ly t i c a l solut ions In this subsection we explain how we arrived at the analytical solutions of the three examples forementioned. Method of eigenfunctions Our objective is a construction of solutions of IVP for P D E iH + v V i i - uxx = 0 (3.2) where v is a constant velocity field, v = (vi,v2). The method is known by the name "ex-pansion of solution in terms of eigenfunctions'''. We use the following differential operators H = at + P2 - a2, K = dt + vxdx - dl One-dimensional eigenvalue problem Let us consider the following eigenvalue problem (- dl + p2)</; = X6, 4>{0) = </>(!) = 0. (3.3) 18 Solutions are <t>k{x) = 2 1 / 2 sin;/;'-,:;. A f c = p2 + fcV. The functions (fik make an orthonormal basis of L2(0,1)- The orthogonality can be easily' checked. Consider first the scalar product (<pk — 4>i) f ° r k K (4>k\4>i) = 2 / sin(/c7ra;) sin(/7ra;)(i.'c Jo = 2 j ^ cos(/c7r.'j; — ITTX) — cos(kiTX + ITTX)^ dx 11 1 - ii 1 = — K sin((fc - l)ivx) - — sin((fc + l)irx) = 0. lo7r(A;-/) u ' ' Io7r(fc + 0 v v Now let's look at the case k = I. /•1 sin(fc7r.x) sin(A;7ra;)o!a; = 2 I'sm2(knx)dx= 2 Z' 1 . 1 ~ cos(2kirx) ^ = 2\'lx + J - l ' sm(2k*x) = 1. Jo Jo 2 lo 2 4/=7rlo In the above calculations we used the trigonometry formulas for the product of two sin functions [21]. In terms of the inner product (• | •) [10] this means {(f>k\4>i) = Ski-It can be proved that any ./^-function can be approximated arbitrarily well by a finite linear combination of basis functions 4>k [16]. Therefore, "any" / on (0,1) can be represented by the series k where ck = (4>k\f) are called Fourier coefficients. Equality of two sides holds point wise for sufficiently smooth / . Therefore, we try to represent a solution of Hu(t, x) = 0 as a series i(t,x) = <*kit)4>k{x), where ak(t) can be considered as the the Fourier coefficients of u(t, •). If we apply H from the left onto u, and demand Hu = 0 we get (dt + A t ) ak{t) = 0, giving us the only possibility a.k(t) = ck exp ( — \kt) where ck 6 R. In a step by step computation we derive this result in the following way dtu{t,x) = Y^dtctk{t)<t>k(x), k p2u(t,x) = ^2p2ak(t)<t>k(x), k -d2:a(t:x) = ~J2Mt)dlMx) = ^^ak(t)k2ir2(j)k(x), 19 Hu = Y, [dtc*k(t) + Xkak(t)^k(x) = 0. k This gives us an ordinary differential equation in a(t) dta(t) + Xkak(t) = 0, which has the solution ak(t) = ck exp(—Xkt). Thus w(£, x) = 5 3 cfc exp ( - A fc t) <t>k{x), fc is a general solution of Hu = 0. Next we consider an extension (- dl + v x d x ^ = Xip, V(0) = V-(l) = 0. (3.4) and try to apply the same construction. Let us make an ansatz tp(x) = exp(vix/2)cj>(x). ff we plug ip(x) into the formulated eigenvalue problem (3.4), we get the problem (3.3) for 4> where p2 — v2/A. Hence, the eigenfunctions are 'ipk(x) = exp(vix/2)4>k(x) and the corresponding eigenvalues are as before. Therefore, we now write u as a Fourier series of the new eigen-basis u(t,x) = Y^(t)Mx)eM~), fc and try to solve Ku = 0. dtu(t, x) = 5 3 dtf3(t) exp{^y )<f>k(x), fc dxu(Lx) = E^W[y'exp(^)^(.x-)+exP(^)0:(.x-)] = £ f t ( t ) e x p ( ^ ) [^k(x)+<p'k(x)], fc fc Lx) = ^ / 3 f c ( £ ) ! l e x p ( ^ ) [ ^ ^ ( . T ) + 4>'k(x)] + pk(t)exp(^) '^fc(x) + 4'k(*)], / f u = ^exp(^)[a t /? f c ( t )^( .xO+/? f c W (^ F C ( A ; ) -^( .T))] =0. (3.5) fc , Since we know that <f>k — — \/2k2'K2sm{kixx) — —k2Tr2cj)k(x), the above calculation (3.5) can be written as Ku = 53exp (^) [c? t / 3 f c ( f ) +pk(t)(?± + k2^)]</>k(x) = 0. fc ff we identify ^ as p2, then the equation (3.5) can be rewritten as Ku = 5 3 e x p ( ^ ) [dt(3k(t) + \kpk(t)]Mx) = 0-Similarily, as we did for ak(t) we now have an ordinary differential equation in Pk(t), with solution pk(t) = c fc exp(-Afei), 20 and thus u(t,x) = ^2exp(-\kt) exp(vlx/2)4>k(x). (3.6) fc The eigenfunctions i/'fc('';) a L e not orthonormal. Yet they define a basis, because the function exp in front of <j>k is bounded from bellow and above by the numbers exp(=p|-t;1 /2). This gives us a possibility to built the above solutions of Ku = 0 in terms of ipk- In order to show that ipk are not orthonormal it is enough to show that either f^ip^dx yt. 1, or Jo i>k'i>i dx 0. Here we show that 'ip2. dx ^ 1 in general. Again, for the integration of a product of exp and sin we consulted [21]. W>A#fc) = 2 / exp{vix) sin2(kirx)dx •Jo , , 1 - COS (2/CTT.T) , = 2 J- exp{vxx) i '-dx l exp('Uia;) — exp(^iX') cos(2£;7r2;)dx o = (— exp^a:)! 1 - [ T^toK (vi cos(2A.-7r.T) + 2knsin(2knx))\1} ivy lo lv( + Akiixi JoJ = | - ( e x p ( , 1 ) - l ) - [ v 2 + A k ^ 2 j}-The construction of solutions in terms of the eigenfunctions 4>K of the boundary value problem (3.4) is general so it does not depend on a particular choice of the boundary conditions. As an alternative to using the boundary conditions of (3.4) one can impose zero derivatives at the end points of the interval (0,1). Also, zero value can be set for a subset of the end-points of the domain, while zero derivative on the remaining end-points. This would, for example, require eigenfunctions <f>k with zero derivative at x = 0 and zero value at x = 1. For this particular case the sequence of orthonormal eigenfunctions is Mx) = V-2 cos ( £ * ± f e ) , . = 0,1,. . and the corresponding eigenvalues are A. . + , = 0,1 ' • 3.1.3 T w o - d i m e n s i o n a l eigenvalue p rob lems To solve the 2-D problem (3.2) we utilize the same technique as for one-dimensional-prob-lems. Thus we start with oo u(t,x,y) = fk(t,y) exp(vix/2)<pk(x), 21 where fk{t,y) are the new Fourier coefficients, and apply from the left the differential oper-ator M = K + v2dy = 8t + Vldx + v2dy - d\ in order to calculate Mu = 0. The following result is obtained: Mu{t,x,y) = (dt + v2dyya(t,x,y) + (yidx - <92) u{t, x, y) = (dt + v2d^u{t,x,y) + Xku{t,x,y) =0 = Efc (°t + v2dy + Xk\ fk{t,y) exp{v1x/2)<pk{x) implying (dt + v2dy + A f e) fk(t,y) = 0. The only non-trivial solution of this problem is fk(t, y) = exp {-Xkt) gk(y - tv2), where gk is a function on R. Hence, oo u{t, x, y) = ^ e x P ( " A<t *) cJk{y ~ tv2) exp(vix/2) (f>k{x). (3.7) fc=i Apparently, the constructed solution is on the set [0, oo) x (0,1) x K. 3.1.4 Solu t ions on squares ff all the functions gk have supports in a set [0, a] C [0,1] then for times 0 < t < (1 — a)/v2 the packages gk travel over the domain (0, f )2 from their initial supports towards 1 and eventully hit the right end-point at T = (1 — a)/v2. Within t 6 [0, T] solution (3.7) has zero values at the boundary of [0, l ] 2 . An example function is u{t,x,y) = exp{-Xkt)x{y ~v2t) Qxp{vxx/2)(j>k{x), where ' f sin(2™) for x G [0,1/2] ' U ' \ 0 for x i [0,1/2]. Therefore, when executing the code for the constant velocity field for Degenerate heat equation the running time of the simulation can be set to at most 0.5. 3.1.5 S o l u t i o n w i t h the vor t ex ve loc i ty field We consider an fVP for the second order P D E of the following general form dt u{t, x) + v{x)Vu{t, x) + Au(t, x) = 0 for t > 0, x G D, u(t,-)\dD = 0 for t > 0, (3.8) u{0,x) = u0(x) for x G D, 22 where D C R 2 is a bounded domain with the boundary dD, x >—> v(x) is a divergence free velocity field on D, UQ is a function on D and A is a second order differential operator (degenerate or non-degenetate) of the following form 2 A = - ]T 0; d,. ij=l with a constant diffusion tensor a^. One of possibilites for the construction of solutions u(t, x) in terms of known func-tions is given by the following two rules 1) Solution has the general form u(t,x) = fi(t)f2(x). 2) The function f2 has the property v(x)^ f2(x) = 0 and Af2(x) = Xf2(x). Let us point out that Af2 = Af2 means that f2 is an eigenfunction of A with the corre-sponding eigenvalue A. When one inserts the assumed form of the solution into the P D E the following result is obtained dtfi(t) + A/ i ( i ) = 0, -u0 (a;) = f2(x). Thus, fi(t) = exp(—At), the function f2 must be an eigenfunction of A and the demanded equality vVf2 = 0 is ensured by choosing the divergence free velocity field of the form vi(x) = d2f2(x), • v2(x) = -dif2(x). We can conclude this brief description by the following conclusion. One has to construct explicitly in closed form an eigenfunction of the problem - E i , = i 0, i)j o[x) = X(p(x), 4>(x) = 0, for x e dD. For simpler domains such as a square or rectangle and diagonal diffusion tensors this problem can be easily solved. Let us consider the case D = (0, l ) 2 . Any such eigenfunction is the product of two functions, <f>(x) = <fii{xi)<fi2(x2), and the functions 4>i must have zero values at 0 and 1. In addition there must hold d24> = number x <j>. Only sin functions satisfy these conditions. Now one can easily utilize this procedure. Given vortex velocity field v = (—d2'ip, dyip), where 'ip(x.y) = — sin2(7ra;) sin 2(7rj/), a solution to (3.8) is the following function u(t,x,y) = exp(—7r2i) sin(7r.x-) sin(7ry). The following simple fact is used in the constuction of a solution to (3.8). If u(t,x,y) = exp(—Xt)x{x,y) is a solution to Af2 — Xf2 defined in 2) then we have vVu = 0 for any velocity field v = (d2'4>, —di'ip) where the stream function is a polynomial of the function x- In our case we chose the second order power of X-, L e - iJ = A'2/71"-23 3.1.6 Findings The list of programs implementing Degenerate heat equation can be found in Appendix. For a detailed explanation of the in language implementation the user can consult the manual of the Toolbox of Level Set Methods 1.1 and visit Section 5.1 of this thesis report. The figures pertaining to our findings are included at the end of this section. The C F L bound computed in the C implementation is the theoretical bound as discussed in Section 2.3. Therefore, we use C implementation to generate errors and figures of analytical and numerical solutions at different time points of the time interval [0,2']. For the constant velocity field examples, the final time T , has to be < 1/2 in order for the boundary conditions to be satisfied. Note that for Problem B the final time can be any in R + . For both problems the number of nodes taken in each spatial direction is 100, and two sets of four images display the analytical and the numerical solutions at different times. Figures 3.1 and 3.2 show the results of the implementation of Problem A while Figures 3.3 and 3.4 show the results of the implementation of Problem B. The errors for the two problems, A and B, are calculated in the norm and are displayed to the console for a user-defined number of time points of the time interval [0,2']. Each line of the program's output contains three numbers. T'he first number is the time, the second number is the absolute error, and the third number is the relative error. For Problem A errors are maximum error at time 0.000000 i s 0.000000 ( in 7.) 0.00000000 maximum error at time 0.166668 i s 0.0II709 ( in 7.) 8.57906057 maximum er ror at time 0.333336 i s 0.002948 ( in 7.) 11. 1912132/ maximum error at time 0.500004 i s 0.000665 ( in 7.) 13.08010965 For the vortex velocity case errors are maximum error at time 0.000000 i s 0.000000 ( in 7,) 0.00000596 maximum error at time 0.166663 i s 0.000349 ( in 7.) 0.18079901 maximum er ror at time 0.333327 i s 0.000139 ( in 7.) 0.37364951 maximum error at time 0.499990 i s 0.000041 ( in 7.) 0.56455006 The solutions go quickly to 0 as time increases, since they contain exp(—tir2) term. T'he time step taken is le—5 using C F L number 0.8, and the initial conditions for the two problems are their respective solutions at time t — 0. Numerical errors shown above, and illustrated by graphs, enable us to say that the monotone scheme of Degenerate heat equation appears to be better suited for an approx-imation of problems with the vortex velocity field than the constant velocity field. The error of 13% for Problem A seems rather large. Looking at this figure alone we could doubt the quality of accuracy of monotone schemes. Especially since taking a smaller -time step produces almost the same error. Still, the error of 0.56% of Problem B leaves us in the need for further investigation of their accuracy. What can be said with certainty is that the maximal time step calculated from the TCFL bound produces a stable method in each of these three problems. 24 Figure 3.1: Analytical solution of Problem A computed for the final time 0.5 , bound-ary condition 0. velocity v = (0,1), grid size [100 x 100], and initial condition sin(27ry)sin(7rri;)/v/2 V.-t £ domain D, and V-y e [0,1/2]. 25 Figure 3.2: Numerical solution of Problem A computed for the final time 0.5 bound-ary condition 0. velocity v = (0. 1). grid size [100 x 100], and initial condition sin(27ri/)sin(7rx)/v /2 Va: £ domain D, and Vi/ € [0,1/2]. 26 Figure 3.3: Analytical solution of Problem B computed for the final time 0.5 , bound-ary condition 0, vortex velocity field, grid size [100 x 100], and initial condition sin(7ra;) sin(7rj/) Vx,y 6 domain D. analytic solution, t = 0 analytic solution, t = 0.16666 analytic solution, t = 0.33333 analytic solution, t = 0.49999 0 0.5 1 0 0.5 1 27 Figure 3.4: Numerical solution of Problem B computed for the final time 0.5 boundary condi-tion 0, vortex velocity field, grid size [100 x 100], and initial condition sin(7ra) sin(7r.x) Va;, y £ domain D. numeric solution, t = 0 numeric solution, t = 0.16666 0 0.5 1 0 0.5 1 • numeric solution, t = 0.33333 numeric solution, t = 0.49999 0 0.5 1 0 0.5 1 28 3.2 General degenerate parabolic equation Next we consider the third equation of interest, u t(x) - uxx(x) - 2uxy(x) - uxx(x) = 0. (3.9) This parabolic equation was obtained by adding time dependency to a particular example, d2u - {uxx + 2uxy + uyy) = - — = 0, v = (1,1) (3.10) of [2]. The example (3.10) was designed to show how choice of a difference scheme has an impact on the convergence of a method. The author of [2] demonstrates that the difference scheme Tjiuix + h,y + h) - 2u{x,y) + u(x - h,y - h)) (3.11) produces a convergent method of (3.10), while the scheme jir (2W(.T + h, y) + 2u(x, y + h) + 2u{x - h, y) + 2u[x, y — h) — 6u(x, y) — u{x + h, y — h) — u(x — h, y + h)) results in a non-convergent method. Both schemes are central difference schemes. The first scheme takes central differences in the direction (1,1), while the second in the direction (-1,1). Oberman in [2] suggests an additional way to derive a convergent scheme for the above problem via coordinate axis rotation. This is an alternative, less elegant than central differences. 3.2.1 G o a l and p r o b l e m set t ing Our goal is to confirm or contradict the findings of [2] by analysing associated parabolic equation (3.9). The findings of [2] state that the direction in which the central differences are taken has impact on the convergence of a method. Again, we'created implementations in both C programming and MATLAB's m languages., The m language implementation is a routine in the Toolbox of Level Set Methods 1.1 and therefore can not be run as a stand-alone program. The starting point is construction of an analytical solution. It is easily obtainable from the boundary conditions given in [2]. Therefore, we can now postulate the initial boundary value problem for General degenerate parabolic equation. Consider P D E (3.9) on the domain D = (0, l ) 2 , with initial condition u(0,x,y) =0, in the interior of the domain, and boundary condition sin(67r (x — y)) \/(x,y) G dp, t > 0 29 Apparently, the approximation converges to the classical solution u(t,,x,y) = sin(67r (x - y)) as 1: —i oo. We must bring into the correspondence the schemes (3.11) and (3.12) with the schemes of Section 2.2.1. Recall that (3.9) coincides with General degenerate parabolic equation, where the degenerate differential operator is A(x) = ~^2ij=iaij(,x)didj, and the elliptic diffusion tensor is atj = 1 for all i, j 6 {1, 2}. In accordance with the proposed discretizations of Section 2.2.1 we have at our disposal four possibilites for dido, defined by dijii. 2/i 7'u(x ± e;/i ± ejh) — u(x ± e;/i) — u(x ± ejh) + u(x), — '«(x ± eih =F ejh) + v/(x ± e;/i) + u(x =F ejh) — u(x). Due to a 12 > 0 we have to use the first possibility, |^ u(x ± e.;/i ± ejh) — u(x ± e;/i) — u(x ± ejh) + tt(x)J, to arrive at a converging method. When we apply the prescribed discretization to the differential operator . aijdidj we get duu(x) + 2 B,y«(x) + 3JJU(X) = = ^ |u(x + eih) - 2u(x) + w(x - e.;/i)] + + 2 ^'«(x + e.;/i + ejh) + u(x — eih — ejh) — u(x + eih) —u(x — eih) — u(x + ejh) — u(x — ejh) + 2-u(x)j + -f-^ y u(x + ejh) — 2u(x) + u(x — ej/i)j = = £ u(x + eih + ejh) — 2ti(x) -I- -u(x — e t/i — ejh) This is precisely the convergent scheme (3.11). If we use intentionally the second possibility to discretize the operator didj, i ^ j, we get a scheme which is not monotone B,,'«(x) + 2 Biju(x) + Ejjii,(x) = = ^ ['«(* + e,/i) - 2-u(x) + -u(x - e ; / i ) j + + 2 2775 ^ — w(x + e.;/i — ejh) — u(x — e;/i. + ejh) + u(x -f eih) + u(x — e^h) + u(x — ejh) + u(x + ejh) — 2u(x)j + u(x + ejh) — 2u(x) + u(x — e,-/i)J = 2 u(x + eth) + 2u(x - e;/?,) + 2 u(x + ejh) + 2 u{x - ejh) + + 6 u(x) + —w(x + eth — ejh) — u(x — eih + e^/i)]. _ 1 ~ h 30 I Actually this is the non-convergent scheme (3.12). The absence of monotonicity results in the method being non-convergent. For' further discussion on monotonicity see the next section. The code in the Toolbox of Level Set Methods 1.1 handles a completely separate discretization of mixed partial derivatives, discussed in Section 2.2.2. It applies operator 3if first in the variable i and then in the variable j. The resulting approximation operator is ^ 2 (u(x- + eih + ejh) — i'.(x — e,/i + ejh) — u(x + e /^i — ejh) + u(x — e-^h — ejh) Using this operator the discretization of the considered J2ij <>•<). is a.„u(x) + 2 r3jj'u(x) + 3JJU(X) = = jiy [u(x + eth) - 2-u(x) + u(x - e.(/i)j + (u(x + e;/?, + ejh) - u(x - eji + ejh) - (3.13) —u(x + eih — ejh) + i/,(x — e.;/i — ejh)) + + jp [u(x + ejh) — 2u(x) + 'u(x — ej/i)j. It is easy to check that the resulting scheme is not monotone. 3.2.2 M o n o t o n i c i t y of the schemes In this section we consider monotonicity of the schemes (3.11), (3.12), and (3.13) and show that only (3.11) is monotone, while (3.12) and (3.13) are not monotone. The analysed elliptic problem is defined by (3.10). For the definition of monotonicity refer to Section 2.1. We first consider scheme (3.11) to approximate the equation (3.10). In order to check the monotonicity of the approximation we write 1-1 i j c lS Si function of its neighbouring nodes, i.e. its arguments, Uij = /('«(,fc | W(,fc £ N(uij)), 2uij = ui+xj+i +u.i-ij-i. Since both coefficients next to and ' U ; _ X J _ I are positive, Uij is nondecreasing in each of its arguments, and therefore scheme (3.11) is monotone. Next we write the scheme (3.12) in the same form 6u-ij = 2iti+ij + 2 « i j + i + 2ui-itj + 2'(tjJ_i — Ui+itj_\ — Since u^j is decreasing in 'u; + i j_ i and t t i _ i J + i , the scheme is not monotone. Finally, we look at (3.13) 16uj,j = 4ui+ij + 4u i_u + u i + i j + i - ui-ij+i - iLi+ij_1 + Ui_itj-i + Auij+i + 4w i j_ 1 . This scheme is non monotone since uitj is decreasing in and u l + i ^ - \ . The following can be said for the mixed partial derivative approximations imple-mented in the Toolbox of Level Set Methods 1.1. Any scheme implemented in the Toolbox 31 which uses Toolbox's mixed partial derivative approximation is not guaranteed to be mono-tone. This comes from the fact that schemes examined are built from schemes of individual terms. For example, Heat equation is build from two individual schemes. One for the u x x term, and one for u,jy term each of which happen to be monotone. Therefore, the scheme for Heat equation which adds these two monotone schemes is monotone. Let us now examine the individual schemes existing in the Toolbox of Level Set Methods 1.1. The Toolbox of Level Set Methods uses upwinding for the approximation of first order derivatives, which is a monotone approximation. It uses the differences approx-imation operator, da, for the non-mixed second order term, which is also monotone. Yet, it uses the following expression 2ai2 -^(ui+ij+i - — + W j - i j - i ) for the approximation of the mixed partial derivatives for the diffusion tensor a = {ay} 2 2-As readily visible, depending on the sign of ayo two of the terms always have positive values, while the other two always have negative values. Therefore, once this approximation is combined with monotone approximations for upwinding and Bu, the result is a non-monotone approximation. 3.2.3 R o t a t i o n of coord ina te axes for G e n e r a l degenerate pa rabo l i c equa t ion As mentioned in Section 2.2.4 a way to obtain a monotone scheme would be to perform a coordinate axes rotation of a problem [4]. Here we consider this approach. . Let us recall the definition of the elliptic differential operator for a two dimensional problem 2 A = - a-ijdid,, (3.14) ij=l where the diffusion tensor a = {ay} 2 2 is positive semidefinite. Thus, either it is definite and has two real positive eigenvalues or it is semidefinite with one eigenvalue equal to zero implying det(a) — 0. In the original coordinates, which are denoted here by a: 1,2:2, the partial differentiation.operators are denoted by d{. We define new coordinates, x[. x'2 as X[ = t U X i + tX2X2; x'2 = £ 2 i xi + t.22 X2, and the corresponding partial differentiation operators d[. The following transformation between differentations in the original and new variables can be established di = t n d [ + t 2 1 d 2 , d2 = d[ + t22 &2-32 This can be rewritten in matrix form 0i d2 TT To obtain the elliptic operator A in new coordinates, we have to use these transfor-mations of partial differential operators. Therefore, we insert expressions 6\ in terms of d[ into (3.14) to obtain 2 • 2 A = - Y *vhid'kthdl = - J2 <ud'kd'h ijkl=\ fc(=l where the new diffusion tensor has entries a'ki = 5Z* f c*a'J'*'J = ( T a T l ) ki-ln the case of det(a) = 0 we have det(a') = det(T)det(a)det(TT) = 0, i.e. A is degenerate in the new coordinates as well. We know from Strang [10] that there exists a rotation defined by [ cos(a) sin(a) 1 ( a) = • / > ( \ [ — sm(a) cos(a) such that the matrix a1 has one of the following forms: " p 0 " " 0 0 ' _ 0 0 _ 0 p . For non-degenarate diffusion tensors there exist difference schemes which discretize the elliptic differential operator (3.14) by matrices An of positive type. The idea is given in [4]. There are many elaborations such as [20]. In order to get discretizations An of positive type one has to use stencils of'different length in the xx and x2 directions. This idea we use for degenerate schemes to develop the method of Chapter 4. 3.2.4 F i n d i n g s The C F L bound on the schemes which are monotone was computed as given in Section 2.3. The C F L bound for the non-monotone schemes cannot be computed using our described procedure. Still, 'we used the same C F L bound as the one for monotone schemes, namely hi 2 maxfci [an(xu) + a 2 2 (x f c ( ) - |a12(ccfc()|] This time step when used with non-monotone schemes produces erroneous results which show non-convergence of these schemes. We must note that for parabolic problems with the zero initial condition and nonho-mogeneous boundary conditions we do not show error values. Instead we show the difference between the numerical solution and the analytical solution at a point in time. In addition, 33 when we say that a solution is analytical, in some cases it is really an asymptotic solution. These cases are the ones for which the initial condition is not equal to the analytical solution at time zero. • ' • Figure 3.5 shows the result of using discretization (3.11) when running the imple-mentation up to time 0.2 for grid size 100 x 100. The solution driven by the boundary conditions is quickly approaching the analytic solution as can be seen in the figure. The differences in norm are maximum difference at time 0 000000 is 0 998867 (in 7.) 100.00000000 maximum difference at time 0 024997 is 0 919144 (in 7.) 92 01866571 maximum difference at time 0 049995 is 0 706837 (in 7.) 70 76385770 maximum difference at time 0 074992 is 0 529321 (in 7.) 52 99214182 maximum difference at time 0 099990 is 0 395380 (in 7.) 39 58280635 maximum difference at time 0 124987 is 0 296519 (in 7.) 29 68551002 maximum difference at time 0 149985 is 0 222830 (in 7.) 22 30829146 maximum difference at time 0 174982 is 0 167454 (in 7.) 16 76441874 maximum difference at time 0 199980 is 0 125842 (in 7.) 12 59842562 Figure 3.6 shows the result at time 1.0 for both the analytic and numeric solutions. The differences are maximum difference at time 0.000000 is 0.998867 (in 7.) 100.00000000 maximum difference at time 1.000026 is 0.000135 (in 7.) 0.01350979 An additional test was run with the initial condition being the solution in order to assert the validity of the implementation, and thus potentially observe the calculational error on a 100 x 100 grid with final time 0.5. maximum error at time maximum error at time maximum error at time maximum error at time maximum error at time maximum error at time maximum error at time maximum error at time maximum error at time 0.000000 is 0.000000 0.124987 is 0.000000 0.249974 is 0.000000 0.374962 is 0.000000 0.499949 is 0.000000 0.624936 is 0.000000 0.749923 is 0.000000 0.874911 is 0.000000 0.999898 is 0.000000 (in %) 0.00000000 (in 7.) 0.00004028 (in 7.) 0.00004028 (in 7.) 0.00004028 (in 7.) 0.00004028 (in 7.) 0.00004028 (in 7.) 0.00004028 (in 7.) 0.00004028 (in 7.) 0.00004028 We conclude that the 4e—5 percent error is due to round off error. Figure 3.7 shows the same setting as in Figure 3.5 using discretization (3.12) on a 100 x 100 nodes grid for final time 1.2. It is readily visible how the non-convergence is materialized and remains with the running time. The differences obtained are maximum difference at time 0.150000 is 0.826923 (in 7.) 82.7861089311 34 maximum difference at time 0.300000 i s 0.906501 ( i n '/,) 90.7528988305 maximum difference at time 0.450000 i s 0.900932 ( i n */,) 90.1953603932 maximum difference at time 0.600000 i s 0.861522 ( in 7.) 86.2498673530 maximum difference at time 0.750000 i s 0.665061 ( i n 7.) 66.5815557683 maximum difference at time 0.900000 i s 0.946972 ( i n 7.) 94.8045391399 maximum difference at time 1.050000 i s 1.081382 ( i n 7.) 108.2608615755 maximum difference at time 1.200000 i s 0.941579 ( i n 7,) 94.2646321863 When compared to the results of discretiazation (3.11) it is plain that discretization (3.11) produces a convergent method, while discretization (3.12) doesn't. Figure 3.8 shows the differences obtained by using the discretization (3.13). The non-convergence happens due to the scheme being non-monotone. The differences for a grid with 100 x 100 nodes and final time 1.2 are maximum difference at time 0.150000 i s 0.492477 ( i n 7.) 49.3035721126 maximum difference at time 0.300000 i s 0.456759 ( i n 7.) 45.7276746742 maximum difference at time 0.450000 i s 0.454310 ( i n 7.) 45.4825515735 maximum difference at time 0.600000 i s 0.451359 ( in 7.) 45.1870767231 maximum difference at time 0.750000 i s 0.450665 ( i n 7.) 45.1176407953 maximum difference at time 0.900000 i s 0.461395 ( i n 7.) 46.1918515911 maximum difference at time 1.050000 i s 0.465639 ( i n 7.) 46.6167256667 maximum difference at time 1.200000 i s 0.546891 ( i n 7.) 54.7510753188 Let us point out once more that schemes (3.12) and (3.13), which are not monotone, result in a non-convergence, while the monotone scheme (3.11) is convergent. We also investigated what would happen on a grid with periodic boundary condi-tions when domain is (—1, l ) 2 using the Toolbox of Level Set Methods 1.1 implementation of the boundary calculation, addGhostPeriodic .m. The initial condition is sin(67r(a; — y)) everywhere within the domain, and the central differences approximation is (3.11). The result for the final time 1.2 on a grid with 100 x 100 nodes using the discretization (3.11) surprisingly produces a non-convergent scheme. This can be seen from the percentage errors computed at nine different time points from 0 to 1.2 maximum error at time 0.150000 i s 0.521210 ( i n 7.) 52.1801004934 maximum error at time 0.300000 i s 0.514805 ( in 7.) 51.5388485854 maximum error at time 0.450000 i s 0.376802 ( in 7.) 37.7229136800 maximum error at time 0.600000 i s 0.377520 ( in 7.) 37.7948549638 maximum error at time 0.750000 i s 0.428266 ( i n 7.) 42.8751591476 maximum error a t ' t ime 0.900000 i s 0.522602 ( i n 7.) 52.3194115335 maximum error at time 1.050000 i s 0.395414 ( i n 7.) 39.5862203809 maximum error at time 1.200000 i s 0.555725 ( i n 7.) 55.6355194335 Figure 3.9 shows the above example at four time points. Further on, we tested the discretization (3.13) with the periodic boundary condi-35 tions of the Toolbox of Level Set Methods 1.1. For a grid with 100 x 100 nodes, and final time 0.2, the result is shown in Figure 3.10. The errors generated are maximum error at time 0.025000 is 0.565357 (in '/.) 56.5998184766 maximum error at time 0.050000 is 0.961642 (in 7.) 96.2732027284 maximum error at time 0.075000 is 0.952116 (in 7.) 95.3195222936 maximum error at time 0.100000 is 1.004376 (in 7.) 100.5514906582 maximum error at time 0.125000 is 1.256236 (in 7.) 125.7660333185 maximum error at time 0.150000 is 1.016355 (in 7.) 101.7507094230 maximum error at time 0.175000 -is 1.133993 (in 7.) 113.5278499173 maximum error at time 0.200000 is 1.076474 (in 7.) 107.7694694336 We conducted a test with boundary conditions sin(67r(.x + y)) on a grid of 100 x 100 nodes for final time 1.0 with discretizations (3.11) and (3.12). We expected that if we were to use discretization (3.11) we would obtain a non-convergent method, while if we were to use discretization (3.12) the method would be convergent. Our results depicted in Figures 3.11, 3.12 show our findings, and the following difference figures confirm our expectation. Differences for the discretization (3.12) are maximum difference at time 0 000000 is 0 998867 (in 7.) 100.00000000 maximum difference at time 0 333333 is 0 027840 (in 7.) 2.78714085 maximum difference at time 0 666667 is 0 000715 (in 7.) 0.07162458 maximum difference at time 1 000000 is 0 000135 (in 7.) 0.01350979 while those foi discretization (3.1.1 ) are maximum difference at time 0 125000 is 0 709033 (in 7.) 70.9837158757 maximum difference at time 0 250000 is 0 894754 (in 7.) 89.5768376666 maximum difference at time 0 375000 is 1 061009 (in 7.) 106.2211725151 maximum difference at time 0 500000 is 0 937179 (in 7.) 93.8241215025 maximum difference at time 0 625000 is 0 697497 (in 7.) 69.8287873810 maximum difference at time 0 750000 is 1 031804 (in 7.) 103.2973533648 maximum difference at time 0 875000 is 0 703944 (in 7.) 70.4741875121 maximum difference at time 1 000000 is 0 864096 (in 7.) 86.5075695024 We conducted convergence studies of the examples with asymptotic solution sin(67r(.x'— v/)), sin(67r(.x —2y.)), and sin(67r(3.x — 2y)). The results obtained in convergence studies show non-convergence of the first two examples, and a jump in the convergence of the third example. The negative values appear due to the fact that the numerical solution at a time point is compared to an asymptotic solution, rather than the analytic solution. In order to obtain true convergence studies' values, we would have to let time go to infinity. Since we execute the implementations up to a finite time t, both convergent, and divergent results are possible. 36 Our findings show that montone schemes are stable and in most cases very accurate approximations of the analytical solutions. 37 Figure 3.5: Numerical solution for the final time 0.2 for the discretization (3.11), boundary condition sin(67r(.T — y)), grid size [100 x 100], and initial condition 0 Vx,y. 38 Figure 3.6: Analytical and numerical solutions for the final time 1.0 for the discretization (3.11), boundary condition sin(67r(.T - y)), grid size [100 x 100], and initial condition zero. 39 Figure 3.7: Final time 1.2 for the discretization (3.12), boundary condition sin(67r(a; — grid size [100 x 100], and initial condition 0 Vx,y. numerict = 0 numeric t = 0.15 numeric t = 0.3 0 0.5 1 0 0.5 1 0 0.5 1 40 Figure 3.8: Final time 1.2 for the discretization (3.13), boundary condition sin(67r(.x -y)), grid size [100 x 100], and initial condition zero for all x,y. numeric t = 0 1 i • 0 ' ' 0 0.5 1 numeric t = 0.45 0 0.5 1 numeric t = 0.9 0 0.5 1 numeric t = 0.15 0 0.5 1 numeric t = 0.6 0 0.5 1 numeric t = 1.05 0 0.5 1 numeric t = 0.3 0 0.5 1 numeric t = 0.75 0 0.5 1 numeric t = 1.2 0 0.5 1 41 Figure 3.9: Final time 0.2 for the discretization (3.11), periodic boundary conditions, size [100 x 100], and initial condition sin(6Tr(x — y)). numerict = 0 numeric t = 0.4 - 1 0 1 _ 1 0 1 numeric t = 0.8 numeric t = 1.2 -1 0 1 -1 • • 0 1 42 Figure 3.10: Final time 0.2 for the discretization (3.13), periodic boundary conditions, size [100 x 100], and initial condition sin(67r(.r — y)). Figure 3.11: Numerical solution for the final time 1.0 for the discretization (3.12). boundary condition sin(67r(.-i; + y)), grid size [100 x 100], and initial condition 0\/x,y 6 D. numeric t = 0 1 | • 0.8 0.6 0.4 0.2 0 L _ . 0 0.5 numeric t = 0.67 numeric t = 1 0 0.5 1 0 0.5 1 numeric t = 0.33 1 0 0.5 1 44 Figure 3.12: Numerical and analytical solutions for the final time 1.2 for the discretization (3.11), boundary condition s'm(6n(x+y)), grid size [100 x 100], and initial condition 0 V.T, y € D. numeric t = 0 1 i • 0.5 0 0.5 1 numeric t = 0.75 0 0.5 1 numeric t = 0.125 0 0.5 1 numeric t = 0.5 0 .. 0.5 1 numeric t = 0.875 numeric t = 0.25 0 0.5 1 45 Chapter 4 Build ing monotone schemes In this chapter we present our general procedure for building monotone schemes from positive semidefinite and positive definite diffusion tensors. We restrict our research to PDEs with purely diffusive terms such as General degenerate parabolic equation. 4.1 Central difference approximations Let a be an m x n matrix and D be the Hessian, i.e. the square matrix of second order derivatives. Then we have 2 trace^aDa1) = trace(oL aD) = ajjdjdj, ij=i where the entries of the matrix a = a1a are denoted by ay, and the matrix a is positive definite or semidefinite. Discretizations of differential operators in Section 2.2.1 are given for the case of equal spatial step size h, in both discretization directions. We now define these discretizations with spatial steps /ij and h2. where h\ is the step size in the x-axis direction, while / i 2 is the step size in the y-axis direction. The corresponding expressions for the forward and backward difference operators are 3 J ( x ) - + l ( / ( x ) - / ( x - e A ) ) -Then r3;;/(x) has the form as before with /ij instead of h. Similarity, the mixed partial differential operator can be approximated by one of the following four possibilities 9 i 3 j / ( x ) - a i j / ( x ) = 1 f / (x ± e A ± ejhj) - / (x ± e.,;/?,) - / (x ± eff^) + /(x), (4.1) hih2 1 - / ( x ± eihi T ejhj) + /(x ± e,^) + /(x T ejhj) - /(x). We construct discretizations of the operator a 1 1(9i) 2 + 2ai2<9i<92 + a.22(c>2)2 by u s ' n g o.ii3n + 2o.i2 3i2 + a22 322 a s i n the previous case. The resulting entries of the system matrix are 46 Au k±u Au. u.±i Aid k±u±i = — 2 \hl 2an + h2 = h~l\h~lau hi h2 \ai2 (4.2) now -2 <-2 a22 -h~l\al2\ h~l[h~la22 - h~l\a,i2\ / i f ' / i2 l\ai2\ (= Auk±u^i) In order to get monotone schemes we must have system matrices of positive type. This implies that the following conditions must be satisfied hi l a n /i 2" 1a 22 h-2l\ai2\ >0, h~l\ai2\ >0. (4.3) Let us analyse these conditions. We know that the matrix a is positive definite or semidefinite. In the former case we have a\2 < a u a 2 2 and in the latter case a\2 = ano.22-This conclusion follows from the conditions det(a) > 0 and 'det(a) = 0, respectively, ff the matrix a is positive definite, there exist plenty of possibilities of choosing the pairs hx,h2 in order to satisfy (4.3). For instance if we fix h\ > 0, then there are as many h2 as there are real numbers. If a is semidefinite. then the choices are very restricted. If we fix hi > 0. then / i 2 is also fixed since it is calculated from the expression h2 hi l f l12| an (4.4) assuming that an ^ 0. The above ratio was computed knowing that in case of a degeneracy the equations (4.3) are one and the same equation for which a strict equality holds as the following computation shows. Take (4.3) and plug into each one |a,i2| = N / a , u a 2 2 , h i v/ 0 ' ll°'22 h o > 0 a22 ~h~2 s/ana22 hi > 0 Now divide the one on the left with v/aTT / a n , and one on the right with Ja22, h hi V°22 h2 > o s/0,22 h2 Win >o, hi - > /022 h2 \^2~2 h2 - hi ' In other words But this means that /Q-22 v a n hi / a n hi /«22 hi l a u - h2 l\ai2\ = 0, h2 l a 2 2 - hx l\an\ = 0. The reason for equation (4.4) is simple. There exists only one rotation such that in the new coordinates the diffusion tensor has the form 1 0 0 0 47 Thus we can draw the following conclusion. Instead of rotations we use the original coordinate system and spatial grid steps of different sizes. Let x = (fc, I) be a grid node. Then the system matrix has non-trivial entries for those grid nodes which are on the cross through x and additional two grid-knots as illus-trated in Figure 4.1. These sets can be called numerical neighbourhoods. 4 . 2 Computational procedure 1. Given a calculate a = o1 a. 2. Let M l be a given number of grid intervals in the x-axis direction, which means that there are M l + 1 grid nodes in the ei direction. Then h\ = 1/M1 is the corresponding grid step size. 3. Calculate K2I K2I ?'i = , 'f'2 = • ° ' l l °'22 If both rk < 1, then the spatial differential operator is of the second order in more than one coordinate directions so take h2 = hi and construct the system matrix. 4. Otherwise, calculate the number and set h 2 = ghi. • The number g is calculated from (4.3) 5. The number p, = l/h2 is not an integer so one needs to take the nearest integer M2. This integer determines the number of grid intervals in the second direction. In this way the domain D is equal to a rectangle of sides h\Ml and h2M2. This rectangle is very close to the unit square for fine grids. 6. If Q 4 2 > 0 the system matrix is constructed by using the numerical neighbourhoods depicted in Figure 4.1(a), and for 'a\2 < 0 one uses the numerical neighbourhoods depicted in Figure 4.1(b). Figure 4.1: The numerical neighbourhoods used in the calculation of the terms of the system matrix A in (4.2). (a): 012 > 0 (b): al2 < 0 48 4.3 Problem definition We would like to test out implementation on three problems. Two degenerate and one non-degenerate. We call them Problem C, Problem D, and Problem E. All three problems are examples of General degenerate parabolic equation. The P D E analysed in Problem C is ut - {Auxx + Auxy + uyy) = 0, (4.5) while the one for Problem D is ut - {4uxx + l2uXtV + Quyy) = 0. (4.6) We define these two problems next. Problem C. Consider domain D = (0, l ) 2 , equation (4.5) with analytical solution u(t,x,y) = sin(67r(a; - 2y)), where boundary condition is given by the function u{0,x,y) = sin(67r(a: - 2y)), and zero initial condition. We analyse numerical solution as t —> co. Problem D. The domain is D = (0, l ) 2 , the equation is (4.6), the analytical solution is u(t,x,y) = sin(67r(3.-r - 2y)), the boundary condition is given by the function ' u(0,x,y) = sin(67r(3rt - 2y)), and the initial condition is zero. Again, we analyse numerical solution as t —> oo. The difference between the above two problems and Problem E is in the fact that Problem E is designed to test our implementation with a positive definite diffusion tensor. The P D E studied in Problem E is ut - {a,lxuxx + 2o i 2 uXiV + a 2 2 U y y ) = 0, (4.7) with the corresponding diffusion tensor a = {ay}?? , "1 0.5 " 0.5 1 We chose to write the P D E in this way to point out that we could have used any other positive definite diffusion tensor of the form a b b a ' 49 for the solution (1-2/) u(t, x, y) = -xK- + yy- J-L a\\ a22 to still hold. Given our choice of the diffusion tensor we define Problem E as follows. Problem E. Consider again the domain D = (0, l ) 2 , equation (4.7) with the solution V i is u(t, x, y) = -x{l - x) + 2,(1 -y). - (4.8) Numerical boundary condition on the top and bottom borders of the domain is — x(l — x), while the one on the left and right borders of the domain is y{\ — y). We analyse numerical solution as t —> oo. 4.4 Findings Al l three problems were tested on the domain discretized by Nl = 100 grid nodes in the x-axis direction. The number of nodes in the y-axis direction is computed from Nl as discussed in Section 4.2. For Problems C and D we generate the results using both non-equally and equally spaced grids. Problem E we will investigate only on an equally spaced grid since its associ-ated diffusion tensor requires an equally spaced grid. We must note that for problems C, D, and E the difference figures do not show the actual numerical error. Instead they show the difference between the numerical solution and the analytical solution at a point in time. This is why the percentage change decreases with the increase in the running time. The differences of Problem C generated for unequal grid step sizes and final time 0.1 are maximum difference at time 0.000000 is 1.000000 (in 7.) 100.00000000 maximum difference at time 0.033330 is 0.341236 (in 7.) 34.12359953 maximum difference at time 0.066660 is 0.091538 (in '/.) 9.15384889 maximum difference at time 0.099990 is 0.024577 (in '/,) 2 .45770812 As the final running time of the simulation increases, the difference decreases. Fig-ures 4.2 and 4.3 depict the numerical and the analytical solutions for Problem C respectively for the final time 0.1. Figure 4.4 shows the result obtained when using equally spaced grid for final time 0.25. We include the following differences generated by the equally spaced grid scheme for the final time 0.25. maximum difference at time 0.000000 is 0.998867 (in 7.) 100.00000000 maximum difference at time 0.031247 is 0.493938 (in '/,) 49.44984886 maximum difference at time 0.062494 is 0.341225 (in 7.) 34.16121397 maximum difference at time 0.093740 is 0.310776 (in 7.) 31.11286388 maximum difference at time 0.124987 is 0.304707 (in 7.) 30.50522153 50 maximum difference at time 0.156234 is 0.303492 (in 7.) 30.38361548 maximum difference at time 0.187481 is 0.303254 (in 7.) 30.35983013 maximum difference at time 0.218728 is 0.303245 (in 7.) 30.35888731 maximum difference at time 0.249975 is 0.303245 (in 7.) 30.35888731 The relative difference remains close to 30.35%, as illustrated in Figure 4.4 and evident from the following differences maximum difference at time 0.000000 is 0.998867 (in 7.) 100.00000000 maximum difference at time 0.333333 is 0.303245 (in 7.) 30.35888731 maximum difference at time 0.666667 is 0.303245 (in 7.) 30.35888731 maximum difference at time 1.000000 is 0.303245 (in 7.) 30.35888731 Problem D on an unequally spaced grid is shown in Figure 4.5. The difference fig-ures for the Problem are shown below for final time 0.25. maximum difference at time 0.000000 is 0.999998 (in 7.) 100.00000000 maximum difference at time 0.083333 is 0.003905 (in '/,) 0.39048283 maximum difference at time 0.166667 is 0.003208 (in */.) 0.32080485 maximum difference at time 0.250000 is 0.003208 (in 7.) 0.32080485 Figure 4.6 illustrates the numerical solution to Problem D computed on an equally spaced grid. Again, the equally spaced grid scheme produces a non-convergent method as can be seen from the following differences. maximum difference at time 0.000000 is 0.998867 (in 7.) 100.00000000 maximum difference at time 0.083333 is 0.689878 (in */,) 69.06607524 maximum difference at time 0.166667 is 0.689878 (in 7.) 69.06607524 maximum difference at time 0.250000 is 0.689878 (in 7.) 69.06607524 This is also visible in Figure 4.6. We could say that a lower bound on difference is 69.06%. Figures 4.7 and 4.8 show the numerical and analytical solution to Problem E for final time 0.05. In Figure 4.7 it can be seen how the solution is being formed. We ran the same setup for final time 0.1 to obtain the following difference values maximum difference at time 0.000000 is 0.239976 (in 7.) 100.00000000 maximum difference at time 0.033330 is 0.006672 (in 7.) 2.78036850 maximum difference at time 0.066660 is 0.000312 (in 7.) 0.12986445 maximum difference at time 0.099990 is 0.000015 (in 7.) 0.00604490 From the presented tests we conclude that methods with unequal step sizes give far better results than methods with equal step sizes. 51 Figure 4.2: Numerical solution for the final time 0.1 using unequal grid step sizes dis-cretization, boundary condition sin(67r(.x — 2y)), grid size [100 x 201], and initial condition O V i . y e D . 52 Figure 4.3: Analytical solution for the final time 0.1 using unequal grid step sizes discretiza-tion, and grid size [100 x 201]. analytic .t = 0 analytic t = 0.033' 0 0.5 1 0 0.5 1 analytic t = 0.067 analytic t = 0.1 0 0.5 1 0 0.5 1 53 Figure 4.4: Numerical solution for the final time 0.25 using equal grid step sizes dis-cretization, boundary condition sin(67r(x — 2-t/)), grid size [100 x 100], and initial condition 0\/x,yeD. 54 Figure 4.5: Numerical and analytical solutions for the final time 0.25 using unequal grid step sizes discretization, boundary condition sin(67r(3a: — 2y)), grid size [lOOx 68], and initial condition 0 Vrt, ye D. 55 Figure 4.6: Numerical solution for the final time 0.25 using equally spaced grid step sizes, boundary condition sin(67r(3.T — 2j/)), grid size [100 x 100], and initial condition 0 V.T, y 6 D. Figure 4.7: Numerical solution for the final time 0.05, with boundary condition —x on the top and bottom ends of the domain, o n the left and right ends of the domain, and grid size [100 x 100]. numeric t = 0 numeric t = 0.017 0 0.5 1 0 0.5 1 numeric t = 0.033 . numeric t = 0.05 0 0.5 1 0 0.5 1 57 Figure 4.8: Analytical solution for the final time 0.05, with boundary condition —x ' on the top and bottom ends of the domain, y on the left and right ends of the domain, and grid size [100 x 100]. analytic t = 0 analytic t = 0.017 0 0.5 1 0 0.5 1 58 4.5 Analysing Problem 4. Now that we have developed our method of unequal step sizes in we can look at the final equation of interest, General degenerate parabolic equation with advection, u t(x) + v v '« - (uxx + 2uxy + uyy) = 0. (4.9) This equation is interesting since it combines convection with diffusion terms discretized using our unequal step size method. The first order term vy i t is discretized using upwinding, while the second order terms using By developed in Section 2.2.1. 4.5.1 P r o b l e m D e f i n i t i o n In accordance with our alphabetical naming we call this example Problem F. Problem F. Consider P D E (4.9) over domain D = (0, l ) 2 . The solution is u(t, x, y) = sin(67r(.T — (y — V2t))) for v\ = 0, where t is time and (x,y) are grid-node coordinates. The boundary condition is computed at every time step from the value of the analytical solution. The velocity is constant v = (vi,vo). We have chosen v = (0,0.5). 4.5 .2 F i n d i n g s We ran the implementation of Problem F until final time 1.0 on D discretized by 100 x 100 nodes. The differences are maximum difference at time 0.000000 is 0.000000 (in 7.) 0.00000000 maximum difference at time 0.333333 is 0.091789 (in 7.) 9.18926738 maximum difference at time 0.666667 is 0.093522 (in 7.) 9.36283003 maximum difference at time 1.000000 is 0.093567 (in 7.) 9.36730545 The difference is noticeable in Figure 4.9 for t = 0.33, for exampje. A closer look at difference values shows that the magnitude of velocity drives the magnitude of percentage difference. I.e. the larger the velocity is, the larger is the percentage difference between the numeric and analytic solutions. To show this we include the following difference values for final running time 0.5 and velocity v = (0,1) on D maximum difference at time 0.000000.is 0.000000 (in 7.) 0.00000000 maximum difference at time 0.166667 is 0.153738 (in 7.) 15.39127710 maximum difference at time 0.333333 is 0.170976 (in 7.) 17.11701003 maximum difference at time 0.500000 is 0.171830 (in 7.) 17.20244873 These differences are depicted in Figure 4.10. The difference figures can appear surprising at first. However, the error magnitude of h coming from the discretization of 59 first order derivatives when multiplied by the velocity v does not render the differences completely unexpected. We also note that the percentage differences decrease with the increase in the number of grid-nodes. Table 4.1 shows the results of the convergence study of Problem F. Since the example does compare the numerical solution at a time point to the actual analytical solution the values in the table confirm the method's convergence. Table 4.1: Example of the IVP sin(67r(.x- - (y - v2t))) for final time 1.0, and v = (0,0.5). JV1 Loo r L i r L 2 r 50 0.172639 0.04204 0.05861 100 0.093567 0.88 0.02254 0.90 0.03122 0.90 200 0.050461 0.89 0.01505 0.58 0.01839 0.76 60 Figure 4.9: Numerical solution for the final time 1.0, with boundary condition sin(67r(x-(y — V2t))), where v2 = 0.5, grid size [100 x 100], and initial condition sin(67r(x- — y)). -numeric t = 0 numeric t = 0.33 0 0.5 1 . 0 0.5 1 numeric t = 0.67 numeric t = 1 0 0.5 1 0 0.5 1 61 Figure 4.10: Numerical solution for the final time 0.5, with boundary condition sin(67r(.7; (y — v2t))), where v2 = 1, grid size [100 x 100], and initial condition sin(67r(.T — y)). 62 Chapter 5 Software implementat ion In this chapter we take a close look at the way we have implemented our programs in MATEAB ' s M language. The implementations are located within files with extension m, for example filename.rn. Therefore, we call such a file m, file,. A list of programs created in this thesis can be found in Appendix where it is also explained how to execute the programs. 5.1 Implementation details of m files The Toolbox of Level Set Methods 1.1. package developed by I. Mitchell [1] is a set of in files that are ready to be executed and edited by the user. The package is distributed from Http://www.cs.ubc.ca/ mitcflell/ToolboxLS/index.Html and comes with a manual explaining its contents and usage. In order to understand the m language implementations of this thesis it is important to understand the simple convective flow example explained on the pages 14 to 27 of the manual. The thesis m file implemen-tations follow the layout of the convective flow example. In this section we explain on one example implementation the most iniortant parts of the code under separate headings. The other implementations have layout very similar to the one of the example implementation. Let our example implementation be the m program GDPdiricHlet .m which imple-ments the monotone discretization (3.11) of General degenerate parabolic equation for the boundary values sin(67r(n; — y)) and initial condition zero on a cartesian grid. The P D E solved by this problem is (3.9). and an example run of this implementation is depicted in Figure 3.5. A reader will notice that we use the name of an m file to name the function the file contains. In effect, the name of the function can be any. Execution of M A T L A B functions is done by a call to the filename containing the function, not by the function name itself. This, however, is true only if a file contains a single function. Function calls within an m file are clone by function names. Most of our in files contain just one function and so when refering to a function we interchange the file name with the function name. When specifying file names below we use asterisk (*) in the usual UNIX wildcard sense. I.e. * means zero or more characters. 63 BOUNDARY CALCULATION The boundary values are computed by using two functions located in two separate m files. The first m file calculates the (x,y) coordinates of the boundary nodes of the grid. It then calls the other m file which contains the function to be used to calculate the value of the numerical solution on the boundary nodes given their (x,y) positions at time t. In GDPdirichlet .m example the function that calculates the (x,y) coordinates of grid-nodes is called computeBoundary.m. The function that computeBoundary .m calls is named boundaryFunction.m. boundaryFunction.m contains function sin(67r(a; — y)) to be used in the calculation of the numerical solution for (x, y) position at time t. boundaryFunction.m is designed to operate on a set of values, rather than indi-vidual (x,y) positions to reduce the computation time. A user should not feel the need to edit-the computeBoundary .m file. This is due to the fact that the Toolbox of Level Set Methods 1.1 works only with rectangular domains without "holes", and changing the computeBoundary .m would change the look of the domain. The boundaryFunction.m in the GDPdirichlet .m implementation is set up to hold the analytical solution to the problem, so it should not be changed either. COMPARISON OF SOLUTIONS The analytical solution to a problem is kept in a separate file. In the case of GDPdirichlet. m the solution is kept in the file solutionGDP .m. The function that actually does the compari-son of the numerical and analytical solutions, compareSols .m, is called by GDPdirichlet .m. Any file containing an analytical solution can be edited by the user. However, the solution for the GDPdirichlet ,m is set to the one solving the P D E (3.9), so the user should not feel the need to change it. DERIVATIVE DISCRETIZATION • The discretization of the first and second order terms, &\,<9y, is computed by a separate function called hessian* .m. All of the hessian* .m functions follow the same layout. In the. case of GDPdirichlet .m the file is called hessianDegenerate .m. hessianDegenerate .m computes the derivative approximations in two directions: bottom-left to top-right, and bottom-right to top-left. Bottom-left to top-right direction was used for P D E (3.9) where the analytical solution is sin{Qix(x - y)). Bottom-right to top-left direction was used to generate the example with the analytical solution .S'm(67r(x + '</))• THE term* FUNCTIONS There are several term*.m functions to be found in the Toolbox of Level Set Methods 1.1, and we have build ours following the same layout. The integrator functions can be found in 7ToolboxLS- l . l /Kernel /Expl ic i t ln tegrat ion/Term directory, where " is the directory into which the Toolbox of Level Level Set Methods 1.1 is installed. A term*.m function can be viewed as a driver of the numerical approximation. It is the file that repeatedly calls the discretization and integration routines. We mention it here since at a first glance its role in the implementation might be confusing to a novice user of the Toolbox of Level Set Methods 1.1. 64 A term function is also responsible for computing the C F L bound on the time step for every iteration. Finally, a few notes on the integrator function used in the thesis even though the implementation of the integrator function was not a part of this thesis but was already provided with the Toolbox of Level Set Methods 1.1. In all the m file implementations of the thesis, the integrator function is odeCFLl.m, which approximates Forward Euler integration. A user could try to use a different integrator (all integrator functions can be found inside of 7ToolboxLS- l . l /Kernel /Expl ic i t ln tegrat ion/ In tegrators directory), but then the resulting system matrix would not be of positive type. It is noticeable that GDPdirichlet. m is implemented in a way so as not to be changeable. So why then do we have separate files containing the analytical solution and the function to be used in the calculation of the boundary values when we could have written the analytical solution and the boundary function in GDPdirichlet .m.? We wanted to keep the same layout for all our implementations. Examples of Chapter 4 can all be tested on one m file implementation by changing the m files that contain the solution and the function to be used in the calculation of the boundary values. In this way we made the implementations easier for a user to understand. Since the Toolbox of Level Set Methods 1.1 is designed to calculate the C F L condi-tion and the derivative approximations at every time step, we had to follow this pattern. The m files contain comments on the implementation of a numerical example, but not comments on the in language syntax. .For questions regarding the m language syntax we refer the user to the MATLAB's online help website http://www.mathworks.com/access/helpdesk/help/techdoc/matlab.html, or the MAT-L A B help provided with a M A T L A B distribution. In addition, a simple google search (http://www.google.com) can return many examples. 5.2 m file v.s. C language implementations The monotone implementations of Problems 1., 2., 3., and 4. have been implemented both in the m files, as a part of the Toolbox of Level Set Methods 1.1, and as. stand-alone C programs. The non-monotone examples exist only as m file implementations. Some reasons for advocating a C implementation vs. an m file implementation stem from the fact that a C implementation is faster and requires less memory storage. In addition, the C F L condition computed in the C implementations is guaranteed to be the theoretical one. This bound is computed once, before the numerical solver is called. In this way, the time step is calculated only once and used in every iteration of the numerical solver. The arrays that contain the values of the system matrix are also calculated only once, as are the values of the boundary for examples where the calculation of the boundary values is time independent. Contrary to this, the implementations of the Toolbox of Level Set Methods 1.1 compute a guess on the C F L bound at every time step. Also, the system matrix and the values on the boundary are computed at every time step. Moreover, the boundary is allocated and deallocated at every time step. This means additional storage space is required 65 to store the values that will be copied into the boundary at every time step once the boundary for that time step gets created. Finally, the passing of structures and arrays to a function written in the MATLAB's m language is a "pass by value" if the passed in argument(s) are changed within the function. This means that every array is copied into a new memory location when being passed to a function that changes it. Thus executing M A T L A B m files can quickly appropriates a substantial amount of available memory, depending on the implementation and the size of the arrays. We include the following comparison of the m files versus C language implemen-tations with respect to memory and speed. An m language implementation of any of the explicit method examples presented in this thesis cannot be run on a 256MB R A M , Pentium 4, machine with 2GB of swap space if the grid has 100 x 100 nodes, due to thrashing. Yet a grid of 100 x 100 .can be used on the same machine using the C implemetations without problems. The domain in the Toolbox of Level Set Methods 1.1 can only be a square or a rectangle without "holes". The same holds for the domain in the current C language implementations. However, the C implementations are made in a way to allow a user-specified domain, but then the computation of the diagonal terms of the system matrix would have to be changed. We note that the implicit Midpoint Runge-Kutta example is only implemented in C language. 5.3 Implementation of C files The C implementations are not a requirement of this thesis. Therefore, their implemen-tation will not be discussed here. They were done by the principal investigator as an aid in overcoming the memory and running time difficulties of the m file implementations. In effect, most of the results shown in this thesis are obtained using C language implementa-tions. If, however, an interested reader would like to obtain a detailed description of the C language implementations as well as the implementations themselves, they should contact the principal investigator at mirnalimScs .ubc . ca. 66 C h a p t e r 6 Future work • Extension to three and n dimensions The work conducted in this thesis involves only two-dimensional problems. One can continue in the reseach of the importance of stencil directions for three and higher dimensional spaces. The three-dimensional case, we expect, would present more difficulty to program, but not to draw given the MATLAB's graphing capabilities. It is uncertain what demands on programming and displaying would higher than three dimensions pose. In order to test the implementation in three-dimensional space one would have to obtain an analytical solution to a three-dimensional problem. One could argue that a two dimensional example could be used by setting diffusion to be zero in one of the spatial dimensions. However, the question remains on whether two-dimensional examples can be utilized to test properly a three-dimensional implementation. Further research could involve the implementation of the method presented by Ober-man in [14], and its comparison to the method of Chapter 4. The idea of [14] is to use an equally spaced grid, but a wider stencil, thus allowing for more choices of the direction of dis-cretization. The discretization direction which is closest to the direction of diffusion should be taken since then the central differences approximations will be aligned with the diffusion as much as possible. We did not implement the method of [14] for a two-dimensional case. Once implemented it could be easily compared to the method of varying step sizes given the examples of this thesis. • Implicit and implicit-explicit methods The Midpoint Runge-Kutta is the only implicit method implemented in this thesis research. We decided not to include it as part of the thesis due to the fact that it was our introduction into the implicit methods that has already been thoroughly researched in published work. Further research of implicit methods could involve implementation of another equation(s) using Midpoint or other implicit Runge-Kutta methods such as Hammer-Hollingsworth and Diagonal implicit Runge-Kutta. Since the implementation of the Midpoint Runge-Kutta was only written in C programming language, the next step would be to write a Toolbox of Level Set Methods 1.1 implementation. [t is known that implementations of implicit schemes take more time to run than 67 do the implementations of explicit schemes. Speed and memory requirements would need to be investigated when considering implicit implementations within the Toolbox of Level Set Methods 1.1. Additional research could involve the implicit-explicit (IMEX) methods. The idea is to solve the second order terms in an equation via an implicit solver, and first order terms using an explicit solver. The reason being that the second order terms require a much smaller time step than do the first order terms if solved using an explicit solver. So to have as big of a time step as possible in an IMEX scheme, the second order terms could be solved by using an implicit method, while the first order terms by using an explicit method. The work of Ascher, et. al. [19] develops several Runge-Kutta based IMEX schemes. The issues remaining to be considered are the running time and memory require-ments of IMEX implementations. Since we have not implemented an IMEX scheme we cannot claim any findings. C language implementations would give insight into the speed of IMEX schemes since in C, unlike in M A T L A B , pointers to arrays can be passed between functions, and so the running time would be that of the underlying solver rather than that of the solver plus time for memory allocations and deallocations. 68 Bibliography [1] author: I. Mitchell, Assistant Professor, Toolbox of Level Set Methods 1.1, Department of Computer Science, University of British Columbia, http:/ /www.cs.ubc.ca/ mitchell /ToolboxLS/index.html [2] A. M. Oberman, Convergent difference schemes for nonlinear elliptic and parabolic equations: Hamilton-Jacobi equations and. free boundary problems, SINUM, 44: No. 2: 879-895, 2006 [3] D. H. Anderson, Comparlmental Modeling and Tracer Kinetics, Lecture Notes in Biomathematics, vol. 50, Springer-Verlag, Berlin, 1983 [4] T. S. Motzkyn, W. Wasov, On the approximation of linear elliptic differential equations by difference equations with positive coefficients, J. Math. Phys. 31: 253-259, 1953 [5] G. Barles, P. E. Souganidis, Convergence of approximation schemes for fully nonlinear second order equations, Asymptotic Anal., 4(3):271-283, 1991 [6] Class notes from the course Numerical methods taught by U. Ascher, University of British Columbia, 2005-2006 [7] D. Enright, R. Fedkiw, J . Ferziger, I. Mitchell, A hybrid particle level set method for improved, interface capturing, J. Comp. Physics, 183:83-116, 2002 [8] J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Wadsworth & Brooks/Cole, 1989 [9] W. E. Boyce, R. C. DiPrima, Elementary Differential Equations and Boundary Value Problems, 7. ed., John Wiley & Sons, 2001 [10] G. Strang, Linear Algebra and its Applications, 3. ed., Brooks/Cole, 1988 [11] U. M. Ascher, L. R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, SIAM, 1998 [12] W. H. Press, W. T. Vetterling, S. A. Teukolsky, B. P. Flannery, Numerical Recipies in C++: The Art of Scientific Computing, 2. ed., 2002 [13] J. Bramble, B. Hubbard, A theorem on error estimation for finite difference analogues of the dirichlet problem for elliptic problems, Contrib. Diff. Eq., 2:319-340, 1963 69 [14] A. M. Oberman, .4 convergent monotone difference scheme for motion of level sets by mean curvature, Numer. Math., 99:365-379, 2004 [15] S. Osher, Ft. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, 2003 [16] M. Spiegel, Fourier analysis with applications to boundary value problems (Shaurn's outline series), McGraw-Hill Companies, 1974 [17] S. M. A. Bruin, seminar talk, www.win.tue.nl/casa/meetings/seminar/previous/runge-kutta/talksem.pdf, 2002 [IS] M. S. Gockenbach, Partial Differential Equations, Analytical and Numerial Methods, SIAM, 2002 [19] U. A. Ascher, S. J. Ruuth, R. J. Spiteri, Implicit- Explicit Runge-Kutta Methods for Time-Dependent Partial Differential Equations, Appl. Numer. Math., 25(2-3):151-167, 1997 [20] N. L l M I C , M. R.OGINA, Explicit stable methods for second order parabolic systems, Math. Comraun. Vol 5, (2000), 97-115 [21] I.N. BRONSTE.IN, K. A. SEMEND.IA.IEV, Matematicki prirucnik za inzenjere i siudenle, Tehnicka knjiga, Zagreb, 1991. 70 Appendix The Toolbox of Level Set Methods 1.1 can be downloaded from http:/ /www.cs.ubc.ca/ mitchell /ToolboxLS/index.html website together with the installation instructions. The Toolboox of Level Set Methods 1.1 is a set of m files to be executed within a M A T L A B process. M A T L A B m files do not require compiling. In order to execute an m file type the name of the file followed by the list of arguments enclosed in parentheses at the M A T L A B prompt. The following is an example call to a function stored in file filename .m which takes arguments argi to arg n : f i lename(argi, a r g 2 , a r g n ) ; The description of the arguments accepted by a function is located at the beginning of the file that contains the function itself, ft is important to note that if the program called does not reside in the current directory (i.e. the directory the user is currently in), the entire path to that program needs to be specified, or the directory in which the program resides must be added to the MATLAB's search path using addpathO. Fo example, to run an m file named filename.m from outside of the directory in which it resides one would execute /path.to_toolbox/ToolboxLS-l .1/Examples/Basic/f i lename(argi, . . . , a r g „ ) ; or one could store in a file, path- f i le , the following command addpath(genpath(' /path_to_toolbox/ToolboxLS-l. l /path_to_desired_directory')); and then from a subdirectory call run( 'path_to_path_f i le /path- f i l e ' ) ; where path_to_toolbox and path_to_pathJ: i l e are the absolute paths from the root direc-tory, / , to the directory containing the Toolbox of Level Set Methods 1.1, and the path_f i l e , respectively. Before we list the programs developed as a part of this thesis we would like to comment on the writing conventions used in order to simplify the table listing of the m file implementations. We will assume that the user has installed the Toolbox of Level Set Methods 1.1 in directory dir. Then, the m file implementation of Heat equation (Problem 1) is stored in dir /ToolboxLS-1.1 /Examples/Basic /pl / directory. The m file implementation of Degenerate heat equation (Problem 2) is contained in dir /ToolboxLS-1. l /Examples/Basic/p2/ 71 directory, and those of General degenerate parabolic equation (Problem 3) in d l r /Too lboxLS- I . I /ExampIes /Bas ic /p3 / directory. The rn file implementation of General degenerate parabolic equation with advec-tion (Problem 4) can be found in directory d i r /Too lboxLS- I . I /ExampIes /Bas ic /p4 / . Having said that, we write program filenames in Table A . l without specifying the path to them. We will also use abreviations of P D E problems • H for Heat equation, • DH for Degenerate heat equation, • G D P for General degenerate parabolic equation, • GDPA~ gor General degenerate parabolic equation with advection. Problems A~-F each fall into one of the above P D E problems • A , B are under D H , • C, D, E are under GDP, • F is under G D P A . When specifying file names in the table, we use asterisk (*) in the usual UNIX wildcard sense. I.e. * means zero or more characters. OSER identifies the user of a M A T L A B session who initiates an execution of a program. When both numerical and analytical solutions to a problem are graphed, then two separate graphs are used where at a certain time ti, one graph shows the numerical solu-tion at ti, while the other graph shows the analytical solution at ti. Otherwise, only the numerical solution is graphed. 72 Table AM list of m ( M A T E A B ) files implemented in the Toolbox of Level Set Methods 1.1. Program Name P D E Problem Is called by heat2d.m H USER heat2dSolution.m H heat2d.m lid. in D H USER hdsConstant.m D H hd.m hdsVortex.m D H hd.ni add GhostPuiictiori2D. m G D P GDP* . in coinputeBoimdary.nl GDP, G D P A GDP* . in b o u n d ar y Fu n c t ion. in G D P GDPdirichlet.ni boundaryPYmctionMono.nl G D P GDPmono.m GDPdirichlet.ni G D P USER GDPmono.m G D P USER. GDPperiodic.in G D P USER hessianDegenerate.nl G D P GDPdirichlet.m, GDPperiodic.in hessianMono.m G D P GDPmoiio.m soliitionGDP.ni G D P GDPdirichlet.m, GDPperiodic.ni solutionGDPA.m G D P GDPdirichletl.m solutionMono.in G D P GDPmono.m terniDegeiierate.nl G D P GDPdirichlet.ni., GDPperiodic.in terniMono.ni G D P • GDPmono.m GDPdirichletl.ru G D P A USER coniputeBoundaryt.nl G D P A GDPdirichletl.m b ou n d ary Function 1. m G D P A GSPdirichletl.m compareSols.m A L L A E E called by U S E R visualizeEevelSetl.m A L E • A E E called by U S E R 73
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Monotone schemes for degenerate and non-degenerate...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Monotone schemes for degenerate and non-degenerate parabolic problems Limic, Mirna 2006
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Monotone schemes for degenerate and non-degenerate parabolic problems |
Creator |
Limic, Mirna |
Publisher | University of British Columbia |
Date Issued | 2006 |
Description | This thesis analyses several monotone schemes for degenerate second order parabolic PDEs over two-dimensional domains which is in most cases equal (0, l )2 . For such problems the degeneracy means that in addition to the axis spanned by the standard basis vectors, there exists another pair of orthogonal cordinate axes such that the spatial difference operator is of the second order in one of the directions, and of the first order in the remaining directions. The direction of one axis along which the spatial difference operator is of the second order we call the direction of diffusion. The thesis considers only constant coefficient PDE's and therefore the degeneracy can be easily determined by solving the eigenvalue problem for diffusion tensor. We analyse the impact of the direction of diffusion on the convergence of a scheme. Previous work on a second order elliptic problem has shown that central differences taken in the direction of diffusion produce a convergent scheme, wherease disalignment between the two result in non-convergent schemes. One of our aims was to check the validity of this finding and possibly improve the construction of schemes. As a result we present a novel approach to building monotone schemes from the diffusion tensor by taking spatial step sizes of different length in order to align the second order central differences with the direction of diffusion. We give a step by step algorithm and supply all of the findings. Our findings are based on MATLAB m file, and C language implementations. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-02-24 |
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.0052071 |
URI | http://hdl.handle.net/2429/31763 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2007-0168.pdf [ 3.24MB ]
- Metadata
- JSON: 831-1.0052071.json
- JSON-LD: 831-1.0052071-ld.json
- RDF/XML (Pretty): 831-1.0052071-rdf.xml
- RDF/JSON: 831-1.0052071-rdf.json
- Turtle: 831-1.0052071-turtle.txt
- N-Triples: 831-1.0052071-rdf-ntriples.txt
- Original Record: 831-1.0052071-source.json
- Full Text
- 831-1.0052071-fulltext.txt
- Citation
- 831-1.0052071.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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052071/manifest