LEAST-SQUARES METHODS FOR COMPUTATIONALELECTROMAGNETICSA DissertationbyTZANIO VALENTINOV KOLEVSubmitted to the Oﬃce of Graduate Studies ofTexas A&M Universityin partial fulﬁllment of the requirements for the degree ofDOCTOR OF PHILOSOPHYAugust 2004Major Subject: MathematicsLEAST-SQUARES METHODS FOR COMPUTATIONALELECTROMAGNETICSA DissertationbyTZANIO VALENTINOV KOLEVSubmitted to Texas A&M Universityin partial fulﬁllment of the requirementsfor the degree ofDOCTOR OF PHILOSOPHYApproved as to style and content by:James H. Bramble(Co-Chair of Committee)Joseph E. Pasciak(Co-Chair of Committee)Raytcho D. Lazarov(Member)Vivek Sarin(Member)Albert Boggess(Head of Department)August 2004Major Subject: MathematicsiiiABSTRACTLeast-squares Methods for Computational Electromagnetics. (August 2004)Tzanio Valentinov Kolev, M.S., Soﬁa University “St. Kliment Ohridski”, BulgariaCo–Chairs of Advisory Committee: Dr. James H. BrambleDr. Joseph E. PasciakThe modeling of electromagnetic phenomena described by the Maxwell’s equationsis of critical importance in many practical applications. The numerical simulationof these equations is challenging and much more involved than initially believed.Consequently, many discretization techniques, most of them quite complicated, havebeen proposed.In this dissertation, we present and analyze a new methodology for approximationof the time-harmonic Maxwell’s equations. It is an extension of the negative-normleast-squares ﬁnite element approach which has been applied successfully to a varietyof other problems.The main advantages of our method are that it uses simple, piecewise polynomial,ﬁnite element spaces, while giving quasi-optimal approximation, even for solutionswith low regularity (such as the ones found in practical applications). The numericalsolution can be eﬃciently computed using standard and well-known tools, such asiterative methods and eigensolvers for symmetric and positive deﬁnite systems (e.g.PCG and LOBPCG) and preconditioners for second-order problems (e.g. Multigrid).Additionally, approximation of varying polynomial degrees is allowed and spuriouseigenmodes are provably avoided.ivWe consider the following problems related to the Maxwell’s equations in the fre-quency domain: the magnetostatic problem, the electrostatic problem, the eigenvalueproblem and the full time-harmonic system. For each of these problems, we present anatural (very) weak variational formulation assuming minimal regularity of the solu-tion. In each case, we prove error estimates for the approximation with two diﬀerentdiscrete least-squares methods. We also show how to deal with problems posed ondomains that are multiply connected or have multiple boundary components.Besides the theoretical analysis of the methods, the dissertation provides variousnumerical results in two and three dimensions that illustrate and support the theory.vTo my grandfather, (2. II. 1936 - 15.VI. 1996)in loving memory.viACKNOWLEDGMENTSI was very lucky that I had the opportunity to study in the Numerical Analysis groupat Texas A&M University. This has been a life-changing experience and I am mostgrateful to my advisors, Professors James Bramble and Joseph Pasciak. They showedme what it means to be a mathematician and taught me how to strive for perfectionin everything I do. Their insight, friendliness and enthusiasm are things I will alwaysremember and try to emulate in my own career.I wish to thank Professor Raytcho Lazarov, who suggested the Ph.D. programat Texas A&M to me. He has been the one constant support from the very beginningof my studies.I acknowledge the ﬁnancial support provided by the Department of Mathematicsand the Institute for Scientiﬁc Computations throughout my studies. I would alsolike to thank Ms. Monique Stewart for the many occasions on which I came to herfor help.I am grateful to Dr. Panayot Vasilevski for his help and mentoring during myvisits to the Center for Applied Scientiﬁc Computing (CASC) at Lawrence LivermoreNational Laboratory. I believe that those internships and the interaction with thegroup at CASC were essential to my education.This work would have not been possible without the help and support of mycolleagues, friends and my family, who encouraged me to study what I enjoy. Theone person however, whose support contributed the most to the completion of thisdissertation, is my ﬁance´e. Catrina, my apologies and love.viiTABLE OF CONTENTSCHAPTER PageI INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . 1II FUNCTION SPACES . . . . . . . . . . . . . . . . . . . . . . . . 8A. Hilbert spaces and operators . . . . . . . . . . . . . . . . . 8B. Sobolev spaces . . . . . . . . . . . . . . . . . . . . . . . . . 12C. Spaces of vector ﬁelds . . . . . . . . . . . . . . . . . . . . . 16D. Finite element subspaces . . . . . . . . . . . . . . . . . . . 23III AN ABSTRACT LEAST-SQUARES METHOD . . . . . . . . . 27A. Operator equations . . . . . . . . . . . . . . . . . . . . . . 28B. Approximation . . . . . . . . . . . . . . . . . . . . . . . . 37C. Implementation . . . . . . . . . . . . . . . . . . . . . . . . 43IV THEMAGNETOSTATIC AND THE ELECTROSTATIC PROB-LEMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46A. Weak formulation of the magnetostatic problem . . . . . . 47B. Weak formulation of the electrostatic problem . . . . . . . 54C. Least-squares approximation . . . . . . . . . . . . . . . . . 581. Approximation based on a discrete inf-sup condition . 59a. Pairs of stable approximation subspaces . . . . . 612. Approximation based on form modiﬁcation . . . . . . 663. Extensions to more general domains . . . . . . . . . . 70a. Domains with curved boundaries . . . . . . . . . 70b. Multiply-connected domains with holes . . . . . . 74V THE EIGENVALUE PROBLEM . . . . . . . . . . . . . . . . . 82A. Reformulation of the eigenvalue problem . . . . . . . . . . 86B. Approximation of the least-squares solution operators . . . 92C. The eigenvalue and eigenvector discretization . . . . . . . . 961. Improved estimate of the eigenvalue convergencerate for smooth eigenfunctions on a convex domain . . 99D. Extensions to more general domains . . . . . . . . . . . . . 103VI THE TIME-HARMONIC MAXWELL SYSTEM . . . . . . . . . 107viiiCHAPTER PageA. Weak formulation . . . . . . . . . . . . . . . . . . . . . . . 110B. Least-squares approximation . . . . . . . . . . . . . . . . . 1151. Approximation based on a discrete inf-sup condition . 1162. Approximation based on form modiﬁcation . . . . . . 1173. Error estimates . . . . . . . . . . . . . . . . . . . . . . 119VII NUMERICAL RESULTS . . . . . . . . . . . . . . . . . . . . . . 120A. Implementation issues . . . . . . . . . . . . . . . . . . . . 1201. Mesh characteristics . . . . . . . . . . . . . . . . . . . 122B. The magnetostatic problem . . . . . . . . . . . . . . . . . 1241. A problem with a known smooth solution . . . . . . . 1252. Magnetostatics in a L-shaped domain . . . . . . . . . 1263. Cross-section of a magnet . . . . . . . . . . . . . . . . 1274. Magnetic ﬁeld in a transformer . . . . . . . . . . . . . 129C. The eigenvalue problem . . . . . . . . . . . . . . . . . . . . 1291. Eigenvalues of the unit cube . . . . . . . . . . . . . . 1312. Eigenvalues of the unit ball . . . . . . . . . . . . . . . 1323. Eigenvalues of the Fichera corner . . . . . . . . . . . . 1354. Eigenvalues of a linear accelerator cell . . . . . . . . . 137D. The time-harmonic problem . . . . . . . . . . . . . . . . . 137VIII CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . 140REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141APPENDIX A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153APPENDIX B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156VITA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165ixLIST OF TABLESTABLE Page7.1 Mesh characteristics after uniform reﬁnement. . . . . . . . . . . . . . 1247.2 Numerical results for magnetostatic problem with a known smoothsolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1257.3 Numerical results for magnetostatics in an L-shaped domain. . . . . . 1267.4 Numerical results for the cross-section of a magnet. . . . . . . . . . . 1287.5 Numerical results for magnetostatics in a transformer. . . . . . . . . 1297.6 Unit ball, exact eigenvalues. . . . . . . . . . . . . . . . . . . . . . . . 1347.7 Unit ball, test meshes and number of LOBPCG iterations. . . . . . . 1347.8 Fichera corner, benchmark results from [100]. . . . . . . . . . . . . . 1357.9 Fichera corner, results for tetrahedral mesh (column 3) and hex-ahedral mesh (column 4). . . . . . . . . . . . . . . . . . . . . . . . . 1367.10 Numerical results for the time-harmonic test using exact solver. . . . 1387.11 Numerical results for the time-harmonic test using Multigrid. . . . . 139xLIST OF FIGURESFIGURE Page2.1 Typical geometry of the domain Ω. . . . . . . . . . . . . . . . . . . . 132.2 Face bubble functions: element of BFh in 2D and the bubbles foreach face of a tetrahedron in 3D. . . . . . . . . . . . . . . . . . . . . 254.1 Curved boundary approximation. . . . . . . . . . . . . . . . . . . . . 705.1 Cross-section of a coaxial cable with an oﬀset center conductor:pollution by spurious modes. After Paulsen and Lynch [84], c©1991 IEEE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 847.1 Comparison of the dimensions of test and solution spaces afteruniform reﬁnement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1237.2 Magnetostatics in an L-shaped domain, computed magnetic ﬁeld. . . 1277.3 Cross-section of a magnet: geometry and coarse mesh. . . . . . . . . 1277.4 Magnetostatic in transformer, geometry and coarse mesh. . . . . . . 1287.5 Unit cube, eigenvalue convergence. . . . . . . . . . . . . . . . . . . . 1317.6 Unit cube, approximation error. . . . . . . . . . . . . . . . . . . . . . 1327.7 Unit ball, initial mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . 1337.8 Unit ball, eigenvalue convergence. . . . . . . . . . . . . . . . . . . . . 1347.9 Fichera corner, initial meshes. . . . . . . . . . . . . . . . . . . . . . . 136B.1 Approximation to the magnetic ﬁeld in the iron core of the transformer.156B.2 Cross-section of the approximate solution ﬁeld for the transformerproblem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156B.3 Approximation to the magnetic ﬁeld in the transformer. . . . . . . . 157xiFIGURE PageB.4 Unit ball, eigenmode 1. . . . . . . . . . . . . . . . . . . . . . . . . . 158B.5 Unit ball, eigenmode 2. . . . . . . . . . . . . . . . . . . . . . . . . . 158B.6 Unit ball, eigenmode 3. . . . . . . . . . . . . . . . . . . . . . . . . . 158B.7 Unit ball, eigenmode 4. . . . . . . . . . . . . . . . . . . . . . . . . . 159B.8 Unit ball, eigenmode 5. . . . . . . . . . . . . . . . . . . . . . . . . . 159B.9 Unit ball, eigenmode 6. . . . . . . . . . . . . . . . . . . . . . . . . . 159B.10 Unit ball, eigenmode 7. . . . . . . . . . . . . . . . . . . . . . . . . . 160B.11 Unit ball, eigenmode 8. . . . . . . . . . . . . . . . . . . . . . . . . . 160B.12 Unit ball, eigenmode 9. . . . . . . . . . . . . . . . . . . . . . . . . . 160B.13 Unit ball, eigenmode 10. . . . . . . . . . . . . . . . . . . . . . . . . . 161B.14 Linear accelerator cell, eigenmode 1. . . . . . . . . . . . . . . . . . . 161B.15 Linear accelerator cell, eigenmode 2. . . . . . . . . . . . . . . . . . . 162B.16 Linear accelerator cell, eigenmode 3. . . . . . . . . . . . . . . . . . . 162B.17 Linear accelerator cell, eigenmode 4. . . . . . . . . . . . . . . . . . . 162B.18 Linear accelerator cell, eigenmode 5. . . . . . . . . . . . . . . . . . . 163B.19 Linear accelerator cell, eigenmode 6. . . . . . . . . . . . . . . . . . . 163B.20 Linear accelerator cell, eigenmode 7. . . . . . . . . . . . . . . . . . . 163B.21 Linear accelerator cell, eigenmode 8. . . . . . . . . . . . . . . . . . . 164B.22 Linear accelerator cell, eigenmode 9. . . . . . . . . . . . . . . . . . . 164B.23 Linear accelerator cell, eigenmode 10. . . . . . . . . . . . . . . . . . . 1641CHAPTER IINTRODUCTIONComputational electromagnetics is the science of applying modern computationaltechniques to numerically simulate the physical interactions and phenomena betweenelectromagnetic waves and material structures. This is of critical importance in manypractical applications, including the design of various devices: antennas, radars, mi-crowaves, waveguides and particle accelerators. Electromagnetic problems appearnaturally in diverse areas such as geophysics, relativity theory and optics. Speciﬁcapplications are discussed in many references, cf. [5, 59, 93, 47, 92]. The importance ofdeveloping advanced methods in computational electromagnetics is illustrated by thefollowing excerpt from the SciDAC project “Advanced Computing for 21st CenturyAccelerator Science & Technology” (see [103]):Particle accelerators have helped enable some of the most remarkable discov-eries of the 20th century. They have also led to substantial advances in appliedscience and technology, many of which greatly beneﬁt society. . . .Given theimportance of particle accelerators, it is imperative that the most advancedhigh performance computing tools be brought to bear on their design, opti-mization, technology development, and operation.Consider an isotropic, linear medium Ω with electric permittivity ε and magneticpermeability µ. Let E be the intensity of the electric ﬁeld generated by charges withvolume density ρ, and B be the intensity of the magnetic ﬁeld generated by currentwith volume density J. Maxwell suggested (see [73] for the original and [18, 90, 10, 57]for a modern presentation) that, when these ﬁelds depend on time, they are coupledThis dissertation follows the format of SIAM Journal of Numerical Analysis.2by the following system of equations: 1⎧⎪⎨⎪⎩∇×E = − ∂∂tB∇·D = ρ,⎧⎪⎨⎪⎩∇×H = ∂∂tD+ J∇·B = 0. (1.1)Here D and H are the densities of the electric and the magnetic ﬂux, which in thelinear case are given byD = εE , H = µ−1 B . (1.2)Theoretically (1.1) should be solved on all of R3. However, one usually computes in asuﬃciently large domain, which is assumed to be surrounded by a perfect conductor.The boundary conditions in this case are:E×n = 0 , B · n = 0 on ∂Ω , (1.3)where n denotes the outward unit normal on the boundary.Even though they will not be considered in this dissertation, we should remarkthat physically more meaningful radiation boundary conditions are possible, see e.g.[79]. A more advanced treatment can be achieved by using absorbing boundaryconditions as the perfectly matched layer technique given in [12, 13], see also [59].We also note that there are more general frameworks in which to understand theabove equations. For example, in [56] the electromagnetic phenomena are describedin the language of diﬀerential geometry and algebraic topology. The discretizationis based on discrete diﬀerential forms, which are a generalization of the Lagrangianﬁnite elements.Commonly in practice, only one or few frequencies of propagations are considered.Based on that, or by applying the Fourier transform, one can reduce the Maxwell’s1The equations involving the curl operator correspond to Faraday’s and Ampere’slaws, while the divergence equations are called Gauss’ electric and magnetic laws.3equations to their time-harmonic form. The assumption that the ﬁelds vary harmoni-cally in time with frequency ω means that E(x, t) = e0(x) cos(ωt+φE) = (e(x) eiωt)and H(x, t) = h0(x) cos(ωt + φH) = (h(x) eiωt), where e(x) = e0(x) eiφE andh(x) = h0(x) eiφH are some complex ﬁelds. Assuming that the data are also time-harmonic, J(x, t) = (j(x) eiωt), the equations (1.1)–(1.3) take the following form,known as the time-harmonic Maxwell system⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩∇×e = −λµh in Ω,∇×h = λ εe+ j in Ω,e×n = 0 on ∂Ω,µh · n = 0 on ∂Ω .(1.4)Here λ = i ω, the current density j is given, and we are looking for the magnetic andelectric ﬁelds h ,e : Ω→ C3.In realistic computations this problem is posed on complicated, three-dimensionaldomains where a natural choice for a discretization technique is the ﬁnite elementmethod. There is extensive literature on the use of ﬁnite elements in computationalelectromagnetics, see [75, 59, 93, 58].In two dimensions, most of the electromagnetic problems can be reduced tosecond-order problems for one of the ﬁelds or for a potential. However, the threedimensional problems are signiﬁcantly more complicated, in particular due to thelarge nullspace of the curl operator. This suggests that a new set of methods isrequired for the problem (1.4). Indeed, the straightforward application of standardpiecewise linear elements to the eigenvalue problem (1.7), related to (1.4), leads tospurious eigenmodes as shown in [16, 93].A considerable amount of research has been targeted speciﬁcally to computa-tional electromagnetics. Many methods have been proposed, each with its advantages4and drawbacks. Some of them are discussed below.A new set of ﬁnite element spaces that seems to ﬁt the Maxwell problem wasgiven by Ne´de´lec in [77]. Their curl-conforming property eliminates the spuriousmodes and leads to optimal convergence. Since their introduction, Ne´de´lec elementshave been considered the natural choice in many electromagnetic problems, and theresearch activity on this topic has been very active (cf. [75]). However, the Ne´de´lec el-ements, especially those of higher order, have the drawback of being relatively diﬃcultto implement. The resulting algebraic system usually needs special, sophisticated so-lution algorithms. There also seems to be a lack of clear theory for general hexahedralmeshes.Some methods use the standard nodal ﬁnite element spaces but modify the bi-linear form to ensure ellipticity. This is the approach taken in [43, 40, 42, 85]. Thedrawback of these methods is that the added complexity in the form evaluation maysurpass the convenience of working with simple ﬁnite elements. Furthermore, whenapplied to the eigenvalue problem, the modiﬁed form may introduce additional familyof eigenpairs as discussed in [41].A diﬀerent set of ideas, which are closest to the one considered in this dissertation,are based on the least-squares ﬁnite element method. The standard functionals usedwidely in the engineering community are L2-based (see [58, 96]). Related second-order problems can be treated by these methods after the introduction of additionalvariables which reduce the system to ﬁrst-order system least-squares (FOSLS). Forexample, in [70] a FOSLS method is applied to the scalar Helmholtz equation withexterior radiation boundary conditions to derive an algorithm uniform with respectto the wave number. This result is obtained under the assumption that the domainis convex or has a smooth boundary.5The least-squares ﬁnite element method is well studied, in particular, for second-order problems. Among the many papers that deal with this subject are [14, 58,31, 32, 33]. The dual, or negative-norm, approach is described in [21, 22, 23]. Itseems that [26] is the ﬁrst time when such a method was applied to electromagneticproblems.Motivated by the previous discussion, in this dissertation we develop and analyzea new methodology in computational electromagnetics—the least-squares methodbased in a dual space. Speciﬁcally, this dissertation deals with the approximation ofthe full time-harmonic system (1.4) and the following related problems:the (generalized) magnetostatic problem⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×h = j in Ω,∇ · (µh) = ρ in Ω,µh · n = σ on ∂Ω,(1.5)which may model the magnetic ﬁelds produced by steady currents;the (generalized) electrostatic problem⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×e = j in Ω,∇ · (εe) = ρ in Ω,e×n = σ on ∂Ω,(1.6)which may describe the electric ﬁelds produced by stationary source charges;6and the eigenvalue problem⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩∇×h = λ εe in Ω,∇×e = −λµh in Ω,e×n = 0 on ∂Ω,µh · n = 0 on ∂Ω ,(1.7)which gives the frequencies of the ﬁelds that will propagate through a given medium.The proposed method is based on natural weak variational formulations of (1.5),(1.6) and (1.4) which assume minimal regularity of the solution. The solution oper-ators for the ﬁrst two problems are further used to obtain an approximation to theeigenvalue problem.The resulting discretization method has the advantages of avoiding potentialsand the use of Ne´de´lec spaces. In fact, the mixing of continuous and discontin-uous approximation spaces of varying polynomial degrees is allowed. The theoryand implementation for general hexahedral meshes is analogous to that on tetrahe-dra. Additional advantages are that the matrix of the discrete system is uniformlyequivalent to the mass matrix and that spurious eigenmodes are completely avoided.Finally, the method can be eﬃciently implemented using preconditioners for standardsecond-order problems (e.g. Multigrid).The outline of the contents of the dissertation is as follows. In Chapter II, wediscuss the needed notation and basic facts from ﬁnite element theory and functionalanalysis. These results are standard and are needed in the subsequent development.Next, we present and analyze an abstract least-squares algorithm in Chapter III.Here we formulate general approximation results and address the question of imple-mentation. The following chapter deals with the theory for the electrostatic and themagnetostatic problems. We mostly follow the theory from [26]. This material is7included for completeness, since it is the basis upon which the rest of the dissertationis built. We give a detailed presentation and provide further results concerning stablepairs of approximation spaces, regularity and extensions to domains with holes andcurved boundaries. In Chapter V, the eigenvalue problem is discussed. We start witha reformulation of the original problem to an eigenvalue problem based on the solu-tion operators for (1.5) and (1.6). Then we show how to approximate those solutionoperators and investigate the convergence of eigenvectors and eigenvalues. The topicof Chapter VI is the least-squares method for the full time-harmonic system. Thedevelopment here is similar to the one in Chapter IV, but it is also naturally con-nected to the results for the eigenvalue problem. In the next chapter we present andcomment on various numerical experiments illustrating the theory. The last chapterof the dissertation contains the conclusions including plans for possible future work.A ﬁnal note on notation: we use the symbol C with or without subscript to denotea generic positive constant, which may be diﬀerent in the diﬀerent occurrences. Thisconstant may depend on explicitly stated quantities, but it will always be independentof the mesh size h.8CHAPTER IIFUNCTION SPACESIn this chapter, we recall a few concepts and results that will be needed later. We col-lect material from various sources and try to present it brieﬂy and with the appropriatereferences. For notation, deﬁnitions, and further details, see [62, 3, 72, 74, 53, 97, 98].A. Hilbert spaces and operatorsLet X be a Hilbert space with an inner product (·, ·)X. In this dissertation, we assumethat X is separable and deﬁned over the ﬁeld K, which is either R or C.The dual space of X is denoted by X∗ and consists of all bounded conjugate-linearfunctionals : X → K. Here, conjugate-linear means that (λx + y) = λ (x) + (y)for any λ ∈ K; x,y ∈ X. Clearly is conjugate-linear if and only if ¯, deﬁned by¯(x) = (x), is linear. The norm on X∗ is deﬁned by‖‖X∗ = supx∈X\{0}|〈, x〉|‖x‖X ,where 〈·, ·〉 ≡ 〈·, ·〉X∗×X denotes the duality pairing between X∗ and X. By the RieszRepresentation Theorem, there exists a linear isometry TX : X∗ → X, satisfying(TX, x)X = 〈, x〉 ∀x ∈ X . (2.1)It follows from the polarization identity(x,y)X =14(‖x+ y‖2X− ‖x− y‖2X + i ‖x+ iy‖2X− i ‖x− iy‖2X) , (2.2)and the fact that a Banach space is a Hilbert space if and only if the parallelogramidentity‖x+ y‖2X + ‖x− y‖2X = 2 ‖x‖2X + 2 ‖y‖2X (2.3)9holds, that X∗ is a Hilbert space with an inner product(, j)X∗ = 〈,TXj〉 = 〈j,TX〉 = (TX,TXj)X . (2.4)If L is a subspace of X, the quotient space X/L consists of all equivalence classesunder the equivalence relation u ∼ v ⇐⇒ u− v ∈ L. The orthogonal complementof L in X is deﬁned as L⊥X ≡ L⊥ = {x ∈ X : (x, l)X = 0 ,∀l ∈ L}. We recall thatwhen L is closed, X = L⊕ L⊥, and X/L is isomorphic to L⊥.Let X and Y be two Hilbert spaces. The set of all bounded linear operators fromX to Y is denoted by L(X,Y). This is a Banach space with respect to the operatornorm‖A‖ ≡ ‖A‖X→Y ≡ ‖A‖L(X,Y) = supx∈X\{0}‖Ax‖Y‖x‖X .When X ⊆ Y and the identity operator is in L(X,Y) we use X ↪→ Y to denote thatX is continuously embedded Y. We say that A ∈ L(X,Y) deﬁnes an isomorphismbetween X and Y if A is bijective, bounded and A−1 is also bounded. The operatorA ∈ L(X,Y) is said to be compact if it maps bounded sets in X into sets with compactclosure in Y. The following sets denote the kernel and the image of A:N(A) = {x ∈ X : Ax = 0} , R(A) = {Ax ∈ Y : x ∈ X} .Remark 2.1 Let X and Y be two Hilbert spaces that are continuously embedded in anormed space Z. Then X ∩ Y is a Hilbert space with an inner product(x,y)X∩Y = (x,y)X + (x,y)Y ∀x ,y ∈ X ∩ Y .Remark 2.2 Let X be a real Hilbert space. Analogous to the construction of C asR × R, we can think of X× X as a complex Hilbert space, denoted with XC. Inparticular ‖x+ iy‖2XC = ‖x‖2X + ‖y‖2X, for any x ,y ∈ X.10The operator A ∈ L(X,Y) can be naturally extended to AC ∈ L(XC,YC) by deﬁn-ing AC(x+ iy) = Ax+ iAy. Note that ‖AC‖XC→YC = ‖A‖X→Y.A form a : X×Y → K is said to be bilinear 1 if it is linear with respect to itsﬁrst argument and conjugate-linear with respect to the second. A bilinear form isbounded, with a bound ‖a‖, if|a(x,y)| ≤ ‖a‖ ‖x‖X‖y‖Y ∀(x,y) ∈ X×Y .We say that a(·, ·) satisﬁes the inf-sup condition, if there exists a constant C ∈ R+such thatC ‖x‖X ≤ supy∈Y\{0}|a(x,y)|‖y‖Y , ∀x ∈ X . (2.5)The following result is well known (see e.g. [7]).Theorem 2.1 (Generalized Lax-Milgram) Suppose that a(·, ·) is a bounded bi-linear form on X×Y satisfying the inf-sup condition (2.5). DeﬁneY0 = {y ∈ Y : a(x,y) = 0, for all x ∈ X} .Then, for any f ∈ Y∗ there exists a unique x ∈ X satisfyinga(x,y) = 〈f,y〉 ∀y ∈ Y , (2.6)if and only if〈f,y〉 = 0 ∀y ∈ Y0 . (2.7)Furthermore, the solution satisﬁesC ‖x‖X ≤ ‖f‖Y∗ ≤ ‖a‖ ‖x‖X . (2.8)1In the case K = C, the bilinear forms are also called sesquilinear.11Next, we discuss the spectral properties of an operator A ∈ L(X,X). For anyλ ∈ C, the resolvent operator Rλ(A) is deﬁned as Rλ(A) = (λI − A)−1. Theresolvent set ρ(A), and the spectrum σ(A), are deﬁned by ρ(A) = {λ ∈ C :Rλ(A) is an isomorphism on X}, and σ(A) = C \ ρ(A). We say that λ ∈ C is aneigenvalue of A if there is x = 0 such that Ax = λx. The set of all such x forms thelinear subspace of eigenvectors corresponding to λ, and is denoted with Vλ.The adjoint operator A∗ ∈ L(X,X) is deﬁned by(Ax,y) = (x,A∗y) ∀x,y ∈ X .The operator A is called symmetric, or Hermitian if A = A∗. When A = −A∗, theoperator is called skew-Hermitian. Clearly A is skew-Hermitian if and only if iA isHermitian. A Hermitian operator is positive semi-deﬁnite if(Ax, x) ≥ 0 ∀x ∈ X .When the equality above is achieved only for x = 0, the operator is called positivedeﬁnite. Using this notation, we can formulate some basic theorems from the spectraltheory of operators on Banach spaces.Theorem 2.2 (Hilbert-Schmidt Theory) Let A ∈ L(X,X) be a compact opera-tor. Let {λn} be the set of its nonzero eigenvalues. We have the following results:1. Each of the spaces Vλn is of ﬁnite dimension, called the multiplicity of λn. Thespectrum of A is {λn} ∪ {0}.2. The nonzero eigenvalues of A∗ are precisely {λn}. Furthermore λn and λn havethe same multiplicity.3. If the number of nonzero eigenvalues is not ﬁnite, it is countable, and they canbe ordered in a sequence λn → 0.124. If A is also Hermitian, all the eigenvalues are real. If A is skew-Hermitian, allthe eigenvalues are purely imaginary. In both cases N(A)⊥ =⊕λnVλn.5. If A is symmetric and positive semi-deﬁnite, the eigenvalues are positive and‖A‖ = λ1 > λ2 > . . . > λn > . . . ≥ 0.Theorem 2.3 (Fredholm Alternative) Let A ∈ L(X,X) be a compact and self-adjoint operator. For λ = 0 and b ∈ X consider the equationAx− λx = b . (2.9)Then:1. If λ ∈ σ(A) then (2.9) has a unique solution x for any b.2. If λ is an eigenvalue, then (2.9) has a solution if and only if b ∈ V⊥λ . Thesolution is unique in the quotient space X/Vλ.B. Sobolev spacesLet Ω be a nonempty, bounded connected open set in Rd, d ∈ {2, 3}. Then Ω ismeasurable and its Lebesgue measure, cf. [3], is denoted by µ(Ω). We assume thatthe boundary ∂Ω is Lipschitz continuous (see [55] for the deﬁnition). In this case, theoutward unit normal n is well deﬁned almost everywhere on ∂Ω.The connected components of ∂Ω are denoted by Γi, i = 0, . . . , n1, where Γ0 isthe exterior boundary, i.e. Γi ⊂ int(Γ0) for i = 1, . . . , n1. As in Hypothesis 3.3 from[4], we assume that there exist a ﬁnite number of cutting surfaces Σj, j = 1, . . . , n2,so that the domain Ω0 = Ω \⋃n2j=1 Σj is simply connected.The simplest example of such a domain is a convex open set Ω ⊂ Rd, in whichcase n1 = n2 = 0. In two dimensions we always have n1 = n2. However, in R3, n113equals the number of connected bounded components of R3 \ Ω, while n2 is equal tothe genus of ∂Ω. Informally, we say that the domain has n1 “holes” and n2 “loops”.Clearly these two numbers are independent.Assumption (AΩ) The domain Ω is nonempty, open, bounded, connected and hasLipschitz continuous boundary with n1 “holes” and n2 “loops”.The above assumption allow us to consider domains as the one shown on Figure2.1.Γ0Γ1Γ2Γ3Σ1ΩFig. 2.1. Typical geometry of the domain Ω.We start with the following spaces of functions deﬁned on Ω: D(Ω) ≡ C∞0 (Ω) isthe set of all inﬁnitely smooth functions with compact support in Ω, and D′(Ω) isthe set of all distributions (the continuous linear functionals on D(Ω) with the weakstar topology). For p ∈ [1,∞), Lp(Ω) is the Banach space of classes of Lebesgue-measurable functions for which the norm‖f‖Lp =(∫Ω|f(x)|pdx) 1pis ﬁnite.Remark 2.3 In this section, we concentrate on spaces of real-valued functions, i.e.the case K = R. The extension to complex-valued functions and vector ﬁelds is14straightforward (see Remark 2.2). When we want to emphasize the ﬁeld of scalars fora given space, we will use a subscript notation like LpC(Ω) and LpR(Ω).Let α = (αi)di=1 ∈ Nd be a multiindex and ∂αf denote the distributional (or weak)derivative of f ∈ D′(Ω) of order |α| = ∑di=1 αi. When |α| = 0, we set ∂αf = f. Fors ∈ N0 and integer p ∈ (1,∞), the Sobolev space Ws,p(Ω) consist of distributions fwhich are in Lp(Ω) together with all their derivatives of order less or equal to s. Thisis a Banach space with respect to the norm‖f‖Ws,p =(s∑k=1|f|pWk,p) 1p, where |f|Ws,p =⎛⎝∑|α|=s‖∂αf‖pLp⎞⎠ 1p .In particular, Ws,2(Ω) is a Hilbert space, traditionally denoted by Hs(Ω). For conve-nience we will use ‖ ·‖s, | · |s, and (·, ·)s for the norm, seminorm and the inner producton Hs(Ω).The deﬁnition of Sobolev spaces can be extended to s ∈ R+ as follows: if s =m + σ, with m ∈ N0 and σ ∈ (0, 1), then ‖f‖Ws,p = (‖f‖pWm,p + |f|pWs,p)1p , where|f|Ws,p =⎛⎝∑|α|=m∫Ω∫Ω|∂αf(x)− ∂αf(y)|p‖x− y‖d+σp dx dy⎞⎠ 1p . (2.10)The spaces Hs(Ω), s ∈ R+ can be alternatively deﬁned by the real method ofinterpolation, see [91, 3] and Appendix A in [29]. This is particularly useful since itallows for obtaining estimates for bounded linear operators in intermediate spaces by“interpolation” (cf. Theorem 1.4 in [54]).Introduce Ws,p0 (Ω) as the closure of D(Ω) in ‖ · ‖Ws,p . The space W−s,p(Ω) ≡W−s,p0 (Ω) is deﬁned as the dual of Ws,q0 (Ω), where1p+ 1q= 1. In particular Hs0(Ω) =Ws,20 (Ω), and H−s(Ω) ≡ H−s0 (Ω) = Hs0(Ω)∗.Denote D(Ω) to be the space of restrictions of functions in D(Rd) to Ω. It is15well known that the trace operator γ0, deﬁned on D(Ω), can be uniquely extended toa bounded linear operator from Hs(Ω) onto Hs−12 (∂Ω), for s ∈ (12, 1]. Moreover, fors = 1 we have N(γ0) = H10(Ω). We recall that the Sobolev spaces on the boundarycan be deﬁned by the use of local charts. In particular, Hs(Γi) is well deﬁned fors ∈ (12, 1], i = 0, . . . , n1 and any u ∈ Hs(Ω) has traces u|Γi ∈ Hs−12 (Γi).For s ∈ R+, let H˜s0(Ω) be the space of functions f for which f˜, the extensions of fby 0 outside of Ω is in Hs(Rd). The dual of H˜s0(Ω) is denoted with H˜−s(Ω) or H˜−s0 (Ω).In the next Theorem, we summarize some results that will be needed later. Formore details, including the Sobolev embedding theorem, an equivalent description ofHs(Rd) in terms of the Fourier transform, and the case p =∞ we refer to [3, 55, 54]and Appendix A in [29].Theorem 2.4 Let [X,Y]s denote the interpolation space between X and Y (with s = 0corresponding to X). The following hold true1. Hs1(Ω) is compactly embedded in Hs2(Ω) for any real s1 > s2. This means thatevery bounded sequence in Hs1(Ω) has a convergent subsequence in Hs2(Ω).2. D(Ω) is dense in Hs(Ω) for any s ≥ 0.3. There exists a bounded linear extension operator E : Hs(Ω) → Hs(Rd), indepen-dent of s > 0 such that Ef|Ω = f.4. For any |α| = 1, and s ∈ R, s − 12∈ Z, the weak derivative ∂α is a boundedlinear operator from Hs(Ω) to Hs−1(Ω). In addition, ∂α is a bounded linearoperator from H12 (Ω) to H˜−12 (Ω).5. For any s ∈ R+ there exists C = C(Ω, s) > 0 such that ‖u‖s ≤ C |u|s for allu ∈ Hs0(Ω) (Poincare´’s inequality).6. There exists C = C(Ω) > 0 such that C ‖u‖0 ≤ ‖u‖−1 + ‖∇u‖−1 for all16u ∈ L2(Ω) (Necˇas inequality, see [80]).7. Hs(Ω) = Hs0(Ω), for |s| ≤ 12 .8. H1+s0 (Ω) = H10(Ω) ∩H1+s(Ω), for 0 ≤ s ≤ 12 .9. [L2(Ω),H10(Ω)]s = H˜s0(Ω) for s ∈ [0, 1].10. H˜s0(Ω) = Hs0(Ω) for s ∈ [−1, 1], |s| = 12 .11. [H10(Ω),H10(Ω) ∩H2(Ω)]s = H10(Ω) ∩H1+s(Ω) for s ∈ [0, 1], see [9].Remark 2.4 The space H˜120 (Ω) is a proper subspace of H120 (Ω), usually denoted byH1200(Ω). It is a Hilbert space with norm,‖f‖H1200(Ω)=(‖f‖2H12 (Ω)+ |f|2H1200(Ω)) 12, where |f|H1200(Ω)=∥∥∥∥ f√ρ∥∥∥∥L2and ρ(x) = infy∈∂Ω ‖x− y‖ denotes the distance from x ∈ Ω to the boundary.C. Spaces of vector ﬁeldsWe adopt the notation of using boldface symbols to denote vector quantities andspaces. In particular, D(Ω) = D(Ω)d, D′(Ω) = D′(Ω)d, Lp(Ω) = Lp(Ω)d, Ws,p(Ω) =Ws,p(Ω)d, Hs(Ω) = Hs(Ω)d and H˜s(Ω) = H˜s(Ω)d. The norm and the inner productsare naturally inherited. For example(f,g)L2(Ω) =d∑i=1(fi,gi)L2(Ω) ,for any f = (f1, . . . , fd) ∈ L2(Ω), g = (g1, . . . ,gd) ∈ L2(Ω).The distributional divergence ∇· ≡ div : D′(Ω) → D′(Ω) is deﬁned by〈∇·f, ϕ〉 = −〈f,∇ϕ〉 ∀ϕ ∈ D(Ω) . (2.11)The space H(div) = {v ∈ L2(Ω) : ∇·v ∈ L2(Ω)} is a Hilbert space (see Remark 2.1)17with respect to the inner product(u, v)H(div) = (u, v)L2(Ω) + (∇·u,∇·v)L2(Ω) .The closure of its subspace D(Ω) is denoted by H0(div).Let ∇× ≡ curl : D′(Ω) → D′(Ω) be the distributional curl operator, deﬁnedby〈∇×f,ϕ〉 = 〈f,∇×ϕ〉 ∀ϕ ∈D(Ω) . (2.12)Depending on the argument, the standard curl operator on the right is given by oneof the matrices∇× =⎛⎜⎜⎜⎜⎝0 −∂z ∂y∂z 0 −∂x−∂y ∂x 0⎞⎟⎟⎟⎟⎠ , ∇× =(−∂y ∂x), or ∇× =⎛⎜⎝ ∂y−∂x⎞⎟⎠ . (2.13)Deﬁne H(curl) = {v ∈ L2(Ω) : ∇×v ∈ L2(Ω)}. This is a Hilbert space withrespect to the inner product(u, v)H(curl) = (u, v)L2(Ω) + (∇×u,∇×v)L2(Ω) .An important subspace is H0(curl) = D(Ω)H(curl). The next theorem summarizessome results from [54].Theorem 2.5 The spaces H(div) and H(curl) have the following properties.1. D(Ω) is dense in H(div) and H(curl).2. Any u ∈ D′(Ω) and v ∈D′(Ω) satisfy∇·(∇×v) = 0 , ∇×(∇u) = 0 , and ∇×(∇×v) = −∆v+ ∇(∇·v) .3. The mapping γn : v ∈D(Ω) → v · n can be extended to a surjective continuous18linear map γn : H(div) → H− 12 (∂Ω).4. For any v ∈ H(div) and φ ∈ H1(Ω) we have the Green’s formula 2:(v,∇φ)L2(Ω) + (∇·v, φ)L2(Ω) = 〈v · n, γ0(φ)〉 . (2.14)5. The mapping3 γτ : v ∈ D(Ω) → v×n can be extended to a continuous linearmap4 γτ : H(curl) → H−12 (∂Ω).6. For any v ∈ H(curl) and u ∈ H1(Ω) we have the Green’s formula:(∇×v,u)L2(Ω) − (v,∇×u)L2(Ω) = 〈v×n, γ0(u)〉 . (2.15)7. N(γn) = H0(div), N(γτ ) = H0(curl) and H10(Ω) = H0(div) ∩H0(curl) 5.8. If Ω is simply connected (n2 = 0) and v ∈ L2(Ω), then ∇×v = 0 if and only ifthere exists a unique p ∈ H1(Ω)/R such that v = ∇p.9. If ∂Ω is connected (n1 = 0) and v ∈ L2(Ω), then ∇·v = 0 if and only if thereexists w ∈ H1(Ω) with ∇·w = 0, such that v = ∇×w.10. If Ω ⊂ R2, v = (v1, v2) and v⊥ = (−v2, v1), then ∇·v = ∇×v⊥. In particular,H(curl) = H(div)⊥ = {v⊥ : v ∈ H(div)}.We will need to work with spaces that depend on a real-valued function γ, whichmay be the electric permittivity ε, the magnetic permeability µ or one of their recip-rocals. In some physical applications these may be complex or nonlinear functions2The case φ ≡ 1 is also known as the Divergence Theorem.3In R2, v×n = v · t, where t = n⊥ is the vector, tangential to the boundary.4γτ is not surjective, see the discussion in [75], pp.58-59.5In fact, this is an isometry since for any v ∈ H10(Ω), we have by density|v|2H1(Ω) = ‖∇×v‖2L2(Ω) + ‖∇·v‖2L2(Ω) .19and may even exhibit hysteresis, depending on the solution and its history. However,we shall only consider the case when ε and µ are piecewise smooth, real functionsthat are bounded and bounded away from zero on Ω. This is formalized below.Assumption (Aµ,ε) The functions ε , µ are in L2(Ω) and there exist constants µ0,µ1, ε0, ε1 satisfying 0 < µ0 ≤ µ(x) ≤ µ1 and 0 < ε0 ≤ ε(x) ≤ ε1, a.e. x ∈Ω. Furthermore, Ω can be split into non-overlapping Lipschitz subdomains {Ωi},satisfying (AΩ), such that ε|Ωi , µ|Ωi ∈ H1(Ωi).In practice, diﬀerent Ωi correspond to diﬀerent materials and the nonempty inter-sections ∂Ωi ∩ ∂Ωj are called (material) interfaces. The vector ﬁelds in H(div) andH(curl) satisfy continuity conditions across the interfaces as described below.Theorem 2.6 Suppose that Ω is split into non-overlapping Lipschitz subdomains{Ωi} which are either polygonal or have boundaries of class C1,1 (see [54]). Letv ∈ L2(Ω) be such that v|Ωi ∈ H1(Ωi). Thenv ∈ H(div) if and only if v · n = 0and similarlyv ∈ H(curl) if and only if v×n = 0 ,where · denotes the jump across the interfaces ∂Ωi ∩ ∂Ωj.Proof Deﬁne w ∈ L2(Ω) by w|Ωi = ∇·(v|Ωi). For any ϕ ∈ D(Ω) we have〈∇·v, ϕ〉 = −∑i(v|Ωi ,∇ϕ)L2(Ωi) = (w, ϕ)L2(Ω) −∑i〈v|Ωi · n, ϕ〉 .The assumption on ∂Ωi implies that v|Ωi · n ∈ L2(∂Ωi). Therefore ∇·v = w in L2(Ω)if and only if v|Ωi · n = v|Ωj · n in L2(∂Ωi ∩ ∂Ωj). The argument for H(curl) issimilar.20Relative to the partition, for any real power s, we deﬁne the family of piecewisespacesPHs(Ω) =⊕Hs(Ωi) , PHs0(Ω) =⊕Hs0(Ωi) . (2.16)Note that for s ∈ [0, 1/2), we have PHs(Ω) = PHs0(Ω) = Hs(Ω)For γ ∈ {ε, µ, ε−1, µ−1}, let L2γ(Ω) be the space L2(Ω) equipped with the weightedinner product (u, v)γ = (γ u, v)L2(Ω). The induced norm on L2γ(Ω) will be denotedwith ‖ · ‖γ . Additionally, deﬁneH(div; γ) = {v ∈ L2(Ω) : ∇·(γv) ∈ L2(Ω)},H0(div; γ) = {v ∈ H(div; γ) : n · v = 0 on ∂Ω} ,X1(µ) = H(curl) ∩H0(div;µ) ,X2(ε) = H0(curl) ∩H(div; ε) .These are Hilbert spaces according to Remark 2.1.Physical, as well as mathematical, considerations imply that the natural spacesfor the electromagnetic ﬁelds are h ∈ X1(µ) and e ∈ X2(ε). Thus, the regularity ofthese spaces is of primary interest.Theorem 2.7 The continuous embeddingsX1(µ),X2(ε) ↪→ Hs(Ω) , (2.17)hold true in the following cases:1. If ε and µ are smooth and the domain is a convex polygon/polyhedron, thens = 1 (proved by Saranen in [86]).2. If ε and µ are smooth and the domain is Lipschitz polygon/polyhedron, thens > 1/2 (proved by Amrouche, Bernardi, Dauge and Girault in [4]).213. If ε and µ are piecewise constants, then s can be arbitrarily close to 0 (provedby Costabel, Dauge and Nicaise in [45]).4. For any ε and µ satisfying (Aµ,ε), the embeddings for s = 0 are compact (provedby Weber in [94]). The boundary conditions are essential since the embeddingof H(curl) ∩H(div) in L2(Ω) is not compact, as shown in [4].For constant ε and µ, the paper [42] gives an explicit representation of the ﬁeldsX1(µ) and X2(ε) as a regular part plus a gradient of the solution of Neumann andDirichlet problems posed on Ω. In this case, s in (2.17) can be chosen as the regularityof the above problems minus 1. For piecewise constant ε and µ, as shown in [45],the regularity of X1(µ),X2(ε) is related to the operators −∆Dirε and −∆Neuµ deﬁnedbelow.For f ∈ H−1(Ω), we set −∆Dirε u = f, where u ∈ H10(Ω) satisﬁes(ε∇u,∇ϕ) = 〈f, ϕ〉 ∀ϕ ∈ H10(Ω) . (2.18)Similarly, for f ∈ (H1(Ω))∗, with 〈f, 1〉 = 0 we set −∆Neuµ u = f, where u ∈ H1(Ω)/Rsatisﬁes(µ∇u,∇ϕ) = 〈f, ϕ〉 ∀ϕ ∈ H1(Ω) . (2.19)Recall that if Ω is convex, the operator −∆ : u → −∆u is an isomorphism ofH10(Ω) ∩ H2(Ω) onto L2(Ω). In general, e.g. on polygonal/polyhedral domains, thepresence of reentrant corners leads to lower regularity. This is characterized by anumber s > 0 such that −∆ is an isomorphism of H10(Ω) ∩ H1+(Ω) onto H˜−1+0 (Ω)for any 0 ≤ ≤ s. The above is a motivation for the next two assumptions.Assumption (AL2∆Dirε ,∆Neuµ ) There exists s ∈ (0, 1] such that when f ∈ L2(Ω), thesolutions of the problems (2.18) and (2.19) are in H1+s(Ω) and ‖u‖1+s ≤ C ‖f‖.22Assumption (A∆Dirε ,∆Neuµ ) There exists s0 ∈ (0, 1] such that when 0 ≤ s ≤ s0 andf ∈ H˜−1+s0 (Ω), the solutions of the problems (2.18) and (2.19) are in H1+s(Ω) and‖u‖1+s ≤ C ‖f‖−1+s.The validity of regularity results related to the above assumptions for piecewisesmooth coeﬃcients was investigated in [45, 46].We ﬁnish this section with a list of Helmholtz-like decomposition results 6. Forsimplicity, we assume that Ω is either simply connected (n2 = 0) or it has one bound-ary component (n1 = 0). The more general cases will be addressed in Chapter IV,§C.3.b.Theorem 2.8 (cf. [26]) Let u be in L2(Ω). Then it can be decomposed asu = ∇×w+ µ∇ψ (2.20)in the following spaces1. For n2 = 0, w ∈ H0(curl) and ψ ∈ H1(Ω) with ∇·w = 0.2. For n2 = 0, w ∈ H10(Ω) and ψ ∈ H1(Ω).3. For n1 = 0, w ∈ H1(Ω) and ψ ∈ H10(Ω) with ∇·w = 0.The decompositions are orthogonal in L2µ−1(Ω). In the last two cases, we additionallyhave‖w‖H1(Ω) ≤ C ‖∇×w‖.Proof Let ψ be the unique element of H1(Ω)/R satisfying(µ∇ψ,∇θ) = (u,∇θ) , (2.21)6i.e. a splitting of a ﬁeld in a solenoidal and irrotational parts, see [57].23for any θ ∈ H1(Ω). Then ∇·(u−µ∇ψ) = 0 and (u− µ∇ψ) · n|∂Ω = 0. By Theorem3.6, 2o) from [54], there exists w ∈ H0(curl) such that u − µ∇ψ = ∇×w. Thisgives the ﬁrst decomposition. Using Lemma 2.2 from [83], proven in the case n1 > 0 asLemma IV.2 in [99], one can decompose w = w˜+ ∇ξ where w˜ ∈ H10(Ω), ξ ∈ H1(Ω)and ‖w˜‖H1(Ω) ≤ C‖∇×w˜‖. This proves Decomposition 2. For the last result, chooseψ to be the unique element of H10(Ω), satisfying (2.21) for any θ ∈ H10(Ω). Again,∇·(u− µ∇ψ) = 0 and the rest follows from the proof of Theorem 3.4 in [54].D. Finite element subspacesLet Ωh ⊆ Ω be a polygonal/polyhedral subdomain satisfying assumption (AΩ) 7.Unless stated otherwise, we assume Ωh = Ω.Let Th be a ﬁnite element mesh on Ωh. This means that Ωh is decomposed inthe non-overlapping set Th ≡ {τ} of closed “elements” τ . For each τ ∈ Th, we denoteby hτ and ρτ its diameter and the radius of the largest inscribed ball. We assumethat the mesh Th is aligned with the jumps of µ and ε and is shape regular (see [37]).Furthermore, we require that Th is locally quasi-uniform, i.e. there exists C ∈ R+such thatC ≥ hτρτ, ∀τ ∈ Th .In particular, we allow for meshes obtained by local reﬁnement.The theory presented in the dissertation is applicable to Th composed of triangles,quadrilaterals, tetrahedra and hexahedra. We assume that there exists a referenceelement τˆ such that each τ ∈ Th is obtained from τˆ by a linear, bilinear or trilineartransformation (depending on the type of the mesh). Below, we consider the case7There are simple polyhedral domains which are not Lipschitz, for example the“crossed bricks” domain shown on Figure 3.1 in [75].24of triangular or tetrahedral mesh. The extension to quadrilateral and hexahedralmeshes is routine.Let Pk(τ) be the space of polynomials on τ of degree k. We will use the followingstandard ﬁnite element spaces (see [37])Ŝh(k) = {vh ∈ L2(Ω) : vh|τ ∈ Pk(τ) , ∀τ ∈ Th} ,Sh(k) = Ŝh(k) ∩H1(Ω) , Sh,0(k) = Sh(k) ∩H10(Ω) .(2.22)For convenience, we set Ŝh = Ŝh(0), Sh = Sh(1) and Sh,0 = Sh,0(1).Remark 2.5 It is possible to consider the case where the order of the polynomialschange from element to element, see Corollary 4.5.By mapping to the reference element one can prove various inequalities asC hτ‖v‖2L2(∂τ) ≤ ‖v‖2L2(τ) + h2τ |v|2H1(τ) ∀v ∈ H1(τ) (2.23)andC hτ(‖v · n‖2L2(∂τ) + ‖v×n‖2L2(∂τ))≤ ‖v‖2L2(τ) + h2sτ |v|2Hs(τ) (2.24)for any v ∈ Hs(τ) with 1 ≥ s > 12. The last inequality follows from the existence ofbounded trace operator from Hs(τ) to L2(∂τ) and from the deﬁnition (2.10).We recall the following approximation property for u ∈ Hs(Ω):infuh∈Ŝh(k){∑τ∈Thh−2sτ ‖u− uh‖2L2(τ)}≤ C‖u‖2Hs(Ω) s ∈ [0, k + 1] , (2.25)and the existence of a stable approximation operator Ihu : L2(Ω) → Sh, such that∑τ∈Th{h−2τ ‖u− Ihu‖2L2(τ) + ‖Ihu‖2H1(τ)}≤ C‖u‖2H1(Ω) . (2.26)For (2.26), one can choose uh = Chu, the Cle´ment interpolation operator (see [38]and [54, pp. 109-111]). In this case, we additionally have Ihu : H10(Ω) → Sh,0.25We next describe the spaces of “bubble” functions associated with the faces.Denote with Fh the set of all faces of Th. Fix F ∈ Fh, and let TF be the unionof all elements τ ∈ Th which have F as a face. Let hF be the diameter of F . Bythe quasiuniformity hτ ≈ hF for any τ ∈ TF . The bubble function βF (x) associatedwith F should be in H1(Ω) with support equal to TF . In particular, βF (x) should benonzero on F and should vanish on all other faces in Fh. The simplest deﬁnition ofsuch face bubble function isβF |τ (x) = cFNF∏i=1i(x) ∀τ ∈ TF , (2.27)where NF is the number of vertices of F , {i(x)}NFi=1 are the barycentric coordinatesfor x ∈ τ corresponding to those vertices, and cF is a scaling parameter. For example,the choice cF = 2dNF guarantees that βF ≥ 0 with a maximum of 1 in the barycenterof F .We deﬁne the space of face bubble functions BFh as the linear span of {βF (x) :F ∈ Fh}. The space with zero boundary conditions, BFh,0, is deﬁned, similarly, byignoring the faces on the boundary of Ω. A typical element of BFh on a triangularmesh and the bubbles for each face of a tetrahedron are shown in Figure 2.2.Fig. 2.2. Face bubble functions: element of BFh in 2D and the bubbles for each faceof a tetrahedron in 3D.26One can construct face bubble functions of higher degree as follows: let Pk(F )be the space of polynomials of degree k on a ﬁxed face F . Let dk be the dimensionof this space and {ρjF}dkj=1 be the usual nodal basis. Each function ρF ∈ Pk(F ) canbe extended to a polynomial ρˆF of degree k on Rd by setting it to be constant in thedirection normal to F . The basis bubble functions are deﬁned byβjF∣∣τ(x) = cF ρjF (x)NF∏i=1i(x) ∀τ ∈ TF , (2.28)for each 1 ≤ j ≤ dk. The linear span of all these functions form the space BkFh . Thespace BkFh,0 is deﬁned, similarly, using only the interior faces.We next describe the spaces of bubble functions associated with the elements.For τ ∈ Th, the bubble function βτ (x) is in H1(Ω) with support equal to τ . Inparticular, βτ (x) should be nonzero on τ and should vanish on all other elements.The simplest deﬁnition isβτ (x) = cτNτ∏i=1i(x) ∀x ∈ τ , (2.29)where Nτ is the number of vertices of τ , {i(x)}Nτi=1 are the barycentric coordinatesfor x ∈ τ , and cτ = 2dNτ is a scaling factor which guarantees that βτ ≥ 0 with amaximum of 1 in the barycenter of τ . The space of element bubble functions BTh , isdeﬁned as the linear span of {βτ (x) : τ ∈ Th}. We note that the restriction of aface bubble function βF to F gives the element bubble function for F . One can alsointroduce the space BkTh of element bubbles of order k analogous to (2.28).27CHAPTER IIIAN ABSTRACT LEAST-SQUARES METHODIn this chapter we present and analyze a least-squares method in abstract settings.The name least-squares can be attached to a variety of approaches including Galerkinleast-squares, stabilized mixed methods and discrete least-squares in which discretiza-tion is performed before the formulation of the least-squares functional, see [34]. How-ever, in this dissertation, we will consider only the standard least-squares approachin which one minimizes a quadratic functional based on some a priori estimate.Methods of this type have been extensively developed and analyzed in recentyears. They have been applied to a variety of problems ranging from standard second-order elliptic equations to ﬁrst-order systems, elasticity, Stokes and Navier-Stokesequations, hyperbolic problems and electromagnetics. Some of the advantages of theleast-squares methods are that they always result in a symmetric and positive deﬁnitediscrete problem, and the essential boundary conditions can be weakly imposed. Weare interested in methdos for which optimal order error estimates can be derived,even if the solutions has low regularity.Least-squares method, where the functional involves only ‖·‖2L2 terms, have beenwell-known and often applied in the engineering community, see [58, 96]. We referto this variant of the method as L2-based. Recent trends in the area have been therecasting of the initial problem into ﬁrst-order systems (FOSLS method) and the useof dual norms in the functional (negative-norm least-squares). Below, we commenton some of these approaches.The naive application of L2-based least-squares to a second-order problem has thedrawbacks of higher requirements on the smoothness of the solution, which does notallow the use of standard ﬁnite element spaces. Additionally, the condition number28of the discrete system is the square of the corresponding system obtained by theGalerkin method.The FOSLS method overcomes this diﬃculty by introducing physically meaning-ful, new dependent variables. Usually, this has to be complemented with additionalcompatability equations. This method has the advantage that it can be implementedin a two-stage scheme where one sequentially minimizes the terms corresponding todiﬀerent unknowns. Additionally, the functional is usually local and therefore can beused for a posteriori error estimation.The consideration of the negative-norm least-squares methods was made possibleby the advances in the multilevel preconditioning theory for second-order problems.The paper [27], for example, constructs eﬃciently computable discrete norms equiv-alent to the norm on Hs(Ω) for |s| < 32.Next, we present the abstract approach, which is convenient for the subsequentdevelopment of the least-squares methods in the next chapters. Here, we will provideonly a few examples to illustrate the theory. For speciﬁc applications we refer to[31, 21, 22, 33, 23, 25, 81, 70, 28] as well as to the survey [14] and the referencestherein.A. Operator equationsLet X and Y be two Hilbert spaces. In our theory, it will be natural to consideroperators A ∈ L(X,Y∗). In this case, the operator A∗ ∈ L(Y,X∗) is uniquely deﬁnedby the equality〈A∗y, x〉X∗×X = 〈Ax,y〉Y∗×Y ∀x ∈ X ,y ∈ Y . (3.1)29Introduce the operators A ∈ L(X,Y) and A∗ ∈ L(Y,X), byA = TYA , A∗ = TXA∗ . (3.2)Then(Ax,y)Y = (x,A∗y)X ∀x ∈ X ,y ∈ Y . (3.3)Note that ‖A‖X→Y∗ = ‖A‖X→Y, and the following diagrams commuteXA YAX∗T−1XTXY∗TYT−1Y,X ﬀA∗YA∗X∗T−1XTXY∗TYT−1Y.For a given b ∈ Y∗, A = 0, we consider the problem: Find x ∈ X such thatA x = b . (3.4)This is the same asA x = b , (3.5)where b = TYb. Clearly (3.4) has a solution if and only if b ∈ R(A). The solution isunique if and only if N(A) = {0}.Assume that the operator is bounded from below, i.e. there exists C1 > 0 suchthatC1 ‖x‖X ≤ ‖Ax‖Y∗ = ‖Ax‖Y ∀x ∈ X . (3.6)When X = Y, this is satisﬁed, for example, if the operator is strongly monotone, i.e.C1 ‖x‖2X ≤ 〈Ax, x〉 = (Ax, x) ∀x ∈ X .The condition (3.6) means that ‖Ax‖Y∗ is a norm on X, equivalent to ‖x‖X. In30particular, N(A) = {0} and R(A) is closed1. Therefore, (3.4) has a unique solution ifand only if b is orthogonal to R(A)⊥Y∗ . We summarize this in the following result.Proposition 3.1 Assume (3.6). The problem (3.4) has a solution if and only if thedata b satisfy the compatability condition〈b,y〉 = 0 ∀y ∈ N(A∗) . (3.7)If it exists, the solution is unique and satisﬁes: C1 ‖x‖X ≤ ‖b‖Y∗ ≤ ‖A‖ ‖x‖.Proof By (3.1), TYy ∈ N(A∗)⇔ y ∈ R(A)⊥Y∗ .When the compatability condition is not satisﬁed, one can still try to solve aproblem that is naturally related to, but weaker than, (3.4). The least-squares ideais to consider the functional F : X → R, deﬁned byF(x) = ‖A x− b‖2Y∗ = ‖A x− b‖2Y , (3.8)and replace (3.4) by the problem: Find x ∈ X such thatF(x) = miny∈XF(y) . (3.9)This is appealing, in particular, because it provides a minimization principle forproblems that may not have naturally associated optimization form.The functional F(·) is convex, and its Frechet derivative is〈F′(x),h〉 = limt→0F(x+ th)− F(x)t= 2{(A x− b,Ah)Y∗} ∀h ∈ X .1If A is bounded from below, then it is injective. The converse is true only inﬁnite dimensional spaces. Indeed, take an inﬁnite dimensional space X and let Y beX equipped with any non-equivalent norm ‖ · ‖Y ‖ · ‖X. Then, the identity operatorfrom X to Y is injective, but not bounded from below.31Therefore, x ∈ X is a solution of (3.9), if and only if(A x,Ah)Y∗ = (b,Ah)Y∗ ∀ h ∈ X . (3.10)This is equivalent to: Find x ∈ X such thatA∗A x = A∗b . (3.11)Introduce the subspacesX0 = N(A) ,X1 = X0⊥ = R(A∗) ,Y0 = N(A∗) ,Y1 = Y0⊥ = R(A) ,and let QY0 : Y → Y0 denote the Y-orthogonal projection onto N(A∗). Consider thefollowing problem:A x = (I− QY0)b . (3.12)Note that the right-hand side of (3.12) is probably the most natural way to obtaincompatible data from any b ∈ Y.Proposition 3.2 Assume (3.6). Then the problems (3.9), (3.10), (3.11) and (3.12)are equivalent and have a unique solution (for any data b). The solution satisﬁes thestability estimateC1 ‖x‖X ≤ ‖b‖Y∗ (3.13)If b satisﬁes the compatability condition (3.7), then the solution of these problemscoincides with the solution of (3.4).Proof Condition (3.6) implies that R(A) is closed in Y∗. Now (3.9) has a uniquesolution by the uniqueness of orthogonal projection onto a closed subspace. Theequivalence of (3.11) and (3.12) follows from the fact that I − QY0 is the orthogonalprojection onto R(A). This, together with (3.6), implies the estimate (3.13).The above proof uses essentially the fact that R(A) is closed. Next, we investigate32when this is true. To that end, instead of (3.6), we consider the weaker condition:there exists C2 > 0 such thatC2 ‖x‖X ≤ ‖Ax‖Y∗ = ‖Ax‖Y ∀x ∈ X1 , (3.14)or equivalentlyC2 dist(x,X0) ≤ ‖Ax‖Y∗ ∀x ∈ X .This is a natural condition which, as will follow from the next result, holds e.g.if either X or Y is ﬁnite dimensional.Proposition 3.3 The condition (3.14) holds if and only if R(A) is closed. It isfurthermore equivalent toC2 ‖y‖Y ≤ ‖A∗y‖X∗ = ‖A∗y‖X ∀y ∈ Y1 . (3.15)Proof Indeed, (3.14) implies that every Cauchy sequence in R(A) is Cauchy in X1,and therefore R(A) is closed by the continuity of A. On the other hand, if R(A) isclosed, A : X1 → R(A) is a bijective linear operator and by the Banach ContinuousInverse Theorem, its inverse is bounded. This is precisely (3.14).To ﬁnish the proof it is enough to show that (3.14) implies (3.15). Indeed, assume(3.14). Then y ∈ Y1 implies that y = Az, for some z ∈ X1. Furthermore, for anyx ∈ X‖A∗Ax‖X = suph∈X\{0}(A∗Ax,h)X‖h‖X = suph∈X1\{0}(Ax,Ah)Y‖h‖X ≥ suph∈X1\{0}(Ax,Ah)YC−12 ‖Ah‖Y= C2 ‖Ax‖Y ,which proves (3.15). The proof in the other direction is analogous.Corollary 3.1 Condition (3.6) holds if and only if R(A∗) = X.Another corollary is that (3.14) implies C22 ‖y‖Y ≤ ‖AA∗y‖Y, which is a motiva-33tion to consider the following problem: Find y ∈ Y1 such thatAA∗ y = b , x = A∗y . (3.16)The equation for y has a solution if and only if the compatability condition (3.7)holds. Therefore (3.16) is equivalent to (3.4).Consider also the related problem: y ∈ Y1,(A∗ y,A∗ h)X∗ = 〈b,h〉Y∗ ∀ h ∈ Y1 , x = A∗y . (3.17)Assume (3.15). Then this problem will always have a unique solution, and for com-patible data it is the same as (3.16). The equations (3.16) and (3.17) are the duals of(3.11) and (3.10), respectively. They can be used to devise least-square type methods(see the FOSLL* method in [33]).Theorem 3.1 Assume (3.14). Then the problems (3.9), (3.10), (3.11), (3.12) and(3.17) are equivalent. Without any restrictions on the data, each of them has exactlyone minimal-norm solution (i.e. the functional ‖ · ‖X achieves a unique minimum onthe set of solutions, or equivalently x ∈ X/N(A) ≈ X1). For this solution we havethe estimate (3.13). All solutions are obtained from this one by adding an arbitraryelement of N(A). If b ∈ Y1, these solutions coincide with the solutions of the problems(3.4) and (3.16).Proof Replace X by X1 and apply Proposition 3.2.Next, we give a series of examples of applications of the least-squares method-ology for approximation (or solution) of (3.4). We start with an algorithm proposedby Gauss, which gives the origin of the name “least-squares.”1. Gauss’ least-squares method.34In this case t ∈ Rn is a ﬁxed vector of n ≥ 2 distinct numbers corresponding toobservations of the quantity b ∈ Rn. We set X = R2, Y = Rn with the Euclidianinner products. The operator A is deﬁned as A(λµ)= λ t + µ e, where e ∈ Rnwith ei = 1, i = 1, . . . , n. Since X1 = {0}, the condition (3.6) holds. Thus, byProposition 3.2, the problem (3.9) has a unique solution corresponding to theline λ t + µ which is “closest” to the points (ti,bi) in the sense thatn∑i=1|λ ti + µ− bi|2 → min .Note that in this case the compatability condition (3.7) is not likely to besatisﬁed. Finally, the solution of (3.11) can be computed eﬃciently, sinceA∗A(λµ)=((t, t) (t, e)(e, t) (e, e))(λµ).2. Systems of linear equations.Let X = Kn, Y = Km with ﬁxed bases. Let A be given by a m×n matrix andb ∈ Y. The original problem A x = b has a unique solution for any right-handside, if and only if X0 = {0} and Y0 = {0}, i.e. if A has full column and rowrank. On the other hand, by Theorem 3.1, the problem A∗ A x = A∗ b alwayshave a (unique minimum-norm) solution. This solution will be unique if A hasfull column rank.In numerical computations, A is often large, sparse and invertible. The stabilityof the iterative procedures for solving A x = b depends on the condition numberof A, deﬁned byκ(A) = ‖A‖ ‖A−1‖ . (3.18)For a nonsymmetric matrix, it might be tempting to use the least-squaresmethod in order to obtain a symmetric and positive deﬁnite problem. However,35this should be avoided since it leads to the eﬀective squaring of the conditionnumber: κ(A∗ A) = κ(A)2.3. Dirichlet problem, posed in L2(Ω).For this example, assume that Ω has full elliptic regularity, i.e. the operator∆ : H2(Ω) ∩ H10(Ω) → L2(Ω), is an isomorphism. Set X = H2(Ω) ∩ H10(Ω) andY = L2(Ω). Then A ∈ L(X,Y∗) satisﬁes the requirements of Proposition 3.2.Fix f ∈ L2(Ω). Then the least-squares problem: Find y ∈ H2(Ω) ∩H10(Ω), suchthat(∆x,∆y)L2(Ω) = (f,∆y)L2(Ω) ∀y ∈ H2(Ω) ∩H10(Ω) ,has a unique solution which satisﬁes ∆x = f.We now turn to the drawbacks of this approach as a method for solving theDirichlet problem. First, we require full regularity. This holds only in somelimited cases, and there are many alternative methods that do not require it.Moreover, the smoothness requirement does not allow for the use of standardﬁnite element spaces which are not in H2(Ω). Finally, the discretization of thebilinear form (∆·,∆·)L2(Ω) leads to a matrix with a signiﬁcantly worse conditionnumber compared to the usual Galerkin method.4. FOSLS for the Dirichlet problem.Instead of the second-order problem −∆x = f, consider the equivalent ﬁrst-order system ∇x = u, ∇·u = −f. Let A be an operator that maps (x,u) to(∇·u,∇x− u).It is proven in [31] that A is bounded and bounded from below as an operatorfrom H10(Ω)×H(div) to L2(Ω)× L2(Ω).Even though this allows for the development of a least-squares method that36avoids the discrete inf-sup condition, it leads to error estimates optimal in thenorm on H10(Ω)×H(div) but not with respect to the regularity of the solution.This result was further improved in [21], where it was proven A is bounded andbounded from below as an operator from H10(Ω)×H(div) to L2(Ω)×H−1(Ω).These estimates use the ‖ · ‖L2(Ω) norm on H(div) and lead to quasi-optimalerror estimates.Note that for the discretization of both of these methods, one can not em-ploy the standard piecewise linear ﬁnite element spaces but needs to work withapproximation spaces for H(div).5. Dirichlet problem, posed in H−1(Ω).Let X = Y = H10(Ω), with inner product (∇x,∇y)L2(Ω).Deﬁne A ≡ −∆ : H10(Ω) → H−1(Ω) by〈Ax,y〉 = (∇x,∇y)L2(Ω) ∀y ∈ H10(Ω) .Clearly A ∈ L(X,Y∗), and Poincare´’s inequality implies that condition (3.6)holds. Fix f ∈ H10(Ω), and consider the Dirichlet problem −∆x = f. ByProposition 3.2, the least-squares formulation of this problem has a uniquesolution. It is also a solution of the original problem since Y0 = {0}. Consider(3.10) and note that TYAx = x for any x ∈ X. Therefore we get(∇x,∇h)L2(Ω) = 〈f,h〉 ∀h ∈ H10(Ω) .Thus, in this case, the least-squares method reduces to the standard Galerkinweak formulation.Let us note that compared to the L2(Ω)-based algorithm, this method does37not posses any of the aforementioned deﬁciencies. This is because the a prioriestimate used is the natural stability result for the problem. A similar idea canbe applied to a much more general second-order elliptic operator, as done in[22].As it is evident from the last few examples, the least-squares approach is not astrictly deﬁned method, but rather, a general methodology which produces diﬀerentmethods depending on interpretation of the original problem.B. ApproximationThe operator A ∈ L(X,Y∗) is closely related to a bounded bilinear form on X×Ydeﬁned bya(x,y) = 〈Ax,y〉 ∀x ∈ X ,y ∈ Y .In this notation, the inf-sup condition (2.5) is the same as (3.6), and the problem(2.6) coincides with (3.4). Furthermore, the result of the Lax-Milgram Theorem 2.1is identical to Proposition 3.1. Finally, (3.11) can be rewritten asa(x,y) = 〈b,y〉 ∀y ∈ Y1 ≡ {y = TYAh : h ∈ X} . (3.19)Let Xh ⊂ X, Yh ⊂ Y be a family of ﬁnite-dimensional subspaces with the inheritedinner products. We will refer to the original problem and spaces as continuous and tothe problem and spaces depending on h as discrete. We assume that Xh approximatesX as h → 0. To make that statement precise, let QXh : X → Xh be the X-orthogonalprojection onto Xh. Furthermore, let Xˆ ⊂ X be another Hilbert space, continuouslyembedded in X, i.e.‖x‖X ≤ Cˆ1 ‖x‖Xˆ ∀x ∈ Xˆ . (3.20)38We assume that there exists a function χ(h) with limh→0 χ(h) = 0, such that‖I− QXh‖Xˆ→X ≤ χ(h) . (3.21)Our goal is to construct a discrete operator Ah : Xh → Y∗h and a discrete problemAhxh = bh that approximates (3.4). Similar to (3.6), it is natural to require thatC3 ‖x‖Xh ≤ ‖Ahx‖Y∗h ∀x ∈ Xh , (3.22)with C3 and ‖Ah‖ independent of h.A straightforward way to deﬁne Ah is by〈Ahx,y〉 = a(x,y) = 〈Ax,y〉 ∀x ∈ Xh ,y ∈ Yh . (3.23)In this case, we get a discrete problem similar to (3.19): x ∈ Xh satisﬁesa(x,y) = 〈b,y〉 ∀y ∈ Yh,1 ≡ {y = TYh Ah h : h ∈ Xh} . (3.24)Note that this is a Petrov-Galerkin approximation, as opposed to the standard Galerkinmethod, where the test and the solution spaces are the same.Furthermore, ‖Ah‖ ≤ ‖A‖, and the condition (3.22), is equivalent to the factthat (Xh,Yh) satisfy the discrete inf-sup conditionC3 ‖x‖X ≤ supy∈Yh\{0}a(x,y)‖y‖Y , ∀x ∈ Xh . (3.25)Often in the least-squares theory, a diﬀerent approach is preferred which avoidsthe discrete inf-sup condition. The idea is to start from (3.6) and use integration byparts over each element. Then Ah corresponds to a form ah(·, ·) that includes theresulting jump terms over the elements’ boundaries. To include this possibility, wemake the following assumption: There exists a function α(h) with limh→0 α(h) = 0,39such that‖A−Ah‖Xh→Y∗h ≤ α(h) . (3.26)In particular, if (3.25) holds, we can deﬁne Ah by (3.23) and set α(h) ≡ 0.Next, we consider an approximation to (3.4). We start with the case when theoriginal problem is well-posed. Then, on the discrete level, we use the least-squaresmethod with the same right-hand side. This is natural because b is not likely tosatisfy the discrete compatability condition (even though it satisﬁes the continuousone).Theorem 3.2 Suppose that (3.6), (3.22) and (3.26) hold. Let b ∈ Y∗ satisfy thecompatability conditions (3.7) and x ∈ X be the unique solution of the problem Ax = b.Let xh be the unique solution of the least-squares method for the equation Ahxh = b.ThenC ‖x− xh‖X ≤ ‖(I− QXh) x‖X + α(h) ‖x‖X . (3.27)Proof The approximation xh ∈ Xh satisﬁes(Ahxh,Ahζh)Y∗h = (Ax,Ahζh)Y∗h ∀ζh ∈ Xh .Fix ζh ∈ Xh. The above, together with (3.25) and (3.23), implyC23 ‖xh− ζh‖2X ≤ (Ah(xh− ζh),Ah(xh− ζh))Y∗h= 〈A(x− ζh),TYhAh(xh− ζh)〉Y∗×Yh +〈Aζh −Ahζh,TYhAh(xh− ζh)〉Y∗h×Yh≤ ‖A‖ ‖Ah‖ ‖x− ζh‖X‖xh− ζh‖X+‖A−Ah‖ ‖Ah‖ ‖ζh‖X‖xh− ζh‖X .40Let C˜ = C−23 ‖Ah‖. By the triangle inequality,‖x− xh‖X ≤ (1 + C˜ ‖A‖) ‖x− ζh‖X + C˜ ‖A−Ah‖ ‖ζh‖X ∀ζh ∈ Xh .The result follows by setting ζh = QXhx.Corollary 3.2 If, additionally, x ∈ Xˆ \ {0}, thenlimh→0‖x− xh‖X‖x‖Xˆ≤ limh→0{(1 + C˜ ‖A‖)χ(h) + C˜ Cˆ1 α(h)}= 0 .Corollary 3.3 Suppose that (2.5) and (3.25) hold. Let b, x and xh be as in thetheorem. Then the least-squares approximation is quasi-optimal, i.e.‖x− xh‖X ≤(1 +‖A‖2C23)infζh∈Xh‖x− ζh‖X . (3.28)Corollary 3.4 Replace (3.26) by‖(A−Ah)QXh‖Xˆ→Y∗h ≤ α(h) . (3.29)Repeating the proof of the theorem for x ∈ Xˆ we getC ‖x− xh‖X ≤ ‖(I− QXh) x‖X + α(h) ‖x‖Xˆ .Next, we consider the case of arbitrary right-hand side, i.e. we have no compata-bility conditions on b. This forces us to use the least-squares method on both thecontinuous and discrete levels.For an operator Q ∈ L(Y,Y), we introduce the notation Q˜ for the operator aselement of L(Y∗,Y∗), i.e. we setQ˜ = T−1Y QTY . (3.30)By (3.12), the least-squares solution operator S : Y∗ → X is deﬁned bySb = x , where A x = (I− Q˜Y0)b . (3.31)41Up to this point there were no requirements for the approximation propertiesof the spaces Yh. Now we assume the following: There exists a sequence of closedsubspaces Yh,0 ⊂ Yh ∩ Y0, such that the Y-orthogonal projectors QYh,0 : Y → Yh,0 areapproximations of QY0 . Speciﬁcally, let Yˆ ⊃ Y be another Hilbert space, such that Yis dense and continuously embedded in it:‖y‖Yˆ ≤ Cˆ2 ‖y‖Y ∀y ∈ Y . (3.32)We assume that there exists a function γ(h) with limh→0 γ(h) = 0, and such that‖Q˜Y0−Q˜Yh,0‖Yˆ∗→Y∗ ≤ γ(h) . (3.33)Since the solution of (3.31) does not change if b is perturbed by an element ofY0, we will need to use a modiﬁed right-hand side in the deﬁnition of the discreteleast-squares solution operator. Speciﬁcally, Sh : Y∗ → Xh is deﬁned byShb = xh , where (Ah xh,Ah ζh)Y∗h = ((I− Q˜Yh,0)b,Ah ζh)Y∗h ∀ ζh ∈ Xh . (3.34)Note that this is the standard least-squares method (3.10), but applied for the right-hand side bh = (I− Q˜Yh,0)b ∈ Y∗ ⊂ Y∗h, instead of b.Theorem 3.3 Assume (3.6), (3.22), (3.26), and the following additional condition:C4 ‖x‖Xˆ ≤ ‖Ax‖Yˆ∗ ∀x ∈ Xˆ . (3.35)Assume also that the spaces (Xh,Yh) approximate (X,Y) in the sense of (3.21) and(3.33). Then we have the estimate‖S− Sh‖Yˆ∗→X ≤ C−14 (1 + C˜ ‖A‖)χ(h) + C˜ γ(h) + C˜ Cˆ1 α(h) , (3.36)where C˜ ≤ C−23 (1+α(h)) ‖A‖. In particular, the least-squares method (3.31) provides42a uniform approximation to (3.34) for any b ∈ Yˆ∗.Proof As in the proof of Theorem 3.2, we haveC23 ‖xh− ζh‖2X ≤ (Ah(xh− ζh),Ah(xh− ζh))Y∗h= 〈(I− Q˜Yh,0)b−Ah ζh,TYhAh(xh− ζh)〉Y∗×Yh= 〈A(x− ζh),TYhAh(xh− ζh)〉Y∗×Yh+ 〈(Q˜Y0−Q˜Yh,0)b,TYhAh(xh− ζh)〉Y∗×Yh+ 〈Aζh−Ahζh,TYhAh(xh− ζh)〉Y∗h×YhLet C˜ = C−23 ‖Ah‖, then‖x− xh‖X ≤ (1 + C˜ ‖A‖)‖x− ζh‖X + C˜ ‖(Q˜Y0−Q˜Yh,0)b‖Y∗ + C˜ ‖A−Ah‖ ‖ζh‖X ,for any ζh ∈ Xh. The result follows by combining (3.33), (3.21) and (3.13).Let Yˆ0 and Yˆh,0 be the closures of Y0 and Yh,0 in Yˆ. Denote with QYˆ0 and QYˆh,0 theYˆ-orthogonal projectors onto these subspaces. Furthermore, deﬁne Q˜Yˆ0 = T−1YˆQYˆ0 TYˆand Q˜Yˆh,0 = T−1YˆQYˆh,0 TYˆ. Next, we consider the case when Q˜Y0 and Q˜Yh,0 are replacedby Q˜Yˆ0 and Q˜Yˆh,0 . This is of interest because the projections in the weaker innerproduct might be easier to implement.Corollary 3.5 Consider the operators S : Yˆ∗ → X deﬁned bySˆb = x , where A x = (I− Q˜Yˆ0)b , (3.37)Sh : Yˆ∗ → Xh deﬁned bySˆhb = xh , where (Ah xh,Ah ζh)Y∗h = ((I− Q˜Yˆh,0)b,Ah ζh)Y∗h ∀ ζh ∈ Xh , (3.38)and the condition‖Q˜Yˆ0−Q˜Yˆh,0‖Yˆ∗→Y∗ ≤ γ(h) . (3.39)43Let b ∈ Yˆ∗, then1. Sˆb = Sb, i.e. the problems (3.31) and (3.37) are equivalent.2. The Theorem 3.3 holds for Sˆ and Sˆh with (3.33) replaced by (3.39).Proof For b ∈ Yˆ∗ we have(TYb,y)Y = 〈b,y〉 = (TYˆb,y)Yˆ ∀y ∈ Y .This implies〈Q˜Y0 b,y〉 = 〈Q˜Yˆ0 b,y〉 ∀y ∈ Y0 .In particular, the problems (3.31) and (3.37) have the same solution. The proof ofthe theorem proceeds exactly as before.To summarize the results from this section: under appropriate conditions on theapproximation space Xh and the discrete operator Ah, the least-squares approxima-tion converges to the solution of the original problem when it is unique. In general,when the solution is not unique, under further conditions on the approximation spaceYh and the operator Ah, we have convergence of the continuous least-square solutionoperator to a discrete least-square solution operator. These results will be applied inthe convergence theory of the next chapters.C. ImplementationWe next consider the implementation of the discrete least-squares method (3.19).Since the evaluation of the operator TYh involves the solution of a linear system, weﬁrst replace it with a spectrally equivalent preconditioner. Speciﬁcally, let T̂Yh ∈L(Y∗h,Yh) and there are constants Ĉ1 ≥ Ĉ0 > 0, independent of h, such thatĈ0 〈y, T̂Yhy〉 ≤ ‖y‖2Y∗h ≤ Ĉ1 〈y, T̂Yhy〉 ∀y ∈ Y∗h .44Furthermore, assume that T̂Yh is symmetric, i.e. 〈x, T̂Yhy〉 = 〈y, T̂Yhx〉 for all x ,y ∈Yh∗. We are interested in solving the problem: Find xh ∈ Xh satisfying〈Ahxh, T̂YhAhyh〉 = 〈b, T̂YhAhyh〉 ∀yh ∈ Xh . (3.40)Here 〈·, ·〉 denotes the duality pairing on Y∗h×Yh. The next result shows that themethod is stable with respect to such a uniform perturbation of the form.Proposition 3.4 Let xh be the solution of (3.40). Then under the additional condi-tions in the corresponding theorems, we have1. Theorem 3.2 holds with C˜ ≤ C−23 Ĉ1 Ĉ−10 (1 + α(h)) ‖A‖.2. Theorem 3.3 holds with C˜ ≤ C−23 Ĉ1 Ĉ−10 (1+α(h)) ‖A‖, if (3.34) is replaced byShb = xh, where〈Ah xh, T̂YhAh yh〉 = 〈(I− Q˜Yh,0)b, T̂YhAh yh〉 ∀ yh ∈ Xh .Next we address the solution of the discrete problem (3.40) in the case when Ahis given by (3.23). Let {xih}ni=1 and {yjh}mj=1 be the bases of Xh and Yh. In all ourapplications those are real, piecewise polynomial ﬁnite element spaces.Every element x ∈ Xh is uniquely determined by its coordinates ∼x ∈ Rn in thebasis {xih}, i.e. x =∑i ∼xi xih. Deﬁne also the vector of dual coordinates∼x ∈ Rn by∼xi = (x, xih)Xh . The vectors ∼y ,∼y ∈ Rn are similarly deﬁned for any y ∈ Yh.Introduce the matrix≈A and the vector∼b by≈Aji = a(xih,yjh),∼bj = 〈b,yjh〉. Then,for example, the problem Ahx = b corresponds to the linear system≈A∼x =∼b. Indeed,〈Ahx,y〉 = ∼yt≈A∼x = ∼yt∼b = 〈b,y〉, for any y ∈ Yh. Similarly, the problem A∗hx = bcorresponds to the linear system≈At ∼x =∼b. This can be summarized asA˜hx =≈A∼x , A˜∗hx =≈At ∼x . (3.41)45The operator T̂Yh takes a dual vector and produces a vector of coordinates. Thisis the typical setup, for example, when T̂Yh is a preconditioner, such as Multigrid.Let ≈T be the matrix corresponding to the action of T̂Yh , i.e. x = T̂Yhb if and only if∼x = ≈T∼b. In other wordsT̂Yhb˜= ≈T∼b . (3.42)By the symmetry of T̂Yh , the problem (3.40) is equivalent toA∗h T̂Yh Ahxh = A∗h T̂Yh b .By (3.41) and (3.42) this reduces to the following linear problem:≈At ≈T≈A∼x =≈At ≈T∼b . (3.43)The matrix of this system is full and should not be assembled. Instead, we solve(3.43) by a preconditioned iterative method for a symmetric and positive deﬁnitematrix. These methods are very well understood, and the preconditioned conjugategradient (PCG) is a popular choice. To implement such a method, we only need tocompute the action of the matrix and that of a preconditioner.Since the form 〈A∗h T̂Yh Ah·, ·〉 is uniformly equivalent to ‖ · ‖2Xh , the conditionnumber of the matrix≈At ≈T≈A is of the same order as the condition number of themass matrix for Xh. In the applications considered in this dissertation, this massmatrix is well conditioned, and therefore, there is no need for a better preconditionerthan a simple diagonal scaling. In general, a necessary and suﬃcient condition for auniform (independent of h) convergence is to choose a preconditioner T̂Xh ∈ L(X∗h,Xh)satisfyingĈ2 〈x, T̂Xhx〉 ≤ ‖x‖2X∗h ≤ Ĉ3 〈x, T̂Xhx〉 ∀x ∈ X∗h ,where Ĉ3 ≥ Ĉ2 > 0 are constants independent of h.46CHAPTER IVTHE MAGNETOSTATIC AND THE ELECTROSTATIC PROBLEMSIn this chapter we consider the generalized magnetostatic and electrostatic problems(1.5) and (1.6). These problems have been studied by many authors as indicatedby the detailed literature review in [26]. Here, we just point out that L2(Ω)-basedleast-squares discretization was considered in [36], and there are other studies basedin L2(Ω), e.g. [6]. These methods have well-known drawbacks, such as requirementfor smoothness of the boundary and restriction to two dimensional problems. Amongthe more standard approaches are the introduction of a scalar or vector potential,complemented with a “gauge” condition, and the mixed ﬁnite element methods basedon the Ne´de´lec approximation spaces. These methods can also be problematic. Inparticular, the implementation of Ne´de´lec elements, especially those of higher order, isquite complicated. The solution methods for the resulting algebraic systems have onlybeen recently developed. Finally, let us note that many authors have demonstrated, cf.[18, 42], that the straightforward application of standard node-based ﬁnite elementscan lead to spurious discrete solutions.In the subsequent development, we will always assume that (AΩ) and (Aµ,ε)hold. By Theorem 2.5, see also Remarks 4.1 and 4.3 below, it is enough to considerhomogeneous boundary conditions, in which case, we are looking for the magneticand electric ﬁelds h ,e : Ω→ C3 satisfying⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×h = j in Ω,∇ · (µh) = ρ in Ω,µh · n = 0 on ∂Ω,and⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×e = j in Ω,∇ · (εe) = ρ in Ω,e×n = 0 on ∂Ω.(4.1)The above systems diﬀer essentially only in the boundary conditions, and we often47refer to them together as div-curl systems. We distinguish between the two problemsby the use of a subscript k which equals 1 for the magnetostatic problem and is 2 forthe electrostatic problem. In particular, we call the problem for h, div-curl system oftype 1 and the problem for e, div-curl system of type 2.The standard interpretation of (4.1) is to assume that j ∈ L2(Ω), ρ ∈ L2(Ω) and tosolve for h ∈ X1(µ) and e ∈ X2(ε). We call this the original form of the magnetostaticand electrostatic problems. The next sections will be devoted to weaker formulationsthat allow us to consider solutions with much lower regularity.Let us note that it is enough to devise the theory for real-valued ﬁelds, since theproblems with complex ﬁelds can be split into problems for their real and imaginaryparts. In this case all the spaces are real and there is no use of complex arithmeticin the implementation.A. Weak formulation of the magnetostatic problemIn this section we assume that (AΩ) is satisﬁed with n2 = 0, i.e., Ω is simply connected.Let h ∈ X1(µ) satisfy the magnetostatic problem (4.1). Integration by parts (seeTheorem 2.5) implies that the problem is equivalent to⎧⎪⎨⎪⎩(h,∇×φ)L2(Ω) = 〈j,φ〉 ∀φ ∈ D(Ω) ,(µh,∇ψ)L2(Ω) = −〈ρ, ψ〉 ∀ψ ∈ D(Ω) .By density (Theorems 2.4 and 2.5), this is the same as⎧⎪⎨⎪⎩(h,∇×φ)L2(Ω) = 〈j,φ〉 ∀φ ∈ H0(curl) ,(µh,∇ψ)L2(Ω) = −〈ρ, ψ〉 ∀ψ ∈ H1(Ω) .Introduce the operators curl1 : L2(Ω) → H0(curl)∗ and div1,µ : L2(Ω) → H1(Ω)∗48deﬁned by〈curl1h,φ〉 = (h,∇×φ) ∀h ∈ L2(Ω), φ ∈ H0(curl),〈div1,µh, ψ〉 = −(µh,∇ψ) ∀h ∈ L2(Ω), ψ ∈ H1(Ω) .(4.2)Note that if h ∈ L2(Ω) then curl1(h) = curl(h), and for h ∈ H(div;µ) we havediv1,µ(h) = div(µh)− γn(h).Consider the mapping A1 : h → (curl1h,div1,µh). We can summarize that theoriginal magnetostatic problem is equivalent toA1h = f , f = (j, ρ) , (4.3)where A1 is considered as an operator from X1(µ) to L2(Ω)× L2(Ω).In order to allow for solutions with very low regularity, we want to replace therequirement h ∈ X1(µ) with h ∈ L2(Ω). A natural way to do that is to work withA1 as an operator from L2(Ω) to H0(curl)∗ × H1(Ω)∗. We claim that, in this case,the problem (4.3) will have a unique solution, provided the right-hand side satisﬁescertain compatibility conditions. Indeed, we only need to show that A1 is boundedfrom below. Recall that by Corollary 3.1 this is equivalent to R(A∗1) = L2(Ω). In ourcase, A∗1 : H0(curl)×H1(Ω) → L2(Ω) is given byA∗1(φ, ψ) = ∇×φ− µ∇ψ , (4.4)and the ﬁrst Helmholtz decomposition from Theorem 2.8 implies that A∗1 is onto.The above approach is attractive, but it will lead to a discretization methodinvolving curl-conforming ﬁnite element spaces which we want to avoid. Therefore,as a compromise, we propose to introduce the spacesX1 = L2µ(Ω) , V1 = H10(Ω) , H1 = H1(Ω) , Y1 = V1 ×H1 . (4.5)49and consider (4.3) with A1 : X1 → Y∗1 . This is stronger than the previous formulationbut still much weaker than (4.1). Moreover, the second Helmholtz decompositionfrom Theorem 2.8 implies that A∗1 (which is again given by (4.4)) is onto.Remark 4.1 Before we continue with the properties of A1, let us remark that themore general problem ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×h = j in Ω,∇ · (µh) = ρ in Ω,µh · n = σ on ∂Ω,(4.6)where σ ∈ H− 12 (∂Ω) can be reduced to the same weak formulation, if we deﬁne f =(j, ρ′), where ρ′ = ρ− σ.By the theory in Chapter III, we need to consider compatibility conditions relatedto the spaceN(A∗) ≡ Y1,0 = {(w, ψ) ∈ Y1 : ∇×w− µ∇ψ = 0} .By orthogonality, the fact that ‖∇ψ‖ is an equivalent norm on H1/R and Theorem2.5, it follows that Y1,0 = V1,0 ×H1,0, whereV1,0 = {w ∈ V1 : w = ∇ψ, ψ ∈ H10(Ω)} , H1,0 = {ψ ∈ H1 : ψ = const} .(4.7)Furthermore, Theorem 2.8 implies that‖∇×w‖+ ‖∇ψ‖ is an equivalent norm on Y1/Y1,0 . (4.8)Proposition 4.1 The operator A1 : X1 → Y∗1 is linear, bounded and bounded frombelow. Speciﬁcally, there exist constants C0 and C1 independent of µ and satisfyingC0 µ˜0 ‖x‖2X1 ≤ ‖curl1x‖2V∗1 + ‖div1,µx‖2H∗1 ≤ C1 µ˜1 ‖x‖2X1, (4.9)50for all x ∈ L2(Ω), where µ˜0 = min{µ0, µ−11 } and µ˜1 = max{µ−10 , µ1}.Proof We ﬁrst prove that C ‖x‖X1 ≤ ‖A1x‖Y∗1 . As indicated before, this holdsfor some C since R(A∗1) = L2(Ω). In fact, see Proposition 3.3, the constant C ischaracterized also as the optimal constant in the inequality C ‖y‖Y⊥1,0 ≤ ‖A∗1y‖X∗1 .Let y = (w, ψ), with w ∈ V⊥1,0 and ψ ∈ H⊥1,0. Then, by (4.8)‖∇×w− µ∇ψ‖2µ−1 ≥ µ−11 ‖∇×w‖2 + µ0 ‖∇ψ‖2 ≥ C µ˜0 (‖w‖2V1 + ‖ψ‖2H1) .This proves the left inequality in (4.9). The right one follows from(x,∇×w) ≤ ‖x‖µ ‖∇×w‖µ−1 and (µ x,∇ψ) ≤ ‖x‖µ ‖∇ψ‖µ .An alternative proof, given in [26], is as follows: ﬁx x ∈ X1, and let µ x =∇×w+ µ∇ψ be the decomposition with w ∈ H10(Ω) and ψ ∈ H1(Ω)/R. Then‖x‖2µ = ‖∇×w‖2µ−1 + ‖∇ψ‖2µ =(x,∇×w)2‖∇×w‖2µ−1+(µ x,∇ψ)2‖∇ψ‖2µ.Therefore‖x‖2µ ≤ Cµ1((x,∇×w)‖w‖V1)2+ Cµ−10((µ x,∇ψ)‖ψ‖H1)2.Thus‖x‖2µ ≤ C max{µ1, µ−10 }{(supw∈V1(x,∇×w)‖w‖V1)2+(supψ∈H1(µ x,∇ψ)‖ψ‖H1)2}.Now we are in a position to characterize the solvability of the magnetostaticproblem in the proposed weak formulation, as well as in its original form. We startwith an application of Proposition 3.1.51Theorem 4.1 Let f = (j, ρ) ∈ Y∗1 . Then the problem⎧⎪⎨⎪⎩curl1h = j in V∗1 ,div1,µh = ρ in H∗1 ,(4.10)has a unique solution h ∈ X1 if and only if〈j,w〉 = 0 and 〈ρ, ψ〉 = 0 ∀w ∈ V1,0 ,∀ψ ∈ H1,0 . (4.11)The solution satisﬁes the estimate C0 µ˜0 ‖h‖2 ≤ ‖j‖2V∗1 + ‖ρ‖2H∗1 .Corollary 4.1 If h ∈ X1(µ) is such that ∇×h = 0 and ∇·µh = 0 then h = 0.For the next results, recall that X1(µ) ↪→ Hs(Ω) denotes continuous embeddingProposition 4.2 Assume that X1(µ) ↪→ Hs(Ω) for some s > 0, then there existsC = C(µ) ∈ R+, such thatC(µ) ‖h‖ ≤ ‖∇×h‖+ ‖∇·µh‖ ∀h ∈ X1(µ) . (4.12)In particular ‖∇×h‖+ ‖∇·µh‖ is an equivalent norm on X1(µ).Proof Assume the converse, then there exist a sequence {hn} ⊂ X1(µ) with ‖hn‖s =1 and ‖∇×hn‖ + ‖∇·µhn‖ ≤ 1n . Since Hs(Ω) is compactly embedded in L2(Ω), bypassing to a subsequence, we have hn → h in L2(Ω), for some h ∈ L2(Ω). Bythe continuous embedding X1(µ) ↪→ Hs(Ω), it follows that {hn} is Cauchy in Hs(Ω).Therefore, hn → h in Hs(Ω) and ‖h‖s = 1. On the other hand, ‖∇×h‖ = ‖∇·µh‖ =0 and therefore, h = 0. This is a contradiction which proves the result.Let V1,0 and H1,0 denote the closures of the spaces V1,0 and H1,0 in L2(Ω) andL2(Ω) respectively, i.e.V1,0 = {∇ψ : ψ ∈ H10(Ω)} , H1,0 = {ψ ∈ L2(Ω) : ψ = const} . (4.13)52Next we give a result for the original magnetostatic problem (4.1). The proof is basedon the estimate (4.12) and Proposition 3.1.Theorem 4.2 Let f = (j, ρ) ∈ L2(Ω)× L2(Ω). Then the problem⎧⎪⎨⎪⎩∇×h = j in L2(Ω) ,∇·µh = ρ in L2(Ω) ,(4.14)has a unique solution h ∈ X1(µ) if and only if(j,w) = 0 and (ρ, ψ) = 0 ∀w ∈ V1,0 ,∀ψ ∈ H1,0 . (4.15)The solution satisﬁes the estimate C(µ) (‖h‖2 +‖∇×h‖2 +‖∇·µh‖2) ≤ ‖j‖2 +‖ρ‖2.Remark 4.2 The conditions (4.15) are the same as ∇·j = 0 and (ρ, 1) = 0.Proposition 4.3 For f = (j, ρ) ∈ L2(Ω)×L2(Ω), the problems (4.10) and (4.14) areequivalent.Proof Clearly any solution of (4.14) satisﬁes (4.10). Furthermore, if f = (j, ρ) ∈L2(Ω)× L2(Ω), then the compatability conditions (4.11) and (4.15) are equivalent.As in Chapter III, we can deﬁne a least-squares method for problem (4.10). Theleast-squares functional isF1(h) = ‖curl1h− j‖2V∗1 + ‖div1,µh− ρ‖2H∗1 ,and its minimization over h ∈ L2(Ω) is equivalent (by Proposition 3.2) to(curl1h, curl1y)V∗1 + (div1,µh,div1,µy)H∗1 = (j, curl1y)V∗1 + (ρ,div1,µy)H∗1 ,for any y ∈ L2(Ω).Let QV1,0 and QH1,0 be the V1 and H1 orthogonal projectors onto V1,0 and H1,0,53respectively. This means that, for w ∈ V1, QV1,0w = ∇ϕ, where ϕ ∈ H20(Ω) satisﬁes(∇ϕ,∇θ)V1,0 = (w,∇θ)V1,0 ∀θ ∈ H20(Ω) .Similarly, for ψ ∈ H1, QH1,0ψ = ψ, where ψ = 1µ(Ω) (ψ, 1) denotes the mean value ofψ over Ω.Deﬁne Q˜V1,0 : V∗1 → V∗1 and Q˜H1,0 : H∗1 → H∗1 similarly to (3.30), i.e.〈Q˜V1,0j,w〉 = 〈j,QV1,0w〉 and 〈Q˜H1,0ρ, ψ〉 = 〈ρ,QH1,0ψ〉 ∀w ∈ V1 ,∀ψ ∈ H1 .Then, by Proposition 3.2, the least-squares method for (4.10) is equivalent to⎧⎪⎨⎪⎩curl1h = (I− Q˜V1,0)j in V∗1 ,div1,µh = (I− Q˜H1,0)ρ in H∗1 .(4.16)By the deﬁnition of Y1,0, this can be rewritten asa1(h; (w, ψ)) ≡ (h,∇×w) + (h, µ∇ψ) = 〈j,w〉 + 〈ρ, ψ〉 (4.17)for any (w, ψ) ∈ Y1,1, whereY1,1 = Y⊥1,0 ≡{((I− QV1,0)w, (I− QH1,0)ψ): (w, ψ) ∈ Y1}.Theorem 4.3 The problem (4.16) has a unique solution h ∈ L2(Ω) for any dataf = (j, ρ) ∈ Y∗1 . When f satisﬁes the compatability conditions (4.11), the problem isequivalent to (4.10). When f ∈ L2(Ω)×L2(Ω) and (4.15) is satisﬁed, the least-squaresproblem is equivalent to the original magnetostatic problem (4.14).Corollary 4.2 Let Q1 : V1 → V1,0 be the L2(Ω)-projection, i.e. Q1w = ∇ϕ, whereϕ ∈ H10(Ω) satisﬁes(∇ϕ,∇θ) = (w,∇θ) ∀θ ∈ H10(Ω) .54Then for any (j, ρ) ∈ L2(Ω) × L2(Ω), the problems (4.16), (4.10) and (4.14) withright-hand side f = ((I− Q1)j, ρ− ρ), have the same unique solution.B. Weak formulation of the electrostatic problemIn this section we assume that(AΩ) is satisﬁed with n1 = 0, i.e., ∂Ω is connected.We proceed analogously to the previous section. Let e ∈ X2(ε) satisfy theelectrostatic equations (4.1). Introduce the operators curl2 : L2(Ω) → H(curl)∗ anddiv2,ε : L2(Ω) → H10(Ω)∗ deﬁned by〈curl2e,φ〉 = (e,∇×φ) ∀e ∈ L2(Ω), φ ∈ H(curl),〈div2,εe, ψ〉 = −(εe,∇ψ) ∀e ∈ L2(Ω), ψ ∈ H10(Ω) .(4.18)Note that if e ∈ H0(curl) then curl2(e) = curl(e), and for e ∈ L2(Ω) we havediv2,ε(e) = div(εe).Consider the mapping A2 : e → (curl2e,div2,εe). Then the original electro-static problem is equivalent toA2e = f , f = (j, ρ) , (4.19)where A2 is considered as an operator from X2(ε) to L2(Ω)× L2(Ω).We introduce the spacesX2 = L2ε(Ω) , V2 = H1(Ω) , H2 = H10(Ω) , Y2 = V2 ×H2 . (4.20)and propose to consider (4.19) with A2 : X2 → Y∗2 . By the last Helmholtz decompo-sition from Theorem 2.8, A∗2 : H(curl)×H10(Ω) → L2(Ω) given byA∗2(φ, ψ) = ∇×φ− ε∇ψ (4.21)is onto.55Remark 4.3 The more general problem⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×e = j in Ω,∇ · (εe) = ρ in Ω,e×n = σ on ∂Ω,(4.22)where σ ∈ H− 12 (∂Ω) can be reduced to the same weak formulation, if we deﬁne f =(j′, ρ), where j′ = j− σ.The compatibility conditions are related to the space Y2,0 = V2,0 ×H2,0, whereV2,0 = {w ∈ V2 : w = ∇ψ, ψ ∈ H1(Ω)} , H2,0 = {0} . (4.23)Moreover, by Theorem 2.8,‖∇×w‖+ ‖∇ψ‖ is an equivalent norm on Y2/Y2,0 . (4.24)Proposition 4.4 The operator A2 : X2 → Y∗2 is linear, bounded and bounded frombelow. Speciﬁcally, there exist constants C0 and C1 independent of ε and satisfyingC0 ε˜0 ‖x‖2X2 ≤ ‖curl2x‖2V∗2 + ‖div2,εx‖2H∗2 ≤ C1 ε˜1 ‖x‖2X2, (4.25)for all x ∈ L2(Ω), where ε˜0 = min{ε0, ε−11 } and ε˜1 = max{ε−10 , ε1}.The proof is analogous to Proposition 4.1. As before, using Proposition 3.1, we getthe following results.Theorem 4.4 Let f = (j, ρ) ∈ Y∗2 . Then the problem⎧⎪⎨⎪⎩curl2e = j in V∗2 ,div2,εe = ρ in H∗2 ,(4.26)56has a unique solution e ∈ X2 if and only if〈j,w〉 = 0 ∀w ∈ V2,0 . (4.27)The solution satisﬁes the estimate C0 ε˜0 ‖e‖2 ≤ ‖j‖2V∗2 + ‖ρ‖2H∗2 .Corollary 4.3 If e ∈ X2(ε) is such that ∇×e = 0 and ∇·εe = 0 then e = 0.Proposition 4.5 Assume that X2(ε) ↪→ Hs(Ω) for some s > 0, then there existsC = C(ε) ∈ R+, such thatC(ε) ‖e‖ ≤ ‖∇×e‖+ ‖∇·εe‖ ∀e ∈ X2(ε) . (4.28)In particular, ‖∇×e‖+ ‖∇·εe‖ is an equivalent norm on X2(ε).Let V2,0 and H2,0 denote the closures of the spaces V2,0 and H2,0 in L2(Ω) andL2(Ω), respectively, i.e.V2,0 = {∇ψ : ψ ∈ H1(Ω)} , H2,0 = {0} . (4.29)Theorem 4.5 Let f = (j, ρ) ∈ L2(Ω)× L2(Ω). Then the problem⎧⎪⎨⎪⎩∇×e = j in L2(Ω) ,∇·εe = ρ in L2(Ω) ,(4.30)has a unique solution e ∈ X2(ε) if and only if(j,w) = 0 ∀w ∈ V2,0 . (4.31)The solution satisﬁes the estimate C(ε) (‖e‖2 + ‖∇×e‖2 + ‖∇·εe‖2) ≤ ‖j‖2 + ‖ρ‖2.Remark 4.4 The condition (4.31) is the same as ∇·j = 0 and j · n = 0.Proposition 4.6 For f = (j, ρ) ∈ L2(Ω)×L2(Ω), the problems (4.26) and (4.30) areequivalent.57We proceed with the least-squares method based on the functionalF2(e) = ‖curl2e− j‖2V∗2 + ‖div2,εe− ρ‖2H∗2 ,which is minimized over e ∈ L2(Ω). This is the same as(curl2e, curl2y)V∗2 + (div2,εe,div2,εy)H∗2 = (j, curl2y)V∗2 + (ρ,div2,εy)H∗2 ,for any y ∈ L2(Ω).Let QV2,0 and QH2,0 be the V2 and H2 orthogonal projectors onto V2,0 and H2,0respectively. This means that for w ∈ V2, QV2,0w = ∇ϕ, where ϕ ∈ H2(Ω)/Rsatisﬁes(∇ϕ,∇θ)V2,0 = (w,∇θ)V2,0 ∀θ ∈ H2(Ω)/R .Deﬁne Q˜V2,0 : V∗2 → V∗2 by〈Q˜V2,0j,w〉 = 〈j,QV2,0w〉 ∀w ∈ V2 .Then, by Proposition 3.2, the least-squares method for (4.26) is equivalent to⎧⎪⎨⎪⎩curl2e = (I− Q˜V2,0)j in V∗2 ,div2,εe = ρ in H∗2 .(4.32)By the deﬁnition of Y2,0, this can be rewritten asa2(e; (w, ψ)) ≡ (e,∇×w) + (e, ε∇ψ) = 〈j,w〉 + 〈ρ, ψ〉 (4.33)for any (w, ψ) ∈ Y2,1, whereY2,1 = Y⊥2,0 ≡{((I− QV2,0)w, ψ): (w, ψ) ∈ Y2}.Theorem 4.6 The problem (4.32) has a unique solution e ∈ L2(Ω) for any dataf = (j, ρ) ∈ Y∗2 . When f satisﬁes the compatability conditions (4.27), the problem is58equivalent to (4.26). When f ∈ L2(Ω)×L2(Ω) and (4.31) is satisﬁed, the least-squaresproblem is equivalent to the original electrostatic problem (4.30).Corollary 4.4 Let Q2 : V2 → V2,0 be the L2(Ω)-projection, i.e. Q2w = ∇ϕ, whereϕ ∈ H1(Ω)/R satisﬁes(∇ϕ,∇θ) = (w,∇θ) ∀θ ∈ H1(Ω)/R .Then for any (j, ρ) ∈ L2(Ω) × L2(Ω), the problems (4.32), (4.26) and (4.30) withright-hand side f = ((I− Q2)j, ρ), have the same unique solution.C. Least-squares approximationIn this section, we consider the approximation of the magnetostatic and electrostaticproblems based on discrete least-squares methods. We concentrate on the magneto-static problem, the results for the electrostatic problem are analogous.Following Section III.B, let Xh,1 ⊂ X1, and Yh,1 = Vh,1 × Hh,1, with Vh,1 ⊂ V1,Hh,1 ⊂ H1 be approximation subspaces. Our main example, using the deﬁnitions in§II.D, is the case Ωh = Ω with Xh,1 = Ŝh, Vh,1 = (Sh,0⊕BFh,0)3 and Hh,1 = Sh⊕BFh .For the notation of operators, spaces and functions, we follow Chapter III with anadded subscript 1.Let s ∈ [0, 1]. Then the space Xˆ1 = Hs(Ω) is continuously embedded in X1.Therefore, (3.20) holds with Cˆ1 = 1. We assume that the estimate (3.21) holds withχ(h) = C hs, i.e.infxh∈Xh,1‖x− xh‖ ≤ C hs ‖x‖s ∀x ∈ Hs(Ω) , s ∈ [0, 1] . (4.34)By (2.25), this is true for our reference choice of Xh,1.The magnetostatic problem (4.3) involves the operator A1, and the objective is59to replace it with a discrete approximation Ah,1 : Xh,1 → Y∗h,1. To insure stability, wewill require that Ah,1 is bounded from below, i.e. (3.22) holds. Next, we consider twodiﬀerent choices for the discrete operator Ah,1, which satisfy this requirement. Theresulting discretizations will be shown to be stable in L2(Ω) and to yield ﬁrst-orderconvergence when the solution is in H1(Ω). By interpolation, we have hs convergencewhen the solution is in Hs(Ω) for any s in [0, 1].1. Approximation based on a discrete inf-sup conditionRecall that the operator A1 induces the bilinear form a1(·, ·) deﬁned in (4.17). In thissubsection we assume that the spaces Xh,1 and Yh,1 are chosen appropriately, so thatthe discrete inf-sup condition (3.25) holds.C ‖h‖X ≤ supw∈Vh,1(h,∇×w)‖w‖Vh,1+ supψ∈Hh,1(h, µ∇ψ)‖ψ‖Hh,1∀h ∈ Xh,1 . (4.35)Such a pair of spaces is called stable, and some examples will be considered later in§a.As discussed in (3.23), the discrete inf-sup condition is equivalent to the fact thatthe operator Ah,1 deﬁned by〈Ah,1h, (w, ψ)〉 = (h,∇×w) + (h, µ∇ψ) ∀h ∈ Xh,1 , (w, ψ) ∈ Yh,1 , (4.36)is bounded from below. In this case, we get the estimate (3.26) with α(h) = 0.For a given f = (j, ρ), consider the magnetostatic problem A1x = f. A naturaldiscrete approximation will be Ah,1xh = f. However, the set N(A∗h,1) is not easilycharacterized, and therefore, the problem for Ah,1 will involve awkward discrete com-patability conditions. To avoid those, we propose to use the least-squares method for60the discrete problem, i.e. solve(Ah,1xh,Ah,1yh)Y∗h,1 = (f,Ah,1yh)Y∗h,1 ∀yh ∈ Xh,1 . (4.37)The implementation of this problem reduces to a system of equations with a sym-metric and positive deﬁnite matrix as discussed in Section III.C.The approximation obtained by this method is further examined below. If fsatisﬁes the compatability conditions (4.15), then by Corollary 3.3, xh is a quasi-optimal approximation of x in Xh,1. The next result gives the rate of approximationwhen f ∈ L2(Ω)× L2(Ω). The proof is a combination of Theorems 3.2 and 4.2.Theorem 4.7 Assume that (Xh,1,Yh,1) is a stable pair, i.e. (4.35) holds. Let s ∈ [0, 1]be such that X1(µ) ↪→ Hs(Ω) and (4.34) hold. Denote with x the solution of theoriginal magnetostatic problem (4.14) with data f ∈ L2(Ω)×L2(Ω) which satisﬁes thecompatability conditions (4.15). Let xh be the least-squares approximation obtainedby solving (4.37). Then we have the error estimate‖x− xh‖ ≤ C(µ)hs ‖f‖ .Next, we give a short summary of the analogous results for the electrostaticproblem. We deﬁne Xh,2 = Xh,1, and Yh,2 = Vh,2 ×Hh,2, with Vh,2 ⊂ V2, Hh,2 ⊂ H2.For example, in the case Ωh = Ω, one can use Xh,2 = Ŝh, Vh,2 = (Sh ⊕ Bh)3 andHh,2 = Sh,0 ⊕ BFh,0. These spaces should be chosen in such a way, thatC ‖e‖X ≤ supw∈Vh,2(e,∇×w)‖w‖Vh,2+ supψ∈Hh,2(e, ε∇ψ)‖ψ‖Hh,2∀h ∈ Xh,2 . (4.38)Then the operator Ah,2 deﬁned by〈Ah,2e, (w, ψ)〉 = (e,∇×w) + (e, ε∇ψ) ∀e ∈ Xh,2 , (w, ψ) ∈ Yh,2 , (4.39)is bounded from below.61For a given f satisfying the compatability conditions (4.31), consider the electro-static problem A2x = f. Let xh be the least-squares approximation satisfying(Ah,2xh,Ah,2yh)Y∗h,2 = (f,Ah,2yh)Y∗h,2 ∀yh ∈ Xh,2 . (4.40)Then xh is a quasi-optimal approximation of x in Xh,2. If furthermore, f ∈ L2(Ω) ×L2(Ω) satisﬁes the compatability conditions (4.31), and X2(ε) ↪→ Hs(Ω) then we have‖x− xh‖ ≤ C(ε)hs ‖f‖ .a. Pairs of stable approximation subspacesIn this subsection, we discuss the construction of stable approximation pairs for thediv-curl systems (4.1). We concentrate on the magnetostatic problem, in which casewe need a pair (Xh,1,Yh,1) satisfyingC ‖x‖X1 ≤ supw∈Vh,1(x,∇×w)‖w‖V1+ supψ∈Hh,1(µ x,∇ψ)‖ψ‖H1∀x ∈ Xh,1 .This is similar to the famous LBB conditionC ‖x‖X1 ≤ supψ∈Hh,1(x,∇·ψ)‖ψ‖H1∀x ∈ Xh,1 (4.41)for an approximation pair Xh,1 ⊂ X1 = L2(Ω)/R, Hh,1 ⊂ H1 = H10(Ω) for the Stokesproblem1.For convenience, we restrict our discussion to tetrahedral partitioning of polyhe-dral domains Ω = Ωh ⊂ R3 and assume that µ is piecewise constant. The constructionextends to other element shapes as well as problems on domains in R2.The lowest order approximation is obtained for Xh,1 = Ŝh, i.e. when Xh,1 consists1The continuous LBB condition is equivalent to the Necˇas inequality from Theorem2.4, see [19]62of piecewise constant vector functions. As shown in [26], a compatible choice for thetest space is Yh,1 = Vh,1 ×Hh,1, with Vh,1 = (Sh,0 ⊕ BFh,0)3 and Hh,1 = Sh ⊕ BFh .The proof is based on the estimatessupψ∈H1(µ x,∇ψ)‖ψ‖H1≤ C supψ∈Hh,1(µ x,∇ψ)‖ψ‖H1∀x ∈ Xh,1 (4.42)andsupw∈V1(x,∇×w)‖w‖V1≤ C supw∈Vh,1(x,∇×w)‖w‖V1∀x ∈ Xh,1 , (4.43)with a constant C independent of h.To get a better approximation, we can choose Xh,1 to be the space of piecewiselinear functions, i.e. Xh,1 = Ŝh(1). Then a compatible choice for the test space, see[26], is Vh,1 = (Sh,0 ⊕ B1Fh,0 ⊕ BTh)3 and Hh,1 = Sh ⊕ B1Fh ⊕ BTh .Instead of dealing with these speciﬁc cases, we present the proof of stability inthe following more general case. To keep the notation uniform, we set B−1Th = ∅.Theorem 4.8 Let k ∈ N0. Then Xh,1 = Ŝh(k), Vh,1 = (Sh,0 ⊕ BkFh,0 ⊕ Bk−1Th )3 andHh,1 = Sh ⊕ BkFh ⊕ Bk−1Th satisfy (4.42)-(4.43). In particular (Xh,1,Vh,1 × Hh,1) is astable pair for the magnetostatic problem.Similarly, (Xh,2,Vh,2 ×Hh,2) is a stable pair for the electrostatic problem, whereXh,2 = Ŝh(k), Vh,2 = (Sh ⊕ BkFh ⊕ Bk−1Th )3 and Hh,2 = Sh,0 ⊕ BkFh,0 ⊕ Bk−1Th .Proof Fix x ∈ Xh,1. Let ψ ∈ H1 = H1(Ω) be arbitrary. To show (4.42) it is enoughto construct ψh ∈ Hh,1 such that(µ x,∇ψ) = (µ x,∇ψh) and ‖ψh‖H1 ≤ C ‖ψ‖H1 . (4.44)Let Ihψ be an approximation operator satisfying (2.26). Using that, the barycentriccoordinate functions are nonnegative. One can choose ψTh ∈ Bk−1Th and ψFh ∈ BkFh63such that(ψFh , q)L2(F ) = (ψ − Ihψ, q)L2(F ) ∀F ∈ Fh ,∀q ∈ Pk(F ) ,(ψTh , p)L2(τ) = (ψ − Ihψ − ψFh , p)L2(τ) ∀τ ∈ Th ,∀p ∈ Pk−1(τ) .(4.45)By Schwartz inequality,‖ψFh‖L2(F ) ≤ ‖ψ − Ihψ‖L2(F ) and ‖ψTh‖L2(τ) ≤ ‖ψ − Ihψ − ψFh‖L2(τ) . (4.46)Set ψh = Ihψ + ψFh + ψTh . Then(∇·x, ψ)L2(τ) = (∇·x, ψh)L2(τ) ∀τ ∈ Th ,(x · n, ψ)L2(F ) = (x · n, ψh)L2(F ) ∀F ∈ Fh ,(4.47)which implies the equality in (4.44).Using mapping to the reference element, the equivalence of norms on ﬁnite ele-ment spaces, (4.46), (2.23) and (2.26) we get‖ψFh‖21 ≤ C∑F∈Fhh−1F ‖ψFh‖2L2(F ) ≤ C∑τ∈Thh−1τ ‖(I− Ih)ψ‖2L2(∂τ)≤ C∑τ∈Th{h−2τ ‖(I− Ih)ψ‖2L2(τ) + ‖(I− Ih)ψ‖2H1(τ)}≤ C ‖ψ‖21 .Similarly,‖ψTh‖21 ≤ C∑τ∈Thh−2τ ‖ψTh‖2L2(τ)≤ C∑τ∈Th{h−2τ ‖(I− Ih)ψ‖2L2(τ) + h−2τ ‖ψFh‖2L2(τ)}≤ C ‖ψ‖21 .This implies that ψh satisﬁes (4.44), and therefore (4.42) holds.Clearly the above construction works for the case of zero boundary conditions,provided that they are preserved by Ih. Furthermore, by applying the constructionto each component, the same proof yields the result (4.43).64Next, we show how the above proof can be extended to more general cases. Tothat end, let τ ∈ Th and F ∈ Fh be any element and face of Ωh. Let Xτ = (Xτ )3 ⊂L2(Ω) be an approximation space for the solution on τ . Associated with this are twoadditional spacesXF =⋃τ∈TF{x|F : x ∈ Xτ} , X′τ =d⋃i=1{∂ix : x ∈ Xτ} . (4.48)Let XΩ ⊂ H1(Ω) be such that there exist an operator Ih : H1(Ω) → XΩ satisfying(2.26). As before, this operator should preserve the homogeneous boundary condi-tions. Finally, we need the spaces of face and element “bubbles”, which are just ﬁnitedimensional spaces BF ⊂ H10(TF ) and Bτ ⊂ H10(τ) . Assume that these spaces aredeﬁned through mappings to the reference element, i.e.‖ψF‖2H1(TF ) ≤ C h−1F ‖ψF‖2L2(F ) and ‖ψτ‖2H1(τ) ≤ C h−2τ ‖ψτ‖2L2(τ) (4.49)for any ψF ∈ BF and ψτ ∈ Bτ .Theorem 4.9 Assume thatXF ∩ B⊥F = {0} , X′τ ∩ B⊥τ = {0} . (4.50)Then the condition (4.42) holds for the following spacesXh,1 =⊕τ∈ThXτ and Hh,1 = XΩ⊕(⊕F∈FhBF )⊕(⊕τ∈ThBτ ) .Proof We follow the proof of the previous theorem by ﬁxing x ∈ Xh,1, ψ ∈ H1 andconsidering problems similar to (4.8): ψτ ∈ Bτ and ψF ∈ BF satisfy(ψF , q)L2(F ) = (ψ − Ihψ, q)L2(F ) ∀F ∈ Fh ,∀q ∈ XF ,(ψτ , p)L2(τ) = (ψ − Ihψ − ψFh , p)L2(τ) ∀τ ∈ Th ,∀p ∈ X′τ ,65where ψFh =∑F∈Fh ψF . By Theorem 3.1, the conditions (4.50) imply that theseproblems have unique minimum norm solutions which satisfy estimates similar to(4.46):‖ψF‖L2(F ) ≤ ‖ψ − Ihψ‖L2(F ) and ‖ψτ‖L2(τ) ≤ ‖ψ − Ihψ − ψFh‖L2(τ) .Deﬁne ψTh =∑τ∈Th ψτ and ψh = Ihψ+ψFh +ψTh . By the deﬁnitions (4.48), it followsthat the equalities (4.45) hold. This, together with (4.49), implies that ψh satisﬁes(4.44) and therefore, (4.42) is satisﬁed.Remark 4.5 As before, we can deﬁne Vh,1, just by taking into account the boundaryconditions in each component, i.e. Vh,1 = (Hh,1 ∩H10(Ω))3. Then (Xh,1,Vh,1 ×Hh,1)is a stable pair for the magnetostatic problem. Similarly (Xh,1,Vh,2×Hh,2) is a stablepair for the electrostatic problem, where Vh,2 = (Hh,1)3 and Hh,2 = Hh,1 ∩H10(Ω).Corollary 4.5 One can approximate the solution with polynomials of varying degree,i.e. Xτ = Pkτ (τ) with kτ ∈ N0 depending on τ . A typical example will be to usehigher order in the interior of each material (where the solution is smooth), andlower order close to the interfaces between diﬀerent materials (where the solution hassingularities). Set kF = max{kτ : τ ∈ TF}, then the conditions of the theorem aresatisﬁed for BF = (BkFFh)∣∣∣Fand Bτ = (Bkτ−1Th )∣∣τ.Corollary 4.6 The introduction of special bubble functions may be avoided if the testspace is deﬁned on a ﬁner mesh. Speciﬁcally, consider the case of triangular mesh,then Xh,1 = Sˆh(k) and Hh,1 = Sh/(2k+2) satisfy condition (4.42).Proof In this case Xτ = Pk(τ), XF = Pk(F ) and X′τ = Pk−1(τ). Deﬁne BF as thespan of the basis functions in Sh/(2k+2) with degrees of freedom in the interior of F .Similarly, Bτ is the span of the basis functions with degrees of freedom in the interiorof τ .66Since a nonzero polynomial in XF has at most k zeros, it follows that there exista basis function in BF such that their inner product in L2(F ) is positive. Similarconsiderations in τ show that the second condition in (4.50) is satisﬁed.Remark 4.6 This result can be extended to hexahedral and quadrilateral meshes. Thecorresponding result for tetrahedral meshes can be obtained for Hh,1 = Sh/(2k+3). Thenumerical results, however, indicate that the method is stable even if we use Sh/(2k+2).2. Approximation based on form modiﬁcationIt is possible to get a stable approximation even when the discrete inf-sup conditiondoes not hold. We illustrate this for the magnetostatic problem in the case Ωh = Ωwith Xh,1 = Ŝh, Vh,1 = (Sh,0)3 and Hh,1 = Sh. For simplicity, we take µ to be piecewiseconstant. The idea is to start with the lower bound for A1 given by inequality (4.9),i.e.C0µ˜0 ‖x‖2µ ≤ ‖curl1x‖2V∗1 + ‖div1,µx‖2H∗1 (4.51)and to strengthen the form using integration by parts and discretely deﬁned operators.To that end, let divh1,µ : Xh,1 → Hh,1 and curlh1 : Xh,1 → Vh,1 be deﬁned by(divh1,µxh, ψh) = −(µ xh,∇ψh) ∀ψh ∈ Hh,1 ,(curlh1xh,wh) = (xh,∇×wh) ∀wh ∈ Vh,1 .These operators are well deﬁned by the Riesz Representation Theorem, and theircomputation can be reduced to the inversion of the mass matrices in Hh,1 and Vh,1.We ﬁrst consider the second term in (4.51). For xh ∈ Xh,1 and any ψh ∈ Hh,1,‖div1,µxh‖2H∗1 = supψ∈H1(µ xh,∇ψ)2‖ψ‖2H1≤ 2 supψ∈H1[(µ xh,∇ψh)2‖ψ‖2H1+(µ xh,∇(ψ − ψh))2‖ψ‖2H1].(4.52)67Taking ψh = Ihψ, with Ih satisfying (2.26) and using integration by parts oneach element gives‖div1,µxh‖2H∗1 ≤ C supψh∈Hh,1(µ xh,∇ψh)2‖ψh‖2H1+ C supψ∈H1∑F∈Fh(µxh · n, ψ − ψh)2L2(F )‖ψ‖2H1,where Fh is the set of all faces of Th and · denotes the jump across a given face 2.Recall that hF denotes the diameter of the face F ∈ Fh and TF is the union ofall elements τ ∈ Th which have F as a face. Combining (2.23) and (2.26) we get‖ψ − ψh‖L2(F ) ≤ Ch1/2F ‖ψ‖H1(TF ).Therefore‖div1,µxh‖2H∗1 ≤ C ‖divh1,µxh‖2H∗h,1 + C∑F∈FhhF‖µ xh · n‖2L2(F ). (4.53)Similar manipulations imply that the ﬁrst term of (4.51) is bounded by‖curl1xh‖2V∗1 ≤ C‖curlh1xh‖V∗h,1 + C∑F∈FhhF‖xh×n‖2L2(F ) . (4.54)Consider the Hilbert spaceL2(Fh) = ⊕F∈FhL2(F ),with the inner product(u, v)L2(Fh) =∑F∈FhhF (u, v)L2(F ) .Furthermore, denote with Fh the subspace of L2(Fh) consisting of functions that areconstants on each face. Set L2(Fh) = (L2(Fh))3 and Fh = (Fh)3.2For a boundary face, the argument is assumed to be zero outside of Ω.68Combining the estimates (4.53) and (4.54) givesCµ˜0 ‖xh‖2µ ≤ ‖Ah,1xh‖2Y∗h,1×F∗h×F∗h ,where Ah,1 : Xh,1 → Y∗h,1 × F∗h × F∗h is deﬁned byAh,1xh = (curlh1xh,divh1,µxh, xh×n, µ xh · n) . (4.55)It is standard to see that Ah,1 is bounded. Moreover‖(A1 −Ah,1)xh‖2Y∗h,1×F∗h×F∗h =∑F∈FhhF{‖xh×n‖2L2(F ) + ‖µ xh · n‖2L2(F )} .The next result shows that we get essentially the same approximation rate as in themethod with bubble functions.Theorem 4.10 Let s ∈ [0, 1], s = 12be such that X1(µ) ↪→ Hs(Ω) and (4.34) hold.Denote with x the solution of the original magnetostatic problem (4.14) with dataf ∈ L2(Ω) × L2(Ω) which satisﬁes the compatability conditions (4.15). Let xh bethe least-squares approximation obtained by solving (4.37). Then we have the errorestimate‖x− xh‖ ≤ C(µ)hs ‖f‖ .Proof First, consider the case s > 12. We will apply Corollary 3.4 with Xˆ = X1(µ).Fix x ∈ X1(µ) and let ζh = QXhx. Then‖(A1 −Ah,1)ζh‖2Y∗h,1×F∗h×F∗h ≤ C∑τ∈Thhτ{‖(x− ζh)×n‖2L2(∂τ) + ‖µ (x− ζh) · n‖2L2(∂τ)}≤ C∑τ∈Th‖x− ζh‖2L2(τ) + h2sτ |x|2Hs(τ) .We used (2.24) and the fact that (2.10) implies |ζh|Hs(τ) = 0. The above inequalitymeans that (3.29) holds with α(h) = hs and therefore we get the result of the theoremby Corollary 3.4.69The case s < 12follows by interpolation as shown in [26]. We recall the detailsof the proof below. First, by Schwarz inequality and the boundedness of A1, Ah,1 weget ‖xh‖ ≤ C ‖x‖ which proves the case s = 0, i.e. we have the stability estimate‖x− xh‖ ≤ C ‖x‖ .On the other hand, recall the deﬁnitions (2.16) and let x ∈ (PH10(Ω))3. As in the cases > 12we can apply Corollary 3.4 with Xˆ = (PH10(Ω))3 and conclude that‖x− xh‖ ≤ C h ‖x‖H1(Ω) .Thus, by interpolation, we get‖x− xh‖ ≤ C hs ‖x‖Hs(Ω)for any x ∈ (PHs0(Ω))3 = Hs(Ω), which completes the proof.Remark 4.7 When s > 12it is straightforward to extend the above proof to the casewhen X1(µ) ↪→ Hs(Ω) is replaced with the weaker regularity assumption X1(µ) ↪→(PHs(Ω))3. This is a signiﬁcant improvement over the result of the theorem, sinceX1(µ) ↪→ Hs(Ω) only when µ ≡ const.In this case, the particular form of (4.37) is: Find xh ∈ Xh,1 satisfyingah,1(xh,yh) = (j, curlh1yh)V∗h,1 + (ρ,divh1,µyh)H∗h,1 , ∀yh ∈ Xh,1 , (4.56)where the corresponding bilinear form isah,1(xh,yh) = (curlh1xh, curlh1yh)V∗h,1 + (divh1,µxh,divh1,µyh)H∗h,1+∑F∈FhhF{(xh×n, yh×n)L2(F ) + (µxh · n, µyh · n)L2(F )} .We emphasize that even though it looks complicated, in some cases, the above form70might be easier to implement. For example, if µ = 1, the sum of the two jump termssimplify to (xh, yh)L2(F ).3. Extensions to more general domainsa. Domains with curved boundariesIn this section we consider domains with piecewise smooth boundaries. We concen-trate on the case when Ω has piecewise smooth boundary and Ωh ⊂ Ω are constructedsuch thatmaxx∈∂Ωdist(x, ∂Ωh) ≤ C h2 . (4.57)Here h is the mesh size of a globally quasiuniform mesh that triangulates Ωh. Theabove inequality, in particular, implies that the approximation of the boundary shouldimprove after reﬁnement, as shown on Figure 4.1. Similar construction can be carriedout in 3D, as shown in [68].Fig. 4.1. Curved boundary approximation.For standard elliptic second-order problems, this situation has been investigated,see [89], and is well understood. In particular, it is known that the “variational crime”of computing on Ωh instead of Ω does not aﬀect the approximation order for standardpiecewise linear ﬁnite elements.On the other hand, as stated in the preface of [75], the eﬀects of the approxi-71mation of smooth boundaries is not well understood for Maxwell’s equations. Somenumerical results from computations on a domain with curved boundaries using aweighted regularization method are reported in [44]. We note that their method in-volves the distance to the reentrant corners and edges of ∂Ω. This means that theimplementation in more general domains will be quite sophisticated.Below, we will show that the application of our method on Ωh gives an approx-imation to the solution on Ω for small enough h. This will be further demonstratedin Chapter VII, where we approximate the eigenvalues in a ball by computations ona sequence of inscribed, non-nested hexahedral meshes.Let Ωh ⊂ Ω and ω = Ω \ Ωh. Then (4.57) implies thatµ(ω) ≤ C h2 .It is well known that for any x ∈ H10(Ω),‖x‖L2(ω) ≤ C h‖x‖H10(Ω) .By interpolation, it follows that for any s ∈ [0, 1],‖x‖L2(ω) ≤ C hs‖x‖Hs0(Ω) . (4.58)Consider the weak formulation of a magnetostatic problem deﬁned on Ω whichinvolves the solution space X1 = L2(Ω) and the test space Y1 = V1 ×H1 = H10(Ω)×H1(Ω). Let (Xˇh,1, Yˇh,1) be an approximation pair for the above spaces on Ωh, i.e.Xˇh,1 ⊂ L2(Ωh) and Yˇh,1 = Vˇh,1 × Hˇh,1 ⊂ H10(Ωh)×H1(Ωh).Using extension by zero in ω, we consider Xˇh,1 and Hˇh,1 as subspaces of X1 andH1. Speciﬁcally if E0 : L2(Ωh) → L2(Ω) denotes the extension by zero operator, thenwe set Xh,1 = E0(Xˇh,1) and Vh,1 = E0(Vˇh,1). Assume that the same can be donefor Hˇh,1 using a bounded extension operator Eh, i.e. we set Hh,1 = Eh(Hˇh,1). The72extension Eh can be chosen to satisfy‖Ehψh‖Hs(ω) ≤ C hs‖ψh‖Hs(Ωh) ∀ψh ∈ Hh,1 s ∈ [0, 1] . (4.59)Such extensions are based on reﬂections of the values of the function in Ωh, and aspeciﬁc example is discussed in detail in [68]. However, in our development we willonly use the fact that Eh : H1(Ωh) → H1(Ω) is bounded.Consider the least-squares approximation based on Xh,1 and Yh,1 = Vh,1 ×Hh,1.For any x ∈ X1 and ξˇh ∈ Xˇh,1, we have‖x− E0ξˇh‖L2(Ω) ≤ ‖x‖L2(ω) + ‖x− ξˇh‖L2(Ωh) .By (4.58), if x ∈ Hs0(Ω) theninfξh∈Xh,1‖x− ξh‖L2(Ω) ≤ hs ‖x‖Hs0(Ω) + infξˇh∈Xˇh,1‖x− ξˇh‖L2(Ωh) .In particular, when Xˇh,1 consists of piecewise constants, we get an order of approx-imation s for any x ∈ Hs(Ω) with s ∈ [0, 1/2). This is the same as the order ofapproximation in L2(Ωh).Next, we look at the construction of a stable approximation operator. For sim-plicity, consider the speciﬁc case when the test spaces on Ωh contain Sh or Sh,0. Denotewith Iˇh an operator from either H1(Ωh) or H10(Ωh) to Sh or Sh,0, respectively, whichsatisﬁes (2.26) on Ωh. An examination of the proofs of discrete stability from theprevious sections shows that we need an operator Ih satisfyingh−2‖u− Ihu‖2L2(Ωh) + ‖Ihu‖2H1(Ωh) ≤ C‖u‖2H1(Ω) . (4.60)For u ∈ H1(Ω) we can set Ihu = EhIˇh(u|Ωh), and the inequality above is satisﬁed.The construction of Ih : H10(Ω) → E0(Sh,0) proceeds as follows: for u ∈ H10(Ω)let u˜ ∈ Sh be the function Iˇh(u|Ωh) modiﬁed by setting all degrees of freedom on ∂Ωh73equal to zero, i.e. if {vi} are the set of vertices of Ωh, we setu˜(vi) =⎧⎪⎨⎪⎩Iˇhu(vi) if vi ∈ ∂Ωh ,0 if vi ∈ ∂Ωh .Deﬁne Ihu = E0u˜. Using that for functions uh ∈ Sh, the norm ‖uh‖ is equivalent tohd∑uh(vi). As well as (2.23), we getC ‖(Ih − Iˇh)u‖2L2(Ωh) ≤∑vi∈∂Ωhh2 Iˇhu(vi)2 ≤ h ‖u− Iˇhu‖2L2(∂Ωh) + h ‖u‖2L2(∂Ωh)≤ ‖u− Iˇhu‖2L2(Ωh) + h2 ‖u− Iˇhu‖2H1(Ωh) + h ‖u‖2L2(∂ω) .Applying (2.23), (4.58) and the deﬁnition of Iˇh, we can conclude that‖(Ih − Iˇh)u‖2L2(Ωh) ≤ C h2 ‖u‖2H1(Ω) .Similarly,‖(Ih − Iˇh)u‖2H1(Ω) ≤ C∑vi∈∂ΩhIˇhu(vi)2 ≤ C ‖u‖2H1(Ω) ,and therefore Ih satisﬁes (4.60).Using the operator Ih, we can repeat the analysis in the case of stabilization byform modiﬁcation. The result is summarized below.Theorem 4.11 Let s ∈ [0, 12) be such that X1(µ) ↪→ Hs(Ω). Denote with x the solu-tion of the original magnetostatic problem (4.14) with data f ∈ L2(Ω)× L2(Ω) whichsatisﬁes the compatability conditions (4.15). Let xh be the least-squares approximationobtained by solving (4.37) with the form (4.55). Then we have the error estimate‖x− xh‖L2(Ω) ≤ C(µ)hs ‖f‖L2(Ω) .In the case of stabilization by adding bubble functions, we can perform theanalysis if we introduce bubbles on ∂Ωh, regardless of the boundary conditions. This74is because, for u ∈ H10(Ω), the diﬀerence u − Ihu is not zero on ∂Ωh. However, weshow below that this can be avoided, provided that the mesh size is small enough.Fix xh ∈ Xh,1, and let ψ ∈ Hh,1. Integration by parts yields(xh,∇ψ)‖ψ‖1 =∑F∈Fh, F ⊂∂Ωh(xh · n, ψ)L2(F )‖ψ‖1 +∑F∈Fh, F⊂∂Ωh(xh · n, ψ)L2(F )‖ψ‖1 .The sum over the boundary faces can be estimated as∑F⊂∂Ωh|(xh · n, ψ)L2(F )| ≤∑F⊂∂Ωh‖xh‖L2(F )‖ψ‖L2(F ) ≤ C h−12 ‖xh‖L2(Ωh) h ‖ψ‖H1(Ωh) .Therefore(1− C h 12 )(xh,∇ψ)‖ψ‖1 ≤∑F∈Fh, F ⊂∂Ωh(xh · n, ψ)L2(F )‖ψ‖1 ,and the result follows from the previous considerations, provided that h < C−2.We can conclude that we get the result of the previous theorem for the least-squares method based on discrete inf-sup condition obtained by enriching the testspaces with bubble functions. Speciﬁcally, we can either add face bubble functionson the whole boundary, regardless of the boundary conditions, or we can skip thebubbles on the boundary, provided that the mesh size is suﬃciently small.b. Multiply-connected domains with holesHere we consider a general domain Ω = Ωh, as the one shown in Figure 2.1, whichsatisﬁes assumption (AΩ). Speciﬁcally, we allow for domains with multiple boundarycomponents {Γi}n1i=0, where Γ0 is the outer boundary and domains that are multiply-connected, depending on the number of cuts {Σj}n2j=1. The extension, especially thecase n2 > 0, is based on the theory developed in [4] and §5 of [20]. Some discussionof the deﬁnition and the implementation of the cuts Σj can be found in §8.3.4 of [18].Recall that Ω0 = Ω \⋃n2j=1 Σj is simply connected. For ψ ∈ H1(Ω0), we denote75with ∇˜ψ its distributional gradient with respect to Ω0, considered as an element ofL2(Ω). This is generally diﬀerent from ∇ψ ∈D′(Ω).In contrast to the special cases considered before, when n1 > 0 and n2 > 0,the magnetostatic and the electrostatic problems (4.1) does no longer have uniquesolutions. Namely, there are nonzero ﬁelds that solve the corresponding homogeneousproblems. These ﬁelds form the spaces K1(µ) and K2(ε) introduced belowK1(µ) = {u ∈ X1(µ) : ∇×u = 0 and ∇·µu = 0} ,K2(ε) = {u ∈ X2(ε) : ∇×u = 0 and ∇·εu = 0} .(4.61)It is shown in [4] that these spaces are fundamentally related to the topologicalcharacteristics n1 and n2. In fact,dim(K1(µ)) = n2 and dim(K2(ε)) = n1 .It is also known, see Theorem 8’ in [35], that the followig graphH1(Ω)∇ H(curl) ∇× H(div) ∇· L2(Ω)K1(1)K2(1)L2(Ω) ﬀ∇·H0(div)ﬀ∇× H0(curl)ﬀ∇ H10(Ω)(4.62)is “exact”, i.e. the kernel of each operator is the (direct) sum of the images of theprevious operators in the sequence.Moreover, K2(ε) has the basis {∇ψi}n1i=1, where ψi ∈ H1(Ω) satisﬁes⎧⎪⎨⎪⎩−∇·ε∇ψi = 0, in Ω ,ψi = δij, on Γj , 0 ≤ j ≤ n1 .(4.63)Here δij denotes the Kronecker Delta. The functions {ψi}n1i=1 form a linear spacewhich we denote with K2(ε).76Similarly, K1(µ) has the basis {∇˜ζj}n2j=1, where ζj ∈ H1(Ω0) satisﬁes⎧⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩−∇·µ∇ζj = 0, in Ω0,µ∂ζj∂n= 0, on ∂Ω,ζji = δij andµ∂ζj∂ni= 0 on Σi , 1 ≤ i ≤ n2 .(4.64)Here ·i denotes the jump across Σi. The functions {ζj}n2j=1 form a linear spacedenoted by K1(µ).These characterizations, together with a compactness argument similar to theone in Proposition 4.2, can be used to prove that‖∇×h‖+ ‖∇·µh‖+n2∑j=1|〈µh · n, 1〉Σj | is an equivalent norm on X1(µ) .Similarly,‖∇×e‖+ ‖∇·εe‖+n1∑i=1|〈εe · n, 1〉Γi| is an equivalent norm on X2(ε) .It is therefore natural to consider the following generalized magnetostatic and elec-trostatic problems⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩∇×h = j in Ω,∇ · (µh) = ρ in Ω,µh · n = σ on ∂Ω,〈µh · n, 1〉Σj = Cj 1 ≤ j≤ n2 ,⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩∇×e = j in Ω,∇ · (εe) = ρ in Ω,e×n = σ on ∂Ω,〈εe · n, 1〉Γi = Ci 1 ≤ i≤ n1 .(4.65)As a straightforward corollary of the norm equivalence we get that for data j ∈ L2(Ω),ρ ∈ L2(Ω), σ ∈ H− 12 (∂Ω), σ ∈ H− 12 (∂Ω), and {Ck} ⊂ R satisfying appropriatecompatability conditions, the above problems have unique solutions h ∈ H(curl) ∩H(div;µ) and e ∈ H(curl) ∩H(div; ε).77Next, we discuss how we can extend our weak formulations to include thesesystems. We start with a generalization of Theorem 2.8. We will need the followingresults proven as Theorem 3.4 in [54] and Theorem 3.17 in [4].Lemma 4.1 Let v ∈ L2(Ω) with ∇·v = 0. Then there exist a vector potentialv = ∇×win the following cases:1. If 〈v · n, 1〉Γi = 0, 0 ≤ i ≤ n1 then w ∈ H1(Ω).2. If v · n = 0, 〈v · n, 1〉Σj = 0, 0 ≤ j ≤ n2 then w ∈ H0(curl).Using the Lemma, we can extend the Helmholtz decomposition results.Theorem 4.12 Let u be in L2(Ω). Then it can be decomposed asu = ∇×w+ µ∇ψ (4.66)in the following spaces 31. w ∈ H0(curl) and ψ ∈ H1(Ω)⊕ K1(µ).2. w ∈ H10(Ω) and ψ ∈ H1(Ω)⊕ K1(µ).3. w ∈ H1(Ω) and ψ ∈ H10(Ω)⊕ K2(µ).In the last two cases, we additionally have‖w‖H1(Ω) ≤ C‖∇×w‖L2(Ω).Proof Let ϕ be the unique element of H1(Ω)/R satisfying(µ∇ϕ,∇θ) = (u,∇θ) , (4.67)3For ψ = φ + ζ ∈ H1(Ω)⊕ K1(µ), ∇ψ is understood as ∇φ + ∇˜ζ.78for any θ ∈ H1(Ω). Then ∇·(u− µ∇ϕ− µ ∇˜ζ) = 0 and (u− µ∇ϕ− µ ∇˜ζ) · n = 0on ∂Ω, for any ζ ∈ K1(µ). Furthermore, for ﬁxed Σj with 1 ≤ j ≤ n2〈µ ∇˜ζ · n, 1〉Σj =n2∑i=1〈µ ∇˜ζ · n, ζj〉Σi = (µ∇ζ,∇ζj)L2(Ω0) ,by the generalization of the Green’s formula given as Lemma 3.10 in [4]. This impliesthat ζ can be chosen appropriately, so that u− µ∇ϕ− µ ∇˜ζ satisﬁes the second setof conditions in Lemma 4.1. This gives the ﬁrst decomposition.Examining the proofs of Lemma IV.1 and Lemma IV.2 in [99], one can concludethat the decomposition w = w˜ + ∇ξ with w˜ ∈ H10(Ω), ξ ∈ H1(Ω) and ‖w˜‖H1(Ω) ≤C‖∇×w˜‖ holds in the general case n1 > 0, n2 > 0. This proves decomposition 2.For the last result, choose ϕ to be the unique element of H10(Ω), satisfying (4.67)for any θ ∈ H10(Ω). Then ∇·(u−µ∇ϕ−µ∇ψ) = 0 for any ψ ∈ K2(µ). Furthermore,for ﬁxed Γi with 1 ≤ i ≤ n1〈µ∇ψ · n, 1〉Γi = 〈µ∇ψ · n, ψi〉Γ = (µ∇ψ,∇ψi)L2(Ω) ,and therefore ψ can be chosen in such a way that 〈(u− µ∇ϕ− µ∇ψ) · n, 1〉Γi = 0.Now, the result follows from the proof of the existence for the ﬁrst vector potentialin Lemma 4.1.Since the weak formulations for the magnetostatic and electrostatic problems pre-sented in the ﬁrst two sections of this chapter were essentially based on the Helmholtzdecompositions, it follows that the theory presented there will work if we increase thespaces H1 and H2 asH1 = H1(Ω)⊕ K1(µ) and H2 = H10(Ω)⊕ K2(ε) . (4.68)79Speciﬁcally, the weak formulation of the magnetostatic problem from (4.65) is⎧⎪⎨⎪⎩curl1h = j′ in V∗1 ,div1,µh = ρ′ in H∗1 ,where div1,µh was naturally extended to a bounded linear functional on K1(µ)∗ by〈div1,µh, ζ〉 = −(µh, ∇˜ζ)L2(Ω) ∀h ∈ L2(Ω), ζ ∈ K1(µ) .Furthermore, j′ = j is in H10(Ω)∗and ρ′ ∈ (H1(Ω) ⊕ K1(µ))∗ is deﬁned as follows: ifψ = φ +∑n2j=1 αjζj ∈ H1(Ω)⊕ K1(µ), then〈ρ′, ψ〉 = 〈ρ, ψ〉 − 〈σ, ψ|∂Ω〉 −n2∑j=1αjCj ,where ρ ∈ H1(Ω)∗ and σ ∈ H− 12 (∂Ω).Similarly, the weak formulation of the electrostatic problem from (4.65) reads⎧⎪⎨⎪⎩curl2e = j′ in V∗2 ,div2,εe = ρ′ in H∗2 ,where, div2,εe was naturally extended to a bounded linear functional on K2(ε)∗,j′ = j − σ with j ∈ H1(Ω)∗, σ ∈ H− 12 (∂Ω) and ρ′ ∈ (H10(Ω) ⊕ K2(ε))∗ is deﬁned asfollows: if ψ = φ +∑n1i=1 αiψi ∈ H10(Ω)⊕ K2(ε), then〈ρ′, ψ〉 = 〈ρ, ψ〉 −n1∑i=1αiCi ,where ρ ∈ H10(Ω)∗.We next consider discretization. Without loss of generality, we may assume thatthe cuts {Σj} align with the mesh. Note that an alternative characterization of (4.68)80isH1 = {ζ ∈ H1(Ω0) : ζj = const , 1 ≤ j ≤ n2} ,H2 = {ψ ∈ H1(Ω) : ψ|Γ0 = 0 , ψ|Γi = const , 1 ≤ i ≤ n1 .}(4.69)Therefore, to deﬁne Hh,1, we start with the usual approximation space for H1(Ω)and append functions which are discontinuous on the cuts. Speciﬁcally, we add basisfunctions which are 1 on the nodes on one side of Σj and vanish on all remainingnodes (including those on the opposite side of Σj).The discrete least-squares methods are still stable. For example, the methodwith bubble functions was based on the fact that for given x ∈ Xh,1 and (v,h) ∈ Y1,one then constructs a pair (vh,hh) ∈ Yh,1 satisfyinga1(x, (vh,hh) = a1(x, (v,h)) (4.70)and‖(vh,hh)‖Y1 ≤ C‖(v,h)‖Y1 .The construction started with a stable approximation operator Ih as an initial ap-proximation and then used the bubble functions to enforce (4.70) on the remainder.Similarly, the method based on form modiﬁcation depended only on integration byparts and the properties of Ih.Thus, to prove stability of the two discrete least-squares methods, we only needto demonstrate the construction of Ih satisfying (2.26). We simply use a modiﬁedapproximation operator for H1. Speciﬁcally, let I˜h be a stable approximation operatorinto the subspace of piecewise linear functions with arbitrary discontinuities acrossthe cuts, and deﬁne Ihh equal to I˜hh on the nodes not on the cut and by a boundaryaveraging operator (on each side of the cut) such as that given in [87]. This resultsin a stable approximation operator. Moreover, since h diﬀers by a constant on each81side of the cut, and the boundary averaging operator preserves constants, Ihh is inHh,1 and has the same jumps as h. Using Ih, the remainder of the proof consideredbefore goes through.For Hh,2, we start with the ﬁnite element approximation of H10(Ω) and appendbasis functions which are one on a given connected component of the boundary andvanish at all remaining nodes. To prove the stability of the discrete least-squaresmethods, we are again left with the construction of a suitable stable approximationoperator. If I˜h denotes a stable approximation operator into the ﬁnite element sub-space with arbitrary boundary values, we set Ihh to be I˜hh at the interior nodes andinterpolate h at the boundary nodes. It is easy to prove, similar to the case of curvedboundary, that Ih is a stable interpolation operator which reproduces h on ∂Ω. Withthis operator, the proof proceeds as before.82CHAPTER VTHE EIGENVALUE PROBLEMIn this chapter we consider the time-harmonic eigenvalue problem (1.7), i.e. we arelooking for the eigenvalues λ ∈ C and their corresponding magnetic and electriceigenfunctions h ,e : Ω→ C3 satisfying1⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩∇×h = λ εe in Ω,∇×e = −λµh in Ω,e×n = 0 on ∂Ω,µh · n = 0 on ∂Ω .(5.1)Clearly λ = 0 is an eigenvalue with an eigenspace consisting of gradients2. Thisdissertation deals only with the physically more interesting case λ = 0. Then astandard interpretation of (5.1) is to look for h ∈ X1(µ) and e ∈ X2(ε). We refer tothis case as the original form of the eigenvalue problem.Note that since λ = 0, we can use Theorem 2.5 to deduce the usual divergenceequations from (5.1): ⎧⎪⎨⎪⎩∇ · (µh) = 0 in Ω,∇ · (εe) = 0 in Ω .(5.2)Even though (5.2) is a corollary of (5.1), we will see that a good approximationmethod should take these equations into account explicitly.One of the more popular approaches to the eigenvalue problem is to eliminate1The following additional terminology is often used: the Maxwell eigenvalues arecalled eigenfrequencies, and the eigenfunctions are called eigenmodes, eigenﬁelds oreigenvectors.2The eigenvectors in this case can be completely characterized using the exactsequence with zero boundary conditions from (4.62).83one of the ﬁelds, e.g. h, and reduce it to a second-order problem for e. The reducedproblem involves the curl-curl operator and reads∇×µ−1∇×e = ω2εe , (5.3)where ∇·εe = 0, e×n = 0 and ω2 = −λ2. This, of course, is understood in the sensethat ω2 ∈ C and e ∈ HC0 (curl) satisfy(µ−1∇×e,∇×w) = ω2 (εe,w) ∀w ∈ HC0 (curl) . (5.4)A straightforward corollary of (5.4) is that ω2 ∈ R, and therefore, the eigenvaluesof the original eigenvalue problem (5.1)-(5.2) are purely imaginary and symmetricwith respect to the origin:λ = ± i ω , ω ∈ R+ . (5.5)Consequently, (e,h, λ) is an eigenpair of (5.1) if and only if ((e), i(h), λ) and((e),−i(h), λ) are eigenpairs of (5.1). This means that to exhibit a basis for theeigenspace corresponding to an eigenvalue of the form (5.5), we can restrict to realelectric and purely imaginary magnetic eigenfunctions.In practical applications ω corresponds to the frequency of propagation, andthe goal is to compute the ﬁrst few minimal positive ω with their correspondingeigenﬁelds. This is critical, for example, in the design of accelerator structures wherethe computed eigenfunctions are used as a “wake ﬁeld”, see [2] as well as [103] andthe reports therein.The importance of the Maxwell eigenvalue problem has led many authors toinvestigate its numerical approximation. A detailed survey of a variety of diﬀerentmethods was published recently in [44]. Early engineering approximations used con-forming ﬁnite element spaces to approximate (5.4), transformed to a vector Helmholtz84equation by Theorem 2.5. It was observed, see [84, 16], that the discrete method con-verges, but to a wrong solution3! Such solutions are called spurious and can beavoided if, e.g., the divergence equations (5.2) are properly taken into account.Fig. 5.1. Cross-section of a coaxial cable with an oﬀset center conductor: pollution byspurious modes. After Paulsen and Lynch [84], c© 1991 IEEE.This is illustrated by Figure 5.1 where we present an example taken from [84].On the left, we have a mesh showing two regions with diﬀerent material properties.In the middle, we have the real part of the reference electric ﬁeld computed by solvingthe eigenvalue problem. On the right, we show the real part of the discrete solutionobtained by straightforward application of node-based ﬁnite elements. Clearly thisapproximation is severely distorted by spurious modes. For more details see ??, or[59], pp. 200–202.Other approaches based on conforming ﬁnite elements are known to have prob-3Speciﬁcally, recall the deﬁnition of PHs(Ω) in (2.16) and setH1(µ) = X1(µ) ∩ (PHs(Ω))3 and H2(ε) = X2(ε) ∩ (PHs(Ω))3 .It is shown in [45], that those spaces are closed in X1(µ) and X2(ε) respectively.However, if Ω is not convex, H1(µ) X1(µ) with inﬁnite co-dimension. In fact,X1(µ) = H1(µ) ⊕ ∇Sµ, where Sµ is a space of singular functions for the operator−∆Neuµ (see [50] for details). We conclude that there are elements in X1(µ), X2(ε)which can not be approximated (in ‖·‖X1(µ) and ‖·‖X2(ε)) by continuous ﬁnite elements.85lems due to low regularity solutions and multiple valued potentials [52, 61, 71].Various alternatives have been proposed in order to avoid spurious solutions.One of the more popular approaches for this problem is based on the curl-conformingspaces, such as those developed by Ne´de´lec (cf. [77, 78]). Analysis of the eigenvalueproblem using these spaces either involves proving collective compactness4 [69, 76] orproving convergence in norm [17, 15]. The discrete eigenvalue problem is then solvedby the use of a shift-and-invert algorithm. A prerequisite for this algorithm is anestimate for the eigenvalue, which may be diﬃcult to obtain.New methods for dealing with these problems have been introduced recently[43, 48, 85]. The methods of [43] depend on a weighted functional with weightsdepending on the strength of the singularities at corners and edges. In [48] the authorsproved discrete compactness in two dimensions for a class of hp ﬁnite elements. Aninterior penalty discontinuous Galerkin method is proposed in [85].The approach which is presented in this chapter is based on [20]. We ﬁrst relatethe problem to a block system involving the solution of two div-curl systems. Thesesystems are formulated as variational problems corresponding to a magnetostatic andan electrostatic problems following the theory developed in Chapter IV. We then showthat the eigenfunctions with non-zero eigenvalues are also eigenfunctions of a compactskew-Hermitian problem and use our div-curl approximation to derive a sequence ofapproximation operators. Note that since the curl-curl operator is not elliptic, itsinverse is not compact which leads to much more complicated analysis. In contrast,our formulation involves the compact “pseudo” inverse mentioned above. To obtaina system which is more amenable to iterative computation, we next show that theoriginal eigenpairs can be computed from those of a compact symmetric real operator.4We say that {Kn} ⊂ L(X,Y) are collectively compact if for any bounded setM ⊂ X, the set ∪nKn(M) has a compact closure in Y. See [75], pp. 32.86This represents a signiﬁcant computational advantage since the iterative techniquesfor computing the eigenvalues of large symmetric problems are more eﬃcient androbust than those developed for non-symmetric and/or indeﬁnite systems.A. Reformulation of the eigenvalue problemFor simplicity, we shall assume that Ω = Ωh is simply connected with a connectedboundary, i.e. n1 = 0, n2 = 0. The case of more general domains will be addressed in§D.We quote the following result for the original eigenvalue problem, in the form(5.4), given as Theorem 4.18 in [75].Theorem 5.1 There is an inﬁnite discrete set of eigenvalues 0 < ω21 ≤ ω22 ≤ . . . withcorresponding eigenfunctions en = 0, such that ωn → ∞ and (en,em)L2ε(Ω) = 0 form = n.Next, we reformulate (5.1)-(5.2) by showing that it is related to an eigenvalueproblem involving a compact Hermitian semideﬁnite operator.Suppose that e ∈ X2(ε), h ∈ X1(µ) is an eigenpair corresponding to a nonzeroeigenvalue λ. The idea is to split the original problem into two independent mag-netostatic and electrostatic systems. Namely, it is natural to consider the followingsource problems ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×h = f1 in Ω,∇ · (µh) = 0 in Ω,µh · n = 0 on ∂Ω,(5.6)87and ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×e = f2 in Ω,∇ · (εe) = 0 in Ω,e×n = 0 on ∂Ω.(5.7)Clearly these are equivalent to (5.1)-(5.2) if we set f1 = λ εe and f2 = −λµh. Belowwe refer to both of these problems by the use of the subscript k which equals 1 for(5.6) and is 2 for (5.7).Recall that the weak formulations introduced in Chapter IV, involved the solutionspaces Xk and the test spaces Yk = Vk×Hk deﬁned by (4.5) and (4.20). Furthermore,recall deﬁnitions (4.13) and (4.29) of the spaces Vk,0 related to the compatabilityconditions. Let Q1 : L2(Ω) → V1,0 be the L2ε(Ω) orthogonal projection onto V1,0, i.e.Q1w = ∇ϕ, where ϕ ∈ H10(Ω) satisﬁes(ε∇ϕ,∇θ) = (εw,∇θ) ∀θ ∈ H10(Ω) .Similarly, Q2 : L2(Ω) → V2,0 is the L2µ(Ω)-projection deﬁned by Q2w = ∇ϕ, whereϕ ∈ H1(Ω)/R satisﬁes(µ∇ϕ,∇θ) = (µw,∇θ) ∀θ ∈ H1(Ω)/R .By Theorem 4.1, for any g1 ∈ L2ε(Ω), the weak formulation of the magnetostaticproblem (5.6) with data f1 = ε (I − Q1)g1 will have a unique solution h ∈ L2µ(Ω).Therefore, we can deﬁne the solution operator S1 : L2ε(Ω) → L2µ(Ω), by x1 = S1g1,where x1 ∈ L2µ(Ω) satisﬁesa1(x1, (v,h)) ≡ (x1,∇×v) + (x1, µ∇h) = (ε(I− Q1)g1, v), ∀(v,h) ∈ Y1 . (5.8)88Furthermore, Corollary 4.2 implies thatS1 : L2ε(Ω) → X1(µ) , ∇×S1g1 = ε(I− Q1)g1 and ∇·(µS1g1) = 0 . (5.9)Similarly, Theorem 4.4, implies that for any g2 ∈ L2µ(Ω), the weak formulation ofthe electrostatic problem (5.7) with data f2 = µ (I− Q2)g2 will have unique solutione ∈ L2ε(Ω). Therefore, we can deﬁne the solution operator S2 : L2µ(Ω) → L2ε(Ω), byx2 = S2g2, where x2 ∈ L2ε(Ω) satisﬁesa2(x2, (v,h)) ≡ (x2,∇×v) + (x2, ε∇h) = (µ(I− Q2)g2, v), ∀(v,h) ∈ Y2 .By Corollary 4.4S2 : L2µ(Ω) → X2(ε) , ∇×S2g2 = µ(I− Q2)g2 and ∇·(εS2g2) = 0 . (5.10)The ﬁrst result of this section is that the solution operators Sk can be used toobtain an equivalent formulation of the eigenvalue problem.Theorem 5.2 Consider the block-matrix operator B : L2ε(Ω) × L2µ(Ω) → L2ε(Ω) ×L2µ(Ω) deﬁned byB =⎛⎜⎝ 0 −S2S1 0⎞⎟⎠ . (5.11)Then, (e,h, λ) with λ = 0 is an eigenpair of the original eigenvalue problem (5.1)-(5.2) if and only ifB⎛⎜⎝eh⎞⎟⎠ = σ⎛⎜⎝eh⎞⎟⎠ , (5.12)with σ = λ−1.89Proof Let e ∈ X2(ε), h ∈ X1(µ), λ = 0 satisfy (5.1). By density (Theorem 2.5),(∇×h,∇ψ) = 0 = λ(e,∇ψ) ∀ψ ∈ H10(Ω) ,(∇×e,∇ψ) = 0 = −λ(µh,∇ψ) ∀ψ ∈ H1(Ω) .It follows that the eigenfunctions deﬁne compatible dataQ1e = 0 and Q2h = 0 .Therefore,S1(λe) = h and S2(−λh) = e .On the other hand, let e ∈ L2ε(Ω), h ∈ L2µ(Ω) and σ = 0 satisfy (5.12). By (5.9)and (5.10), it follows that e ∈ X2(ε), h ∈ X1(µ) and∇×σh = ε(I− Q1)e , ∇×σe = −µ(I− Q2)h .Furthermore, the divergence equations from (5.9) and (5.10) plus the boundary con-ditions in X1(µ) implyQ2S1e = 0 and Q1S2h = 0 ∀e ∈ L2ε(Ω) , h ∈ L2µ(Ω) . (5.13)Therefore, (e,h, σ−1) is an eigenpair of the original eigenvalue problem.Next we investigate the properties of B. Theorem 2.7 implies that Sk, andtherefore B, is a compact operator. We claim that B is also skew-Hermitian onL2ε(Ω)× L2µ(Ω). Indeed,⎛⎜⎝B⎛⎜⎝eh⎞⎟⎠ ,⎛⎜⎝e˜h˜⎞⎟⎠⎞⎟⎠ = −(S2h, εe˜) + (S1e, µh˜).90Using (5.9), (5.10) and (5.13) gives(S2h, εe˜) = (S2h, ε(I− Q1)e˜) = (S2h,∇×S1e˜)= (∇×S2h, S1e˜) = (µ(I− Q2)h, S1e˜) = (µh, S1e˜)from which it follows that B is skew-Hermitian. Note that the above identity is justS∗1 = S2 , (5.14)where S1 is considered as an operator from L2ε(Ω) to L2µ(Ω). When S1 is consideredas an operator on L2(Ω), we will denote its adjoint by St1. Clearly S∗1 = ε−1St1µ.The eigenvectors and eigenvalues ofB are related to the compact positive semidef-inite operator−B2 =⎛⎜⎝S2S1 00 S1S2⎞⎟⎠ . (5.15)This operator is Hermitian relative to the inner product on L2ε(Ω) × L2µ(Ω). Thenonzero eigenvalues and the corresponding eigenvectors for B can be recovered fromthose of either diagonal block above. For example, S2S1 : L2ε(Ω)→ L2ε(Ω) is Hermitianand if the real function e satisﬁesS2S1e = τ2e (5.16)then, ⎛⎜⎝ eiτS1e⎞⎟⎠ and⎛⎜⎝ e−iτS1e⎞⎟⎠are eigenvectors for B with eigenvalues −iτ and iτ , respectively. We get all nonzeroeigenvalues and their corresponding eigenvectors this way. Indeed, (e,h, σ) satisfy5.12 if and only if S1e = σh and S2h = −σe. By elimination of h, this is clearlyequivalent to S2S1e = −σ2e and h = σ−1S1e. The rest follows from the fact that σ91is purely imaginary (since B is skew-Hermitian).Remark 5.1 Parallel to B, consider the real compact Hermitian operatorB˜ =⎛⎜⎝ 0 S2S1 0⎞⎟⎠ . (5.17)Clearly (e,∓ ih) is an eigenfunction of B corresponding to the eigenvalue σ = ± i τif and only if (e,h) is an eigenfunction of B˜, corresponding to the eigenvalue τ .We summarize the above considerations, plus the Hilbert-Schmidt theory fromTheorem 2.2, in the main result of this section:Theorem 5.3 The operator S2S1 : L2ε(Ω) → L2ε(Ω) is compact, Hermitian and pos-itive semi-deﬁnite. It has a countable sequence of nonzero positive eigenvalues τ 2n ∈R+, τ 21 > τ22 > . . . > 0, each of ﬁnite multiplicity and such that τn → 0. Theeigenvectors corresponding to diﬀerent eigenvalues are orthogonal in L2ε(Ω).Furthermore, the eigenvalues of the original eigenvalue problem (5.1)-(5.2) aregiven by λ = ±i τ−1. Therefore, the problem of computing the few minimal (in modu-lus) eigenvalues λ translates into computing the few maximal eigenvalues of S2S1. Abasis for the eigenspace corresponding to λ is given by (e, λS1e), where e is a basisof the real eigenfunctions of S2S1 corresponding to τ2.Finally, the problems (5.3) and (5.16) are equivalent, with ω2 = τ−2. In partic-ular, this theorem is equivalent to Theorem 5.1.Remark 5.2 Analogous result holds for the lower right diagonal block of −B2, i.e.S1S2 : L2µ(Ω) → L2µ(Ω) is compact, Hermitian (in L2µ(Ω)) and positive semi-deﬁnite.The diﬀerence is that, this way, we exhibit a diﬀerent basis for the eigenspace con-sisting of real magnetic and purely imaginary electric ﬁelds.92B. Approximation of the least-squares solution operatorsWe next consider approximation to the eigenvalue problem for the operator S2S1. Ourﬁrst goal is to deﬁne discrete approximations to each of S2 and S1.Let Xh,1 = Xh,2 = Xh ⊂ L2(Ω) and Yh,k ⊂ Yk be approximation subspaces asdiscussed in §IV.C. We consider both of the discrete least-squares methods (basedon a discrete inf-sup condition and based on form modiﬁcation) presented there. Thesimplest discretization involved setting Xh to be the space of piecewise constant vectorﬁelds with respect to the mesh, with companion spaces Yh,k consisting of continuouspiecewise linear functions (satisfying the appropriate boundary conditions). For themethod based on a discrete inf-sup condition, Yh,k had to be enriched with bubblefunctions on the faces.As before, we assume that (AΩ) and (Aµ,ε) hold, and that there exists s ∈ [0, 1],such that the estimate (3.21) holds with χ(h) = C hs, i.e.infxh,k∈Xh,k‖xk − xh,k‖ ≤ C hs ‖xk‖s ∀xk ∈ Hs(Ω) . (5.18)Additionally, we assume the continuous embeddings (see Theorem 2.7)X1(µ) ,X2(ε) ↪→ Hs(Ω) . (5.19)Below, we brieﬂy recall the setup of the discrete least-squares approximations tothe magnetostatic and electrostatic problems as presented in Chapter IV. We ﬁrstset Th,k ≡ TYh,k : Y∗h,k → Yh,k by(Th,k, v)1 = 〈, v〉 ∀v ∈ Yh,k .The approximation xh,k is then deﬁned to be the unique function in Xh satisfying〈Ah,kxh,k,Th,kAh,ky〉 = 〈fk,Th,kAh,ky〉, ∀y ∈ Xh, (5.20)93where Ah,k is a map of Xh into Y∗h,k. The deﬁnition of Ah,k is given by (4.36)-(4.39)or (4.55) depending on the method.To deﬁne our approximation for Sk, we ﬁx gk ∈ L2(Ω) and set fk in (5.20) by〈f1, (v,h)〉 = (ε(I− Qh,1)g1, v) , 〈f2, (v,h)〉 = (µ(I− Qh,2)g2, v)and deﬁne Sh,kgk = xh,k. The operators Qh,k, k ∈ {1, 2}, are deﬁned in terms of theapproximation subspace for Vk,0. For example, if Hh,k is the approximation subspaceassociated with Yh,k, then we deﬁne Qh,1v = ∇φ where φ ∈ Hh,2 satisﬁes(ε∇φ,∇θ) = (ε v,∇θ) ∀θ ∈ Hh,2. (5.21)Similarly, we deﬁne Qh,2v = ∇φ where φ ∈ Hh,1 satisﬁes(µ∇φ,∇θ) = (µ v,∇θ) ∀θ ∈ Hh,1. (5.22)Remark 5.3 Actually, as will become clear later, the bubble functions are not neededfor Qh,k. For example, for the case when Xh is piecewise constant, it suﬃces to usethe subspaces of piecewise linear functions with appropriate boundary conditions.To analyze the approximation properties of the above operators, we shall needregularity results for second-order problems with piecewise smooth coeﬃcients. Specif-ically, we will assume that either (AL2∆Dirε ,∆Neuµ ) or (A∆Dirε ,∆Neuµ ) holds.Additionally, we will need the followingAssumption (Aε,µ×) There exists γ0 ∈ (0, 1] such that the operators of multiplica-tion by ε and µ are bounded from H1(Ω) to Hγ(Ω) for any 0 ≤ γ < γ0.This assumption holds for 0 < γ < 12when the coeﬃcients are piecewise smooth withrespect to the polygonal subdomains {Ωi} from (Aµ,ε). Indeed, for 0 < γ < 12 andLipschitz continuous domains D, Hγ(D) = Hγ0(D) by Theorem 2.4, from which it94follows, by interpolation between∑i H10(Ωi) and L2(Ω), that Hγ(Ω) is isomorphic to∑i Hγ(Ωi). Since ε is piecewise smooth, multiplication by ε is a bounded operatoron∑i Hγ(Ωi). For smooth coeﬃcients, one can take γ = 1.Remark 5.4 The boundedness of multiplication by a function ε in general Sobolevspaces is characterized in Corollary 1.1 from [54]. For example, if ε ∈ Hs(Ω), s > 12then we can take γ0 = s, if d = 2, and γ0 = s− 12 , if d = 3.In the next result, we characterize the rate with which Qh,k approximates Qk inthe space L(L2(Ω),H−γ(Ω)).Lemma 5.1 Let s ∈ [0, 1] be such that (5.18) and (AL2∆Dirε ,∆Neuµ ) hold. For any γ ∈[0, 1] and k ∈ {1, 2} there exists C = C(s) independent of h, such that‖(Qk − Qh,k)gk‖−γ ≤ C hsγ‖gk‖ , ∀gk ∈ L2(Ω) . (5.23)If we assume the stronger shift theorem (A∆Dirε ,∆Neuµ ), then‖(Qk − Qh,k)gk‖−γ ≤ C hmin(γ,s)‖gk‖ , ∀gk ∈ L2(Ω) . (5.24)Proof We concentrate on the case k = 1. Let w = (Qh,1 − Q1)g1. Recall thatQ1g1 = ∇u where u ∈ H10(Ω) satisﬁes (2.18) with 〈f, θ〉 = (εg1,∇θ). In addition,Qh,1g1 = ∇uh where uh is the elliptic projection of u into Hh,2, i.e., w = ∇(u−uh).Now,‖∇(u− uh)‖ ≤ C‖g1‖.Furthermore, by ﬁnite element duality and (AL2∆Dirε ,∆Neuµ ),‖∇(u− uh)‖−1 ≤ ‖u− uh‖ ≤ Chs‖u‖1 ≤ Chs‖g1‖.By interpolation,‖w‖−γ = ‖∇(u− uh)‖−γ ≤ Chsγ‖g1‖.95If we assume (A∆Dirε ,∆Neuµ ), then by ﬁnite element duality‖∇(u− uh)‖−s ≤ ‖u− uh‖1−s ≤ Chs‖u‖1.This implies (5.24) in each of the cases γ ≤ s and γ > s.We now can formulate the main result of this section, which states that Sh,kapproximates Sk in norm.Theorem 5.4 Let γ and s be two numbers in [0, 1], such that (Aε,µ×), (5.18), (5.19)and (AL2∆Dirε ,∆Neuµ ) hold. Then, there is a positive constant C = C(γ, s) independent ofh such that for k = 1, 2,‖Sk − Sh,k‖ ≤ Chsγ.Here ‖ · ‖ is the operator norm in L(L2(Ω),L2(Ω)).If the stronger shift theorem (A∆Dirε ,∆Neuµ ) holds, then we get the improved estimate‖Sk − Sk,h‖ ≤ Chmin(γ,s).Proof We consider k = 1. The case of k = 2 is similar.We are going to apply Corollary 3.5 of Theorem 3.3 with Xˆ = Hs(Ω) and Yˆ =L2(Ω). In this case, the condition (3.35) follows by combining (4.12) and (5.19). Wealso have χ(h) ≈ hs, and α(h) = 0 for both least-squares approximation methods.Due to the weight in the projectors, the proof of the theorem should be slightlymodiﬁed. Speciﬁcally, for g1 ∈ L2(Ω) we need the following additional estimate basedon (Aε,µ×):|(ε(Q1 − Qh,1)g1,Th,1Ah,1(xh− ζh))| ≤ ‖(Q1 − Qh,1)g1‖−γ‖Th,1Ah,1(xh− ζh)‖1 .This allows us to obtain the estimate (3.36):‖S− Sh‖L2(Ω)→L2(Ω) ≤ C (χ(h) + γ(h) + α(h))96with γ(h) = ‖Qk − Qh,k‖L2(Ω)→H−γ(Ω)). Now the result follows from Lemma 5.1.An alternative presentation of the proof, given in [20], proceeds as follows: ﬁxg1 ∈ L2(Ω), let x1 and xh,1 be the solutions of (5.8) and (5.20), respectively, withdata〈f1, (v, h)〉 = (ε(I− Q1)g1, v).Then,S1g1 − Sh,1g1 = x1 − xh,1 + Rh,1(Qh,1 − Q1)g1 .Here Rh,1 denotes the operator on L2(Ω) deﬁned by Rh,1w = xh,1 where xh,1 solves(5.20) with data 〈f1, (v,h)〉 = (εw, v). By Theorems 4.7 and 4.10‖x1 − xh,1‖ ≤ C hs ‖ε(I− Q1)g1‖ ≤ C hs ‖g1‖.Let w = (Qh,1 − Q1)g1. To complete the proof we only need to estimate ‖Rh,1w‖.Recall that ‖Ah,1x‖Y∗h,1 is equivalent to the norm ‖x‖Xh , uniformly in h. Thus,‖Rh,1w‖2 ≤ C (εw,Th,1Ah,1Rh,1w) ≤ C ‖w‖−γ‖εTh,1Ah,1Rh,1w‖γ≤ C‖w‖−γ‖Th,1Ah,1Rh,1w‖1.Moreover,‖Th,1Ah,1Rh,1w‖1 ≤ ‖Ah,1Rh,1w‖Y∗h,1 ≤ C‖Rh,1w‖.The result follows from Lemma 5.1.C. The eigenvalue and eigenvector discretizationIn this section, we deﬁne and analyze an approximation to the original Maxwelleigenvalue problem (1.7). As previously observed, this reduces to approximating theeigenvalues and eigenvectors for either of the symmetric semi-deﬁnite operators S2S1or S1S2. We could directly use the discrete operators Sh,k, k = 1, 2. However, this97will be avoided for two reasons. First, one would have to code both Sh,1 and Sh,2.In addition, even though the product of the continuous operators is symmetric, theproduct of their discrete counterparts is not likely to be symmetric.We circumvent the above mentioned problems by implementing only one of thediscrete operators, e.g., Sh,1. Then, instead of implementing Sh,2, we implement theadjoint S∗h,1 of Sh,1 considered as an operator of L2ε(Ω) into L2µ(Ω). The implementa-tion of S∗h,1 = ε−1Sth,1µ is relatively straightforward given the implementation of Sh,1.Indeed, Sh,1 is implemented as a sequence of matrix operations and the implementa-tion of Sth,1 just reduces to transposing the matrix operations, and running them inreverse order. Note that S∗h,1Sh,1 is symmetric by deﬁnition.The symmetry of the approximation is an important property. This is because re-alistic computations for three dimensional electromagnetic devices necessarily involveminimal problem sizes on the order of 106 unknowns. The eigenvalues and eigenvec-tors of such systems cannot be computed by direct methods. As we mentioned, itis often of interest to compute a block of the smallest eigenvalues and eigenvectorsof (5.1) [60, 88]. This means that we are required to iteratively compute the largesteigenvalues and their corresponding eigenvectors for the problem S∗h,1Sh,1x = τ2x.The problem of iteratively computing the largest eigenvalues of a symmetric positivesemi-deﬁnite problem has been well studied, see, for example, [63, 51, 65]. Even blockversions of the power method work, although not as well as other iterative strategies.A survey of iterative methods for eigenvalue problems can be found in [64].By Theorem 5.4, Sh,1 converges to S1 in norm. It immediately follows that S∗h,1converges to S∗1 = S2. It follows from the identityS∗h,1Sh,1 − S2S1 = (S∗h,1 − S2)Sh,1 + S2(Sh,1 − S1)that S∗h,1Sh,1 converges to S2S1 in norm. By standard perturbation theory, see Theo-98rem 3.16, as well as IV-§3.5 in [62], one can conclude that if τ 2 > 0 is an eigenvalueof S2S1 of multiplicity k and ν > 0 is given, such that there are no other eigenvaluesin the interval δ = (τ 2 − ν, τ 2 + ν), then for h small enough there will be exactly kdiscrete eigenvalues {τ 2i (h)}ki=1 (counted up to multiplicity) in δ. Thus, there will beno spurious discrete eigenvalues.Alternatively, we can use Sh,1S∗h,1 to approximate S1S2, Sh,2S∗h,2 to approximateS2S1, and S∗h,2Sh,2 to approximate S1S2. The analogous results for eigenvalue/eigenvectorconvergence follow for these operators as well.Using the general results for spectral approximation of compact operators (seee.g. [24, 82, 8]), we get that there is a constant C = C(τ) > 0, such that if V isthe eigenspace corresponding to τ 2, and Vh is the eigenspace corresponding to theeigenvalues of S∗h,1Sh,1 in δ, then for small enough hδˆ(V, Vh) ≡ supv∈V,‖v‖=1dist(v, Vh) ≤ C ‖S2S1 − S∗h,1Sh,1‖ . (5.25)The quantity δˆ(V, Vh) is called the “gap” between V and Vh. It is a measure forcloseness of subspaces which, in this case, is related the angle between them 5. Fur-ther details and results concerning δˆ can be found in [62], pp. 197–198. Relatedestimates demonstrating that each orthonormal basis of V can be approximated byan orthonormal basis of Vh, with the same rate, are given in [24], pp. 532–533.Combining (5.25) with Theorem 5.4, we obtain the following convergence resultfor the eigenvectors.Theorem 5.5 Let ω > 0 be ﬁxed, such that λ = iω is an eigenvalue of (5.1). Let5In fact δˆ(V, Vh) = sin(θ), where θ is the (acute) “angle” between the two spaces,i.e. the maximum of the angles between elements of V and their orthogonal projec-tions on Vh. This relation can be used to compute the rate of approximation of theeigenspaces, see [66].99τ = ω−1, and V , Vh are the eigenspaces deﬁned above. Then, for small enough h,there is a positive constant C = C(ω) independent of h such that,δˆ(V, Vh) ≤ Chsγ .Regarding the eigenvalues, the general theory states that there exists a constantC = C(τ) > 0, such that if h is small enough|τ 2 − τ 2i (h)| ≤ C ‖S2S1 − S∗h,1Sh,1‖ ,for all i = 1, . . . , k. Thus, in general, we get the following convergence result for theeigenvalues.Theorem 5.6 Let ω > 0 be ﬁxed, such that λ = iω is an eigenvalue of (5.1). Letτ = ω−1, and {τ 2i (h)}ki=1 are the eigenvalues deﬁned above. Then, for small enough h,there is a positive constant C = C(ω) independent of h such that for all i = 1, . . . , k,|τ 2 − τ 2i (h)| ≤ Chsγ .1. Improved estimate of the eigenvalue convergence rate for smooth eigenfunctionson a convex domainThe Theorems 5.5 and 5.6 imply that the convergence rate of the eigenvectors andeigenvalues are the same. Our numerical results however, indicate that sometimesthe rate of convergence of the eigenvalues is signiﬁcantly better than the rate ofconvergence of the eigenvectors. Below, we outline a proof of this fact in the case of“smooth” eigenvectors.For the remainder of this subsection, we assume that Ω is a convex polyhedron,ε = µ = 1, Th,k corresponds to a direct solve (not a preconditioner) and the eigen-100vectors are such that e · n ∈ H 32 (F ) on each face F of ∂Ω 6. This is the case, forexample, if the domain is the unit cube. By Theorem 5.4 we have ‖Sk − Sh,k‖ ≤ Ch,for k = 1, 2.Fix an eigenvector e of S2S1 corresponding to an eigenvalue τ2 and let 0 < < 12.We will prove that the approximation of τ 2 converges at rate at least h2−.Consider the biharmonic problem∆2ψ = 0 in Ω,ψ = 0 on ∂Ω,∂ψ∂n= θ on ∂Ω.(5.26)with data θ = e · n. By our assumptions, θ ∈ H 32 (F ) and θ = 0 on ∂F on every faceF of ∂Ω. By examination of the proof of the regularity result from [55], one can showthat this implies ψ ∈ H3−(Ω).Set w = τ−2e+ ∇∆ψ and consider the div-curl system∇×v = w in Ω,∇·v = 0 in Ω,v · n = 0 on ∂Ω.(5.27)By construction, w is in H−(Ω) and satisﬁes the compatability conditions, so theabove problem is well-posed. Moreover, we show in Appendix A that the solution isin H1−(Ω) and there exist C > 0 such that ‖v‖1− ≤ C‖w‖.Deﬁne T1 : H−1(Ω) → H10(Ω) by(∇T1,∇z) = 〈, z〉 ∀z ∈ H10(Ω) .6By Theorem 3.10 from [75] this implies that e · n can be extended to a functionin Hs(Ω) for any 32< s < 2.101We claim that∇×T1∇×v = ∇×e . (5.28)Indeed, e−∇ψ ∈ H10(Ω) by (5.26), and thereforee−∇ψ = T1(−∆(e−∇ψ)) = T1(τ−2 e+ ∇∆ψ)) .The result follows by applying the curl operator to both sides.Let τ 2h and eh be the eigenvalue and eigenvector approximations to τ2 and e,respectively. Set u = (e, τ−1S1e)t and uh = (eh, τ−1h Sh,1eh)t. We assume that u anduh are scaled so that ‖u‖ = ‖uh‖ = 1 where ‖ · ‖ denotes the square root of the sumof the squares of the L2(Ω)-norms on the two components. We then have, see Remark5.1, B˜u = τu and B˜huh = τhuh whereB˜ ≡⎛⎜⎝ 0 S2S1 0⎞⎟⎠ and B˜h ≡⎛⎜⎝ 0 S∗h,1Sh,1 0⎞⎟⎠ .Simple algebraic manipulations show thatτ − τh =((τI− B˜)(u− uh),u− uh)− ((B˜− B˜h)(u+ uh),u− uh)+ ((B˜− B˜h)u,u) .Note that eigenvector convergence implies that‖u− uh‖ ≤ Ch.In addition, ‖B˜− B˜h‖ ≤ Ch, so it will be enough to get a higher order bound for theterm ((B˜− B˜h)u,u).Let x1 = S1e and xh,1 = Sh,1e. Then ((B˜ − B˜h)u,u) = 2(x1 − xh,1, h˜) whereh˜ = τ−1S1e = τ∇×e.102Introduce AVh,1 as the map of Xh into V∗h,1 deﬁned by〈AVh,1x, vh〉 = (x,∇×vh) ∀vh ∈ Vh,1.Similarly, let TVh,1 : V∗h,1 → Vh,1 be deﬁned by(∇TVh,1,∇vh) = 〈, vh〉 ∀vh ∈ Vh,1 . (5.29)We assume that TVh,1 is used to deﬁne the Vh,1 component in the deﬁnition of Sh,1. Itfollows that for any vh ∈ Xh, Th,1Ah,1vh consists of two components TVh,1AVh,1vh andTHh,1AHh,1vh where THh,1 is the H1 part of Th,1 and AHh,1 is the H1 part or Ah,1, i.e.,〈AHh,1x, ψh〉 = (x,∇ψh) ∀ψh ∈ Hh,1.The deﬁnition of xh,1 states that(xh,1,∇×TVh,1AVh,1vh) + (xh,1,∇THh,1BHh,1vh) = (e,TVh,1AVh,1vh), ∀vh ∈ Xh. (5.30)Note that we used the fact that Qh,1e = 0 above.Using (5.30), the deﬁnition of x1 and (5.28) gives(x1 − xh,1, h˜) = τ(x1 − xh,1,∇×T1∇×v−∇×TVh,1AVh,1vh)− τ(x1 − xh,1,∇TVh,1AVh,1vh),(5.31)for any vh ∈ Xh. The ﬁrst term in (5.31) can be estimated byC h ‖e‖{‖TVh,1 (∇×v−AVh,1vh) ‖1 + ‖ (T1 − TVh,1)∇×v‖1} .We then have‖TVh,1(∇×v−AVh,1vh) ‖1 ≤ supφ∈Vh,1〈∇×v−AVh,1vh,φ〉‖φ‖1 = supφ∈Vh,1(v− vh,∇×φ)‖φ‖1≤ infvh∈Xh‖v− vh‖ ≤ Ch1−‖v‖1−103and‖(T1 − TVh,1)∇×v‖1 ≤ Ch1−‖∇×v‖− ≤ Ch1−‖v‖1−.Finally, the second term in (5.31) is the same asτ(x1 − xh,1,∇TVh,1(∇ · v−AVh,1vh)) ≤ Ch supφ∈Hh,1((∇ · v−AVh,1vh), φ)‖φ‖1= Ch supφ∈Hh,1(v− vh,∇φ)‖φ‖1 ≤ Ch2−‖v‖1−.Combining the above results we conclude that |(x1 − xh,1,h)| ≤ C(τ)h2− for any0 < < 12, and therefore, we proved the following improved convergence estimate.Theorem 5.7 Assume that Ω is a convex polyhedron, ε = µ = 1, Th,1 is deﬁned interms of the direct solve (5.29), and the eigenvectors are such that e · n ∈ H 32 (Γ) foreach face Γ of ∂Ω.Let λ = iω be a ﬁxed eigenvalue of (5.1), τ 2 = ω−2, and {τ 2i (h)}ki=1 be theeigenvalues of S∗h,1Sh,1 that are approximation of τ2. Fix 0 < < 12. Then there existsa positive constant C = C(λ) independent of h such that for all i = 1, . . . , k,|τ 2 − τ 2i (h)| ≤ Ch2− .D. Extensions to more general domainsIn this section, we discuss the modiﬁcations necessary to deal with curved domainsor non-simply connected domains with holes.We ﬁrst consider the case Ωh ⊂ Ω, in the settings of §IV.C.3.a. By the theorydeveloped there, and presented in Theorem 4.11, we know that the discrete solutionsof the div-curl systems on Ωh approximate the solutions on Ω with order hs, s < 12,provided that the right-hand side is compatible. To extend this result to the eigenvalueproblem, it is enough to obtain an upper bound for ‖Qk − Qh,k‖L2(Ω)→H−γ(Ω)), where104Qk is deﬁned on Ω, while Qh,k is deﬁned on Ωh.Fix gk ∈ L2(Ω). For simplicity, we assume that Qh,k is deﬁned by using onlypiecewise linear functions, see Remark 5.3, in which case the extension of Qh,kgk toΩ is trivial. Furthermore, we still have that Qh,kgk is the elliptic projection of Qkgk,and therefore, the term ‖(Qk − Qh,k)gk‖H−γ(Ωh)) can be estimated by Lemma 5.1.Thus, it remains to estimate the term ‖Qkgk‖H−γ(ω) which is done below.Lemma 5.2 Let 0 ≤ γ < 12. There exists C ∈ R+ independent of h, such that‖Qkgk‖H−γ(ω) ≤ Chγ‖gk‖ .for any gk ∈ L2(Ω).Proof Let E0 be the extension by zero from L2(ω) to L2(Ω). By Theorem 2.4, this isa bounded operator from Hγ(ω) to Hγ(Ω). First, consider the case ε = µ = 1. Usingthe deﬁnitions we get‖Qkgk‖H−γ(ω) = supw∈Hγ(ω)(Qkgk,w)L2(ω)‖w‖Hγ(ω) = supw∈Hγ(ω)(Qkgk, E0w)L2(Ω)‖E0w‖Hγ(Ω)≤ Chγ supw∈Hγ(ω)(gk,QkE0w)L2(Ω)‖w‖L2(ω) ≤ Chγ‖gk‖ ,where we applied the estimate (4.58) in the form ‖E0w‖Hγ(Ω) ≤ C hγ‖w‖L2(ω). Thecase of piecewise smooth ε and µ presents no additional diﬃculties, since the operatorsof multiplication by ε, µ, ε−1 and µ−1 are bounded from Hγ(ω) to Hγ(ω).We summarize the above considerations in the next result.Theorem 5.8 Let s ∈ [0, 12] be such that (5.18), (5.19) and (A∆Dirε ,∆Neuµ ) hold. Letω > 0 be ﬁxed, such that λ = iω is an eigenvalue of (5.1). Let τ 2 = ω−2 be aneigenvalue of S2S1 of multiplicity k and ν > 0 is given, such that there are no othereigenvalues in the interval δ = (τ 2 − ν, τ 2 + ν). Then, for small enough h, there will105be exactly k discrete eigenvalues {τ 2i (h)}ki=1 of S∗h,1Sh,1 (counted up to multiplicity) inδ. Furthermore, there exists C = C(ω) independent of h such that for all i = 1, . . . , k,|τ 2 − τ 2i (h)| ≤ Chs .Let V be the eigenspace corresponding to τ 2 and Vh be the eigenspace correspond-ing to the eigenvalues of S∗h,1Sh,1 in δ. Then, there is a positive constant C = C(ω)independent of h such thatδˆ(V, Vh) ≤ Chs .This completes the analysis for domains with curved boundary. Next we considerthe case Ω = Ωh with n1 > 0, n2 > 0. As in §IV.C.3.b, the only essential diﬀerenceis that we shall have to increase the spaces H1 and H2 with an analogous increase intheir discrete counterparts.Speciﬁcally, we deﬁne H1 and H2 byH1 = H1(Ω)⊕ K1(µ) and H2 = H10(Ω)⊕ K2(ε) , (5.32)where K1(µ) and K2(ε) have bases deﬁned in (4.63) and (4.64). Note that if (e,h)are eigenfunctions of (5.1) corresponding to λ = 0, then by density (Theorem 2.5)and (4.61) we have(∇×h,∇ψi) = 0 = λ(e,∇ψi) ∀ψi ∈ K2(ε) ,(∇×e, ∇˜ζj) = 0 = −λ(µh, ∇˜ζj) ∀ζj ∈ K1(µ) .Therefore, the eigenvectors of (5.1) with nonzero eigenvalues still satisfy (5.12) withSk deﬁned using (5.32). Furthermore, Sk still satisfy (5.9) and (5.10), so the rest ofthe proof of Theorem 5.2 goes through as well.The rest of the analysis in Sections A and B does not need to be changed. As aresult, we get the same theorems for convergence of the eigenvectors and eigenvalues106as the one presented in Section C.107CHAPTER VITHE TIME-HARMONIC MAXWELL SYSTEMIn this chapter, we consider the full time-harmonic system (1.4). We will analyze amore general version in which, given the data j, m and the number λ = −i ω, ω ∈ R,we are looking for the magnetic and electric ﬁelds h ,e : Ω→ R3 satisfying⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩∇×h = λ εe+ j in Ω,∇×e = −λµh+m in Ω,µh · n = 0 on ∂Ω ,e×n = 0 on ∂Ω.(6.1)The ﬁeld m represents the source magnetic current density.For simplicity, we consider only the case Ω = Ωh which is simply connected witha connected boundary, i.e. n1 = 0, n2 = 0. The extension to the more general casesfollows the previously presented theory in §IV.C.b and §V.D.As in Chapter V, we assume λ = 0 and formally obtain the following divergenceequations1 ⎧⎪⎨⎪⎩∇ · (µh) = λ−1∇·m in Ω,∇ · (εe) = −λ−1∇·j in Ω .(6.2)Now, a standard interpretation of (6.1)-(6.2) is to assume j ,m ∈ H(div) and to lookfor h ∈ X1(µ) and e ∈ X2(ε). We call this the original form of the time-harmonicproblem.In some problems, the electric and magnetic charge densities vanish in the whole1When λ = 0, the problem splits into independent magnetostatic and electrostaticproblems which have been already analyzed.108domain, see [50, 49]. This means that the following compatability conditions hold∇·j = 0 , ∇·m = 0 , m · n = 0 . (6.3)Usually in practice m = 0. Then, a popular reformulation of the above systemsis to reduce them to a curl-curl problem for the electric ﬁeld e. Speciﬁcally, byeliminating the magnetic ﬁeld, we get∇×µ−1∇×e = ω2 εe+ j˜ , (6.4)where j˜ = −λ j, e×n = 0 and ω2 = −λ2. The introduction of j˜ is in accordance with(1.10a)–(1.11b) from [75]. The weak formulation of (6.4) reads: Find e ∈ HC0 (curl)such that(µ−1∇×e,∇×w) = ω2 (εe,w) + (˜j,w) ∀w ∈ HC0 (curl) . (6.5)Similar to the eigenvalue problem, since ω ∈ R, we can split the above problem intotwo real problems.Note that (h,e, j,m) is solution to (6.1) with λ = −i ω, ω ∈ R if and onlyif ((h),(e),(j),(m)) and (−(h),(e),−(j),(m)) satisfy the related realproblem ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩∇×h = ω εe+ j in Ω,∇×e = ω µh+m in Ω,µh · n = 0 on ∂Ω ,e×n = 0 on ∂Ω ,(6.6)with corresponding divergence equations⎧⎪⎨⎪⎩∇ · (µh) = −ω−1∇·m in Ω,∇ · (εe) = −ω−1∇·j in Ω,(6.7)109In the implementation, one may prefer to restrict to real ﬁelds and consider (6.6)-(6.7) instead of (6.1)-(6.2). In doing so, one would be able to avoid using any complexarithmetic.There have been a variety of methods for approximation of the time-harmonicMaxwell’s equations. See, for example, [39, 40, 85, 50, 49, 70, 11, 95, 75]. Most ofthese discretization methods are based on the curl-conforming Ne´de´lec spaces (cf.[77, 78]). In this case, the resulting matrix problem is indeﬁnite and it is well knownthat the eﬃcient solution of such systems presents a serious challenge.Recently, an interior penalty discontinuous Galerkin method was analyzed in[85]. This method allows for diﬀerent orders of approximation in the diﬀerent regionsof the grid. Error estimates were proven in the case of smooth coeﬃcients. However,we remark that the resulting bilinear form is quite complicated.A diﬀerent approach proposed in [50] (where they also address the more generalproblem of regions with screens) is based on the singular function method. Themethod uses the splitting of the solution into regular part, which can be approximatedby nodal ﬁnite elements, and a singular part which is treated explicitly by adding thesingular functions to the space. The implementation of this procedure may be quiteinvolved, since one needs to deal explicitly with the singular functions.In some applications, e.g. in the mixed digital and analog signal packages, oneneeds a methodology that is independent of the frequency ω. In [49], such resultswere obtained for the case of ω in a neighborhood of zero. The approach there usesa mixed formulation with a “dummy” Lagrange multiplier (one that is identicallyzero). Further results in this direction are given in [70], where a FOSLS method isapplied to the scalar Helmholtz equation with exterior radiation boundary conditions.The convergence of the resulting Multigrid algorithm is uniform with respect to thewave number ω, under the assumption that the domain is convex or has a smooth110boundary.We ﬁnish the overview of the available literature on the subject by mentioningsome other approaches as the mortar and the FETI methods applied to the Maxwell’sequations, see [11, 95]. Adaptive hp solvers have also been implemented, see [49].In this chapter, we describe a new approximation technique for the time-harmonicsystem, which is a natural extension of the ideas from the previous chapters. The mainadvantage of this approach is that it directly approximates the variables of interest,avoiding potentials and “gauge conditions”. In fact, we approximate simultaneously hand e. This is in contrast to many methods where one of the unknowns is eliminatedand is later computed by diﬀerentiation of the approximate solution. Furthermore,the resulting numerical algorithm can be eﬃciently implemented and is convergentfor problems with low regularity which appear in practical applications.In the next sections, we present the weak variational formulation of (6.1) andrelate it to the eigenvalue problem from Chapter V. We then discuss discretizationbased on a least-squares method similar to the ones developed for the magnetostaticand electrostatic problems in Chapter IV.A. Weak formulationThe following result for the original time-harmonic problem in the form (6.5) is givenas Corollary 4.19 in [75].Theorem 6.1 Suppose that ω2 = 0 is not a Maxwell eigenvalue, i.e. does not satisfy(5.4). Then the curl-curl problem (6.5) has a unique solution e for any data j˜ ∈L2(Ω), and we have the stability estimate‖e‖H0(curl) ≤ C ‖j˜‖ . (6.8)Remark 6.1 Assume that j ∈ H(div) and m ∈ H0(div) satisfy (6.3). Let f and g111be the unique solutions to the div-curl problems⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×f = j in Ω,∇ · (µf) = 0 in Ω,f · n = 0 on ∂Ω,⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×g = m in Ω,∇ · (εg) = 0 in Ω,g×n = 0 on ∂Ω,(6.9)i.e.B⎛⎜⎝ jm⎞⎟⎠ =⎛⎜⎝−gf⎞⎟⎠ .Then, analogous to §V.A, it can be shown that (e,h) is a solution to (6.1)-(6.2) withdata (j,m) and λ = 0 if and only if h ∈ L2µ(Ω), e ∈ L2ε(Ω), f, g and τ = λ−1 satisfyB⎛⎜⎝eh⎞⎟⎠− τ⎛⎜⎝eh⎞⎟⎠ = −τ⎛⎜⎝gf⎞⎟⎠ . (6.10)By the Fredholm Alternative (Theorem 2.3), when τ is not an eigenvalue of B, theabove problem has unique solution.As a natural extension of the weak formulations in Chapter IV, we proposeto replace the diﬀerential operators in (6.1)-(6.2) with the weaker operators curl1,curl2, div1,µ and div2,ε deﬁned in (4.2) and (4.18). Then we can consider an evenmore general problem: given j ∈ V∗1, m ∈ V∗2, ρ ∈ H∗1 and q ∈ H∗2 ﬁnd h ,e ∈ L2(Ω)satisfying ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩curl1h− λεe = jcurl2e+ λµh = mdiv1,µh = qdiv2,εe = ρ(6.11)112Let Aλ : L2(Ω) × L2(Ω) → Y∗1 × Y∗2 denote the operator on the left. Then our weakformulation of the time-harmonic problem isAλ⎛⎜⎝eh⎞⎟⎠ ≡⎛⎜⎜⎜⎜⎜⎜⎜⎝curl1h− λεecurl2e+ λµhdiv1,µhdiv2,εe⎞⎟⎟⎟⎟⎟⎟⎟⎠=⎛⎜⎜⎜⎜⎜⎜⎜⎝jmqρ⎞⎟⎟⎟⎟⎟⎟⎟⎠. (6.12)First, we note that any solution to the original time-harmonic problem satisﬁes(6.12) with q = λ−1∇·m and ρ = −λ−1∇·j. On the other hand, if j ,m ∈ L2(Ω),then (6.12) implies (6.1). In particular,Aλ(e,h) = 0 if and only if (h,e, λ) satisfy (5.1) . (6.13)It follows that if λ is not a Maxwell eigenvalue, then the kernel of Aλ is trivial. Infact, a stronger statement is true.Lemma 6.1 Assume that (Aε,µ×) holds for some γ < 12 . Then Aλ is bounded frombelow, i.e.‖e‖+ ‖h‖ ≤ C ‖Aλ(e,h)‖Y∗1×Y∗2 , (6.14)provided that λ = 0 is not a Maxwell eigenvalue, i.e. does not satisfy (1.7).Proof Assume that (6.14) does not hold. Then there exist a sequence {xn =(en,hn)} ⊂ L2(Ω)×L2(Ω) such that ‖xn‖2 ≡ ‖en‖2+‖hn‖2 = 1, while ‖Aλxn‖2Y∗1×Y∗2 ≤1n. Using the compact embedding L2(Ω) ↪→ H−γ(Ω) and passing to a subsequence,we get hnH−γ−→ h and en H−γ−→ e for some h ,e ∈ H−γ(Ω). Since γ < 12, we also havethe continuous embeddings ‖v‖V∗k ≤ C ‖v‖−γ for k = 1, 2. In particular, h ∈ V∗1 ande ∈ V∗2. Note that (Aε,µ×) implies‖µhn‖V∗1 ≤ C supw∈Hγ|(hn, µw)|‖w‖γ ≤ Cµ supv∈Hγ|(hn, v)|‖v‖γ = ‖hn‖−γ (6.15)113for any hn ∈ L2(Ω). By (4.9) and (4.25), for m,n ∈ N,‖xm−xn‖2 ≤ Cε,µ{‖Aλ(xm − xn)‖2Y∗1×Y∗2 + ω2‖µhm − µhn‖2V∗1 + ω2‖εem − εen‖2V∗2}.Using (6.15), ‖Aλ(xm − xn)‖2Y∗1×Y∗2 ≤ 2(1n+ 1m)and the above inequality we get that{xn} is Cauchy in L2(Ω)×L2(Ω). Therefore, we can conclude that x = (e,h) ∈ L2(Ω)2and xnL2(Ω)2−→ x.After passing to a limit, we get ‖x‖ = 1, while Aλx = 0. Since ω is not aneigenvalue, (6.13) implies x = 0, which is a contradiction.To characterize the solvability of the weak formulation (6.11), we need to deter-mine what are the compatability conditions on the data. To that end, let us considerthe bilinear form aλ(·, ·), corresponding to Aλ, which is deﬁned on L2(Ω)2× (Y1×Y2)byaλ(h,e;v1, v2,h1,h2) =a1(h; v1,h1)− λ(εe, v1) + a2(e; v2,h2) + λ(µh, v2) ,(6.16)where the forms ak(·, ·) were introduced in (4.17) and (4.33).Let (v1, v2,h1,h2) belong to the compatability spaceN(A∗λ) = {(v1, v2,h1,h2) ∈ Y1×Y2 : aλ(h,e; v1, v2,h1,h2) = 0 ∀(h,e) ∈ L2(Ω)2} .Then ⎧⎪⎨⎪⎩∇×v1 + λµ v2 + µ∇h1 = 0∇×v2 − λ ε v1 + ε∇h2 = 0 .Set v˜1 = −v1 + λ−1 ∇h2 and v˜2 = v2 + λ−1 ∇h1. It follows that v˜1 ∈ X2(ε) andv˜2 ∈ X1(µ) satisfy the eigenvalue problem (5.1). This impliesv1 = λ−1 ∇h2 and v2 = −λ−1 ∇h1 . (6.17)114Conversely, if (6.17) holds then (v1, v2,h1,h2) ∈ N(A∗λ). The above considerationsprove the following result.Lemma 6.2 Assume that λ is not a Maxwell eigenvalue. Then the compatabilityspace for (6.12) is given byVλ,0 ≡ N(A∗λ) = {(∇h2,−∇h1, λh1, λh2) : h1 ∈ H2(Ω) ,h2 ∈ H20(Ω)} . (6.18)Consequently, the data (j,m, q, ρ) are compatible if and only if〈ρ,h2〉 = λ−1〈j,∇h2〉 , 〈q,h1〉 = −λ−1〈m,∇h1〉 , (6.19)for all h2 ∈ H20(Ω), h1 ∈ H2(Ω). When j ∈ H(div), m ∈ H0(div), the aboveconditions simplify toρ = −λ−1∇·j , q = λ−1∇·m . (6.20)We combine the results of Lemma 6.1, Lemma 6.2 and Proposition 3.1 in themain result of this section.Theorem 6.2 Assume that (Aε,µ×) holds for some γ < 12 , and λ = 0 is not aMaxwell eigenvalue. Then the weak formulation (6.12) has unique solution for anydata satisfying the compatability conditions (6.19), and the following stability estimateholds‖h‖+ ‖e‖ ≤ C ‖j‖V∗1 + ‖w‖V∗2 + ‖q‖H∗1 + ‖ρ‖H∗2 .When j ∈ H(div), m ∈ H0(div) and ρ and q are deﬁned by (6.20), the weak solutioncoincides with the solution of the original time-harmonic problem (6.1)-(6.2). Inaddition, if X1(µ) and X2(ε) are compactly embedded in Hs(Ω) for some s > 0, thenwe have the stability estimate‖h‖X1(µ) + ‖e‖X2(ε) ≤ C ‖j‖H(div) + ‖w‖H0(div) .115Proof We only need to show how to get the second stability estimate. Using (4.12)and (4.28) for any h ∈ X1(µ) and e ∈ X2(ε) we getC (‖e‖2X1(µ) + ‖h‖2X2(ε)) ≤ ‖∇×h− λεe‖2 + ‖∇·µh‖2+ ‖∇×e+ λµh‖2 + ‖∇·εe‖2+ ω2 ‖µh‖2 + ω2 ‖εe‖2 .By a compactness argument, identical to the one in the proof of Lemma 6.1, we canconclude thatC (‖e‖X1(µ) + ‖h‖X2(ε)) ≤ ‖∇×h− λεe‖+ ‖∇·µh‖+ ‖∇×e+ λµh‖+ ‖∇·εe‖for any λ = 0 which is not a Maxwell eigenvalue. This implies that the operatorcorresponding to the original time-harmonic problem is bounded from below, andhence, the stability estimate follows from Proposition 3.1.Remark 6.2 Theorem 6.1 is a special case of the above result.B. Least-squares approximationWe next consider the discrete approximation to our weak formulation. As in theprevious chapters, see §IV.C, we choose approximation subspaces Xh,1 = Xh,2 =Xh ⊂ L2(Ω) and Yh,k ⊂ Yk.We will consider extensions of both of the discrete least-squares methods—theone based on a discrete inf-sup condition and the one based on form modiﬁcation. Forsimplicity, we concentrate on the case of real ﬁelds, i.e. we approximate (6.6)-(6.7).The weak form, in this case, is based on the operatorARλ(e,h) = Aω(e,h) = (curl1h− ωεe, curl2e− ωµh,div1,µh,div2,εe) , (6.21)116where all ﬁelds, spaces and operators are real. The corresponding bilinear form isaRλ(·, ·) = aω(·, ·) given bya1(h; v1,h1)− ω(εe, v1) + a2(e; v2,h2)− ω(µh, v2) .1. Approximation based on a discrete inf-sup conditionFirst, consider the method from §IV.C.1.a. There, for ﬁxed xk ∈ Xh,k and anyψk ∈ Hk, vk ∈ Vk, we showed that one can choose ψh,k ∈ Hh,k and vh,k ∈ Vh,k suchthat (µx1,∇ψ1) = (µx1,∇ψh,1) , (εx1,∇ψ2) = (εx2,∇ψh,2) and(xk,∇×vk) = (xk,∇×vh,k) (6.22)with ‖ψh,k‖1 ≤ C ‖ψk‖1 and ‖vh,k‖1 ≤ C ‖vk‖1. The idea was to use a stable ap-proximation operator, plus face and element bubble functions, in order to satisfy theequalities. For the time-harmonic problem, (6.22) should be replaced by(h,∇×v1)− ω (εe, v1) = (h,∇×vh,1)− ω (εe, vh,1) (6.23)for h ,e ∈ Xh.Below, we illustrate the needed modiﬁcations in the case Xh = Ŝh(k) and piece-wise constant ε and µ. Let v1 = (vc)dc=1 and vh,1 = (vch)dc=1. Set vch = Ihvc +vcFh +vcTh ,where Ih is an approximation operator satisfying (2.26) and(vcFh , q)L2(F ) = (vc − Ihvc, q)L2(F ) ∀F ∈ Fh ,∀q ∈ Pk(F ) ,(vcTh , p)L2(τ) = (vc − Ihvc − vcFh , p)L2(τ) ∀τ ∈ Th ,∀p ∈ Pk(τ) .This implies(∇×h− ω εe, v1)L2(τ) = (∇×h− ω εe, vh,1)L2(τ) ∀τ ∈ Th ,(h×n, v1)L2(F ) = (h×n, vh,1)L2(F ) ∀F ∈ Fh ,117and therefore (6.23) follows. Clearly the diﬀerence with Theorem 4.8 is that we needelement bubbles of degree k instead of k − 1. This is summarized in the next result.Theorem 6.3 Let k ∈ N0. Then the least-squares method for the time-harmonicproblem based on the spaces Xh,1 = Xh,2 = Ŝh(k), Vh,1 = (Sh,0 ⊕ BkFh,0 ⊕ BkTh)3,Hh,1 = Sh ⊕ BkFh ⊕ Bk−1Th , Vh,2 = (Sh ⊕ BkFh ⊕ BkTh)3 and Hh,2 = Sh,0 ⊕ BkFh,0 ⊕ Bk−1Th isstable (i.e. the discrete inf-sup condition for the form aω(·, ·) holds).As a corollary, the discrete least-squares method based on the form aω(·, ·) will haveunique solution, provided that λ = −i ω is not a Maxwell eigenvalue.2. Approximation based on form modiﬁcationNext, we consider the least-squares approach based on form modiﬁcation and pre-sented in §IV.C.2. The only diﬀerence here is that we have to estimate two additionalterms:supv1∈V1(ω εe, v1)‖v1‖V1and supv2∈V2(ω µh, v2)‖v2‖V2for e ,h ∈ Xh. For any v1 ∈ V1, let vh,1 ∈ Vh,1 be obtained by applying the stableapproximation operator Ih to each component of v1. Thensupv1∈V1(εe, v1)2‖v1‖2V1≤ C supv1∈V1(εe, v1 − vh,1)2‖v1‖2V1+ C supvh,1∈Vh,1(εe, vh,1)2‖vh,1‖2V1.Deﬁne projh1,ε : Xh → Vh,1 by(projh1,εe, vh,1) = (εe, vh,1) ∀vh,1 ∈ Vh,1 . (6.24)Using this deﬁnition and (2.26), we obtainsupv1∈V1(εe, v1)2‖v1‖2V1≤ C∑τ∈Thh2τ ‖εe‖2L2(τ) + C ‖projh1,εe‖2V∗h,1 .118Similarly,supv2∈V2(µh, v2)2‖v2‖2V2≤ C∑τ∈Thh2τ ‖µh‖2L2(τ) + C ‖projh2,µh‖2V∗h,2 ,where projh2,µ : Xh → Vh,2 is deﬁned by(projh2,µh, vh,2) = (µh, vh,2) ∀vh,2 ∈ Vh,2 . (6.25)Combining this with the results for the magnetostatic and electrostatic problems, weget that the least-squares problemah,ω(h,e; h˜, e˜) = (j, curlh1h˜)V∗h,1 + (m, curlh2 e˜)V∗h,2+ (q,divh1,µh˜)H∗h,1 + (ρ,divh2,εe˜)H∗h,2 ∀h˜ , e˜ ∈ Xh .will have a unique solution, provided the conditions of Theorem 6.2 are met. Thebilinear form on the left is given byah,ω(h,e; h˜, e˜) = (curlh1h, curlh1h˜)V∗h,1 + (divh1,µh,divh1,µh˜)H∗h,1+ (projh2,µh,projh2,µh˜)V∗h,2 +∑τ∈Thh2τ (µh, µ h˜)L2(τ)+∑F∈FhhF{(h×n, h˜×n)L2(F ) + (µh · n, µh˜ · n)L2(F )}+ (curlh2e, curlh2 e˜)V∗h,2 + (divh2,εe,divh2,εe˜)H∗h,2+ (projh1,εe,projh1,εe˜)V∗h,1 +∑τ∈Thh2τ (εe, ε e˜)L2(τ)+∑F∈FhhF{(e×n, e˜×n)L2(F ) + (εe · n, εe˜ · n)L2(F )} .1193. Error estimatesIn this subsection we assume that (AΩ) and (Aµ,ε) hold, and that there exists s ∈[0, 1], such that the estimate (3.21) holds with χ(h) = C hs, i.e.infxh,k∈Xh,k‖xk − xh,k‖ ≤ C hs ‖xk‖s ∀xk ∈ Hs(Ω) . (6.26)Additionally, we assume the continuous embeddings (see Theorem 2.7)X1(µ) ,X2(ε) ↪→ Hs(Ω) . (6.27)By combining the results of Theorem 3.2 and Theorem 6.2, we get the followingestimate for the approximation error of each of the methods.Theorem 6.4 Let s, γ ∈ [0, 1] be such that (6.26), (6.27) and (Aε,µ×) hold. As-sume that λ = 0 is not an eigenvalue, and let (h,e) be the solution of the time-harmonic problem with data j ,m ∈ L2(Ω), q and ρ satisfying the compatability con-ditions (6.20). Let (hh,eh) be the least-squares approximation obtained by either ofthe methods presented in the previous two subsections (for the method based on formmodiﬁcation s = 12). Then we have the error estimate‖h− hh‖+ ‖e− eh‖ ≤ Cµ,ε,ω hs(‖j‖H(div) + ‖m‖H0(div)) .120CHAPTER VIINUMERICAL RESULTSIn this chapter, we discuss some results from computer simulations with a programwhich implements the dual least-squares methods described in the previous chap-ters (see [102]). It is written in C++, in the framework of the AggieFEM ﬁniteelement library, which supports complex geometries, local reﬁnement, Multigrid pre-conditioning and OpenGL visualization. The code is based on the solvers for themagnetostatic and electrostatic problems. It works on triangular, tetrahedral andhexahedral meshes. It provides an eigenvalue solver (based on [65]), which allows forcomputations of blocks of eigenvalues, and a solver for the full time-harmonic system.A. Implementation issuesWe concentrate on the case of a simply connected domain Ωh, which is either polygonalor polyhedral. As mentioned before, the theory of the previous chapters extends totwo-dimensional problems without diﬃculty. We give speciﬁc details in the followingsections.In all of our examples, we partition the domain Ωh into a shape regular mesh,which is triangular in 2D and either tetrahedral or hexahedral in 3D. We mainlyfocus on the least-squares method with bubbles, i.e. we use piecewise constant vectorfunctions for the space Xh and piecewise linear, plus face bubble vector functions, foreach component of Yh,k. Instead of dealing with complex arithmetic, we recast eachof the problems as an equivalent real problem as discussed in the introductions of theprevious chapters.The implementation basically follows the description in §III.C. We note thefollowing speciﬁcs:121• The actions of TYh,k for k = 1, 2, are implemented using a two level algorithminvolving a Gauss-Seidel sweep over the bubble functions and V-cycle Multi-grid preconditioner for the remaining piecewise linear functions. A comparisonbetween this operator and the exact solver is given in §D.• The operators Qh,k were deﬁned using only the piecewise linear part of the testspaces, i.e. disregarding the bubbles as discussed in Remark 5.3.In the numerical tests, we start with a coarse mesh and apply few levels of uniformreﬁnement. On each mesh level, we either compute the solution by PCG or computea number of the maximal eigenvalues and eigenfunctions of S∗h,1Sh,1 a modiﬁed versionof LOBPCG. The results are reported using the following notation: level denotes thereﬁnement level, h is the mesh size, ‖e‖0 denotes the error in L2(Ω), ratio is the ratiobetween the errors on two consecutive levels, nit equals the number of iterations ofPCG/LOBPCG and N denotes the total number of unknowns.Here are some highlights for the reminder of the chapter. The connection betweenthe memory requirements of our algorithm and the geometrical characteristics of themesh is discussed in Subsection 1. An optimal conversance rate for a low regularitymagnetostatic problem is presented in §B.2. A problem with jumping coeﬃcientsand fairly anisotropic mesh is solved in §B.3. The case Ω = Ωh is considered in§C.2. The singular eigenvalue problem on a Fichera corner is compared in §C.3 withother previously available results. Finally, some of the approximate solutions fromthe problems with unknown exact solution are visualized in Appendix B.Our general conclusion from the experiments is that the new method performsquite well in a variety of applications. The eigenvalue approximation deals well withmultiple eigenvalues. Let us stress, again, that spurious eigenmodes are completelyavoided. We also conclude that LOBPCG seems to be a good choice for an eigensolver,122yielding a constant number of iterations in the tests presented.Finally, let us mention that further numerical experiments seem to suggest thatthe use of the projectors Qh,k and the stabilizing face bubble functions are essentialfor the convergence and cannot be avoided.1. Mesh characteristicsConsider a ﬁnite element mesh Th on Ωh ⊂ R3. LetV, E, F and T denote, respectively,the number of vertices, edges, faces and elements in the mesh. The Euler-Poincare´formula (see e.g §5.3 in [18]) states thatV − E + F − T = 1 + n1 − n2.The memory requirements for the discrete least-squares algorithm are directlyrelated to the above quantities. Consider, for example, the simplest method for themagnetostatic problem in three dimensions. Recall that the solution space consists ofpiecewise constants, while the test space is build of piecewise linears plus face bubblefunctions. Let M be the dimension of the test space, and N denotes the dimension ofthe solution space. Then,M = 4(V+F) and N = 3E .Suppose we have a sequence of meshes obtained by uniform reﬁnement. Then,the mesh characteristics are transformed as follows:• for Th consisting of tetrahedra(V,E,F,T) → (V+ E, 4F + 8T, 2E+ 3F+ T, 8T) .123• for Th consisting of hexahedra(V,E,F,T) → (V+ E+F+ T, 4F + 12T, 2E+ 4F+ 6T, 8T) .A typical example is shown in Figure 7.1, where we plot the ratio M/N afteruniform reﬁnement of hexahedral and tetrahedral meshes.1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 1812.95.31020level of refinmenet (l+1)test vs. solution space dimension (M/N)hexahedraltetrahedralFig. 7.1. Comparison of the dimensions of test and solution spaces after uniform re-ﬁnement.Observe that, in this model case, the dimension of the test space is always bigger,but the ratio eventually stabilizes. In particular, we can conclude that the memoryrequirements in both cases are proportional to N.This is illustrated further in Table 7.1, where we examine the case of comparableinitial tetrahedral and hexahedral meshes on the unit cube, which are subject to 6levels of uniform reﬁnement. On each reﬁnement level l, the mesh size in both casesis of order 2−l, and therefore we get comparable order of approximation. However,the above data indicate that the discretization with tetrahedral elements will requiresigniﬁcantly more memory.124Table 7.1. Mesh characteristics after uniform reﬁnement.l V F T V F Ttetrahedral mesh hexahedral mesh0 9 30 12 8 6 11 35 216 96 27 36 82 189 1632 768 125 240 643 1241 12672 6144 729 1728 5124 9009 99840 49152 4913 13056 40965 68705 792576 393216 35937 101376 327686 536769 6316032 3145728 274625 798720 262144B. The magnetostatic problemIn this section, we report the results of numerical experiments for the magnetostaticproblem (1.5). Some of the problems considered are two-dimensional, and below, wesummarize this special case of our theory. Speciﬁcally, for a polygonal domain Ω, themagnetostatic problem is: Find h ∈ X1 ≡ (L2(Ω))2 satisfying⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×h = j in Ω,∇·µh = ρ in Ω,(µh) · n = σ on ∂Ω.Here, we used the scalar curl deﬁned in (2.13).For this problem, both test spaces are scalar. In fact, we take Y1 = V1 × H1where V1 ≡ H10(Ω) and H1 = H1(Ω). The least-squares approximation satisﬁes(h,∇×w) + (µh,∇ψ) = 〈j,w〉 + (σ, ψ)∂Ω − (ρ, ψ) (7.1)for all (w, ψ) ∈ Y1. Here we used the deﬁnition of vector curl from (2.13).125As in three dimensions (see Theorem 3.2 of [54]), we have that each functionu ∈ X1 can be decomposedu = ∇×w+ µ∇ψ with (w, ψ) ∈ Y1 .Consequently (7.1) is well-posed.1. A problem with a known smooth solutionThe ﬁrst test problem is posed on the unit square and involves known smooth solution.We take µ = 1, j = 0, ρ = cos(πx) cos(πy) and σ = 0. Then the solution ish =12π(sin(πx) cos(πy), cos(πx) sin(πy)) .The numerical results on a uniform triangular mesh are presented in Table 7.2.The error behavior in (L2(Ω))2 clearly illustrates the expected ﬁrst-order convergenceTable 7.2. Numerical results for magnetostatic problem with a known smooth solution.h ‖e‖0 ratio nit N1/8 0.576961 6 2561/16 0.290813 1.98396 6 10241/32 0.145741 1.99541 6 40961/64 0.072897 1.99926 5 163841/128 0.036451 1.99984 5 655361/256 0.018226 1.99997 4 262144rate. Note that the number of iterations required to reduce the residual by a factorof 10−6 remains bounded independently of the number of unknowns.1262. Magnetostatics in a L-shaped domainFor the second example, we consider a problem on the L-shaped domain [−1, 1]2 \[0, 1] × [−1, 0]. Solutions of problems on this domain are not smooth in general. Toillustrate the typical singularity, we take j, ρ, and σ so that the solution in polarcoordinates is given byh = ∇(rβ cos(β θ)) with β = 2/3 .Note that h is only in (Hs(Ω))2 for s < 23. Therefore, we expect that a mesh reductionof a factor of two should result in an error reduction of 22/3 ≈ 1.587.This is clearly illustrated by the convergence results in Table 7.3. Again, we seeTable 7.3. Numerical results for magnetostatics in an L-shaped domain.h ‖e‖0 ratio nit N0.176777 0.223524 11 5120.0883883 0.143219 1.56072 11 20480.0441942 0.091108 1.57196 11 81920.0220971 0.057727 1.57826 11 327680.0110485 0.036492 1.58188 11 1310720.00552427 0.023038 1.58483 11 524288that the number of iterations remains bounded as the mesh size is decreased.The components of the computed approximation to the magnetic ﬁeld are shownon Figure 7.2.127Fig. 7.2. Magnetostatics in an L-shaped domain, computed magnetic ﬁeld.3. Cross-section of a magnetWe next report numerical results for a problem with jumps in the coeﬃcient µ. Weconsider the geometry given in Figure 7.3, which models the cross-section of a mag-net. This consists of a iron segment with ﬁxed magnetic permeability µ1 = 1000surrounded by an air region with permeability µ0 = 1. A uniform current of j and−j (shaded regions) is applied in the z direction. There is also a small air gap ofsize d = .01. For this problem, we do not report the error behavior as the analyticsolution is not available.µµ01-j10.3 0.8010.80.20djFig. 7.3. Cross-section of a magnet: geometry and coarse mesh.128Our goal was to illustrate the iterative convergence rate. The numerical exper-iments reported in Table 7.4 show that, even though there are large jumps in thepermeability, the iterative process still converges in relatively few iterations. It alsoshows that the method performs well, even in the case of a fairly anisotropic mesh(see Figure 7.3).Table 7.4. Numerical results for the cross-section of a magnet.hmin hmax nit N0.0316111 0.316228 9 1520.0158055 0.158114 11 6080.0079027 0.079056 12 24320.0039513 0.039528 11 97280.0019756 0.019764 13 38912Fig. 7.4. Magnetostatic in transformer, geometry and coarse mesh.1294. Magnetic ﬁeld in a transformerOur last magnetostatic example models a three-dimensional transformer. The geom-etry and the initial mesh are given in Figure 7.4. Speciﬁcally, we have an iron core,where µ = 103, and three coils, on the exterior two of which a rotational current f isapplied. We set µ = 0 and f = 0 in the rest of the region.Numerical experiments were performed on three tetrahedral meshes obtained byuniform reﬁnement. Their characteristics are listed in Table 7.5Table 7.5. Numerical results for magnetostatics in a transformer.hmin hmax V F T0.632805 4.32786 784 8302 40940.316402 2.16393 5775 65960 3275200.158201 1.08197 44757 525856 262016Diﬀerent views of the computed approximate solution are shown in Appendix Bon pages 156 and 157. As expected, we observe a magnetic ﬁeld following the ironcore.C. The eigenvalue problemIn this section, we report results from some numerical experiments with the least-squares method for the problem (5.16).We report computations involving both tetrahedral and hexahedral meshes. Al-though there are many analyses available for tetrahedral meshes using methods basedon curl conforming ﬁnite element approximations [17, 69, 75], very little has been donefor general hexahedral meshes. In contrast, our analysis easily extends to general hex-130ahedral meshes.The eigensolver that we use is based on the Locally Optimal Block Precondi-tioned Conjugate Gradient Method (LOBPCG), introduced in [65]. A very detaileddescription of LOBPCG from implementation point of view is given in [67], §8. Orig-inally, LOBPCG was designed to compute a block of few minimal eigenvalues of asymmetric and positive deﬁnite matrix with their corresponding eigenvectors. Thealgorithm uses only the action of the matrix and is based on a local optimizationof a three term recurrence, similar to the one from the Conjugate Gradient method.This produces a sequence of discrete approximation subspaces for the eigenvectors.The Rayleigh-Ritz procedure, combined with the soft-locking1 of the converged eigen-vectors, is then used to determine the approximate eigenvalues on each step. Letus recall, that the Rayleigh-Ritz method computes optimal approximation to theeigenvalues and eigenvectors of the matrix, given a trial subspace. It employs thesolution of generalized eigenvalue problem of dimension k, where k is the number ofeigenvalues we wish to compute (typically 10-20).As it was shown in Chapter V, the Maxwell eigenvalue problem reduces to com-putation of a block of few maximal eigenvalues of a symmetric and positive deﬁnitematrix. LOBPCG can be applied to that problem after a simple modiﬁcation in thegeneralized eigenvalue problem solver mentioned above. Our experience is that withthis modiﬁcation, LOBPCG is a very robust eigensolver. The number of iterationsfor our (well-conditioned) problems is usually independent of the mesh parameter h.1This means that even if an approximate eigenvector has already converged, it stillparticipates in the Rayleigh-Ritz procedure (which, in particular, can change it). Formore details, see §7 in [67]1311. Eigenvalues of the unit cubeThe ﬁrst test problem is posed on the unit cube partitioned into a uniform tetrahedralmesh. The eigenvalues and eigenfunctions of this problem can be computed exactly,see [2]. Speciﬁcally, the eigenfunctions are tensor products of trigonometric functions,and the eigenvalues are of the form {τ 2i } ={1k π2}, where k = k21 +k22 +k23 and {ki}3i=1are non-negative integers satisfying k1k2 + k2k3 + k3k1 > 0. Triplets with k1k2k3 > 0generate two linearly independent eigenfunctions.Figure 7.5 gives the eigenvalue approximation error (S∗h,1Sh,1 approximating S2S1)as a function of the number of reﬁnement levels. Observe that the method performswell with multiple eigenvalues. In addition, the eigenvalue convergence appears to bemonotone.1 2 3 4 5 6 7 8 90.010.0150.020.0250.030.0350.040.0450.050.055eigenvalue numbereigenvalue approximationexactlevel 1level 2level 3level 4Fig. 7.5. Unit cube, eigenvalue convergence.Figure 7.6 presents the same results in diﬀerent formats. On the left side, weshow the approximation in the error for each {τ 2i }. We note that the approximationbecomes slightly worse with the increase of the eigenvalue number. This is furtherexamined on the right, where we are looking at the error in three representative1321 2 3 4 5 6 7 8 910−510−410−310−2eigenvalue numberapproximation errorlevel 1level 2level 3level 41 2 3 410−510−410−310−2levelapproximation errorτ2 (1st)τ2 (5th)τ2 (9th)4−lFig. 7.6. Unit cube, approximation error.eigenvalues, τ 21 , τ25 and τ29 , on the diﬀerent levels of approximation. As expected from§V.C.1, we have almost quadratic convergence of the eigenvalues, twice the order ofapproximation of the eigenfunctions.2. Eigenvalues of the unit ballOur second example is the computation of the eigenmodes of the unit ball. Theeigenvalues and eigenfunctions are known, see §10.4 in [10], but they are not assimple as in the previous test.Speciﬁcally, the eigenvalues {ω2i } = {ω2mn, ωˆ2mn : m,n = 1, 2, ...} are split intotwo groups:• Transverse Electric (TE), which satisfyjm(ω2mn) = 0 ,and133• Transverse Magnetic (TM), which satisfyjm(ωˆ2mn) + ωˆ2mn j′m(ωˆ2mn) = 0 .Here jm is the m-th order spherical Bessel function and j′m is its derivative. They areobtained by the formulasj0(x) =sin(x)x, j1(x) =sin(x)x2−cos(x)x, . . . jn(x) = (−x)n(1xddx)n(sin(x)x).They are also related to the Bessel functions of ﬁrst kind byjn(x) =√π2xJn+ 12(x) .There are tables with the zeros of jn (e.g. in §10.1 of [1]), but there are no simpleformulas for the zeros of j′n.The numerical values for the ﬁrst few eigenvalues {ω2i }, together with their mul-tiplicities, are given in Table 7.6. We used a set of hexahedral meshes, starting withthe coarse mesh shown in Figure 7.7. Their characteristics, together with the numberof iterations of the eigensolver, are given in Table 7.7.Fig. 7.7. Unit ball, initial mesh.134Table 7.6. Unit ball, exact eigenvalues.i ω2i type multiplicity1 7.5279e+00 TM (ωˆ211) 32 1.4979e+01 TM (ωˆ221) 53 2.0191e+01 TE (ω211 ) 34 2.4735e+01 TM (ωˆ231) 75 3.3217e+01 TE (ω221 ) 56 3.6747e+01 TM (ωˆ241) 97 3.7415e+01 TM (ωˆ212) 3Table 7.7. Unit ball, test meshes and number of LOBPCG iterations.level hmin hmax V F T nit1 0.109665 0.255241 976 2700 875 222 0.046295 0.124278 9736 28314 9317 133 0.023515 0.066545 66256 195804 64827 131 2 3 4 5 6 7 8 9 1000.020.040.060.080.10.120.14eigenvalue numbereigenvalue approximationexactlevel 1level 2level 31 2 310−410−310−2levelapproximation errorτ2 (1st)τ2 (5th)τ2 (9th)4−lFig. 7.8. Unit ball, eigenvalue convergence.135We proceed to compute the ﬁrst ten eigenfunctions. The approximation errorsfor the eigenvalues of (1.7) and S∗h,1Sh,1 are presented in Figure 7.8. The results aresimilar to the previous test problem.Each of the ﬁrst ten computed electric eigenﬁelds, both as a magnitude plot onthe surface and as a vector ﬁeld in the interior, are shown in Appendix B on pages158 to 161.3. Eigenvalues of the Fichera cornerOur third example is the computation of the eigenvalues in the Fichera corner [−1, 1]3\[−1, 0]3. The exact eigenfunctions are not known, but some of them have singularitiesat the origin which makes the problem diﬃcult to approximate. We will compare ourresults with the ones from Table 7.8. These are taken from M. Dauge’s benchmarkwebsite [100], see also the survey [44].Table 7.8. Fichera corner, benchmark results from [100].i ω2i reliable digits conjectured eigenvalue1 3.31381e+00 1 3.2???e+002 5.88635e+00 3 5.88??e+003 5.88635e+00 3 5.88??e+004 1.06945e+01 4 1.0694e+015 1.06945e+01 4 1.0694e+016 1.07006e+01 2 1.07??e+017 1.23345e+01 3 1.232?e+018 1.23345e+01 3 1.232?e+01Two tests were performed for this problem using unstructured tetrahedral anduniform hexahedral meshes. The initial meshes are shown in Figure 7.9.136Fig. 7.9. Fichera corner, initial meshes.The computations were performed on reﬁned grids consisting of 28489 vertices,323072 faces and 159744 tetrahedra and 31841 vertices, 89088 faces and 28672 hexa-hedra, respectively.The results of the eigenvalue approximations for the ﬁrst eight eigenfunctions ofS∗h,1Sh,1, in each case, are reported in Table 7.9.Table 7.9. Fichera corner, results for tetrahedral mesh (column 3) and hexahedralmesh (column 4).i ω2h,i |ω2i − ω2h,i| |ω2i − ω2h,i|1 3.23432e+00 7.94855e-02 2.63062e-022 5.88267e+00 3.67742e-03 1.69117e-023 5.88371e+00 2.64462e-03 1.69511e-024 1.06789e+01 1.55709e-02 6.22111e-025 1.06832e+01 1.12777e-02 6.22377e-026 1.06945e+01 6.08114e-03 1.03244e-017 1.23653e+01 3.07189e-02 1.20678e-018 1.23723e+01 3.77137e-02 1.22141e-01We note that the hexahedral mesh oﬀers better approximation with signiﬁcantly137less memory usage. This can be explained by the fact that the mesh is uniform andthat the dimensions of Xh and Yh,k are balanced better in this case.4. Eigenvalues of a linear accelerator cellOur ﬁnal problem involves complicated geometry modeled with ﬁne hexahedral mesh.It is a linear accelerator induction cell taken from Lawrence Livermore National Lab-oratory’s EMSolve project, see [101]. The mesh has 46382 vertices, 128992 faces and41344 elements and comes from a real-world application.Our code successfully computed the ﬁrst ten eigenvalues of this diﬃcult problem.The magnitudes of the ﬁrst ten electric eigenmodes are visualized in Appendix B onpages 161 to 164.D. The time-harmonic problemIn this section, we report the results of computation for the full time-harmonic system.For ease of implementation, we report results in two dimensions. We also, assumethat the ﬁelds are real, i.e. we are approximating the problem (6.6)-(6.7).Speciﬁcally, the weak formulation is⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩(e,∇×v) + ω (µh, v) = 0 ∀v ∈ V1 = H1(Ω) ,(h,∇×w) + ω (εe,w) = 〈j,w〉 ∀w ∈ V2 = H10(Ω) ,(µh,∇ψ) = 0 ∀ψ ∈ H2 := H1(Ω) .Here, we used the scalar and vector curls deﬁned in (2.13).As in the three-dimensional case, we get that when ω is not an eigenvalue, theleast-squares method is well posed (i.e. has unique solution for compatible data). Asan illustration, we consider an application involving a known smooth solution. We138Table 7.10. Numerical results for the time-harmonic test using exact solver.h ||e||0 ratio nit N time0.125 0.57812 9 384 0.030.0625 0.29133 1.9844 9 1536 0.150.03125 0.14599 1.9955 9 6144 0.930.015625 0.07302 1.9992 9 24576 5.490.0078125 0.03645 1.9998 9 98304 46.70.00390625 0.01826 1.9999 9 393216 425.let Ω be the unit square and take ω = µ = ε = 1. This problem has solutione = x (1− x) sin(πy) , h = ∇(cos(πx) cos(πy)). (7.2)The numerical results on a uniform triuangular mesh with two diﬀerent choices forTh,k are given in Table 7.10 and Table 7.11.The error behavior in L2(Ω) clearly illustrates the expected ﬁrst-order conver-gence rate. Note that, in both cases, the number of iterations required to reducethe residual by a factor of 10−6 remains bounded independently of the number ofunknowns.We also remark that using Multigrid, instead of the exact solver, leads to amodest increase of the number of iterations, while signiﬁcantly reducing the overallsolution time.139Table 7.11. Numerical results for the time-harmonic test using Multigrid.h ||e||0 ratio nit N time0.125 0.57812 9 384 0.030.0625 0.29134 1.9844 11 1536 0.110.03125 0.14600 1.9956 12 6144 0.480.015625 0.07302 1.9993 12 24576 2.380.0078125 0.03645 1.9999 13 98304 11.50.00390625 0.01826 2.0000 13 393216 58.9140CHAPTER VIIICONCLUSIONSThis dissertation introduces the dual least-squares technique for approximation of avariety of problems related to the time-harmonic Maxwell’s equations. We presentedtheoretical results concerning the stability of the discrete problems, as well as resultscharacterizing the error of approximation. We also showed that the methods can beeﬃciently implemented.The abstract least-squares framework was introduced and analyzed in ChapterIII. The results concerning the approximation of the solution and the least-squaressolution operator can possibly be applied to other problems. The theory for the mag-netostatic and electrostatic problems from [26] was expanded in Chapter IV to moregeneral stable spaces and to domains with curved boundaries. Improved regularityresults were obtained in Appendix A. The eigenvalue problem was treated in ChapterV by introducing a new reformulation based on a compact skew-Hermitian operator.We gave estimates for the convergence of the eigenvalues and eigenvectors and showedthat spurious modes are avoided. This method is new and appears to be quite diﬀer-ent from the previously available algorithms for this problem. Finally, we extendedthe method to the full time-harmonic system and characterized its solvability andapproximation error in Chapter VI.As with any introduction of new methodology in a well-developed ﬁeld, thereare many interesting open questions related to the topic of the dissertation. Forexample, it will be interesting to investigate the existence of other pairs of stablespaces. Also, the application of similar ideas to the more general equations describingphotonic crystals looks promising. Additional models that might be of interest arethose involving perfectly matched layers and eddy current problems.141REFERENCES[1] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover,New York, 1965.[2] S. Adam, P. Arbenz, and R. Geus, Eigenvalue solvers for electromagnetic ﬁeldsin cavities, Tech. Report TR 275, Swiss Federal Institute of Technology, 1997.[3] R. A. Adams and J. F. Fournier, Sobolev Spaces, vol. 140 of Pure and AppliedMathematics, Academic Press, Boston, 2003.[4] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault, Vector potentials inthree-dimensional nonsmooth domains, Math. Methods Appl. Sci., 21 (1998),pp. 823–864.[5] D. A. Aruliah, Fast Solvers for Time-Harmonic Maxwell’s Equations in3D, Ph.D. dissertation, Department of Mathematics, University of BritishColumbia, Vancouver, Canada, 2001.[6] G. Auchmuty and J. C. Alexander, L2 well-posedness of planar div-curl systems,Arch. Ration. Mech. Anal., 160 (2001), pp. 91–134.[7] A. Aziz and I. M. Babuˇska, Part I, survey lectures on the mathematical foun-dations of the ﬁnite element method, in The Mathematical Foundations of theFinite Element Method with Applications to Partial Diﬀerential Equations,A. Aziz, ed., Academic Press, New York, 1972, pp. 1–362.[8] I. M. Babuˇska and J. E. Osborn, Eigenvalue problems, in Finite Element Meth-ods (Part 1), vol. II of Handbook of Numerical Analysis, North-Holland, Ams-terdam, The Netherlands, 1991, pp. 641–787.142[9] C. Bacuta, J. H. Bramble, and J. E. Pasciak, New interpolation results andapplications to ﬁnite element methods for elliptic boundary value problems, East-West J. Numer. Math., 3 (2001), pp. 179–198.[10] C. A. Balanis, Advanced Engeneering Electromagnetics, John Wiley & Sons,New York, 1989.[11] F. Ben Belgacem, A. Buﬀa, and Y. Maday, The mortar method for the Maxwell’sequations in 3D, C.R. Acad. Sci. Paris, 329 (1999), pp. 903–908.[12] J.-P. Berenger, A perfectly matched layer for the absorption of electromagneticwaves, J. Comp. Phys., 114 (1994), pp. 185–200.[13] J.-P. Berenger, Three-dimensional perfectly matched layer for the absorption ofelectromagnetic waves, J. Comp. Phys., 127 (1994), pp. 363–379.[14] P. B. Bochev and M. D. Gunzburger, Finite element methods of least-squarestype, SIAM Review, 40 (1998), pp. 789–837.[15] D. Boﬃ, F. Brezzi, and L. Gastaldi, On the convergence of eigenvalues for mixedformulations, Ann. Scuola Norm. Sup. Pisa Cl. Sci, XXV (1997), pp. 131–154.[16] D. Boﬃ, F. Brezzi, and L. Gastaldi, Mixed ﬁnite elements for Maxwell’s eigen-problem: The question of spurious modes, in ENUMATH 97, 2nd EuropeanConference on Numerical Mathematics and Advanced Applications, World Sci-entiﬁc, Singapore, 1998, pp. 180–187. Symposium held at Heidelberg, Germany,September 28 to October 3, 1997.[17] D. Boﬃ, P. Fernandes, L. Gastaldi, and I. Perugia, Computational models ofelectromagnetic resonators: Analysis of edge element approximation, SIAM J.Numer. Anal., 36 (1999), pp. 1264–1290.143[18] A. Bossavit, Computational Electromagnetism: Variational Formulations,Complementarity, Edge Elements, Academic Press, San Diego, 1998.[19] J. H. Bramble, A proof of the inf-sup condition for the Stokes equations onLipschitz domains., Math. Model. and Meth. in Appl. Sci., 13 (2003), pp. 361–371.[20] J. H. Bramble, T. V. Kolev, and J. E. Pasciak, The approximation of theMaxwell eigenvalue problem using a least-squares method, (2004). preprint.[21] J. H. Bramble, R. D. Lazarov, and J. E. Pasciak, A least-squares approach basedon a discrete minus one inner product for ﬁrst order systems, Math. Comp., 66(1997), pp. 935–955.[22] J. H. Bramble, R. D. Lazarov, and J. E. Pasciak, Least-squares for second-orderelliptic problems, Comp. Meth. Appl. Mech. Eng., 152 (1998), pp. 195–210.[23] J. H. Bramble, R. D. Lazarov, and J. E. Pasciak, Least-squares methods forlinear elasticity based on a discrete minus one inner product, Comp. Meth.Appl. Mech. Eng., 191 (2001), pp. 727–744.[24] J. H. Bramble and J. E. Osborn, Rate of convergence estimates for nonselfad-joint eigenvalue approximations, Math. Comp., 27 (1973), pp. 525–549.[25] J. H. Bramble and J. E. Pasciak, Least-squares methods for Stokes equationsbased on a discrete minus one inner product, J. Comp. App. Math., 74 (1996),pp. 155–173.[26] J. H. Bramble and J. E. Pasciak, A new approximation technique fordiv-curl systems, Math. Comp., (2003). Posted electronically on Au-144gust 26, 2003. Accessed online at http://www.ams.org/mcom/0000-000-00/S0025-5718-03-01616-8/S0025-5718-03%-01616-8.pdf on May 21, 2004.[27] J. H. Bramble, J. E. Pasciak, and P. S. Vassilevski, Computational scales ofSobolev norms with application to preconditioning, Math. Comp., 69 (200),pp. 463–480.[28] J. H. Bramble and T. Sun, A negative-norm least squares method for Reissner-Mindlin plates, Math. Comp., 67 (1998), pp. 901–916.[29] J. H. Bramble and X. Zhang, The analysis of multigrid methods, in Techniquesof Scientiﬁc Computing (Part 3), vol. VII of Handbook of Numerical Analysis,North-Holland, Amsterdam, The Netherlands, 2000, pp. 173–415.[30] S. C. Brenner, Convergence of the multigrid V-cycle algorithm for second-orderboundary value problems without full elliptic regularity, Math. Comp., 71 (2001),pp. 507–525.[31] Z. Cai, R. D. Lazarov, T. A. Manteuﬀel, and S. F. McCormick, First-ordersystem least squares for second-order partial diﬀerential equations: Part I, SIAMJ. Num. Anal., 31 (1994), pp. 1785–1802.[32] Z. Cai, T. A. Manteuﬀel, and S. F. McCormick, First-order system least squaresfor second-order partial diﬀerential equations: Part II, SIAM J. Num. Anal., 34(1997), pp. 425–454.[33] Z. Cai, T. A. Manteuﬀel, S. F. McCormick, and J. Ruge, First-order system LL*(FOSLL*): Scalar elliptic partial diﬀerential equations, SIAM J. Num. Anal.,39 (2001), pp. 1418–1445.145[34] Z. Cai and B. C. Shin, The discrete ﬁrst-order system least squares: The second-order elliptic boundary value problem, SIAM J. Num. Anal., 40 (2002), pp. 307–318.[35] M. Cessenat, Mathematical Methods in Electromagnetism: Linear Theory andApplications, vol. 41 of Series on Advances in Mathematics for Applied Sciences,World Scientiﬁc, Singapore, 1996.[36] C. L. Chang, Finite element approximation for grad-div type of systems in theplane, SIAM J. Numerical Analysis, 29 (1992), pp. 590–601.[37] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, vol. 40 ofClassics in Applied Mathematics, SIAM, Philadelphia, 2002.[38] P. Cle´ment, Approximation by ﬁnite element functions using local regularization,Rev. Franc¸aise Automat. Informat. Recherche Ope´rationnelle Se´r. Rouge Anal.Nume´r., 9 (1975), pp. 77–84.[39] M. Costabel, A remark on the regularity of solutions of Maxwell’s equations onLipschitz domains, MMAS, 12 (1990), pp. 365–368.[40] M. Costabel, A coercive bilinear form for Maxwell’s equations, Jour. Math.Anal. Appl., 157 (1991), pp. 527–541.[41] M. Costabel and M. Dauge, Maxwell and Lame´ eigenvalues on polyhedra, Math.Methods Appl. Sci., 22 (1999), pp. 243–258.[42] M. Costabel and M. Dauge, Singularities of electromagnetic ﬁelds in polyhedraldomains, Arch. Rational Mech. Anal., 151 (2000), pp. 221–276.[43] M. Costabel and M. Dauge, Weighted regularization of Maxwell equations inpolyhedral domains, Numer. Math., 93 (2002), pp. 239–277.146[44] M. Costabel and M. Dauge, Computation of resonance frequencies for Maxwellequations in nonsmooth domains, in Topics in Computational Wave PropagationDirect and Inverse Problems, vol. 31 of Lecture Notes in Computational Scienceand Engineering, Springer-Verlag, New York, 2003.[45] M. Costabel, M. Dauge, and S. Nicaise, Singularities of Maxwell interface prob-lems, Mode´l. Math. Anal. Nume´r, 33 (1999), pp. 627–649.[46] M. Dauge, Elliptic Boundary Value Problems on Corner Domains, LectureNotes in Mathematics, 1341, Springer-Verlag, New York, 1988.[47] L. A. Da´valos and D. Zanette, Fundamentals of Electromagnetism: VacuumElectrodynamics, Media and Relativity, Springer-Verlag, New York, 1999.[48] L. Demkowicz, P. Monk, C. Schwab, and L. Vardapetyan, Maxwell eigenvaluesand discrete compactness in two dimensions, Tech. Report 99-12, TICAM, 1999.[49] L. Demkowicz and L. Vardapetyan, Modeling of electromagnetic absorp-tion/scattering problems using hp–adaptive ﬁnite elements, Comp. Meth. Appl.Mech. Eng., 152 (1998), pp. 103–124.[50] A.-S. B. Dhia, C. Hazard, and S. Lohrengel, A singular ﬁeld method for thesolution of Maxwell’s equations in polyhedral domains, SIAM J. Appl. Math.,59 (1999), pp. 2028–2044.[51] E. D’yakonov and M. Orekhov, Minimization of the computational labor in de-termining the ﬁrst eigenvalues of diﬀerential operators, Math. Notes, 27 (1980),pp. 382–391.[52] C. Emson, J. Simkin, and C. Trowbridge, Further developments in three-dimensional eddy current analysis, IEEE Trans. on Magnetics, MAG-21 (1985),147pp. 2231–2234.[53] M. Fabian, P. Habala, P. Ha´jek, V. M. Santalucia, J. Pelant, and V. Zizler,Functional Analysis and Inﬁnite Dimensional Geometry, vol. 8 of CMS Booksin Mathematics, Springer-Verlag, New York, 2001.[54] V. Girault and P. Raviart, Finite Element Approximation of the Navier-StokesEquations, vol. 749 of Lecture Notes in Mathematics, Springer-Verlag, NewYork, 1981.[55] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, Boston, 1985.[56] R. Hiptmair, Finite elements in computational electromagnetism, Acta Numer-ica, 11 (2002), pp. 237–339.[57] J. D. Jackson, Classical Electrodynamics, John Wiley & Sons, New York, 1999.[58] B.-n. Jiang, The Least-Squares Finite Element Method: Theory and Applica-tions in Computational Fluid Dynamics and Electromagnetics, Springer-Verlag,New York, 1998.[59] J.-M. Jin, The Finite Element Method in Electromagnetics, John Wiley & Sons,New York, 2002.[60] J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals, PrincetonUniversity Press, Princeton NJ, 1995.[61] A. Kameari, Three-dimensional eddy current calculation using ﬁnite elementmethod with a-v in conductor and ω in vacuum, IEEE Trans. on Magnetics, 24(1988), pp. 118–121.148[62] T. Kato, Perturbation Theory for Linear Operators, vol. 132 of Lecture Notesin Mathematics, Springer-Verlag, Berlin, 1966.[63] A. V. Knyazev, Convergence rate estimates for iterative methods for a meshsymmetric eigenvalue problem., Sov. J. Num. Anal. Math. Modeling,, 2 (1987),pp. 371–396.[64] A. V. Knyazev, Preconditioned eigensolvers—an oxymoron?, Electron. Trans.Numer. Anal., 7 (1998), pp. 104–123.[65] A. V. Knyazev, Toward the optimal preconditioned eigensolver: Locally optimalblock preconditioned conjugate gradient method, SIAM J. Sci. Comp., 23 (2001),pp. 517–541.[66] A. V. Knyazev and M. E. Argentati, Principal angles between subspaces in anA-based scalar product: Algorithms and perturbation estimates, SIAM J. Sci.Comp., 23 (2002), pp. 2008–2040.[67] A. V. Knyazev and M. E. Argentati, Implementation of a preconditioned eigen-solver using Hypre, Numerical Linear Algebra with Applications, (2004). sub-mitted.[68] S. Korotov and M. Krizek, Finite element analysis of variational crimes for aquasilinear elliptic problem in 3D, Numer. Math., 84 (2000), pp. 549–576.[69] F. Kukuchi, On a discrete compactness property for the Ne´de´lec ﬁnite elements,J. Fac. Sci. Univ. Tokyo, Sect. 1A, Math, 36 (1989), pp. 479–490.[70] B. Lee, T. A. Manteuﬀel, S. F. McCormick, and J. Ruge, First-order sys-tem least-squares for the Helmholtz equation, SIAM J. Sci. Comp., 21 (2000),pp. 1927–1949.149[71] P. Leonard and D. Rodger, Finite element scheme for transient 3D eddy cur-rents, IEEE Trans. on Magnetics, 24 (1988), pp. 58–66.[72] J.-L. Lions and E. Magenes, Non-homogeneous Boundary Value Problems andApplications, Springer-Verlag, New York, 1972.[73] J. C. Maxwell, A Treatise on Electricity and Magnetism, Clarendon Press, Ox-ford, UK, 1873.[74] V. G. Maz’ya and S. V. Poborchi, Diﬀerentiable Functions on Bad Domains,World Scientiﬁc, Singapore, 1997.[75] P. Monk, Finite Element Methods for Maxwell’s Equations, Numerical Mathe-matics and Scientiﬁc Computation, Oxford University Press, Oxford, UK, 2003.[76] P. Monk and L. Demkowicz, Discrete compactness and approximation ofMaxwell’s equations in R3, Math. Comp., 70 (2001), pp. 507–523.[77] J.-C. Ne´de´lec, Mixed ﬁnite elements in R3, Numer. Math., 35 (1980), pp. 315–341.[78] J.-C. Ne´de´lec, A new family of mixed ﬁnite elements in R3, Numer. Math., 50(1986), pp. 57–81.[79] J.-C. Ne´de´lec, Acoustic and Electromagnetic Equations: Integral Representa-tions for Harmonic Problems, Springer-Verlag, New York, 2001.[80] J. Necˇas, Les Me´thodes Directes en The´orie des E´quations Elliptiques, Mason,Paris, 1967.150[81] L. Olson, Multilevel Least-Squares Finite Element Methods for Hyperbolic Par-tial Diﬀerential Equations, Ph.D. dissertation, Department of Mathematics,University of Colorado, Boulder, 2003.[82] J. E. Osborn, Spectral approximation for compact operators, Math. Comp., 29(1975), pp. 712–725.[83] J. E. Pasciak and J. Zhao, Overlapping Schwarz Methods in H(curl) on Non-convex Domains, East West J. Num. Anal., 10 (2002), pp. 221–234.[84] K. D. Paulsen and D. R. Lynch, Elimination of vector parisites in ﬁnite ele-ment Maxwell solutions, IEEE Trans. Microwave Theory Tech, MTT-39 (1991),pp. 395–404.[85] I. Perugia, D. Scho¨tzau, and P. Monk, Stabilized interior penalty methods for thetime-harmonic Maxwell equations, Comp. Meth. Appl. Mech. Eng., 191 (2002),pp. 4675–4697.[86] J. Saranen, On an inequality of Friedrichs, Math. Scand., 51 (1982), pp. 310–322.[87] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functionssatisfying boundary conditions, Math. Comp., 54 (1990), pp. 483–493.[88] C. M. Soukoulis, Photonic Band Gap Materials, Kluwer, Dordrecht, TheNetherlands, 1996. editor.[89] G. Strang and A. Berger, The change in solution due to change in domain,in AMS Symposium on Partial Diﬀerential Equations, American MathematicalSociety, Providence, 1971, pp. 199–206.151[90] H. F. Tiersten, A Development of the Equations of Electromagnetism in Ma-terial Continua, vol. 36 of Tracts in Natural Philosophy, Springer-Verlag, NewYork, 1990.[91] H. Triebel, Interpolation Theory, Function Spaces, Diﬀerential Operators,vol. 18 of North-Holland Mathematical Library, North-Holland, Amsterdam,The Netherlands, 1978.[92] F. T. Ulaby, Fundamentals of Applied Electromagnetics, Prentice Hall, UpperSaddle River, NJ, 2000.[93] J. L. Volakis, A. Chatterjee, and L. C. Kempel, Finite Element Method forElectromagnetics : Antennas, Microwave Circuits, and Scattering Applications,IEEE Press, New York, 1998.[94] C. Weber, A local compactness theorem for Maxwell’s equations, Math. Meth.Appl. Sci., 2 (1980), pp. 12–25.[95] C. Wolfe, U. Navsariwala, and S. Gedney, A parallel ﬁnite-element tearing andinterconnecting algorithm for solution of the vector wave equation with PMLabsorbing medium, IEEE Transactions on Antennas and Propagation, 47 (2000),pp. 278–284.[96] J. Wu and B.-n. Jiang, A least-squares ﬁnite element method for electromagneticscattering problems, Tech. Report CR-202313, NASA, 1996.[97] E. Zeidler, Applied Functional Analysis: Applications to Mathematical Physics,vol. 108 of Applied Mathematical Sciences, Springer-Verlag, New York, 1995.[98] E. Zeidler, Applied Functional Analysis: Main Principles and Their Applica-tions, vol. 109 of Applied Mathematical Sciences, Springer-Verlag, New York,1521995.[99] J. Zhao, Analysis of Finite Element Approximation and Iterative Methods forTime-Dependent Maxwell Problems, Ph.D. dissertation, Department of Mathe-matics, Texas A&M University, College Station, 2002.[100] Benchmark Computations for Maxwell Equations for the Approximation ofHighly Singular Solutions, M. Dauge, Universite´ de Rennes 1, Rennes,France. Accessed online at http://perso.univ-rennes1.fr/monique.dauge/benchmax.html on May 21, 2004.[101] EMSolve Project: Unstructured Grid Computational Electromagnetics UsingMixed Finite Element Methods, Center for Applied Scientiﬁc Computing,Lawrence Livermore National Laboratory. Accessed online at http://www.llnl.gov/CASC/emsolve/ on May 21, 2004.[102] Least-Squares Methods for Computational Electromagnetics Project, Depart-ment of Mathematics, Texas A&M University. Accessed online at http://www.math.tamu.edu/research/numerical analysis/least square/ on May 21, 2004.[103] SciDAC Accelerator Modeling Project, National Energy Research ScientiﬁcComputing Center, Lawrence Berkeley National Laboratory. Accessed onlineat http://scidac.nersc.gov/accelerator/ on May 21, 2004.153APPENDIX AREGULARITY RESULT ON A CONVEX DOMAINIn this appendix, we will provide a regularity estimate for the solution of amagnetostatic div-curl system of a special type. Assume that the domain Ω is convex.Let f ∈ H−1(Ω), and ∇·f = 0 in the sense of distributions. As discussed earlier, thesystem ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩∇×x = f in Ω,∇·x = 0 in Ω,x · n = 0 on ∂Ω,(A.1)has a unique “weak” solution x ∈ L2(Ω) satisfying(x,∇×v) + (x,∇h) = 〈f, v〉 ∀(v,h) ∈ Y1 . (A.2)For g ∈ (H (Ω))∗, consider the problem of ﬁnding ψ ∈ H˜1−0 (Ω) such that〈−∆θ, ψ〉 = 〈g,∇θ〉 ∀θ ∈ H10(Ω) ∩H1+(Ω). (A.3)This problem has a unique solution since the operator −∆ deﬁnes an isomorphism ofH10(Ω) ∩H1+(Ω) onto H˜−1+(Ω), see [30]. Additionally, we have‖∇ψ‖H˜−(Ω)≤ ‖ψ‖H˜1−0 (Ω) ≤ C ‖g‖(H(Ω))∗ ,and therefore, by setting Q1g = ∇ψ we get the unique continuous extension of Q1 asan operator from (H (Ω))∗ to H˜−(Ω). Thus, in particular,‖Q1g‖H˜−1(Ω) ≤ C ‖g‖(H1(Ω))∗ .Analogous to the deﬁnition in §V.A, let S1 be the solution operator deﬁned as154S1g = x, where x solves (A.2) with data f = (I− Q1)g.Lemma A.1 For any ∈ [0, 1], S1 : (H (Ω))∗ → H1−(Ω) is a bounded linear oper-ator.Proof For any g ∈ (H (Ω))∗ we have that (I− Q1)g ∈ H˜−1(Ω) and〈(I− Q1)g,∇θ〉 = 0 (A.4)for arbitrary θ ∈ D(Ω). Therefore, the compatability condition (4.11) is satisﬁed,and S1g is well deﬁned. Moreover, the convexity of Ω and the inf-sup condition(Proposition 4.1) imply that‖S1g‖1 ≤ C ‖(I− Q1)g‖0 , and ‖S1g‖0 ≤ C ‖(I− Q1)g‖H˜−1 ,for g ∈ L2(Ω) and g ∈ (H 1(Ω))∗, respectively. Using the boundedness of Q1 we get‖S1g‖1 ≤ C ‖g‖0 , and ‖S1g‖0 ≤ C ‖g‖(H1(Ω))∗ .Thus, by interpolation,‖S1g‖1− ≤ C ‖g‖(H)∗ , ∀g ∈ (H (Ω))∗ .Corollary A.1 Let ∈ (0, 12). There exists C = C() > 0, such that for dataf ∈ H−(Ω), with ∇·f = 0, the solution of (A.2) is in H 1−(Ω) and we have thestability estimate‖x‖1− ≤ C‖f‖− . (A.5)Proof Since < 12, we have H−(Ω) = H˜−(Ω) = (H (Ω))∗, see Theorem 2.4. By(A.3) and the fact that D(Ω) is dense in H10(Ω) ∩ H1+(Ω), it follows that Q1f = 0155when f ∈ H−(Ω) and ∇·f = 0. For such f, S1f coincides with the solution x of (A.2).The corollary follows from Lemma A.1.Remark A.1 When f ∈ (H(Ω))∗ with ∈ [12, 1], the condition Q1f = 0 implies∇·f = 0 by (A.4). The converse is false. Indeed, for example, let = 1 and φ ∈ L2(Ω)be a non-constant harmonic function. Deﬁne f ∈ (H1(Ω))∗ by〈f, v〉 = −(φ,∇·v) ∀v ∈ H1(Ω) .Clearly ∇·f = 0. On the other hand, (A.3) implies that Q1f = f|H−1 = ∇φ = 0.156APPENDIX BVISUALIZATION OF SOME APPROXIMATE SOLUTIONSFig. B.1. Approximation to the magnetic ﬁeld in the iron core of the transformer.Fig. B.2. Cross-section of the approximate solution ﬁeld for the transformer problem.157Fig. B.3. Approximation to the magnetic ﬁeld in the transformer.158Fig. B.4. Unit ball, eigenmode 1.Fig. B.5. Unit ball, eigenmode 2.Fig. B.6. Unit ball, eigenmode 3.159Fig. B.7. Unit ball, eigenmode 4.Fig. B.8. Unit ball, eigenmode 5.Fig. B.9. Unit ball, eigenmode 6.160Fig. B.10. Unit ball, eigenmode 7.Fig. B.11. Unit ball, eigenmode 8.Fig. B.12. Unit ball, eigenmode 9.161Fig. B.13. Unit ball, eigenmode 10.Fig. B.14. Linear accelerator cell, eigenmode 1.162Fig. B.15. Linear accelerator cell, eigenmode 2.Fig. B.16. Linear accelerator cell, eigenmode 3.Fig. B.17. Linear accelerator cell, eigenmode 4.163Fig. B.18. Linear accelerator cell, eigenmode 5.Fig. B.19. Linear accelerator cell, eigenmode 6.Fig. B.20. Linear accelerator cell, eigenmode 7.164Fig. B.21. Linear accelerator cell, eigenmode 8.Fig. B.22. Linear accelerator cell, eigenmode 9.Fig. B.23. Linear accelerator cell, eigenmode 10.165VITATzanio Valentinov Kolev was born in Plovdiv, Bulgaria on December 3, 1975 toSonia and Valentin Kolev. He graduated summa cum laude with a M.S. degree in ap-plied mathematics from Soﬁa University “St. Kliment Ohridski” in September 1998.The topic of his thesis was Two-level Preconditioning of Elasticity Non-conformingFEM Systems. After graduation, he worked as a researcher at the Central Laboratoryfor Parallel Processing at the Bulgarian Academy of Sciences. Following his militaryservice, from January to September 1999, Mr. Kolev began his graduate studies inmathematics at Texas A&M University in January 2000. He spent the summers of2001, 2002 and 2003 as an intern at the Center for Applied Scientiﬁc Computingat Lawrence Livermore National Laboratory. He defended his dissertation on Least-squares Methods for Computational Electromagnetics in May 2004 and received hisPh.D. in August 2004. His academic advisors were Prof. James H. Bramble and Prof.Joseph E. Pasciak.Tzanio Kolev can be contacted by writing to: Department of Mathematics,Texas A&M University, College Station, TX 77843-3368, or at the e-mail address:tkolev@math.tamu.edu .

LEAST-SQUARES METHODS FOR COMPUTATIONAL

Abstract

Extracted data

LEAST-SQUARES METHODS FOR COMPUTATIONAL

Abstract

Extracted data

Related items

Related items