@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Computer Science, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Huber, Sarah"@en ; dcterms:issued "2015-12-01T03:19:13"@en, "2015"@en ; vivo:relatedDegree "Master of Science - MSc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description "A variety of finite difference schemes are explored for the numerical solution of elliptic partial differential equations, specifically the Poisson and convection-diffusion equations. Problems are investigated that require the use of a non-uniform or non-square mesh. This may be due to a non-square domain or a problem with a singularity. We explore the properties of the linear operators in the resulting systems of linear equations. In particular, we investigate the conditioning and eigenvalues of these operators, both numerically and in search of an approximation of these eigenvalues. We also investigate the choice of finite difference scheme with respect to accuracy and cost."@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/55607?expand=metadata"@en ; skos:note "Finite Difference Schemes for EllipticPartial Differential Equationsrequiring a Non-Uniform MeshbySarah HuberB.Sc., The University of British Columbia, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Computer Science)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)November 2015c© Sarah Huber 2015AbstractA variety of finite difference schemes are explored for the numerical solu-tion of elliptic partial differential equations, specifically the Poisson andconvection-diffusion equations. Problems are investigated that require theuse of a non-uniform or non-square mesh. This may be due to a non-squaredomain or a problem with a singularity. We explore the properties of thelinear operators in the resulting systems of linear equations. In particular,we investigate the conditioning and eigenvalues of these operators, both nu-merically and in search of an approximation of these eigenvalues. We alsoinvestigate the choice of finite difference scheme with respect to accuracyand cost.iiPrefaceThis thesis is the result of work done by Chen Greif and myself. I wasresponsible for the numerical results, analysis and writeup, with input frommy coauthor.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . ix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Mathematical Tools . . . . . . . . . . . . . . . . . . . . . . . 22 Centered Differences for a Singularly Perturbed BVP . . 32.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Matrix Structure . . . . . . . . . . . . . . . . . . . . . . . . . 52.3 Discrete Solution for Uniform Mesh . . . . . . . . . . . . . . 72.4 Discrete Solution for Non-Uniform Mesh . . . . . . . . . . . 83 On the Eigenvalues of a Finite Difference Operator on aNon-Uniform Mesh . . . . . . . . . . . . . . . . . . . . . . . . . 123.1 Matrix Transformation . . . . . . . . . . . . . . . . . . . . . 123.2 Finding the Secular Equation . . . . . . . . . . . . . . . . . . 143.2.1 Approximating the Roots of the Secular Equation . . 154 Finite Differences on a Non-Square Domain . . . . . . . . . 194.1 Convection-Diffusion on Non-Square Domain . . . . . . . . . 194.2 Motivating Problems . . . . . . . . . . . . . . . . . . . . . . 204.3 Computational Schemes . . . . . . . . . . . . . . . . . . . . . 214.4 Domain with a Curved Boundary . . . . . . . . . . . . . . . 234.5 Reentrant Corner . . . . . . . . . . . . . . . . . . . . . . . . 27ivTable of Contents4.6 Alternative Discretizations . . . . . . . . . . . . . . . . . . . 315 Cell Centered Finite Differences . . . . . . . . . . . . . . . . 395.1 A Symmetric Scheme . . . . . . . . . . . . . . . . . . . . . . 395.2 Data Structure . . . . . . . . . . . . . . . . . . . . . . . . . . 405.3 Domain Definition and Refinement Strategies . . . . . . . . . 415.4 The Discretization Scheme . . . . . . . . . . . . . . . . . . . 425.4.1 Boundary Conditions and Operator Modifications . . 455.5 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . 466 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55vList of Tables4.1 Convergence of FD scheme for (4.1) on (4.6) σ = 2, τ = 0,r = 0.25, s = 1. E2 = ||U − uexact||2 and E∞ = ||U − uexact||∞. 274.2 Convergence of FD scheme for (5.2) on (4.12). E2 = ||U −uexact||2 and E∞ = ||U − uexact||∞. . . . . . . . . . . . . . . . 325.1 Convergence of FD scheme for (5.2) with (5.14) on circleshaped domain (5.9). E2 = ||U2h − Uh||2 and E∞ = ||U2h −Uh||∞. k is the largest factor of refinement. . . . . . . . . . . 485.2 Convergence of FD scheme for (5.2) with (5.14) on L-shapeddomain (5.10). E2 = ||U2h − Uh||2 and E∞ = ||U2h − Uh||∞.k is the largest factor of refinement. . . . . . . . . . . . . . . 485.3 Convergence of FD scheme for (5.2) with (5.14) on L-shapeddomain (5.10). E2 = ||U2h − Uh||2 and E∞ = ||U2h − Uh||∞.k is the largest factor of refinement. . . . . . . . . . . . . . . 51viList of Figures2.1 Exact solution to (2.1) . . . . . . . . . . . . . . . . . . . . . . 42.2 Uniform mesh with N + 1 points. . . . . . . . . . . . . . . . . 52.3 Non-uniform mesh with N+1 mesh points, two distinct spac-ings h1, h2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.4 Discrete solution for uniform and non-uniform mesh of (2.1)with  = .01, N = 22, σ = 1. . . . . . . . . . . . . . . . . . . . 102.5 Size of coefficients a, e in discrete solution for  = .03, N = 10,β1 =53 , β2 =512 . . . . . . . . . . . . . . . . . . . . . . . . . . 113.1 Eigenvalues of finite difference operator for (2.1) with uniformand non-uniform mesh. A1 and A3 are the submatrices of Adefined in (3.3). . . . . . . . . . . . . . . . . . . . . . . . . . . 163.2 Graph of f(λ) and h(λ) for N = 20,  = 0.1, σ = 1 and k = 4. 184.1 FD solution values on a uniform mesh. . . . . . . . . . . . . . 214.2 Lexicographic labelling of a uniform mesh. . . . . . . . . . . . 224.3 Lexicographic indices for a locally refined grid. Spaces cor-responding to a missing index are in yellow, and grid pointsaffected by the missing indices are in green. . . . . . . . . . . 234.4 Row and column indices for a locally refined grid. Spacescorresponding to a missing index are in yellow, and grid pointsaffected by the missing indices are in green. . . . . . . . . . . 244.5 Grid points on (4.6) with h = 0.0625, N = 176. Grid pointslying on a uniform mesh are shown in blue, and additionalgrid points on the boundary of the domain shown as red circles. 254.6 Modified FD scheme for a curved boundary with Dirichletboundary conditions. . . . . . . . . . . . . . . . . . . . . . . . 264.7 FD solution of (4.1) on the domain (4.6) with σ = 2, τ = 0,r = 0.25, s = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . 284.8 Error in FD solution of (4.1) on the domain (4.6) with σ = 2,τ = 0, r = 0.25, s = 1. . . . . . . . . . . . . . . . . . . . . . . 294.9 FD solution of (5.2) on L-shaped domain. . . . . . . . . . . . 30viiList of Figures4.10 Error in FD solution for (5.2) on L-shaped domain with auniform grid, h = 0.0156, N = 2945. . . . . . . . . . . . . . . 314.11 FD solution of Poisson problem on L-shaped mesh with Dirich-let boundary conditions. . . . . . . . . . . . . . . . . . . . . . 324.12 Locally refined grid for (4.12), hmax = 0.0625, k = 64. . . . . 334.13 Locally refined grid for (4.12), hmax = 0.03125, k = 256. . . . 344.14 FD solution values on non-uniform mesh, non-uniform meshsize. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.15 FD solution values on non-uniform mesh, grid point missingnorthern neighbour. . . . . . . . . . . . . . . . . . . . . . . . 364.16 FD operator for (5.2) on L-shaped domain with locally refinedgrid, N = 205. Grid points are labelled lexicographically. . . 375.1 Coordinates of grid points in cell centred grid. . . . . . . . . . 405.2 Entries in the Depth Map array before and after the cell inleft hand corner is refined. Cells undergoing refinement arein yellow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 415.3 Quadtree refinement according to (5.1). . . . . . . . . . . . . 425.4 Grid point and side labelling of quadtree grid. . . . . . . . . . 435.5 Flux across sides of quadtree grid. . . . . . . . . . . . . . . . 445.6 Domains considered in Section 5.5. Boundary points for thedomain are red circles, and point inside the domain greendiamonds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 475.7 Error U2h − Uh on cross shaped domain (5.7). . . . . . . . . 495.8 L-shaped domain refinement for Div-Grad scheme. . . . . . . 505.9 Condition numbers for operatorR in Div-Grad and FD schemevs. finest mesh size. . . . . . . . . . . . . . . . . . . . . . . . 515.10 Operator for Div-Grad scheme, hmax = 0.25, k = 2, N = 208. 52viiiAcknowledgementsMy thanks to my supervisor, Chen Greif, who has given good advice andmuch encouragement throughout my program. I am particularly thankfulfor his well timed and practical feedback, which was often able to set meback on the right track without feeling disheartened about having gone downthe wrong one. This work was supported financially by the Natural Sciencesand Engineering Research Council of Canada as well as the University ofBritish Columbia Affiliated Fellowships program.I am also grateful for my undergraduate supervisors, Noham Weinbergand Manfred Trummer, who gave me the chance to learn about the researchprocess and who encouraged me in my efforts.I am thankful for the support, encouragement and graciousness of myfriends, family, and church family. In particular, my husband Liam has beenpatient and supportive throughout this process. His curiosity and daring arean inspiration for my own research.Finally, I am grateful to all those who have taught me that the best wayto learn is to teach.ixChapter 1IntroductionElliptic partial differential equations (PDEs), [20] such as the convection-diffusion [16] and Poisson [20] equations, describe physical phenomena suchas the movement of a pollutant in a stream or electrostatic charge. Since theanalytical solution to these problems may be difficult or impossible to find,they are generally solved numerically. A wide range of different solutiontechniques exist, both in the discretization scheme for the differential equa-tion and the solution technique for the resulting system of linear equations.Our focus in this work will be on the choice of discretization scheme, andon finite difference schemes [20] in particular.These solution techniques are usually introduced in the context of anidealized problem or domain, usually a square domain and a uniform mesh.We are interested in the choice and complexity of the necessary changes inthe discretization scheme when the situation changes to a more complexdomain or mesh. We are also interested in the resulting changes in thesystem of linear equations, and whether we can still perform analysis on theoperator as in the simplest of cases. As most real-world applications are fornon-trivial domains, this is an important topic. We consider starting froma standard second order centred finite difference scheme [16] with a uniformmesh on a square domain, and consider various adaptations.In 1D, we investigate problems with boundary layers, which require anarea of fine resolution near the boundary for good resolution. In 2D, weexplore domains that are non-polymoino [8], that is, not made up of somecombination of squares. We also investigate solving an elliptic PDE on adomain for which the exact solution contains a singularity; e.g., the Poissonequation on an L-shaped domain.These problems all require the use of a non-uniform grid. We explorethe use of bi-trees and quadtrees [7] (the reduction of an Octree to 1D and2D respectively) for these problems and for those with non-trivial domains.These meshes have all mesh sizes as a factor of a power of 2 of the smallestmesh size. This allows for ease when performing adaptive mesh refinement,as well as in storage and labelling.When a grid in two dimensions or higher is non-uniform, the neighbour-11.1. Mathematical Toolsing points for a given finite difference scheme may not be present. In thiscase, modifications must be made to the scheme, or a different scheme mustbe chosen. We discuss the use of several discretization schemes with respectto accuracy, ease of implementation, and the properties of the operator inthe resulting system of linear equations. In particular, one property of theoperator we are interested in is symmetry; we will investigate how to choosea scheme to maintain this feature.Symmetry is a desirable attribute for a matrix, as the solvers available fora symmetric matrix are generally computationally stable, less computation-ally costly, and require less storage than those for nonsymmetric matrices.Furthermore, we are also interested in eigenvalue estimates, which may beeasier to obtain for a symmetric matrix.1.1 Mathematical ToolsCode was developed and run using Matlab [14] in Chapters 2 3, and 4. InChapter 5, much of the code used was adapted for use in Julia [2] from aquadtree implementation in Matlab by Eldad Haber [10].2Chapter 2Centered Differences for aSingularly Perturbed BVP2.1 IntroductionConsider, for 0 <  1, the convection-diffusion problem in 1D−u′′ + σu′ = S u ∈ (0, 1),u(0) = 0, u(1) = 1.(2.1)Problem (2.1) has an exact solution in the homogeneous case of S = 0, asshown in 2.1,eσx/ − 1eσ/ − 1 .For  as above and σ = O(1), this problem has a boundary layer near 1.Since this boundary layer is dependant on a potentially small value of , itis called a singularly perturbed boundary value problem (BVP).A uniform mesh with N+1 points and origin x0, as shown in Figure 2.2,includes the uniformly spaced points xj = x0 + jh, where h =1N . We defineUj as our numerical approximation to u(xj) on this grid. The standardcentred difference approximation to (2.1) is then:LUj := −Uj+1 − 2Uj + Uj−1h2+ σUj+1 − Uj−12h= Sj , j = 1, ..., N − 1 (2.2)with U0 = 0, UN = 1.As  −→ 0, the boundary layer grows thinner and more difficult toresolve, and a fine mesh is increasingly necessary for the stability of thediscrete solution. A measure of the fineness of the mesh relative to theproblem is the Pe´clet number, defined as β = σh2 . We require β < 1 forstability [15]. Without this constraint, non-physical oscillations appear andgrow large as the Pe´clet number increases.However, in the area away from the boundary layer, the solution issmooth and a fine mesh is not so necessary. For this reason, we consider32.1. IntroductionFigure 2.1: Exact solution to (2.1)discretizing the problem with a piecewise uniform mesh, with all mesh sizes2−k×h1, where h1 is the largest mesh size and k some non-negative integer.We call this discretization a “bi-tree”, in relation to the quadtree in 2D andthe octree in 3D.Specifically, we consider a mesh with two different mesh sizes; a coarsemesh with spacing h1 outside of the boundary layer and a fine mesh withspacing h2 = 2−k × h1 in the boundary layer. Each of the fine mesh andcoarse mesh contain N/2 grid points; we define the index of the N/2thgridpoint as M .We must also consider the choice of size for the two regions. Segal [19]suggests dividing the interval into [0, 1−8], [1−8, 1], with an equal numberof points in each region. For the values of h1 ∈ [0.1, 0.001],  ∈ [0.1, 0.01]and σ ≈ 1 considered, however, this does not appear to be a large enoughregion. In practice, oscillations were observable outside the boundary regionwith this division. Heuristically we find the division [0, .9], [.9, 1] to workwell for the values of , σ and h1 considered. In order to accommodate themesh size restrictions of the bi-tree, we allow the size of the two regions to42.2. Matrix StructureFigure 2.2: Uniform mesh with N + 1 points.vary slightly.2.2 Matrix StructureConsider a mesh with two distinct mesh spacings, as in Figure 2.3. Giventhe previous definition of β, we haveβ1 =σh12, β2 =σh22.Since there are two distinct mesh spacings, we have three different com-ponents in the centered difference approximation. Below, (2.3) and (2.5)correspond to the fine and coarse mesh spacings, and (2.4) to the grid point52.2. Matrix StructureFigure 2.3: Non-uniform mesh with N+1 mesh points, two distinct spacingsh1, h2.with a different left and right mesh spacing. As before, U0 = 0, UN = 1 and(β1 − 1)Uj+1 + 2Uj − (β1 + 1)Uj−1 =(h21)Sj , j = 1, 2, ...,M − 1; (2.3)(β1 − 1h2)Uj+1 +(β2 + 1h1− β1 − 1h2)Uj −(β2 + 1h1)Uj−1 =(h1 + h2)Sj , j = M ; (2.4)(β2 − 1)Uj+1 + 2Uj − (β2 + 1)Uj−1 =(h22)Sj , j = M + 1, ..., N − 1. (2.5)Therefore, our operator A isA =2 −1 + β1−1 + β1 2 −1− β1. . .. . .. . .−1 + β1 2 −1− β1−1+β2h1β2+1h1− β1−1h2−β2−1h1−1 + β2 2 −1− β2. . .. . .. . .−1 + β2 2(2.6)62.3. Discrete Solution for Uniform MeshWhen β1, β2 < 1 this is a diagonally dominant M-matrix. This meansthat its eigenvalues will be positive and real. The conclusions drawn inSection 3.1 for the eigenvalues of A will therefore focus on systems whereβ1, β2 < 1 to ensure real eigenvalues. However, we also wish to considernon-uniform meshes in order to avoid limiting β outside of the boundarylayer. Many of the possible cases considered in Section 2.4, therefore, willhave operators that are not M-matrices, meaning that their eigenvalues maybe complex.2.3 Discrete Solution for Uniform MeshConsider the one-dimensional convection-diffusion problem (2.1). Followingthe example of [15, p. 11], for a uniform mesh with N + 1 points, we havethe centered difference approximation for (2.1) as given in (2.2).We can write (2.2) as(β − 1)Uj+1 + 2Uj − (β + 1)Uj−1 =(h2)Sj , j = 1, 2, ..., N − 1. (2.7)We assume that the solution is of the formUj = ξj ,which allows us to write (2.7) as(β − 1) ξj+1 + 2ξj − (β + 1) ξj−1 = 0.This gives us the resulting characteristic equation(β − 1) ξ2 + 2ξ − (β + 1) = 0,by which we findξ− = 1, ξ+ =1 + β1− β . (2.8)This means our solution isUj = aξj− + bξj+.where a, b are determined by our boundary conditions. For U0 = 0, UN = 1,Uj =ξj− + bξj+ξN+ − 1.72.4. Discrete Solution for Non-Uniform MeshBy the definition in (2.8), when β > 1, ξ+ < 0, and terms of the form ξj+will be oscillatory. These non-physical oscillations in the discrete solutionmay be large and lead to a poor approximation of the exact solution. Whenβ ≤ 1 we have no such difficulties. If non-physical oscillations do occur, bydecreasing h we may decrease β to a value ≤ 1 and restore accuracy to thenumerical solution.2.4 Discrete Solution for Non-Uniform MeshConsider now a mesh divided into two distinct mesh spacings, as in Figure2.3. Assuming again that our solution is of the form Uj = ξj , from thecentered difference approximation in (2.3), (2.4), (2.5) we get the followingset of characteristic equations(β1 − 2) ξ21 + 2ξ1 − (β1 + 2) = 0, j = 1, 2, ...,M − 1; (2.9)(β1 − 2h2)ξ2M +(β2 + 2h1− β1 − 2h2)ξM −(β2 + 2h1)= 0, j =M ; (2.10)(β2 − 2) ξ22 + 2ξ2 − (β2 + 2) = 0, j =M + 1, ..., N − 1, (2.11)with solutionξ1− = 1, ξ1+ =1 + β11− β1 , j = 1, 2, ...,M − 1;ξM− = 1, ξM+ =h2 (1 + β2)h1 (1− β1) , j = M ;ξ2− = 1, ξ2+ =1 + β21− β2 , j = M + 1, ..., N − 1.Thus, our solution isUj = aξj1+ + bξj1−, j = 1, ..,M − 1;Uj = cξjM+ + dξjM−, j = M ;Uj = eξj2+ + fξj2−, j = M + 1, .., N − 1.When continuity of the solution is enforced at the transition point, theequation for j = M is unnecessary and we haveUj = aξj1+ + bξj1−, j = 1, ..,M ; (2.12)Uj = eξj2+ + fξj2−, j = M + 1, .., N − 1. (2.13)82.4. Discrete Solution for Non-Uniform MeshWe now need to find a, b, e, f . If we have boundary conditions U0 = 0,UN = 1, then the following two equations must be satisfied:aξ01+ + bξ01− = 0,eξN2+ + fξN2− = 1.By enforcing continuity at M we obtain the equationaξM1+ + bξM1− = eξM2+ + fξM2−.We may also specify a continuous first derivative. We enforce this byrequiring the equality of a first order finite difference (FD) at M :aξM−11+ (ξ1+ − 1) + bξM−11− (ξ1− − 1)h1=eξM2+ (ξ2+ − 1) + fξM2− (ξ2− − 1)h2.We now have four equations and four unknowns; a solvable system ofequations. Therefore our solution is (2.12), (2.13) with coefficientsa =αα(ξM1+ − 1)− ξM2+ + ξN2+;b = − αα(ξM1+ − 1)− ξM2+ + ξN2+;e =1α(ξM1+ − 1)− ξM2+ + ξN2+;f = 1− ξN2+α(ξM1+ − 1)− ξM2+ + ξN2+;defined in terms of a parameter αα =ξM2+ (ξ2+ − 1)h1ξM−11+ (ξ1+ − 1)h2.As in the case of the uniform mesh, a value of βi > 1 on any area of themesh means that ξi+ < 0. This corresponds to a solution where the termsξji+ give rise to non-physical oscillations.We now consider again the case of a piecewise uniform mesh with twodistinct mesh sizes, as in Figure 2.3. Suppose that β2 < 1, meaning thatno oscillations are present inside the boundary layer. We would hope thatoscillations outside of the boundary layer would be reduced even if β1 > 1.92.4. Discrete Solution for Non-Uniform MeshFigure 2.4: Discrete solution for uniform and non-uniform mesh of (2.1)with  = .01, N = 22, σ = 1.Farrell et. al. [6] investigate a similar situation using a piecewise-uniformmesh with N/2 points in each of the fine and coarse regions, though withoutthe restriction that the mesh must be a bi-tree. They find that non-physicaloscillations outside of the boundary layer are bounded, and greatly reducedin comparison with the unbounded oscillations of the uniform mesh. Inparticular, if the solution inside of the boundary layer is oscillation free,numerical oscillations in the coarse region will decrease in magnitude furtherfrom the fine mesh. Though it is impossible to remove these oscillationswhen β1 > 1, they may be small to the point of invisibility. We observea similar effect, as shown in Figure 2.4. Here, for the non-uniform coarsemesh β1 ≈ 4 and for the uniform mesh β ≈ 2. However, although the Pe´cletnumber is larger for the coarse mesh, the oscillations are not visible as theyare for the uniform mesh. This is due to the oscillation free solution in theboundary layer on the non-uniform mesh.We consider the effect on the coefficients of varying the location of thetransition point t as shown in Figure 2.5. Note that t is chosen such that t ∈[0.5, 1]. The fine and coarse mesh sizes are kept constant, so as the transition102.4. Discrete Solution for Non-Uniform Meshpoint moves right, less of the overall mesh is finely resolved. It follows thata transition point of t = 1 corresponds to a uniform coarse mesh. Thecoefficients a and e generally decrease as more and more of the mesh is finelyresolved. In particular, e will always decrease, since as the transition pointt −→ 0.5, the two terms in the denominator, α and−ξM2++ξN2+, both increase.The coefficient a may increase or decrease. As t −→ 1, −ξM2 + ξN2 −→ 0 anda −→ αα(ξM1+−1). Whether a increases or decreases depends on whether thenumerator α or the denominator α(ξM1+ − 1) − ξM2+ + ξN2+ is approaching 0more quickly. Regardless of its change in size, however, a will be very smallfor the choices of N and  we investigate. For some choices of M , N , and, the coefficients a and e are below machine precision, visually removingthe oscillatory part of the solution. In summary, the small values of thecoefficients a and e such as those shown in Figure 2.5 will control the size ofnon-physical oscillations arising from large powers of ξ+ < 0 in the solution.Figure 2.5: Size of coefficients a, e in discrete solution for  = .03, N = 10,β1 =53 , β2 =512 .11Chapter 3On the Eigenvalues of aFinite Difference Operatoron a Non-Uniform Mesh3.1 Matrix TransformationConsider the (N − 1) × (N − 1) operator A corresponding to the centreddifference approximation on a uniform mesh in (2.2). A is a tridiagonalToepliz matrix so its eigenvalues are known analytically. A general tridiag-onal Toeplitz matrix with entries [b, a, c] will have eigenvalues as shown in[4]:λj = a+ sign(c)2√bc cos(pijN), j = 1, ...N − 1. (3.1)The matrix A defined on the uniform mesh will have eigenvalues2− sign(β + 1)2√(β − 1)(−1− β) cos(pijN), j = 1, ...N − 1. (3.2)We show the eigenvalues calculated experimentally in Matlab for a uni-form and non-uniform mesh in Figure 3.1. We also show the eigenvalues ofthe two Toeplitz principal sub-matrices of the non-uniform mesh.We now turn our attention to investigating the analytical form for theeigenvalues of the operator for a non-uniform mesh. With the loss of aconstant mesh size, the Toeplitz structure of the matrix disappears. Anoperator A corresponding to a non-uniform mesh has the structure123.1. Matrix TransformationA =a1 c1b1 a1 c1. . .. . .. . .b1 a1 c1b2 a2 c2b3 a3 c3. . .. . .. . .b3 a3 c3b3 a3.The entries b2 and c2 limit this from being a block diagonal matrix withToeplitz submatrices, for which a closed form solution for the eigenvaluesexists. We attempt to regain a closed form solution.This matrix can be transformed by the orthogonal similarity transfor-mation shown in [4] into a symmetric matrix.A˜ =α1 γ1,1γ1,1 α1 γ1,1. . .. . .. . .γ1,1 α1 γ2,1γ2,1 α2 γ3,2γ3,2 α3 γ3,3. . .. . .. . .γ3,3 α3 γ3,3γ3,3 α3.where γi,j =√bicj and αi = ai.A˜ has formA˜ = A˜1 γ2,1 0γ2,1 α2 γ3,20 γ3,2 A˜3 , (3.3)where A˜1, A˜3 are symmetric tridiagonal Toeplitz matrices.We follow the algorithm shown in [9] for approximating the eigenvaluesof A˜. Since the matrices are related by a similarity transformation, thisequates to approximating the eigenvalues of A.133.2. Finding the Secular Equation3.2 Finding the Secular EquationThe algorithm shown in [9] involves a recursive division of a symmetrictridiagonal matrix into submatrices. However, since each A˜i is a tridiagonalToeplitz matrix, we can find its spectral decomposition as in [4] and onlyrequire one such recursive step, the proceedings of which are detailed below.Suppose that QiDiQTi is such a decomposition of A˜i. We can thentransform A˜ through an orthogonal similarity transformation asA˜ =Q1D1QT1 γ2,1ek 0γ2,1eTk α2 γ3,2eT10 γ3,2e1 Q3D3QT3 ,=0 Q1 01 0 00 0 Q2 α2 γ2,1lT1 γ3,2fT2γ2,1l1 D1 0γ3,2f2 0 D3 0 1 0QT1 0 00 0 QT3 ,= QHQT ,where lT1 is the last row of Q1 and fT2 is the first row of Q2. Since theyare related by an orthogonal similarity transformation, H has the sameeigenvalues as A˜.H is a symmetric arrowhead matrix (SAM). A closed form solution forirreducible symmetric arrowhead matrices exists [17]; that is, a SAM whereall elements of the first column z = [γ2,1l1, γ3,2f2]T are non-zero, and thediagonal elements di obey d2 >, . . . , > dN−1.A symmetric arrowhead matrix may also be transformed into an irre-ducible symmetric arrowhead matrix. Suppose an entry zi = 0. Since thecolumn corresponding to di is a scalar multiple of the unit vector ei, (di, ei)is an eigenvalue eigenvector pair for H. We can therefore remove row andcolumn i of H and solve the reduced problem to find the remaining eigen-values of H.Furthermore, we know that the diagonal entries di are the distinct eigen-values of A˜1, A˜3 given by (3.1). If the entries are not ordered d2 >, . . . , >dN−1, the matrix may be symmetrically permuted to order them withoutaffecting the spectral decomposition.By the Cauchy interlacing theorem, [22], we know that the sorted eigen-values of H, λi interlace the sorted di:λ1 ≥ d2 ≥ λ2 ≥ d3 ≥ · · · ≥ λN−2 ≥ dN−1 ≥ λN−1. (3.4)143.2. Finding the Secular EquationH is now a irreducible symmetric arrowhead matrix, and we can isolatethe remaining eigenvalues of H as the roots of the secular equation [17]f(λ) = λ− α2 +N−1∑j=2z2jdj − λ = 0, (3.5)which may be solved by a root-finding method such as bisection.3.2.1 Approximating the Roots of the Secular EquationEquation (3.5) above gives us the roots of the secular equation. When N > 5the exact roots are impossible to find. We therefore consider approximateroot finding methods to find the eigenvalues λk.We use a modification of Newton’s method in the fashion of Demmel [3].We approximate f(λ) by a simple function with easily computable roots.We call our approximate function h(λ).Since f(λ) has poles at di and di+1 for i = 2, ..., N − 2, when seeking theroot in (di, di+1) we choose to have h(λ) have these poles as well:h(λ) =c1di − λ +c2di+1 − λ + λ+ c3.If c1, c2 and c3 are known, we can solve h(λ) = 0 for λ by solving thecubic equationc1(di+1 − λ) + c2(di − λ) + (λ+ c3)(di − λ)(di+1 − λ) = 0.We now consider how to compute c1, c2, and c3 so that h(λ) ≈ f(λ) nearλ = λk:c1di − λ +c2di+1 − λ + c3 = h(λ) ≈ f(λ) = λ − α2 +N−1∑j=2z2jdj − λ.We considerf(λ) = λ− α2 +i∑j=2z2jdj − λ +N−1∑j=i+1z2jdj − λ,= λ− α2 + φ1(λ) + φ2(λ).153.2. Finding the Secular Equation(a) N = 20,  = .1, σ = 1.(b) N = 120,  = .01, σ = 1.Figure 3.1: Eigenvalues of finite difference operator for (2.1) with uniformand non-uniform mesh. A1 and A3 are the submatrices of A defined in (3.3).163.2. Finding the Secular EquationFor λ ∈ (di, di+1), φ1(λ) is a sum of positive terms, and φ2(λ) a sum ofnegative terms. Each of φ1(λ), φ2(λ) can therefore be computed withoutfear of cancellation error. We choose c1, cˆ1 to satisfy h1(λk) = φ1(λk),h′1(λk) = φ′1(λk), whereh1(λ) = cˆ1 +c1di − λ.These conditions are the same as in Newton’s method, except that wespecify a hyperbolic approximation to the graph of φ1(λ) at λ = λk insteadof a tangent line approximation.This specification gives us c1 = φ′1(λk)(di − λk)2 and cˆ1 = φ1(λk) −φ′1(λk)(di − λk).We similarly choose c2, cˆ2 to satisfy h2(λk) = φ2(λk), h′2(λk) = φ′2(λk),whereh2(λ) = cˆ2 +c2di − λ.This gives us c2 = φ′2(λk)(di+1−λk)2 and cˆ2 = φ2(λk)−φ′1(λk)(di+1−λk).We then seth(λ) = λ− α2 + h1(λ) + h2(λ),= λ− (cˆ1 + cˆ2 − α2) + c1di − λ +c2di+1 − λ,= λ+ c3 +c1di − λ +c2di+1 − λ,Example 3.2.1. We consider using h(λ) to approximate f(λ) in as in Figure3.2. We are solving for λ4 inside the range (d3, d4) = (0.4708, 0.6678), andwe observe that the approximation, shown as a solid blue line, is very goodbetween these two poles, shown as vertical black lines. It is between thesetwo poles that λ4 lies, and so our approximation need only be be accuratein this region. We use an initial approximation of λ¯0 =d3+d42 , to constructh(λ) = −0.7499 + λ + 0.16590.4708−λ + 0.03710.6678−λ = 0, which we solve to obtainλ¯1 = 0.6348. This is within 10−2 of the exact value of λ4, 0.6438181494679.At the next iteration, we use λ¯1 to construct h(λ) = −0.4349+λ+ 0.16900.4708−λ+0.01860.6678−λ = 0, and solve the equation for a value of λ¯2 = 0.6436, which haserror |λ¯2 − λ4| < 10−3. A third iteration has error |λ¯3 − λ4| < 10−6.We now have both a secular equation for the eigenvalues of a non-uniformmesh as well as an efficient method for determining the roots of that equa-173.2. Finding the Secular EquationFigure 3.2: Graph of f(λ) and h(λ) for N = 20,  = 0.1, σ = 1 and k = 4.tion. This is a promising step towards approximating the eigenvalues of anon-uniform mesh in a higher dimension.18Chapter 4Finite Differences on aNon-Square Domain4.1 Convection-Diffusion on Non-Square DomainConsider the convection-diffusion equation in two dimensions on some do-main Ω with Dirichlet boundary conditions:−uxx − uyy + σux + τuy = f u ∈ Ω,u = g u ∈ δΩ.(4.1)The problem, for small values of σ and τ , is well posed. The simplest caseto consider is when Ω is a polyomino [8], that is, made up of a combinationof orthogonal squares. This property implies that (4.1) can be solved numer-ically using the second order centered finite difference (FD) method illus-trated below with a uniform grid. A grid with mesh size h = 1n+1 and origin(x0, y0) is made up of the uniformly spaced points (xi, yj) = (x0+ih, y0+jh)for 1 ≤ i, j,≤ n. We will define Ui,j as our numerical approximation tou(xi, yj). Similar notation may be used to refer to other functions definedon the grid, such as f or u. A representation of a subset of the grid is shownin Figure 4.1. We also assume that our grid points are labelled lexicograph-ically, which for Figure 4.1 would result in the labelling shown in Figure4.2.Given a uniform grid with mesh size h and points as shown in Figure4.1, we will use the following centred finite difference discretizations in ourfinite difference scheme.D0xUi,j =Ui+1,j − Ui−1,j2h, (4.2)D0xxUi,j =Ui+1,j − 2Ui,j + Ui−1,jh2, (4.3)with analogous definitions of D0y, D0yy. These definitions give us the following194.2. Motivating Problemsfinite difference scheme to approximate (4.1) at the point (xi, yj) ∈ Ω.−D0xxUi,j −D0yyUi,j + σD0xUi,j + τD0yUi,j ≈ fi,j . (4.4)We use Dirichlet boundary conditions, meaning thatUi,j = gi,j (xi, yj) ∈ δΩ. (4.5)The linear system resulting from this scheme is of the form Ax = b,where A is a block tridiagonal Toeplitz matrix. If the problem is reduced tothe Poisson problem (σ, τ = 0), then A is also a symmetric positive definitematrix. These properties allow for the linear system to be solved efficiently,and also allow further analysis of the matrix A. For instance, the eigenvaluesof A are known analytically.If either of σ or τ increase such that the Pe´clet number, max{σ,τ}h2 isgreater than 1, then large non-physical oscillations may arise in the numeri-cal solution, as in the 1D case examined previously. Furthermore the matrixA is no longer an M-matrix and the system no longer has real eigenvalues. Inthis situation, a scheme such as upwinding [16] is advised to restore stabilityto the numerical solution. We consider situations where max{σ,τ}h2 < 1. Ourfocus will be on overcoming the difficulties in using finite difference methodsto solve (4.1) with a non-square domain.4.2 Motivating ProblemsThere are several problems which we have chosen to motivate the study ofthe non-uniform mesh in 2D. A domain that is not a polyomino may requirea refinement strategy other than a non-uniform grid. We have chosen twodistinct non-square problems to highlight the different techniques involvedin using FD methods on a non-square domain.The first case is that of a square grid with an ellipse removed. A practicalexample where this may occur is the case of a square plate with a round holein it, as discussed by Morton and Meyers [16]. In the initial situation wherethis domain is discretized with a square uniform grid, there are variousdifficulties. Mathematically, the grid points may not fall on the curvedboundary. Computationally, the labelling of grid points is more complicated.Another situation that we will consider is the numerical solution of thePoisson problem on an L-shaped domain. Here, the singular nature of thesolution at the origin forces a highly refined mesh in the neighbourhood ofthe origin. This non-uniform mesh presents difficulties when finite differ-ences are used.204.3. Computational SchemesFigure 4.1: FD solution values on a uniform mesh.4.3 Computational SchemesWe will first discuss the various computational difficulties associated withusing a FD scheme on a non-uniform mesh.In order to construct the finite difference operator, coefficients must beassigned for the finite difference scheme at each grid point. This requireshaving a list of the neighbours for each grid point.For a uniform grid, orthogonal points are required for a centred FDscheme. Consider the lth grid point at the location (xi, yj) in a uniformmesh with total number of grid points N = n2. If the grid is labelledlexicographically, the orthogonal neighbours required for the centred FDscheme given in (4.4) will be located at the grid points labelled l ± n and214.3. Computational SchemesFigure 4.2: Lexicographic labelling of a uniform mesh.l ± 1.When the grid is non-uniform and non square, the adjacent grid pointsare no longer labelled ±n or ±1. An example of this can be seen in Figure4.3.Since we no longer have consistent spacing of labels, we must use analternate scheme to find the neighbours of each grid point. We find the rowand column index, i, j of each grid point, then add or subtract 1 from theseindices to find the respective neighbour. For example, to find the northernneighbour of the point 15 in Figure 4.3, we find its row and column indicesof (3, 3) , as shown in Figure 4.4. We then find the lexicographic label (14)of the grid point with row and column indices of (2, 3). The lexicographiclabels are stored in an array to allow the described indexing.224.4. Domain with a Curved BoundaryHowever, the situation is complicated when the grid is refined locally.This creates extra row and column indices, and “missing” row and columnindices. Following the method above, when we try to find the western neigh-bour of the point 27 by subtracting 1 from the column index, we find thatthere is no point at (6, 4). These “missing indices” are represented by thepoints in yellow with dashed borders in Figures 4.4 and 4.3. The grid pointsaffected are coloured green in these figures. To avoid repeated and exten-sive searches, we instead add to our array additional points at the missingindices, with the lexicographic label of the next neighbour in the missingdirection. Multiple arrays may be required for the case of multiple possiblelabels, as represented by the multiple labelings for the yellow “missing in-dices” in Figure 4.3. We use the array associated with the direction we aresearching for a neighbour in order to find the correct index.Figure 4.3: Lexicographic indices for a locally refined grid. Spaces corre-sponding to a missing index are in yellow, and grid points affected by themissing indices are in green.4.4 Domain with a Curved BoundaryConsider the domain of a square with an ellipse-shaped hole at its center.The radius and eccentricity of the ellipse are defined by the factors r, s and234.4. Domain with a Curved BoundaryFigure 4.4: Row and column indices for a locally refined grid. Spaces cor-responding to a missing index are in yellow, and grid points affected by themissing indices are in green.the domain has formΩ = {(x, y) ∈ [0, 1]× [0, 1]} ∪ {(x, y)|s(x− 0.5)2 + (y − 0.5)2s> r}. (4.6)This domain is no longer rectangular, meaning that the grid points of arectangular uniform grid may not lie inside the domain. In particular, theboundary of the domain may lie between grid points. We will instead use amodified finite difference scheme, as in [16].This scheme is shown in Figure 4.6. Here, h is the mesh size of the grid.Due to the non-square domain, additional grid points have been definedon the boundary of the domain, at (xA, yj) and (xi, yB). These points areshown as red circles in Figure 4.5.We must use a modified discretization in order to construct an approx-imations of uxx, uyy, ux and uy at (xi, yj). We will begin by discussingFD schemes for for the partial derivatives with respect to x. We find|xA − xi| = αh for some α. Then we find the following Taylor expansionsfor u(xA, yj) and u(xi−1, yj):u(xA, yj) = u(xi, yj) + αhux(xi, yj) +12(αh)2uxx(xi, yj) +O(h3),u(xi−1, yj) = u(xi, yj)− hux(xi, yj) + 12h2uxx(xi, yj)−O(h3).244.4. Domain with a Curved BoundaryFigure 4.5: Grid points on (4.6) with h = 0.0625, N = 176. Grid pointslying on a uniform mesh are shown in blue, and additional grid points onthe boundary of the domain shown as red circles.From these we obtainux(xi, yj) =u(xA, yj)− α2u(xi−1, yj)− (1− α2)u(xi, yj)α(1 + α)h+O(h),uxx(xi, yj) =u(xA, yj) + αu(xi−1, yj)− (1 + α)u(xi, yj)12α(1 + α)h2+O(h).From this we get the finite difference approximations of ux, uxx (respec-tively):DαxUi,j =UA,j − α2Ui−1,j − (1− α2)Ui,jα(1 + α)h, (4.7)DαxxUi,j =UA,j + αUi−1,j − (1 + α)Ui,j12α(1 + α)h2. (4.8)Since the value of u(xA, yj) = g(xA, yj) is known from our boundaryconditions, this is a valid scheme. We can obtain in a similar way the finite254.4. Domain with a Curved BoundaryFigure 4.6: Modified FD scheme for a curved boundary with Dirichletboundary conditions.difference approximations for uyy and uy:DβyUi,j =Ui,B − β2Ui,j−1 − (1− β2)Ui,jβ(1 + β)h, (4.9)DβyyUi,jUi,B + βUi,j−1 − (1 + β)Ui,j12β(1 + β)h2. (4.10)Now that we have finite difference equations for all grid points, we canconstruct and solve the full linear system Ax = b for the discretized problem.We do so for the problem with exact solution:264.5. Reentrant CornerTable 4.1: Convergence of FD scheme for (4.1) on (4.6) σ = 2, τ = 0,r = 0.25, s = 1. E2 = ||U − uexact||2 and E∞ = ||U − uexact||∞.N h ||U − uexact||2 ||U − uexact||∞ E2h2 /Eh2 E2h∞ /Eh∞176 6.250× 10−2 3.999969× 10−5 9.326048× 10−5 4.52 4.25764 3.125× 10−2 8.849087× 10−6 2.193087× 10−5 4.19 4.063172 1.563× 10−2 2.109649× 10−6 5.402625× 10−6 4.12 4.0512920 7.813× 10−3 5.115872× 10−7 1.332949× 10−6 4.06 4.0352172 3.906× 10−3 1.258631× 10−7 3.311233× 10−7 4.12 4.06209688 1.953× 10−3 3.055653× 10−8 8.164156× 10−8 4.80 4.39840668 9.766× 10−4 6.371100× 10−9 1.858327× 10−8 − −uexact(x, y) = exyy(1− y)x(1− x)(s(x− 0.5)2 + 1s(y − 0.5)2 − r2) + 1. (4.11)By definition, this function will have boundary value g = 1.It can easily be seen that A will lose the symmetry it had on a uniformmesh in a square domain, even when σ and τ are zero. Even if the i, jthentry of A, ai,j is nonzero, it is possible that aj,i = 0. This lack of structuralsymmetry means that this matrix is not symmetrizable.It may also be observed that the FD equation for the points on theinterior boundary is only O(h), since it results from a combination of (4.7),(4.8), (4.9), (4.10). However, this does not appear to affect the global secondorder accuracy, as seen in Table 4.1. As discussed in [16, p. 172], this gainin accuracy at the boundary is to be expected. By performing error analysisusing the maximum principle, it can be shown that the second order accuracyof the finite difference scheme for the interior mesh will be retained even witha first order scheme at the boundary.4.5 Reentrant CornerWe now consider a second example of where a non-uniform mesh may berequired, as introduced in Section 4.2. Consider the L-shaped domain:Ω = [−0.5, 0.5]× [−0.5, 0.5]\\{(0, 0.5)× (−0.5, 0)}. (4.12)The 90◦ angle formed by this domain at the origin is called a “reentrantcorner”.274.5. Reentrant CornerFigure 4.7: FD solution of (4.1) on the domain (4.6) with σ = 2, τ = 0,r = 0.25, s = 1.We consider solving Poisson’s equation on Ω with Dirichlet boundaryconditions:−uxx − uyy = f u ∈ Ω,u = g u ∈ δΩ.(4.13)In the case of the homogeneous boundary conditions g = 0 and a smoothright hand side f , (5.2) has solution: [20]uexact(x, y) = (√x2 + y2)1/3 sin(23)arctan(yx). (4.14)The solution and its first derivative are in L2(Ω), but since the secondderivative near the origin is unbounded, the second derivative is singularand an accurate discrete solution will prove difficult to find numerically.An initial plan is to discretize Ω with a uniform grid. Since Ω is apolyomino a finite difference method with a uniform discretization will leadto a valid discrete solution.284.5. Reentrant CornerFigure 4.8: Error in FD solution of (4.1) on the domain (4.6) with σ = 2,τ = 0, r = 0.25, s = 1.Specifically, we can discretize (5.2) on a uniform grid as in Figure 4.1.We can then use the second order centered finite difference scheme:−D0xxUi,j −D0yyUi,j ≈ fi,j . (4.15)A plot of the approximate solution U resulting from this scheme is shownin Figure 4.9.We also observe the error in U for a uniform grid in Figure 4.10. Notein particular that the error at the origin is very large.In fact, at a fixed distance from the singularity, we expect convergenceof O(h4/3), and near the singularity we expect the even worse convergencebehaviour of O(h2/3) [1, 13]. Here, the singularity in the reentrant cornerleads to an increase of the discretization error in all of Ω. This effect is knownas pollution. It follows, especially considering Figure 4.10, that refinementof the mesh is not needed everywhere to improve the accuracy of solution;but concentrated refinement is needed at the reentrant corner.294.5. Reentrant CornerFigure 4.9: FD solution of (5.2) on L-shaped domain.In fact, an “optimal” general formula for returning to second order ac-curacy is given in [1] which references [12]. The scheme, which we follow, isat each refinement step to refine each mesh size by a factor of 2 and thenadd two additional layers of refinement near the singularity. Each refinedregion is also an L-shape and has half the side length of the next largestlocally refined segment. For example, given the mesh in Figure 4.12, thenext refinement step would produce the mesh in Figure 4.13. The resultingmesh will have mesh sizes hmax, hmax/2, hmax/4, ..., hmax/k, where k is thelargest factor of refinement. k will grow by a factor of 22 each time themesh is refined. This scheme follows from an understanding that, near thesingularity, we see convergence of O(h2/3). Therefore, in order to see theerror decrease by a factor of 4 when we halve h, that is(12)2= 14 , we mustrefine the grid by a factor of(12) 22/3 = 18 at the origin. By adding two lay-ers of refinement in addition to halving h, we maintain a graded grid whileachieving the required refinement near the singularity. This restores secondorder accuracy, as shown in Table 5.3.As can be seen in Figure 4.12, some mesh points do not have the four or-304.6. Alternative DiscretizationsFigure 4.10: Error in FD solution for (5.2) on L-shaped domain with auniform grid, h = 0.0156, N = 2945.thogonal neighbours required for the centred finite difference scheme (4.16).Even if a grid point does have all four neighbors, they may not be uni-formly distanced from the central point. We will now discuss the necessaryalterations to the finite difference scheme to accommodate these difficulties.4.6 Alternative DiscretizationsThe first situation we must deal with is when not all mesh sizes are equal;i.e., when some neighbours of Ui,j are closer than others. A visual exampleis given in Figure 4.14. Note that h1 is half the size of h2.With this discretization, we must use alternative finite difference schemesto approximate uxx, uyy. The following is a second order centered finitedifference approximation of uxx:bUi−1,j + aUi,j + cUi+1,jh2, (4.16)314.6. Alternative DiscretizationsTable 4.2: Convergence of FD scheme for (5.2) on (4.12). E2 = ||U−uexact||2and E∞ = ||U − uexact||∞.N hmax k ||U − uexact||2 ||U − uexact||∞ E2h2 /Eh2 E2h∞ /Eh∞205 6.250× 10−2 16 2.469249× 10−4 1.690082× 10−3 5.50 3.90945 3.125× 10−2 64 4.488844× 10−5 4.337456× 10−4 4.44 3.984161 1.563× 10−2 256 1.010690× 10−5 1.088949× 10−4 3.65 3.9916583 7.813× 10−3 1024 2.771935× 10−6 2.726217× 10−5 3.73 4.0067205 3.906× 10−3 4096 7.428256× 10−7 6.820151× 10−6 3.99 4.00281489 1.953× 10−3 16384 1.862977× 10−7 1.705491× 10−6 - -Figure 4.11: FD solution of Poisson problem on L-shaped mesh with Dirich-let boundary conditions.wherea = − 2h1(h1 + h2)− 2h2(h1 + h2),b =2h1(h1 + h2),c =2h2(h1 + h2).324.6. Alternative DiscretizationsFigure 4.12: Locally refined grid for (4.12), hmax = 0.0625, k = 64.Given a situation where one or more neighbours is missing, we must usean alternate finite difference scheme to approximate uxx or uyy.Suppose that we are searching for an alternative discretization for uyywhere there is no northern neighbour, Ui,j+1, in the mesh. An example ofthis situation is given in Figure 4.15.Given the missing grid point, we wish to develop a scheme of the formDunionsqyyUi,j = aUi,j + bUi+1,j+1 + cUi−1,j+1 + dUi−1,j + eUi+1,j + fUi,j−1.A Taylor expansion of Dunionsqyyu(xj , yj) isDunionsqyyu ≈au+ b(u+ h2ux + h3uy +h222uxx + h2h3uxy +h232uyy)+ c(u+ h1ux − h3uy + h212uxx − h1h3uxy + h232uyy)+ d(u− h1ux + h212uxx) + e(u+ h2ux +h222uxx)+ f(u− h4uy + h242uyy)334.6. Alternative DiscretizationsFigure 4.13: Locally refined grid for (4.12), hmax = 0.03125, k = 256.We solve the system of equations so that Dunionsqyyu = uyy + O(h2). Thismeans that we must find coefficients to satisfy(a+ b+ c+ d+ e+ f)u = 0;(h2b− h1c)ux = 0;(h3b+ h3c− h4f)uy = 0;(h22b+ h21c+ h21d+ h22e)uxx = 0;(h23b+ h23c+ h24f)uyy = uyy;(h2h3b+ h1h3c)uxy = 0.We have six equations and six unknowns. We solve the resulting system344.6. Alternative DiscretizationsFigure 4.14: FD solution values on non-uniform mesh, non-uniform meshsize.of linear equations to obtaina =−2h3(h3 + h4)− 2h4(h3 + h4)+2h3(h3 + h4);b =2h2h3(h3 + h4)(h1 + h2);c =2h1h3(h3 + h4)(h1 + h2);d =−2h3(h3 + h4)h2h1h1(h2 + h1);e =−2h3(h3 + h4)h2h1h2(h2 + h1);f =2h4(h3 + h4).354.6. Alternative DiscretizationsFigure 4.15: FD solution values on non-uniform mesh, grid point missingnorthern neighbour.We can obtain a similar operator for a point missing a southern neigh-bour, or to approximate uxx when an eastern or western neighbour is miss-ing.We also wish to expand our methods to solve (4.1). For this, we requiresimilar FD operators to approximate ux, uy. We can use a similar methodto derive these. For example, we obtain the following approximation for uy:Dunionsqy Ui,j = aUi,j + bUi+1,j+1 + cUi−1,j+1 + dUi−1,j + eUi+1,j + fUi,j−1.(4.17)364.6. Alternative DiscretizationsFigure 4.16: FD operator for (5.2) on L-shaped domain with locally refinedgrid, N = 205. Grid points are labelled lexicographically.wherea =h3h24 + h3h4;b =h1h4h3(h1 + h2)(h3 + h4);c =h2h4h3(h1 + h2)(h3 + h4);d = − h2h4h3(h1 + h2)(h3 + h4);e = − h1h4h3(h1 + h2)(h3 + h4);f = − h3h24 + h3h4.374.6. Alternative DiscretizationsThe resulting linear operator, shown in Figure 4.16, has a distinctivesparsity pattern, due to the necessary modified FD schemes. The form ofthis operator is also dependant on the ordering scheme used; as discussedpreviously, lexicographic ordering is used here. Again, the lack of symmetryin the location of non-zeros makes this matrix unsymmetrizable.38Chapter 5Cell Centered FiniteDifferences5.1 A Symmetric SchemeThe finite difference discretization scheme in Chapter 4 resulted in an un-symmetrizable operator. There are many reasons to pursue a scheme result-ing in a symmetric operator including the fast solvers, such as the conjugategradient method [18], only available for symmetric positive definite matricesand the relative ease of analysis for symmetric matrices. Furthermore weare interested in extending the eigenvalue analysis of Chapter 3 to higherdimensions, and a symmetric matrix is at least helpful, if not essential forthis.This pursuit of a symmetric scheme for differentiating elliptic equationshas led us to the finite difference (FD) scheme based on a finite volume(FV) writing of the Poisson equation as in Ewing et. al. [5] and Haber andHeldman [10] which results in a symmetric operator. We will now discussthis scheme and its use for the Poisson problem on a quadtree grid in anon-regular domain.In order to use the approach of [5, 10], we use a cell centred grid, asshown in Figure 5.1, as opposed to the vertex centred grids used previously.This is an important component in considering the differential equation froma FV framework, since a FV scheme involves integrating over a cell’s volumeto find the discrete solution at the center.We also continue to use a quadtree grid, with all mesh sizes restrictedto 2kh0, where h0 is the finest mesh size. We stipulate that the mesh mustbe graded ; that is, the mesh size of a cell must be at most a factor of two ofthe mesh size of all neighbouring cells.Furthermore, we are interested in non-regular domains. In order to con-tinue in a true quadtree style while solving on non square domains, we usea full, square grid with an immersed boundary method as in Rycroft et. al.[21].395.2. Data StructureFigure 5.1: Coordinates of grid points in cell centred grid.5.2 Data StructureA critical component of a quadtree is its capacity for adaptive refinement.For this reason, we require a data structure to store information about thegrid that is cheap to store and to update.We adopt the strategy of [10] of storing grid information in a large,sparse array, which we will call a Depth Map. This array will be of sizeN × 2ML−1, where ML is specified as the maximum number times the gridcan be refined. Each cell in the grid is represented by a value in the upperleft hand corner of its corresponding location in the Depth Map, as shownin Figure 5.2. This value represents the refinement level of the cell, with 1being the finest possible and ML being the coarsest possible. The smallestcells are therefore represented by one entry in the Depth Map array, and thelargest by a sparse square of size 2ML−1 × 2ML−1.When a cell in the grid is refined, its value in the Depth Map is reducedby 1, and new entries in the three new child cells are introduced, as shownin Figure 5.2b.405.3. Domain Definition and Refinement Strategies(a) Before refinement. (b) After refinement.Figure 5.2: Entries in the Depth Map array before and after the cell in lefthand corner is refined. Cells undergoing refinement are in yellow.5.3 Domain Definition and Refinement StrategiesWe continue to be interested in solving problems on non-square domains.Maintaining a quadtree style grid requires a square grid by definition inorder for each parent cell to be divided into four child cells. We maintainthis convention in order to use a quadtree style grid. In order to solveproblems on non-square domains, we will use a modified version of the levelset method shown in [21].These methods allows us to describe our domain in terms of a level-set function, φ(x, y). A true level set method requires this function to beLipschitz continuous in order to characterize the change in the location of theboundary of the domain through time. Since our boundary will not move,we can remove the restriction of Lipschitz continuity from φ. We define ourdomain such that a square grid, say Ω = [−1, 1]× [−1, 1], is divided into twosubsets by the sign of φ(x, y).The domain of our problem is then specified by Ω+, defined whereφ(x, y) > 0. The boundary of the domain, which we will call δΩ+, is definedwhere φ(x, y) = 0. We define all areas where φ(x, y) <= 0 as Ω−. Note thatthis includes the boundary of the domain, δΩ+.As in Rycroft [21] we refine in a uniform band along the boundary,refining a cell c if:minv∈vertices(c)|φ(v)| ≤√2hc, (5.1)415.4. The Discretization Scheme(a) Before refinement. (b) After one refinement.(c) After two refinements. (d) After three refinements.Figure 5.3: Quadtree refinement according to (5.1).where hc is the mesh size of cell c.The implementation in [21] also includes a term related to the Lipschitzconstant of φ; this means that a domain that has more curves will be refinedin a larger uniform band. However, some of the domains investigated arenot Lipschitz continuous and we therefore ignore this term.We also refine cells as necessary to maintain a graded grid. The refine-ment process is shown in Figure 5.3.5.4 The Discretization SchemeConsider again the Poisson equation, as shown in previous chapters, for ageneric domain Ω:425.4. The Discretization Scheme−uxx − uyy = f u ∈ Ω,u = g u ∈ δΩ.(5.2)Figure 5.4: Grid point and side labelling of quadtree grid.We can rewrite this equation as the divergence of the gradient of u:∇ · ∇u = f, u ∈ Ω (5.3)u = g u ∈ δΩ (5.4)We now turn our attention to the discretized problem. Our grid is madeup of square cells but is not necessarily uniform; these cells may be of dif-ferent sizes. Recall, however, that the grid is graded; a cell’s size must beat most a factor of 2 different from any adjacent neighbour.As in [5], we define wx+i,j , wx−i,j as approximations of the flux across thesides of the cell sx+i,j and sx−i,j respectively. Similarly, we define wy+i,j , wy−i,j asapproximations of the flux across sy+i,j and sy−i,j .435.4. The Discretization SchemeFigure 5.5: Flux across sides of quadtree grid.Consider a cell, ci,j , with volume vi,j , as shown in Figure 5.4. By a massbalance approach, we get the divergence equation, as in [10, 11],1vi,j(wx+i,j − wx−i,j + wy+i,j − wy−i,j ) =1vi,j∫∫ci,jf(x, y)dxdy (5.5)We need to define the fluxes wi,j in terms of Ui,j , the discrete solutionvalue at the cell center. We define the gradient operator on the cell faces.For regular cells, the simple difference:wx−i,j =1hi,j(Ui,j − Ui−1,j) (5.6)is a second order approximation. A similar definition exists for the threeother cell faces.At non-regular points, many different possible discretizations exist. Asstated previously, we wish to use a symmetric discretization. As proposed445.4. The Discretization Schemein [5] and implemented in [10, 11], we use the same two point differenceas in the regular cell case, (5.6). This scheme has the advantage of beingsymmetric. More accurate schemes exist, but they are more expensive withrespect to the number of neighbouring points required, and may not besymmetric.Mass conservation gives us thatwx+i,j = wx−i+1,j+1 + wx−i+1,j−1,and therefore, we have the resulting equality, which we substitute in (5.5)wx+i,j =1hi+1,j+1(Ui+1,j+l − Ui,j) + 1hi+1,j−1(Ui+1,j−1 − Ui,j),wx+i,j =2hi,j(Ui+1,j+l − Ui,j) + 2hi,j(Ui+1,j−1 − Ui,j).A similar result follows for irregular grid points in the y direction.Note that this irregular gradient discretization scheme has the lowereddiscretization error of O(h). This result is shown in Horesh and Haber [11].The scaled transpose of the divergence (Div) operator arising from (5.5)is the gradient (Grad) operator, from (5.6). When multiplied together, thesetwo operators are a discrete Laplacian operator for (5.2). We therefore willrefer to this discretization scheme as the Div-Grad scheme.Observe also that we approximate 1vi,j∫∫ci,jf(x, y)dxdy with a midpointquadrature rule; this equates to approximating the integral by multiplyingthe value of f at the cell center, (xi, yj), by the volume of the cell vi,j , andscaling by the factor of 1vi,j . However, in order to construct a symmetricoperator, we multiply both sides of each linear equation by a factor of vi,j .In constructing the system of linear equations, our right hand side valuesare therefore multiplied by a diagonal matrix with entries corresponding tothe volume of each cell.5.4.1 Boundary Conditions and Operator ModificationsTo allow for the non-polyomino [8] domains defined below, we use a level-set function, as discussed in Section 5.3 to separate Ω+, the true domainof our problem, from its enclosing square, Ω = [−1, 1] × [−1, 1]. We defineUi,j = u(xi, yj), the exact solution, for cells ci,j on the boundary of thedomain. Recall that the boundary is defined where φ(x, y) = 0 at some455.5. Numerical Resultspoint inside ci,j . Note that the boundary of the domain might not passexactly through the center of a cell; we extend the value of u(xi, yj) at theboundary over the cell.Cells in the domain where φ(x, y) < 0 are set to have a solution valueof 0. Then, since all points in Ω− have a defined value, we can removeequations corresponding to these points from the system of linear equations.We modify the rows corresponding to points inside the operator that containboundary terms such that the matrix entry is replaced by subtracting theboundary value at that point from the right hand side of the equation.Therefore, our operator contains equations corresponding only to pointsinside Ω+.To summarize, we have found a symmetric scheme that allows for adap-tive domain refinement. This scheme may be used for non-polyomino do-mains and for non-uniform (graded) quadtree meshes. However, it should benoted that the price we have paid for symmetry is a less accurate operator;O(h) instead of the O(h2) accurate scheme considered in Chapter 4.5.5 Numerical ResultsNumerical results for this section were generated using the Julia Language[2]. Much of the code used was adapted from a quadtree implementation inMatlab by Eldad Haber [10].We now look at the results of solving the Poisson equation with theDiv-Grad discretization introduced on one of the domains discussed below.We consider solving (5.2) on the following domain shapes:• A cross shape,Ω+ = {(x, y) ∈ [−0.5, 0.5]× [−1, 1]} (5.7)∪{(x, y) ∈ [−1, 1]× [−0.5, 0.5]}, (5.8)• a circle,Ω+ = {x2 + y2 < 0.5}, (5.9)• an L-shape,Ω+ = [−1, 1]× [−1, 1]\\{(0, 1)× (−1, 0)}, (5.10)465.5. Numerical Results(a) Circle, (5.9). (b) Star, (5.11).(c) Cross, (5.7). (d) L-shape (5.10).Figure 5.6: Domains considered in Section 5.5. Boundary points for thedomain are red circles, and point inside the domain green diamonds.• a star,Ω+ = {(x, y) ∈ [−0.5, 0.5]× [−0.5, 0.5]} (5.11)∪{(x, y)||y| < (2|x| − 1), |x| > 0.5} (5.12)∪{(x, y)||x| < (2|y| − 1), |y| > 0.5}. (5.13)These domains are depicted in Figure 5.6.We consider solving (5.2) with the following right hand side and homo-geneous Dirichlet boundary conditions:475.5. Numerical ResultsTable 5.1: Convergence of FD scheme for (5.2) with (5.14) on circle shapeddomain (5.9). E2 = ||U2h − Uh||2 and E∞ = ||U2h − Uh||∞. k is the largestfactor of refinement.N hmax k ||U2h − Uh||2 ||U2h − Uh||∞ E2h2 /Eh2 E2h∞ /Eh∞2748 0.0625 2 0.004168 0.0112 1.97 2.0111364 0.0313 2 0.002116 0.0056 2.19 2.0746108 0.0156 2 0.000968 0.0027 1.80 1.88186044 0.0078 2 0.000539 0.0014 2.11 2.04747140 0.0039 2 0.000256 0.0007 − −Table 5.2: Convergence of FD scheme for (5.2) with (5.14) on L-shapeddomain (5.10). E2 = ||U2h−Uh||2 and E∞ = ||U2h−Uh||∞. k is the largestfactor of refinement.N hmax k ||U2h − Uh||2 ||U2h − Uh||∞ E2h2 /Eh2 E2h∞ /Eh∞4448 0.1250 2 0.006323 0.0292 2.02 1.5618804 0.0625 2 0.003135 0.0187 2.01 1.5777252 0.0313 2 0.001559 0.0119 2.01 1.58313092 0.0156 2 0.000777 0.0075 2.00 1.581260548 0.0078 2 0.000388 0.0048 − −f(x, y) = 1 (x, y) ∈ Ω+ (5.14)u(x, y) = 0 (x, y) ∈ δΩ+ (5.15)We first analyze the error when Ω+ is refined uniformly. One initialrefinement is performed according to (5.1), resulting in a uniform band ofrefinement around the boundary of the domain. All following refinementsare a uniform halving of the mesh size in each cell in Ω+. Note that Ω−is not refined, except to maintain a graded grid. As stated in Section 5.4,the discretization scheme for the gradient at irregular points (5.6) has adiscretization error of O(h). Upon uniform refinement as described, weobserve this expected accuracy in Table 5.1.We investigate the accuracy of the Div-Grad scheme on an L-shapeddomain (5.10) when the mesh is uniformly refined, again after one initialrefinement along the boundary alone. Since the domain contains a reentrantcorner, we expect the accuracy to drop from O(h2) to O(h43 ) in the L2 normand O(h23 ) in the inf norm. This convergence can be observed in Table 5.2.485.5. Numerical ResultsNote that we observe the lowered O(h) convergence in the L2 norm.We observe similar reentrant corners upon uniform refinement for thecross shaped domain (5.7) as shown in Figure 5.7. For this domain, theeffect of four reentrant corners is observable.Figure 5.7: Error U2h − Uh on cross shaped domain (5.7).We expect that we can return to the first order accuracy of the Div-Grad scheme by refining near a reentrant corner, as in Chapter 4. We wishto refine in such a way that when the mesh size is halved everywhere, theerror is reduced by a factor of 2. In Chapter 4, an additional two layers ofrefinement were necessary to regain O(h2) accuracy; here one layer shouldbe sufficient to return to O(h) accuracy. Since we expect O(h23 ) accuractynear the origin, an additional layer of refinement will mean that the meshsize is reduced by a factor of 4 there. We would then expect the error to bereduced by a factor of (14)23 < 12 . This error reduction is sufficient for O(h)accuracy. This refinement scheme is shown in Figure 5.8 for the L-shapeddomain (5.10). We see that the error converges at the expected rate.We consider the conditioning of the operators for the two different dis-cretization schemes that have been used in Chapters 4 and 5 to solve (5.2) onan L-shaped mesh (5.10). In both schemes, the mesh is refined strategicallynear the origin. Recall that for the FD scheme of Chapter 4, we add two ad-ditional refinement layers near the origin to return to O(h2) accuracy, and forthe Div-Grad scheme, only one layer to return to O(h) accuracy. In Figure5.9, we observe the condition numbers of the operators in the two differentschemes. We see that κ(R) grows approximately at a rate κ(R) ∝ h−2.00190495.5. Numerical ResultsFigure 5.8: L-shaped domain refinement for Div-Grad scheme.for the Finite Difference scheme, and at a rate κ(R) ∝ h−1.17750 for theDiv-Grad scheme. We can currently only observe these values; drawing ex-tensive conclusions may be unwise as the discretization schemes they are theproduct of have different orders of accuracy.We also consider the structure of the Div-Grad operator, as shown inFigure 5.10. In comparison with the FD operator of Chapter 4, this op-erator is much more compact. In particular, we note that the long trainsof connectivity of the FD operator are not present here; instead we haveonly the occasional entry outside of a roughly banded operator. However,the operator is still far from being perfectly banded. This relative lack ofstructure would make an extension of the eigenvalue decomposition schemein Chapter 3 quite difficult. In 1D, an operator for a piecewise uniform meshremains symmetric, though it may no longer be Toeplitz as in the case ofthe uniform mesh. However, the loss of a uniform mesh results in no loss ofstructure in 1D; the operator is still tridiagonal and structurally symmetric,505.5. Numerical ResultsTable 5.3: Convergence of FD scheme for (5.2) with (5.14) on L-shapeddomain (5.10). E2 = ||U2h−Uh||2 and E∞ = ||U2h−Uh||∞. k is the largestfactor of refinement.N hmax k ||U2h − Uh||2 ||U2h − Uh||∞ E2h2 /Eh2 E2h∞ /Eh∞154 0.1250 2 0.041941 0.1066 2.34 1.89850 0.0625 4 0.017956 0.0564 2.27 2.334178 0.0313 8 0.007911 0.0242 2.23 2.4419458 0.0156 16 0.003551 0.0099 − −Figure 5.9: Condition numbers for operator R in Div-Grad and FD schemevs. finest mesh size.and therefore symmetrizable. In 2D, the loss of a uniform mesh means theresulting operator is no longer pentadiagonal, and may no longer be evenstructurally symmetric. This is the case for the operator shown in Chapter4. Structural symmetry may be regained, as shown here, but the matrix isstill not banded. As the banded nature of the operator was critical in findingthe secular equation in Chapter 3, the lack of structure in these operatorsis troubling. Approximation of the eigenvalues may still be possible, for ex-ample based on a decomposition of the operator approximated by its closest515.5. Numerical Resultsbanded counterpart, but these experiments have not yet proved fruitful.Figure 5.10: Operator for Div-Grad scheme, hmax = 0.25, k = 2, N = 208.52Chapter 6ConclusionsWe have discussed and compared various discretization schemes for the solu-tion of elliptic equations on a non-square domain. We explored two strategiesfor dealing with boundaries that did not pass through grid points. A modi-fied discretization scheme can be used, as in Chapter 4, or the value at theboundary can be extended to the grid point, as in Chapter 5. In the formercase, a suitable choice of discretization may not affect the overall accuracy ofthe scheme, even if the truncation error appears to suggest lower accuracy,as discussed in Chapter 4.We also encounter potential accuracy loss for the case of a reentrantcorner. Here, the loss in accuracy is due to the problem itself and not thediscretization scheme. We can restore the scheme to its expected accuracyby refining strategically near the reentrant corner.Alongside our pursuit of accuracy we have sought symmetry for our oper-ator. In 1D, a non-uniform mesh still results in an operator that is at leaststructurally symmetric and thus symmetrizable. However, in 2D the twogoals of accuracy and symmetry seem somewhat at odds with one another.For a given “cost” for example, the number of points in a given scheme, onecan “purchase” either a symmetric and less accurate or nonsymmetric andmore accurate scheme. We have considered both a first order accurate, sym-metric scheme and a second order accurate nonsymmetric scheme. Futurework could include pursuing a second order symmetric scheme.A significant reason for pursuing a symmetric, or at least symmetrizablescheme is the goal of finding a generalized form or approximation for theeigenvalues of the operator in the system of linear equations. In 1D, thishas proven achievable, as discussed in Chapter 3. A significant feature ofthe 1D operator in finding an approximation of its eigenvalues is its bandedstructure, even for a non-uniform mesh. This banded structure allows us todecompose our matrix and obtain the secular equation. In higher dimensionshowever the loss of a uniform, square mesh also means the loss of structurein our matrix. Approximation of the eigenvalues may still be possible, buthas not yet been successfully attempted.We also consider the different finite difference schemes from a computa-53Chapter 6. Conclusionstional point of view. We found implementing adaptivity in Matlab with thefinite difference scheme in Chapter 4 to be expensive as it required searchingfor neighbours and updating list entries as well as storage of several largedense arrays. The large sparse array of the second scheme was much morecapable for the purposes of adaptivity. We also found Julia a useful tool forthis work, though its true capabilities would likely have been more evidentfor larger problems.The use of a quadtree or bitree mesh, while restraining mesh size, provedvery helpful in an efficient implementation, especially with an adaptive mesh.With the added restriction in 2D of a graded mash, adapting and searchingfor possible neighbours was also made more efficient.This graded restriction may also be important in maintaining the form ofresulting operator; with a graded mesh a grid point can have no more thaneight neighbours. Without this restriction, some rows in the matrix couldhave many more entries and the operator would have even less structure. Ifthere is to be hope of finding or approximating eigenvalues in two or greaterdimensions, a highly structured operator is likely to be a requirement.54Bibliography[1][2] J. Bezanson, S. Karpinski, V. B. Shah, and A. Edelman. Julia: Afast dynamic language for technical computing. CoRR, http://arxiv.org/abs/1209.5145, 2012.[3] J. W. Demmel. Applied numerical linear algebra. Society for Industrialand Applied Mathematics (SIAM), Philadelphia, PA, 1997.[4] H. C. Elman and G. H. Golub. Iterative methods for cyclically re-duced non-self-adjoint linear systems. Mathematics of Computation,54(190):671–700, 1990.[5] R. E. Ewing, R. Lazarov, and P. Vassilevski. Local refinement tech-niques for elliptic problems on cell-centered grids. i. error analysis.Mathematics of Computation, 56(194):437–461, 1991.[6] P. Farrell, A. Hegarty, J. M. Miller, E. O’Riordan, and G. I. Shishkin.Robust Computational Techniques for Boundary Layers. CRC Press,2000.[7] R. A. Finkel and J. L. Bentley. Quad trees a data structure for retrievalon composite keys. Acta informatica, 4(1):1–9, 1974.[8] S. W. Golomb. Polyominoes: Puzzles, Patterns, Problems, and Pack-ings. Princeton University Press, 1996.[9] M. Gu and S. C. Eisenstat. A divide-and-conquer algorithm for thesymmetric tridiagonal eigenproblem. SIAM J. Matrix Anal. Appl.,16(1):172–191, Jan. 1995.[10] E. Haber and S. Heldmann. An octree multigrid method for quasi-staticmaxwells equations with highly discontinuous coefficients. Journal ofComputational Physics, 223(2):783–796, 2007.55Bibliography[11] L. Horesh and E. Haber. A second order discretization of maxwell’sequations in the quasi-static regime on octree grids. SIAM Journal onScientific Computing, 33(5):2805–2822, 2011.[12] W. Kaspar and R. Remke. Die numerische behandlung der poisson-gleichung auf einem gebiet mit einspringenden ecken. Computing,22(2):141–151, 1979.[13] P. Laasonen. On the discretization error of the Dirichlet problem ina plane region with corners. Annales AcademiæScientiarum Fennicæ,408:16, 1967.[14] MATLAB. version 8.3 (R2014a). The MathWorks Inc., Natick, Mas-sachusetts, 2014.[15] K. W. Morton. Numerical Solution of Convection-Diffusion Problems.Chapman and Hall, Oxford, UK, 1996.[16] K. W. Morton and D. F. Mayers. Numerical Solution of Partial Differ-ential Equations: An Introduction. Cambridge University Press, NewYork, NY, USA, 2005.[17] D. P. O’Leary and G. W. Stewart. Computing the eigenvalues andeigenvectors of symmetric arrowhead matrices. J. Comput. Phys.,90(2):497–505, Sept. 1990.[18] Y. Saad. Iterative methods for sparse linear systems. Siam, 2003.[19] A. Segal. Aspects of numerical methods for elliptic singular perturba-tion problems. SIAM Journal on Scientific and Statistical Computing,3(3):327–349, 1982.[20] J. Strikwerda. Finite Difference Schemes and Partial Differential Equa-tions, Second Edition. Society for Industrial and Applied Mathematics,2004.[21] M. Theillard, C. Rycroft, and F. Gibou. A multigrid method on non-graded adaptive octree and quadtree cartesian grids. Journal of Scien-tific Computing, 55(1):1–15, 2013.[22] J. H. Wilkinson, editor. The Algebraic Eigenvalue Problem. OxfordUniversity Press, Inc., New York, NY, USA, 1988.56"@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2016-02"@en ; edm:isShownAt "10.14288/1.0220517"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Computer Science"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivs 2.5 Canada"@* ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/2.5/ca/"@* ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Finite difference schemes for elliptic partial differential equations requiring a non-uniform mesh"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/55607"@en .