A mixed method for time-transient acoustic wave propagation in metamaterials
AA MIXED METHOD FOR TIME-TRANSIENT ACOUSTIC WAVEPROPAGATION IN METAMATERIALS
JEONGHUN J. LEE
Abstract.
In this paper we develop a finite element method for acoustic wavepropagation in Drude-type metamaterials. The governing equation is writtenas a symmetrizable hyperbolic system with auxiliary variables. The standardmixed finite elements and discontinuous finite elements are used for spatialdiscretization, and the Crank–Nicolson scheme is used for time discretization.The a priori error analysis of fully discrete scheme is carried out in details.Numerical experiments illustrating the theoretical results and metamaterialwave propagation, are included. Introduction
Metamaterials usually mean the materials with artificial micro/nano-scale struc-tures which show unconventional macro-scale material properties which are notobserved in natural materials. The unconventional material properties of meta-materials have many potential applications in wave propagation. For example,cloaking devices, which hide internal objects from external detection using waverefection, can be made by an appropriate design of metamaterial device. Thereforedevising metamaterials and its numerical simulations are research topics of greatinterest nowadays.There are three major classes of metamaterials, which are for acoustic, elec-tromagnetic, and elastodynamic wave propagation. In this paper we only con-sider acoustic wave propagation in metamaterials. In time-harmonic cases some ofthese wave propagation equations coincide under special circumstances but theyare all different in time transient wave propagation. For the theory of electro-magnetic metamaterials and time-domain finite element methods we refer to, e.g.,[9, 11, 12, 19] and the references in [10] for more comprehensive list of previousstudies. There are also previous studies on elastodynamic metamaterials in, e.g.,[13, 14, 17].To the best of our knowledge there are very limited number of previous studieson numerical methods for time transient acoustic wave propagation in metamateri-als. In [3] some acoustic metamaterial models, the acoustic counterpart of doublynegative index materials in electromagnetics [12, 19], are studied. In the paperthe authors proposed a form of symmetrizable hyperbolic system as the governingequations of acoustic wave propagation in metamaterials. In addition, they provedexistence of weak solutions and showed numerical experiments with the finite dif-ference method.
Date : June, 2, 2020.2000
Mathematics Subject Classification.
Primary: 65N15, 65N30.
Key words and phrases. mixed method, wave propagation, metamaterial. a r X i v : . [ m a t h . NA ] J u l JEONGHUN J. LEE
In this paper we develop a finite element method for the system proposed in[3] and prove the a priori error analysis. For spatial discretization we use themixed finite element for the Poisson equation and some discontinuous finite elementspaces. To circumvent lower convergence rate of the pressure variable in some mixedfinite element pairs, we propose a novel local post-processing which gives numericalpressure with higher order approximation properties (See Subsection 3.3).The paper is organized as follows. In Section 2 we first introduce symbols andnotation in the paper, and then present the governing equations for the acousticwave propagation in metamaterials as well as the energy estimate. In Section 3we introduce finite element discretization for the system and prove the a priorierror analysis. In particular, we show that a local post-processing for the pressurevariable can be used to obtain a numerical pressure which has better approximationproperty than the original numerical pressure. In Section 4 we present the resultsof numerical experiments which illustrate our theoretical results and exotic wavepropagation in metamaterials. 2.
Preliminaries
Notation.
Let Ω ⊂ R d , d = 2 or 3, be a polygonal/polyhedral domain withLipschitz boundary. Throughout this paper we assume that T h is a triangulationof Ω without hanging nodes.We use L r (Ω) to denote the Lebesgue space with the norm (cid:107) v (cid:107) L r = (cid:40)(cid:0)(cid:82) Ω | v | r dx (cid:1) /r , if 1 ≤ r < ∞ , esssup x ∈ Ω {| v ( x ) |} , if r = ∞ . For a domain D ⊂ Ω, L ( D ) and L ( D ; R d ) be the sets of R - and R d -valued squareintegrable functions with inner products ( v, v (cid:48) ) D := (cid:82) D vv (cid:48) dx and ( v , v (cid:48) ) D := (cid:82) D v · v (cid:48) dx . We will use ( · , · ) instead of ( · , · ) D if D = Ω. For an integer l ≥ P l ( D ) and P l ( D ; R d ) are the spaces of R - and R d -valued polynomials of degree ≤ l on D .In the paper H s ( D ), s ≥
0, denotes the Sobolev space based on the L -normwith s -differentiability on the domain D . We refer to [7] for a rigorous definition ofthis space. The norm on H s ( D ) is denoted by (cid:107)·(cid:107) s,D and D is omitted if D = Ω. If ρ is a nonnegative function in L ∞ (Ω), then (cid:107) v (cid:107) ρ and (cid:107) v (cid:107) ρ denotes the ρ -weighted L -norms (cid:0)(cid:82) Ω ρ | v | dx (cid:1) / and (cid:0)(cid:82) Ω ρ v · v dx (cid:1) / .For T > X , let C ([0 , T ]; X ) denote the set offunctions f : [0 , T ] → X that are continuous in t ∈ [0 , T ]. For an integer m ≥
1, wedefine C m ([0 , T ]; X ) = { f | ∂ i f /∂t i ∈ C ([0 , T ]; X ) , ≤ i ≤ m } , where ∂ i f /∂t i is the i -th time derivative in the sense of the Fr´echet derivative in X (cf. [20]). For a function f : [0 , T ] → X , the Bochner norm is defined as (cid:107) f (cid:107) L r (0 ,T ; X ) = (cid:16)(cid:82) T (cid:107) f ( s ) (cid:107) rX ds (cid:17) /r , ≤ r < ∞ , esssup t ∈ (0 ,T ) (cid:107) f ( t ) (cid:107) X , r = ∞ . We define W k,r (0 , T ; X ) for a non-negative integer k and 1 ≤ r ≤ ∞ as the closureof C k ([0 , T ]; X ) with the norm (cid:107) f (cid:107) W k,r (0 ,T ; X ) = (cid:80) ki =0 (cid:107) ∂ i f /∂t i (cid:107) L r (0 ,T ; X ) . The semi-norm (cid:107) f (cid:107) ˙ W k,r (0 ,T ; X ) is defined by (cid:107) f (cid:107) ˙ W k,r (0 ,T ; X ) = (cid:107) ∂ k f /∂t k (cid:107) L r (0 ,T ; X ) . IXED METHOD FOR METAMATERIALS 3
Finally, for a normed space X with its norm (cid:107)·(cid:107) X and functions f , f ∈ X , (cid:107) f , f (cid:107) X will be used to denote (cid:107) f (cid:107) X + (cid:107) f (cid:107) X , and (cid:107) f , f , f (cid:107) X is defined simi-larly.2.2. A metamaterial model of acoustic wave propagation.
A system of equa-tions for the acoustic wave propagation with velocity and pressure unknowns is ρ ∂ v ∂t + grad p = f ,κ − ∂p∂t + div v = g with the density ρ and the bulk modulus κ . In conventional material models thecoefficients ρ and κ are fixed uniformly positive functions in Ω. In this paper we areinterested in metamaterial models such that ρ and κ − are frequency-dependent,more precisely, the temporal Fourier transform of the equations with frequency ω satisfy − iω ˆ ρ ( ω )ˆ v ( ω ) + grad ˆ p ( ω ) = ˆ f ( ω ) , − iω ˆ κ − ( ω )ˆ p ( ω ) + div ˆ v ( ω ) = ˆ g ( ω )withˆ ρ ( ω ) = ρ a (cid:32) − Ω ρ ω − ω ρ (cid:33) ˆ κ − ( ω ) = κ − a (cid:18) − Ω κ ω − ω κ + iγω (cid:19) , γ ≥ ρ a , κ a > ρ ≥ κ ≥ ω ρ > ω κ > γ ≥ v , ˆ p , ˆ f ,ˆ g are the temporal Fourier transforms of v , p , f , g , respectively.To obtain a system of time-dependent equations we introduce new variables u , w , q , r satisfying iω ˆ v ( ω ) = ( ω − ω ρ )ˆ u ( ω ) , − iω ˆ w ( ω ) = ˆ u ( ω ) ,iω ˆ p ( ω ) = ( ω + iγω − ω κ )ˆ q ( ω ) , − iω ˆ r ( ω ) = ˆ q ( ω ) . The system of time-dependent equations are ρ a ∂ v ∂t + grad p + ρ a Ω ρ u = f , (2.1a) κ − a ∂p∂t + div v + κ − a Ω κ q = g, (2.1b) ∂ u ∂t − v + ω ρ w = 0 , (2.1c) ∂ w ∂t − u = 0 , (2.1d) ∂q∂t − p + γq + ω κ r = 0 , (2.1e) ∂r∂t − q = 0 . (2.1f)We assume that the material of wave propagation in Ω consists of a conventionalpositive index material (PIM) on a subdomain Ω P ⊂ Ω and a negative index mate-rial (NIM) on Ω \ Ω P . The PIM and NIM materials are mathematically modeled bythe values of Ω ρ and Ω κ , i.e., Ω ρ = Ω κ = 0 on Ω P and Ω ρ , Ω κ > \ Ω P . Notethat the first two equations in (2.1) are decoupled from the other equations on the JEONGHUN J. LEE domain Ω P because Ω ρ = Ω κ = 0 on Ω P . From this observation we may developnumerical methods which solve different sets of equations on the PIM and NIMdomains. However, we will focus on a monolithic numerical method for the systembecause monolithic approaches can cover problems with varying interfaces betweenPIM and NIM in a unified manner. They can be used for shape optimizationproblems for metamaterial device design, which is our future research interest.For boundary conditions of (2.1) let Γ D , Γ N be the subsets of ∂ Ω such thatΓ D ∩ Γ N = ∅ , Γ D ∪ Γ N = ∂ Ω. Then imposed boundary conditions are p ( t ) = p D ( t ) on Γ D , v ( t ) · n = v N ( t ) on Γ N (2.2)with given functions p D on (0 , T ] × Γ D and v N on (0 , T ] × Γ N , where n is the unitoutward normal vector field on Γ N .To write a variational form of the system let us define function spaces W = L (Ω; R d ) , Q = L (Ω) , V = { v ∈ W : div v ∈ L (Ω) } where div v is defined in the sense of distributions. For future reference we define X := V × Q × W × W × Q × Q with the norm induced by the L norms of thefunction spaces. We also define ρ u , ρ w , ρ q , ρ r to denote (nonnegative) weights ρ u = ρ a Ω ρ , ρ w = ρ a ω ρ Ω ρ , ρ q = κ − a Ω κ , ρ r = κ − a ω κ Ω κ in the rest of this paper.For well-posedness of (2.1) we recall the following result from [3, Theorem 3.1]. Theorem 2.1.
For (2.1) suppose that initial data ( v (0) , p (0) , u (0) , w (0) , q (0) , r (0)) ∈X satisfy v (0) , u (0) , w (0) ∈ H (Ω; R d ) , p (0) , q (0) , r (0) ∈ H (Ω) . In addition, sup-pose that f ∈ C ([0 , T ]; W ) , g ∈ C ([0 , T ]; Q ) hold. Then there exists a uniquesolution ( v , p, u , w , q, r ) ∈ C ([0 , T ]; X ) ∩ C ([0 , T ]; X )(2.3) for the given initial data and f , g . For finite element discretization we need to consider a variational form of (2.1).For simplicity of presentation we assume the homogeneous boundary condition p ( t ) = 0 on ∂ Ω(2.4)for all t ∈ (0 , T ].For simplicity of presentation we will use ˙ v instead of ∂v/∂t in the rest of paper. Definition 2.2.
For f ∈ L ((0 , T ); W ) , g ∈ L ((0 , T ); Q ) , we say ( v , p, u , w , q, r ) ∈ H ([0 , T ] , X ) ∩ L ((0 , T ); X ) a weak solution of (2.1) if it satisfies ( ρ a ˙ v , v (cid:48) ) − ( p, div v (cid:48) ) + ( ρ u u , v (cid:48) ) = ( f , v (cid:48) ) , (2.5a) (cid:0) κ − a ˙ p, p (cid:48) (cid:1) + (div v , p (cid:48) ) + ( ρ q q, p (cid:48) ) = ( g, p (cid:48) ) , (2.5b) ( ˙ u , u (cid:48) ) − ( v , u (cid:48) ) + (cid:0) ω ρ w , u (cid:48) (cid:1) = 0 , (2.5c) ( ˙ w , w (cid:48) ) − ( u , w (cid:48) ) = 0 , (2.5d) ( ˙ q, q (cid:48) ) − ( p, q (cid:48) ) + ( γq, q (cid:48) ) + (cid:0) ω κ r, q (cid:48) (cid:1) = 0 , (2.5e) ( ˙ r, r (cid:48) ) − ( q, r (cid:48) ) = 0(2.5f) for ( v (cid:48) , p (cid:48) , u (cid:48) , w (cid:48) , q (cid:48) , r (cid:48) ) ∈ X and for almost every t ∈ (0 , T ) . IXED METHOD FOR METAMATERIALS 5
One can easily check by the integration by parts that the solution in Theorem 2.1with the boundary condition (2.4) is a weak solution satisfying (2.5).We remark that the above variational form can cover general boundary con-ditions with some necessary modifications. For the boundary condition (2.2) wereplace V by V N = { v ∈ W : div v ∈ L (Ω) , v · n = v N on Γ N } and replace (2.5a) by( ρ a ˙ v , v (cid:48) ) − ( p, div v (cid:48) ) + ( ρ u u , v (cid:48) ) = ( f , v (cid:48) ) − (cid:90) Γ D p D v (cid:48) · n ds for the test function v (cid:48) in V N := { v ∈ W : div v ∈ L (Ω) , v · n = 0 on Γ N } . Theorem 2.3. If ( v , p, u , w , q, r ) is a solution of (2.1) satisfying (2.3) , then (2.6) (cid:107) v , p, u , w , q, r (cid:107) L ∞ ((0 ,T ); X ) ≤ C (cid:107) v (0) , p (0) , u (0) , w (0) , q (0) , r (0) (cid:107) X + C (cid:107) f , g (cid:107) L ((0 ,T ); W × Q ) holds with C which may depend on T . Moreover, if we define E ( t ) = (cid:107) u ( t ) (cid:107) ρ u + (cid:107) v ( t ) (cid:107) ρ a + (cid:107) w ( t ) (cid:107) ρ w + (cid:107) p ( t ) (cid:107) κ − a + (cid:107) q ( t ) (cid:107) ρ q + (cid:107) r ( t ) (cid:107) ρ r (2.7) with the weighted (semi)-norm (cid:107)·(cid:107) ρ , ρ = ρ u , ρ a , ρ w , κ − a , ρ q , ρ r , then E ( t ) ≤ E (0) + C (cid:90) t ( (cid:107) f ( s ) (cid:107) + (cid:107) g ( s ) (cid:107) ) ds. (2.8) with C > depending only on (cid:107) ρ − a (cid:107) L ∞ and (cid:107) κ a (cid:107) L ∞ .Proof. Recall that ( v , p, u , w , q, r ) satisfies (2.5). If we choose ( v (cid:48) , p (cid:48) , u (cid:48) , w (cid:48) , q (cid:48) , r (cid:48) ) =( v , p, u , w , q, r ) in (2.5) and add all the equations, then we get12 ddt (cid:16) (cid:107) u (cid:107) + (cid:107) v (cid:107) ρ a + (cid:107) w (cid:107) + (cid:107) p (cid:107) κ − a + (cid:107) q (cid:107) + (cid:107) r (cid:107) (cid:17) + ( γq, q )(2.9) + (( ρ u − u , v ) + (( ρ q − q, p )+ (cid:0) ( ω ρ − w , u (cid:1) + (cid:0) ( ω κ − r, q (cid:1) = ( f , v ) + ( g, p ) . Let E ( t ) = (cid:107) u ( t ) (cid:107) + (cid:107) v ( t ) (cid:107) ρ a + (cid:107) w ( t ) (cid:107) + (cid:107) p ( t ) (cid:107) κ − a + (cid:107) q ( t ) (cid:107) + (cid:107) r ( t ) (cid:107) . The Cauchy–Schwarz inequality with the above identity gives ddt E ( t ) ≤ CE ( t ) + ( (cid:107) f ( t ) (cid:107) + (cid:107) g ( t ) (cid:107) ) E ( t )with C depending on ρ u , ρ q , ω ρ , ω κ . By Gronwall lemma one can obtain E ( t ) ≤ E (0) + C ( t ) (cid:90) t ( (cid:107) f ( s ) (cid:107) + (cid:107) g ( s ) (cid:107) ) ds. Then (2.6) follows from the equivalence of (cid:112) E ( t ) and (cid:107) ( u ( t ) , v ( t ) , w ( t ) , p ( t ) , q ( t ) , r ( t )) (cid:107) X . To prove (2.8) we choose ( v (cid:48) , p (cid:48) , u (cid:48) , w (cid:48) , q (cid:48) , r (cid:48) ) = ( v , p, ρ u u , ρ w w , ρ q q, ρ r r ) in (2.5)and add all the equations. Then we get12 ddt E ( t ) + ( γρ q q, q ) = ( f , v ) + ( g, p ) . JEONGHUN J. LEE If E ( t ) ≤ E (0), then there is nothing to prove, so we assume E ( t ) > E (0)and will prove (2.8) in the rest of the proof.First, we prove (2.8) assuming that E ( t ) = esssup s ∈ [0 ,t ] E ( s ). Since ( γρ q q, q ) ≥
0, integration of the above identity from 0 to t gives E ( t ) − E (0) ≤ {(cid:107) ρ − a (cid:107) L ∞ , (cid:107) κ a (cid:107) L ∞ } (cid:90) t ( (cid:107) f ( s ) (cid:107) + (cid:107) g ( s ) (cid:107) ) E ( s ) ds ≤ {(cid:107) ρ − a (cid:107) L ∞ , (cid:107) κ a (cid:107) L ∞ } (cid:90) t ( (cid:107) f ( s ) (cid:107) + (cid:107) g ( s ) (cid:107) ) dsE ( t ) . Then E ( t ) ≤ E (0) E ( t ) + 2 max {(cid:107) ρ − a (cid:107) L ∞ , (cid:107) κ a (cid:107) L ∞ } (cid:90) t ( (cid:107) f ( s ) (cid:107) + (cid:107) g (cid:107) ) ds ≤ E (0) + 2 max {(cid:107) ρ − a (cid:107) L ∞ , (cid:107) κ a (cid:107) L ∞ } (cid:90) t ( (cid:107) f ( s ) (cid:107) + (cid:107) g (cid:107) ) ds (2.10)which proves (2.8).If 0 < E ( t ) < esssup s ∈ [0 ,t ] E ( s ), then there exists 0 ≤ t < t such that E ( t ) = esssup s ∈ [0 ,t ] E ( s ) and E ( t ) < E ( t ). By the same argument as above,we can obtain the inequality (2.10) for E ( t ). Then E ( t ) < E ( t ) ≤ E (0) + 2 max {(cid:107) ρ − a (cid:107) L ∞ , (cid:107) κ a (cid:107) L ∞ } (cid:90) t ( (cid:107) f ( s ) (cid:107) + (cid:107) g ( s ) (cid:107) ) ds ≤ E (0) + 2 max {(cid:107) ρ − a (cid:107) L ∞ , (cid:107) κ a (cid:107) L ∞ } (cid:90) t ( (cid:107) f ( s ) (cid:107) + (cid:107) g ( s ) (cid:107) ) ds, so (2.8) is proved. (cid:3) By observing (2.1a) and (2.1b), the auxiliary variables ( u , w , q, r ) interact with( v , p ) only on the NIM domain on which Ω ρ and Ω κ are strictly positive. In fact,the physical meaning of ( u , w , q, r ) on Ω P is not clear, so there is no natural wayto determine the initial data of ( u , w , q, r ) on Ω P . In the following theorem weshow that ( v , p ) in (2.5) is independent on the initial data of ( u , w , q, r ) on Ω P . Asa consequence, any choice of initial data ( u , w , q, r ) on Ω P is allowed to obtain aunique ( v , p ). This argument can be extended to our numerical scheme, so there isno concern in the choice of numerical initial data of ( u , w , q, r ) on Ω P . Theorem 2.4.
Given f ∈ L ((0 , T ); W ) , g ∈ L ((0 , T ); Q ) , and initial data U (0) ∈ X , (2.1) has a unique weak solution. In addition, suppose that U i ∈ H ([0 , T ]; X ) ∩ L ((0 , T ); X ) , i = 1 , are the weak solutions for the two sets ofinitial data U i (0) ∈ X , i = 1 , . For U i ( t ) := ( v i ( t ) , p i ( t ) , u i ( t ) , w i ( t ) , q i ( t ) , r i ( t )) with i = 1 , , if v (0) = v (0) , p (0) = p (0) on Ω , (2.11) u (0) = u (0) , w (0) = w (0) , q (0) = q (0) , r (0) = r (0) on Ω \ Ω P , (2.12) then the same identities hold for the weak solutions U ( t ) , U ( t ) for t ∈ (0 , T ] .Proof. If U i , i = 1 , f , g , and initial data U (0) ∈ X .Then U − U is a weak solution for f = , g = 0, and zero initial data. Then U = U follows by (2.6). IXED METHOD FOR METAMATERIALS 7
To prove the second part of the assertion, let U i , i = 1 , v i (0) , p i (0) , u i (0) , w i (0) , q i (0) , r i (0)) , i = 1 , v , p, u , w , q, r ) be the difference U − U anddefine E ( t ) as in (2.7). Then E ( t ) satisfies (2.8) with f = and g = 0. Moreover, E (0) = 0 because of (2.11), (2.12), so E ( t ) = 0 holds for all t ∈ (0 , T ] and theassertion follows. (cid:3) Finite elements discretization and error analysis
Finite elements for spatial discretization.
Recall that T h is a triangu-lation of Ω without hanging nodes. In the rest of this paper we assume that thedensity functions ρ σ with σ = u, w, q, r are in W , ∞ h ( T h ) where W , ∞ h ( T h ) := { ρ ∈ L (Ω) : ρ | K ∈ L ∞ ( K ) , grad( ρ | K ) ∈ L ∞ ( K ; R d ) ∀ K ∈ T h } with the norm (cid:107) ρ (cid:107) W , ∞ h := sup K ∈T h ( (cid:107) ρ | K (cid:107) L ∞ ( K ) + (cid:107) grad( ρ | K ) (cid:107) L ∞ ( K ) ).Finite element discretization of the first order differential equation form of acous-tic wave equations, is studied in [8]. We extend the approach in [8] to include theauxiliary variables. For discretization of (2.5) with finite elements we use finiteelement spaces V h ⊂ V , W h ⊂ W , Q h ⊂ Q which are defined below. First,BDM l ( K ) for l ≥ l ( K ) for l ≥ l ( K ) = { v ∈ P l ( K ; R d ) } , RTN l ( K ) = { v ∈ P l ( K ; R d ) + x P l ( K ) } where x = ( x , . . . , x d ) T .In the rest of this paper k ≥ k ≥ S k ( K )as either BDM k +1 ( K ) or RTN k ( K ), and define V h as the finite element space V h = { v ∈ V : v | K ∈ S k ( K ) , K ∈ T h } . (3.1)By this definition V h is the Brezzi–Douglas–Marini or the N´ed´elec element of thesecond kind if S k ( K ) = BDM k +1 ( K ) [6, 16], and is the Raviart–Thomas or theN´ed´elec element of the first kind if S k ( K ) = RTN k ( K ) [15, 18]. For the details onthe definition of V h with S k ( K ) = BDM k +1 ( K ) or S k ( K ) = RTN k ( K ), we referto [5, 4] and the original articles [6, 15, 16, 18]. For W h , let W h ( K ) be W h ( K ) = (cid:40) P k +1 ( K ; R d ) if S k ( K ) = BDM k +1 ( K ) , P k ( K ; R d ) if S k ( K ) = RTN k ( K ) , and define W h , Q h as W h = { v ∈ L (Ω; R d ) : v | K ∈ W h ( K ) } , (3.2) Q h = { q ∈ L (Ω) : q | K ∈ P k ( K ) } . (3.3)Let X h be V h × Q h × W h × W h × Q h × Q h . For f ∈ C ([0 , T ]; W ), g ∈ C ([0 , T ]; Q ) the semidiscrete problem is to find ( v h , p h , u h , w h , q h , r h ) ∈ C ((0 , T ] , X h ) JEONGHUN J. LEE which satisfies ( ρ a ˙ v h , v (cid:48) ) − ( p h , div v (cid:48) ) + ( ρ u u h , v (cid:48) ) = ( f , v (cid:48) ) , (3.4a) (cid:0) κ − a ˙ p h , p (cid:48) (cid:1) + (div v h , p (cid:48) ) + ( ρ q q h , p (cid:48) ) = ( g, p (cid:48) ) , (3.4b) ( ˙ u h , u (cid:48) ) − ( v h , u (cid:48) ) + (cid:0) ω ρ w h , u (cid:48) (cid:1) = 0 , (3.4c) ( ˙ w h , w (cid:48) ) − ( u h , w (cid:48) ) = 0 , (3.4d) ( ˙ q h , q (cid:48) ) − ( p h , q (cid:48) ) + ( γq h , q (cid:48) ) + (cid:0) ω κ r h , q (cid:48) (cid:1) = 0 , (3.4e) ( ˙ r h , r (cid:48) ) − ( q h , r (cid:48) ) = 0(3.4f)for ( v (cid:48) , p (cid:48) , u (cid:48) , w (cid:48) , q (cid:48) , r (cid:48) ) ∈ X h and for all t ∈ (0 , T ].We will not discuss an error analysis for semidiscrete solutions in the paper.Instead, we will show a detailed error analysis of fully discrete solutions in thesubsection below.3.2. Error analysis of fully discrete solutions.
In this subsection we considerfully discrete solutions of (3.4) with the Crank–Nicolson scheme and show the apriori error estimates.For
T > t = T /N for a natural number N and define { t n } Nn =0 by t n = n ∆ t . For a variable σ : [0 , T ] → X for a Hilbert space X , we will use σ nh and σ n for the numerical solution of σ at t n and σ ( t n ), respectively. The variable σ can be u , v , w , p, q, r in the problem. As such, f n and g n will denote f ( t n ) and g ( t n ) for a time-dependent functions f ∈ C ([0 , T ]; W ) and g ∈ C ([0 , T ]; Q ). Forsimplicity we will also use the definitions¯ ∂ t v n + := 1∆ t (cid:0) v n +1 − v n (cid:1) , v n + := 12 (cid:0) v n + v n +1 (cid:1) for any sequence { v n } Nn =0 with the upper index n .The Crank–Nicolson scheme of (2.5) is the following: For given U n := ( v nh , p nh , u nh , w nh , q nh , r nh ) , f n , f n +1 , g n , g n +1 , we find U n +1 h := ( v n +1 h , p n +1 h , u n +1 h , w n +1 h , q n +1 h , r n +1 h ) ∈ X h such that (cid:16) ρ a ¯ ∂ t v n + h , v (cid:48) (cid:17) − (cid:16) p n + h , div v (cid:48) (cid:17) + (cid:16) ρ u u n + h , v (cid:48) (cid:17) = (cid:16) f n + , v (cid:48) (cid:17) , (3.5a) (cid:16) κ − a ¯ ∂ t p n + h , p (cid:48) (cid:17) + (cid:16) div v n + h , p (cid:48) (cid:17) + (cid:16) ρ q q n + h , p (cid:48) (cid:17) = (cid:16) g n + , p (cid:48) (cid:17) , (3.5b) (cid:16) ¯ ∂ t u n + h , u (cid:48) (cid:17) − (cid:16) v n + h , u (cid:48) (cid:17) + (cid:16) ω ρ w n + h , u (cid:48) (cid:17) = 0 , (3.5c) (cid:16) ¯ ∂ t w n + h , w (cid:48) (cid:17) − (cid:16) u n + h , w (cid:48) (cid:17) = 0 , (3.5d) (cid:16) ¯ ∂ t q n + h , q (cid:48) (cid:17) − (cid:16) p n + h , q (cid:48) (cid:17) + (cid:16) γq n + h , q (cid:48) (cid:17) + (cid:16) ω κ r n + h , q (cid:48) (cid:17) = 0 , (3.5e) (cid:16) ¯ ∂ t r n + h , r (cid:48) (cid:17) − (cid:16) q n + h , r (cid:48) (cid:17) = 0(3.5f)for all ( v (cid:48) , p (cid:48) , u (cid:48) , w (cid:48) , q (cid:48) , r (cid:48) ) ∈ X h .For the well-definedness of this fully discrete scheme, we show that U n +1 h = if U nh = , f n = f n +1 = , g n = g n +1 = 0 . (3.6) IXED METHOD FOR METAMATERIALS 9
To show it, assume that (3.6) is true. Then (3.5) becomes1∆ t (cid:0) ρ a v n +1 h , v (cid:48) (cid:1) − (cid:0) p n +1 h , div v (cid:48) (cid:1) + 12 (cid:0) ρ u u n +1 h , v (cid:48) (cid:1) = 0 , (3.7a) 1∆ t (cid:0) κ − a p n +1 h , p (cid:48) (cid:1) + 12 (cid:0) div v n +1 h , p (cid:48) (cid:1) + 12 (cid:0) ρ q q n +1 h , p (cid:48) (cid:1) = 0 , (3.7b) 1∆ t (cid:0) u n +1 h , u (cid:48) (cid:1) − (cid:0) v n +1 h , u (cid:48) (cid:1) + 12 (cid:0) ω ρ w n +1 h , u (cid:48) (cid:1) = 0 , (3.7c) 1∆ t (cid:0) w n +1 h , w (cid:48) (cid:1) − (cid:0) u n +1 h , w (cid:48) (cid:1) = 0 , (3.7d) 1∆ t (cid:0) q n +1 h , q (cid:48) (cid:1) − (cid:0) p n +1 h , q (cid:48) (cid:1) + 12 (cid:0) γq n +1 h , q (cid:48) (cid:1) + 12 (cid:0) ω κ r n +1 h , q (cid:48) (cid:1) = 0 , (3.7e) 1∆ t (cid:0) r n +1 h , r (cid:48) (cid:1) − (cid:0) q n +1 h , r (cid:48) (cid:1) = 0 . (3.7f)Let P h and P h be the standard L projections into Q h and W h . If we take v (cid:48) = v n +1 h , p (cid:48) = p n +1 h , u (cid:48) = P h ( ρ u u n +1 h ) , w (cid:48) = ω ρ P h ( ρ u w n +1 h ) , q (cid:48) = P h ( ρ q q n +1 h ) , r (cid:48) = ω κ P h ( ρ q r n +1 h ) , in (3.7a) and add all the equations, then we get1∆ t (cid:16) (cid:107) v n +1 h (cid:107) ρ a + (cid:107) p n +1 h (cid:107) κ − a + (cid:107) u n +1 h (cid:107) ρ u + (cid:107) w n +1 h (cid:107) ρ w + (cid:107) q n +1 h (cid:107) ρ q + (cid:107) r n +1 h (cid:107) ρ r (cid:17) + (cid:0) γω κ ρ q r n +1 h , r n +1 h (cid:1) = 0 , so v n +1 h = , p n +1 h = 0. From these, u n +1 h = w n +1 h = follows by taking u (cid:48) = u n +1 h and w (cid:48) = ω ρ w n +1 h in (3.7c) and (3.7d), and then by adding them. Finally, q n +1 h = r n +1 h = 0 follows by taking q (cid:48) = q n +1 h and r (cid:48) = ω κ r n +1 h in (3.7e) and (3.7f),and then by adding them. Therefore, U n +1 h = .For the error analysis we use e nσ = σ n − σ nh for the error of variable σ ( σ = v , u , w , p, q, r ) at t = t n . For error equations we consider the difference of theaverage of (2.5) at t n and t n +1 , and the fully discrete scheme (3.5). Then the errorequations are (cid:16) ρ a ( ˙ v n + − ¯ ∂ t v n + h ) , v (cid:48) (cid:17) − (cid:16) e n + p , div v (cid:48) (cid:17) + (cid:16) ρ u e n + u , v (cid:48) (cid:17) = 0 , (cid:16) κ − a ( ˙ p n + − ¯ ∂ t p n + h ) , p (cid:48) (cid:17) + (cid:16) div e n + v , p (cid:48) (cid:17) + (cid:16) ρ q e n + q , p (cid:48) (cid:17) = 0 , (cid:16) ˙ u n + − ¯ ∂ t u n + h , u (cid:48) (cid:17) − (cid:16) e n + v , u (cid:48) (cid:17) + (cid:16) ω ρ e n + w , u (cid:48) (cid:17) = 0 , (cid:16) ˙ w n + − ¯ ∂ t w n + h , w (cid:48) (cid:17) − (cid:16) e n + u , w (cid:48) (cid:17) = 0 , (cid:16) ˙ q n + − ¯ ∂ t q n + h , q (cid:48) (cid:17) − (cid:16) e n + p , q (cid:48) (cid:17) + (cid:16) γe n + q , q (cid:48) (cid:17) + (cid:16) ω κ e n + r , q (cid:48) (cid:17) = 0 , (cid:16) ˙ r n + − ¯ ∂ t r n + h , r (cid:48) (cid:17) − (cid:16) e n + q , r (cid:48) (cid:17) = 0 . For v (cid:48) ∈ H s (Ω; R d ) , s > , we define Π h as the canonical interpolation operatorsof RTN or BDM element which satisfydiv Π h v (cid:48) = P h div v (cid:48) , (cid:107) v (cid:48) − Π h v (cid:48) (cid:107) ≤ Ch m (cid:107) v (cid:48) (cid:107) m (3.9)with m := max { s, k + 1 + δ } where δ = 1 if V h is a BDM element and δ = 0 if V h is an RTN element. Using Π h , P h , P h , we can define the decomposition of errors e n v = e I,n v + e h,n v := ( v n − Π h v n ) + (Π h v n − v nh ) , (3.10) e n u = e I,n u + e h,n u := ( u n − P h u n ) + ( P h u n − u nh ) , (3.11) e n w = e I,n w + e h,n w := ( w n − P h w n ) + ( P h w n − w nh ) , (3.12) e np = e I,np + e h,np := ( p n − P h p n ) + ( P h p n − p nh ) , (3.13) e nq = e I,nq + e h,nq := ( q n − P h q n ) + ( P h q n − q nh ) , (3.14) e nr = e I,nr + e h,nr := ( r n − P h r n ) + ( P h r n − r nh ) . (3.15)For estimates of the interpolation errors denoted by e I,nσ for a variable σ , let us usea generic symbol I h σ n to denote the interpolation of the exact solution σ n into thecorresponding finite element space. More specifically, I h = Π h if σ = v , I h = P h if σ = u , w , and I h = P h if σ = p, q, r . Then it holds that (cid:107) e I,nσ (cid:107) = (cid:107) σ n − I h σ n (cid:107) ≤ Ch s (cid:107) σ n (cid:107) s , (3.16)with 12 < s ≤ k + 1 + δ if σ = v , (3.17) 0 ≤ s ≤ k + 1 + δ if σ = u , w , (3.18) 0 ≤ s ≤ k + 1 if σ = p, q, r. (3.19)By (3.9) we can obtain (cid:0) e h,np , div v (cid:48) (cid:1) = 0 ∀ v (cid:48) ∈ V h , (cid:0) div e h,n v , q (cid:48) (cid:1) = 0 ∀ q (cid:48) ∈ Q h . By these identities, the orthogonality of L projections, and some algebraic manip-ulations, the previous error equations are reduced to (cid:16) ρ a ¯ ∂ t e h,n + v , v (cid:48) (cid:17) − (cid:16) e h,n + p , div v (cid:48) (cid:17) + (cid:16) ρ u e h,n + u , v (cid:48) (cid:17) = F nv ( v (cid:48) ) , (3.20a) (cid:16) κ − a ¯ ∂ t e h,n + p , p (cid:48) (cid:17) + (cid:16) div e h,n + v , p (cid:48) (cid:17) + (cid:16) ρ q e h,n + q , p (cid:48) (cid:17) = F np ( p (cid:48) ) , (3.20b) (cid:16) ¯ ∂ t e h,n + u , u (cid:48) (cid:17) − (cid:16) e h,n + v , u (cid:48) (cid:17) + (cid:16) ω ρ e h,n + w , u (cid:48) (cid:17) = F nu ( u (cid:48) ) , (3.20c) (cid:16) ¯ ∂ t e h,n + w , w (cid:48) (cid:17) − (cid:16) e h,n + u , w (cid:48) (cid:17) = F nw ( w (cid:48) ) , (3.20d) (cid:16) ¯ ∂ t e h,n + q , q (cid:48) (cid:17) − (cid:16) e h,n + p , q (cid:48) (cid:17) + (cid:16) γe h,n + q , q (cid:48) (cid:17) + (cid:16) ω κ e h,n + r , q (cid:48) (cid:17) = F nq ( q (cid:48) ) , (3.20e) (cid:16) ¯ ∂ t e h,n + r , r (cid:48) (cid:17) − (cid:16) e h,n + q , r (cid:48) (cid:17) = F nr ( r (cid:48) )(3.20f) IXED METHOD FOR METAMATERIALS 11 where F nv ( v (cid:48) ) = − (cid:16) ρ a (cid:16) Π h ¯ ∂ t v n + − ˙ v n + (cid:17) , v (cid:48) (cid:17) − (cid:16) ρ u e I,n + u , v (cid:48) (cid:17) , (3.21a) F np ( p (cid:48) ) = − (cid:16) κ − a (cid:16) ¯ ∂ t P h p n + − ˙ p n + (cid:17) , p (cid:48) (cid:17) − (cid:16) ρ q e I,n + q , p (cid:48) (cid:17) , (3.21b) F nu ( u (cid:48) ) = − (cid:16)(cid:16) ¯ ∂ t P h u n + − ˙ u n + (cid:17) , u (cid:48) (cid:17) , (3.21c) F nw ( w (cid:48) ) = − (cid:16)(cid:16) ¯ ∂ t P h w n + − ˙ w n + (cid:17) , w (cid:48) (cid:17) , (3.21d) F nq ( q (cid:48) ) = − (cid:16)(cid:16) ¯ ∂ t P h q n + − ˙ q n + (cid:17) , q (cid:48) (cid:17) , (3.21e) F nr ( r (cid:48) ) = − (cid:16)(cid:16) ¯ ∂ t P h r n + − ˙ r n + (cid:17) , r (cid:48) (cid:17) . (3.21f)In the discussions below we will use E n defined as( E n ) = (cid:107) e h,n u (cid:107) ρ u + (cid:107) e h,n v (cid:107) ρ a + (cid:107) e h,n w (cid:107) ρ w + (cid:107) e h,np (cid:107) κ − a + (cid:107) e h,nq (cid:107) ρ q + (cid:107) e h,nr (cid:107) ρ r . (3.22) Proposition 3.1.
For given f ∈ C ([0 , T ]; W ) , g ∈ C ([0 , T ]; Q ) and initial data ( v (0) , u (0) , w (0) , p (0) , q (0) , r (0)) ∈ X suppose that ( v , u , w , p, q, r ) is a weak solu-tion of (2.5) . Assume that numerical initial data ( v h (0) , p h (0) , u h (0) , w h (0) , q h (0) , r h (0)) ∈ X h satisfy (cid:107) u (0) − u h (0) (cid:107) ρ u + (cid:107) v (0) − v h (0) (cid:107) ρ a + (cid:107) w (0) − w h (0) (cid:107) ρ w + (cid:107) P h p (0) − p h (0) (cid:107) ρ p + (cid:107) P h q (0) − q h (0) (cid:107) ρ q + (cid:107) P h r (0) − r h (0) (cid:107) ρ r (3.23) ≤ C (cid:48) h s , < s ≤ k + 1 + δ with C (cid:48) independent of h . We also assume that the exact solution ( v , u , w , p, q, r ) and ρ a , κ − a , ρ u , ρ w , ρ q , ρ r satisfy the regularity assumptions (3.25) , (3.26) , (3.27) below.If { ( v nh , u nh , w nh , p nh , q nh , r nh ) } Nn =1 is a numerical solution obtained by the fully dis-crete scheme (3.5) , then E n ≤ C h s + C (∆ t ) + C h s , < s ≤ k + 1 + δ (3.24) with constants C , C C such that C depends on C (cid:48) in (3.23) and (cid:107) ρ u , ρ a , ρ w (cid:107) L ∞ < ∞ , (3.25) C depends on (cid:107) v , u , w , p, q, r (cid:107) ˙ W , (0 ,T ; L ) , (cid:107) ρ a , ρ u , ρ w , ρ q , ρ r (cid:107) L ∞ , (3.26) and C depends on (cid:107) ˙ v (cid:107) L (0 ,T ; H s ) , (cid:107) ˙ p (cid:107) L (0 ,T ; H s ) , (cid:107) ρ u , κ − a (cid:107) W , ∞ h , (cid:107) κ a , ρ a , ρ w , ρ q , ρ r (cid:107) L ∞ (3.27) with s = max { , s − } .Proof. Let us take ( v (cid:48) , p (cid:48) , u (cid:48) , w (cid:48) , q (cid:48) , r (cid:48) ) ∈ X h as( e h,n + v , e h,n + p , P h ( ρ u e h,n + u ) , P h ( ρ w e h,n + w ) , P h ( ρ q e h,n + q ) , P h ( ρ r e h,n + r )) in (3.20) and add all the equations, then we can get12 (cid:0) ( E n +1 ) − ( E n ) (cid:1) + ∆ t (cid:16) γρ q e h,n + q , e h,n + q (cid:17) (3.28) = ∆ t (cid:16) F nv ( e h,n + v ) + F np ( e h,n + p ) + F nu ( P h ( ρ u e h,n + u )) (cid:17) + ∆ t (cid:16) F nw ( P h ( ρ w e h,n + w )) + F nq ( P h ( ρ q e h,n + q )) + F nr ( P h ( ρ r e h,n + r ))) (cid:17) =: R n . The proof of (3.24) has three steps. In the first step, we shall prove (cid:107) R n (cid:107) ≤ ( C ,n (∆ t ) + C ,n h s )( E n + E n +1 ) , < s ≤ k + 1 + δ (3.29)with C ,n depending on (cid:107) v , u , w , p, q, r (cid:107) ˙ W , ( t n ,t n +1 ; L ) , (cid:107) ρ a , ρ u , ρ w , ρ q , ρ r (cid:107) L ∞ , (3.30)and with C ,n depending on (cid:107) ˙ v (cid:107) L ( t n ,t n +1 ; H s ) , (cid:107) ˙ p (cid:107) L ( t n ,t n +1 ; H s ) , (cid:107) ρ u , κ − a (cid:107) W , ∞ h , (cid:107) κ a , ρ a , ρ w , ρ q , ρ r (cid:107) L ∞ (3.31)with s = max { , s − } . The more detailed dependence of C ,n , C ,n will beclarified in the proof of (3.29). In the second step, we prove E n ≤ E + 2 n − (cid:88) i =0 ( C ,i (∆ t ) + C ,i h s ) , < s ≤ k + 1 + δ. (3.32)In the third step, we prove E ≤ C h s , < s ≤ k + 1 + δ (3.33)with C depending on C (cid:48) , the shape regularity of T h , and (cid:107) ρ u , ρ a , ρ w (cid:107) L ∞ < ∞ .Note that the conclusion (3.24) follows from (3.28), (3.29), (3.32), (3.33) bytaking C = (cid:80) ≤ i ≤ n C ,i , C = (cid:80) ≤ i ≤ n C ,i . Therefore we will devote the rest ofproof to prove (3.29), (3.32), and (3.33).Since the proof of (3.29) is long, we show (3.32) and (3.33) first assuming that(3.29) is proved. For (3.32), note that E n +1 − E n ≤ C ,n (∆ t ) + C ,n h s )is obtained by (3.28) and (3.29). Then (3.32) follows by induction. For (3.33), thetriangle inequality and (cid:107) ρ u , ρ a , ρ w (cid:107) L ∞ < ∞ give E ≤ (cid:107) e I, u (cid:107) ρ u + (cid:107) e I, v (cid:107) ρ a + (cid:107) e I, w (cid:107) ρ w + C (cid:48) h s (3.34) ≤ Ch s (cid:107) u (0) , v (0) , w (0) (cid:107) s + C (cid:48) h s , < s ≤ k + 1 + δ with C > T h and (cid:107) ρ u , ρ a , ρ w (cid:107) L ∞ .Before we prove (3.29), let us review some interpolation error estimates fromtime discretization schemes. Since the interpolation operator I h (= Π h , P h , P h ) isindependent in time, the time derivative of I h σ is same as I h ˙ σ as long as they IXED METHOD FOR METAMATERIALS 13 are well-defined pointwisely in time. Then, assuming that a general variable σ ∈ L ( t n , t n +1 ; L ) is sufficiently regular, we can obtain∆ t (cid:107) ˙ σ n + − ¯ ∂ t σ n + (cid:107) = (cid:107) ∆ t ˙ σ n + − ( σ n +1 − σ n ) (cid:107) (3.35) ≤ C ∆ t (cid:107) σ (cid:107) ˙ W , ( t n ,t n +1 ; L ) , ∆ t (cid:107) ¯ ∂ t σ n + − ¯ ∂ t I h σ n + (cid:107) = (cid:107) σ n +1 − I h σ n +1 − ( σ n − I h σ n ) (cid:107) (3.36) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:90) t n +1 t n ( ˙ σ ( s ) − I h ˙ σ ( s )) ds (cid:13)(cid:13)(cid:13)(cid:13) ≤ Ch s (cid:107) ˙ σ (cid:107) L ( t n ,t n +1 ; H s ) with s satisfying the range conditions in (3.17), (3.18), (3.19).For the proof of (3.29), it suffices to estimate the terms in (3.21) by the definitionof R n in (3.28).By (3.16), (3.35), (3.36), and the triangle inequality, and by assuming that theexact solution v is sufficiently regular, one can show(3.37) ∆ t | F nv ( e h,n + v ) |≤ C (cid:16) (∆ t ) (cid:107) v (cid:107) ˙ W , ( t n ,t n +1 ; L ) + h s (cid:107) ˙ v (cid:107) L ( t n ,t n +1 ; H s ) (cid:17) (cid:107) e h,n + v (cid:107) ρ a with C > (cid:107) ρ a , ρ u (cid:107) L ∞ , (cid:107) ρ a (cid:107) − L ∞ . Noting the identity F nu ( P h ( ρ u e h,n + u )) = − (cid:16) ¯ ∂ t P h u n + − ˙ u n + , P h ( ρ u e h,n + u ) (cid:17) = − (cid:16) ¯ ∂ t u n + − ˙ u n + , P h ( ρ u e h,n + u ) (cid:17) and the inequality (cid:107) P h ( ρ u e h,n + u ) (cid:107) ≤ (cid:107) ρ u e h,n + u (cid:107) ≤ (cid:107)√ ρ u (cid:107) L ∞ (cid:107) e h,n + u (cid:107) ρ u , we ob-tain ∆ t | F nu ( P h ( ρ u e h,n + u )) |≤ C (∆ t ) (cid:107) u (cid:107) ˙ W , ( t n ,t n +1 ; L ) (cid:107) e h,n + u (cid:107) ρ u (3.38)with C > (cid:107) ρ u (cid:107) L ∞ by (3.35). A completely similar argument gives∆ t | F nw ( P h ( ρ w e h,n + w )) | ≤ C (∆ t ) (cid:107) w (cid:107) ˙ W , ( t n ,t n +1 ; L ) (cid:107) e h,n + w (cid:107) ρ w , (3.39) ∆ t | F nq ( P h ( ρ q e h,n + q )) | ≤ C (∆ t ) (cid:107) q (cid:107) ˙ W , ( t n ,t n +1 ; L ) (cid:107) e h,n + q (cid:107) ρ q , (3.40) ∆ t | F nr ( P h ( ρ r e h,n + r )) | ≤ C (∆ t ) (cid:107) r (cid:107) ˙ W , ( t n ,t n +1 ; L ) (cid:107) e h,n + r (cid:107) ρ r (3.41)with C > (cid:107) ρ w (cid:107) L ∞ , (cid:107) ρ q (cid:107) L ∞ , (cid:107) ρ r (cid:107) L ∞ , respectively.Now we only need to estimate the F np ( e h,n + p )-involved term but it needs anadditional discussion because the standard approximation theory with Q h givesonly a bound of O ( h k +1 ). To obtain an estimate of (cid:107) e h,np (cid:107) κ − a with a bound of O ( h s ), < s ≤ k + 1 + δ , we will use (cid:12)(cid:12)(cid:12)(cid:16) κ − a ¯ ∂ t e I,n + p , p (cid:48) (cid:17)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:16) ( κ − a − P κ − a ) ¯ ∂ t e I,n + p , p (cid:48) (cid:17)(cid:12)(cid:12)(cid:12) (3.42) ≤ Ch (cid:107) κ − a (cid:107) W , ∞ h (cid:107)√ κ a (cid:107) L ∞ (cid:107) ¯ ∂ t e I,n + p (cid:107) (cid:107) p (cid:48) (cid:107) κ − a , (cid:12)(cid:12)(cid:12)(cid:16) ρ u e I,n + p , p (cid:48) (cid:17)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:16) ( ρ u − P ρ u ) e I,n + p , p (cid:48) (cid:17)(cid:12)(cid:12)(cid:12) (3.43) ≤ Ch (cid:107) ρ u (cid:107) W , ∞ h (cid:107)√ κ a (cid:107) L ∞ (cid:107) e I,n + p (cid:107) (cid:107) p (cid:48) (cid:107) κ − a . By (3.16), (3.42), (3.43), (3.35), (3.36), assuming that the exact solutions are suf-ficiently regular, one can obtain(3.44) ∆ t | F np ( e h,n + p ) |≤ C (cid:16) (∆ t ) (cid:107) p (cid:107) ˙ W , ( t n ,t n +1 ; L ) + h s (cid:107) ˙ p (cid:107) L ( t n ,t n +1 ; H s ) (cid:17) (cid:107) e h,n + p (cid:107) κ − a for 0 ≤ s ≤ k + 1 + δ with C > (cid:107) κ a (cid:107) L ∞ , (cid:107) ρ u (cid:107) W , ∞ h , (cid:107) κ − a (cid:107) W , ∞ h .Combining (3.37), (3.38), (3.39), (3.40), (3.41), (3.44), and the triangle inequalitywith the definition of R n , we can obtain (3.29) with C ,n , C ,n with the dependencedescribed in (3.30), (3.31). (cid:3) Theorem 3.2.
Suppose that the assumptions of Proposition 3.1 hold and δ isdefined in the same way. Then (cid:107) u n − u nh (cid:107) ρ u + (cid:107) v n − v nh (cid:107) ρ a + (cid:107) w n − w nh (cid:107) ρ w ≤ E n + Ch s (cid:107) u , v , w (cid:107) C ([0 ,T ]; H s ) (3.45) with < s ≤ k + 1 + δ where C depends on the shape regularity of T h and the degree k . Similarly, (cid:107) p n − p nh (cid:107) κ − a + (cid:107) q n − q nh (cid:107) ρ q + (cid:107) r n − r nh (cid:107) ρ r ≤ E n + Ch s (cid:107) p, q, r (cid:107) C ([0 ,T ]; H s ) (3.46) with ≤ s ≤ k + 1 .Proof. By the triangle inequality and the definition of E n ,(3.47) (cid:107) u n − u nh (cid:107) ρ u + (cid:107) v n − v nh (cid:107) ρ a + (cid:107) w n − w nh (cid:107) ρ w ≤ E n + (cid:107) e I,n u (cid:107) ρ u + (cid:107) e I,n v (cid:107) ρ a + (cid:107) e I,n w (cid:107) ρ w . By the approximation properties of Π h and P h ,(3.48) (cid:107) e I,n u (cid:107) ρ u + (cid:107) e I,n v (cid:107) ρ a + (cid:107) e I,n w (cid:107) ρ w ≤ Ch s (cid:107) u , v , w (cid:107) C ([0 ,T ]; H s ) , < s ≤ k + 1 + δ holds with C > T h and (cid:107) ρ u , ρ a , ρ w (cid:107) L ∞ .Then (3.45) follows.Similarly, the triangle inequality gives (cid:107) p n − p nh (cid:107) κ − a + (cid:107) q n − q nh (cid:107) ρ q + (cid:107) r n − r nh (cid:107) ρ r ≤ E n + (cid:107) e I,np (cid:107) κ − a + (cid:107) e I,nq (cid:107) ρ q + (cid:107) e I,nr (cid:107) ρ r . Then (cid:107) e I,np (cid:107) κ − a + (cid:107) e I,nq (cid:107) ρ q + (cid:107) e I,nr (cid:107) ρ r ≤ Ch s (cid:107) p, q, r (cid:107) C ([0 ,T ]; H s ) , ≤ s ≤ k + 1holds with C > T h and (cid:107) κ − a , ρ q , ρ r (cid:107) L ∞ , so(3.46) follows. (cid:3) Error analysis for post-processed solutions. If V h is a BDM element,then δ = 1 and the optimal convergence rate in (3.46) is one order lower thanthe one in (3.45). This lower convergence rate can be circumvented by a localpost-processing which we introduce below.Throughout this subsection we assume that the exact solutions are sufficientlyregular and we will not concern about low regularity of exact solutions. In our localpost-processing, we first find { p n, ∗ h } N − n =0 , a new numerical solution approximating IXED METHOD FOR METAMATERIALS 15 p ( t n + ∆ t/ p ( t n ). The goal is to show that (cid:107) p n, ∗ h − p ( t n + ∆ t/ (cid:107) can have O ( h k +2 ) convergence rate.To define the local post-processing let us define Q ∗ h = { q ∈ L (Ω) : q | K ∈ P k +2 ( K ) , K ∈ T h } . We define p n, ∗ h ∈ Q ∗ h as (cid:90) K p n, ∗ h dx = (cid:90) K p n + h dx, (3.49) (cid:0) grad p n, ∗ h , grad p (cid:48) (cid:1) K = − (cid:16) ρ a ¯ ∂ t v n + , grad p (cid:48) (cid:17) K − (cid:16) ρ u u n + h , grad p (cid:48) (cid:17) K (3.50) + (cid:16) f n + , grad p (cid:48) (cid:17) K , ∀ p (cid:48) ∈ Q ∗ h for all K ∈ T h . Note that this post-processing is solving a system with a blockdiagonal matrix such that the size of each matrix block is the number of DOFsof Q ∗ h on one simplex K . Therefore, the computational costs are negligibly smallcompared to the computational costs of the original linear system. Lemma 3.3.
Suppose that the assumptions of Proposition 3.1 hold and V h is aBDM element. If { p n, ∗ h } N − n =0 is defined as in (3.49) , (3.50) , and the exact solution ( v , p, u , w , q, r ) is sufficiently regular, then (cid:107) p ( t n + ∆ t/ − p ∗ ,n + h (cid:107) ≤ C ((∆ t ) + h k +2 ) with a constant C > which depends on the constants C , C , C in Proposi-tion 3.1, (cid:107) u (cid:107) C ([0 ,T ]; H k +1 ) , (cid:107) p (cid:107) C ([0 ,T ]; H k +2 ) , and (cid:107) p (cid:107) ˙ W , ∞ ( t n ,t n +1 ; L ) .Proof. Recall that p n + = ( p n + p n +1 ) = ( p ( t n ) + p ( t n +1 )). We first note that (cid:107) p ( t n + ∆ t/ − p n + (cid:107) ≤ C (∆ t ) (cid:107) p (cid:107) ˙ W , ∞ ( t n ,t n +1 ; L ) by an argument with the Taylor expansion. By the triangle inequality it suffices toestimate (cid:107) p n + − p ∗ ,n + h (cid:107) .To show an error estimate of (cid:107) p n + − p n, ∗ h (cid:107) consider the error equation (cid:16) grad( p n + − p n, ∗ h ) , grad p (cid:48) (cid:17) = − (cid:16) ρ a ( ˙ v n + − ¯ ∂ t v n + h ) , grad p (cid:48) (cid:17) − (cid:16) ρ u ( u n + − u n + h ) , grad p (cid:48) (cid:17) ∀ p (cid:48) ∈ Q ∗ h from the definition of p n, ∗ h and (2.1a).For P ∗ h , the L projection to Q ∗ h , we can rewrite this equation as (cid:16) grad( P ∗ h p n + − p n, ∗ h ) , grad p (cid:48) (cid:17) = − (cid:16) grad( p n + − P ∗ h p n + ) , grad p (cid:48) (cid:17) − (cid:16) ρ a ( ˙ v n + − ¯ ∂ t v n + h ) , grad p (cid:48) (cid:17) − (cid:16) ρ u ( u n + − u n + h ) , grad p (cid:48) (cid:17) . If we take p (cid:48) = P ∗ h p n + − p n, ∗ h and use the Cauchy–Schwarz inequality, then we get (cid:107) grad( P ∗ h p n + − p n, ∗ h ) (cid:107) ≤ (cid:107) grad( p n + − P ∗ h p n + ) (cid:107) + C (cid:16) (cid:107) ˙ v n + − ¯ ∂ t v n + h (cid:107) + (cid:107) u n + − u n + h (cid:107) ρ u (cid:17) (3.51) =: I + I + I with C > ρ a and Ω ρ .To estimate I , assuming that p is sufficiently regular, we use the Bramble–Hilbert lemma and get (cid:107) grad( p n + − P ∗ h p n + ) (cid:107) ,K ≤ Ch k +1 K (cid:107) p n , p n +1 (cid:107) k +2 ,K , ∀ K ∈ T h . (3.52)An estimate of I is obtained by Theorem 3.2 as (cid:107) u n + − u n + h (cid:107) ρ u ≤ C ((∆ t ) + h k +2 )(3.53)under the assumption that u is sufficiently regular.We estimate I by estimating I a := (cid:107) ˙ v n + − ¯ ∂ t v n + (cid:107) , I b := (cid:107) ¯ ∂ t e I,n + v (cid:107) , I c := (cid:107) ¯ ∂ t e h,n + v (cid:107) . (3.54)By (3.35) and (3.36), I a ≤ C (∆ t ) (cid:107) v (cid:107) ˙ W , ∞ ( t n ,t n +1 ; L ) , (3.55) I b ≤ Ch k +1 (cid:107) ˙ v (cid:107) L ∞ ( t n ,t n +1 ; H k +1 ) . (3.56)The estimate of I c is more technical. First, note that it is enough to estimate (cid:107) ¯ ∂ t e h,n + v (cid:107) ρ a since ρ a > ∂ t e h,n + v = v + v with v , v ∈ V h such thatdiv v = 0 , ( ρ a v , v ) = 0 , (cid:107) div v (cid:107) ≤ C (cid:107) v (cid:107) with C > h . Using this decomposition, we can rewrite (3.20a) as( ρ a ( v + v ) , v (cid:48) ) − (cid:16) e h,n + p , div v (cid:48) (cid:17) + (cid:16) ρ u e h,n + u , v (cid:48) (cid:17) = F nv ( v (cid:48) ) . If v (cid:48) = v , then (cid:107) v (cid:107) ρ a ≤ ( C (cid:107) e h,n + p (cid:107) + (cid:107) e h,n + u (cid:107) ρ u ) (cid:107) v (cid:107) + | F nv ( v ) | . Recall that (cid:107) e h,n + p (cid:107) , (cid:107) e h,n + u (cid:107) ρ u are estimated in Proposition 3.1 and | F nv ( v ) | isestimated by (3.37). As a consequence, (cid:107) v (cid:107) ≤ C ( h k +1 + (∆ t ) )holds with C > (cid:107) v (cid:107) ˙ W , ∞ ( t n ,t n +1 ; L ) , (cid:107) ˙ v (cid:107) L ∞ ( t n ,t n +1 ; H k +1 ) , and theconstants C , C , C in Proposition 3.1. If v (cid:48) = v , then we get (cid:107) v (cid:107) ρ a ≤ (cid:107) e h,n + u (cid:107) (cid:107) v (cid:107) + | F nv ( v ) | . An estimate of (cid:107) v (cid:107) can be obtained by a completely similar argument for theestimate of (cid:107) v (cid:107) . Therefore, by combining the estimates of (cid:107) v (cid:107) and (cid:107) v (cid:107) , wehave I ≤ C ((∆ t ) + h k +1 )(3.57)with C > (cid:107) v (cid:107) ˙ W , ∞ ( t n ,t n +1 ; L ) , (cid:107) ˙ v (cid:107) L ∞ ( t n ,t n +1 ; H k +1 ) , and the con-stants C , C , C in Proposition 3.1. IXED METHOD FOR METAMATERIALS 17
By combining (3.51), (3.52), (3.53), (3.55), (3.56), (3.57), we obtained (cid:107) grad( P ∗ h p n + − p ∗ ,n + h ) (cid:107) ≤ C ((∆ t ) + h k +1 ) . An element-wise Poincar´e inequality and the above estimate give (cid:107) P ∗ h p n + − p n, ∗ h (cid:107) ≤ Ch (cid:107) grad( P ∗ h p n + − p n, ∗ h ) (cid:107) ≤ Ch ((∆ t ) + h k +1 ) . From this we can obtain (cid:107) p n + − p n, ∗ h (cid:107) ≤ (cid:107) p n + − P ∗ h p n, ∗ (cid:107) + (cid:107) P ∗ h p n + − p n, ∗ h (cid:107) ≤ Ch k +2 (cid:107) p n + (cid:107) k +2 + Ch ((∆ t ) + h k +1 ) , which is the desired estimate. (cid:3) Numerical results
In this section we present the results of numerical experiments to illustrate thevalidity of our theoretical analysis. In the first set of experiments we use a manufactured solution and show theconvergence rates of errors with various finite element discretizations. Specifically,let Ω = [0 , × [0 ,
1] with the subdomain Ω = [3 / , / × [0 , ρ a = κ a = ω ρ = ω κ = 1, γ = 0 on Ω whereas Ω ρ , Ω κ are defined asΩ ρ = Ω κ = (cid:40) \ Ω A manufactured solution is constructed with w ( x, y ) = (cid:18) (1 + sin t )( x y + xy )cos(2 t )( x + y + cos x ) (cid:19) , r ( x, y ) = cos(3 t ) xy, (4.1)and the other functions u ( x, y ), v ( x, y ), q ( x, y ), p ( x, y ), f ( x, y ), g ( x, y ) are definedby (2.1). Table 1.
Errors and convergence rates with V h ( K ) × Q h ( K ) × W h ( K ) = BDM ( K ) × P ( K ) × P ( K ; R d ) for the exact solutionin (4.1). h (cid:107) v − v h (cid:107) (cid:107) u − u h (cid:107) (cid:107) w − w h (cid:107) (cid:107) p − p h (cid:107) error rate error rate error rate error rate8 5.53e-03 – 9.51e-03 – 3.65e-03 – 1.65e-01 –16 1.34e-03 2.04 2.39e-03 1.99 9.15e-04 2.00 8.26e-02 1.0032 3.49e-04 1.94 5.98e-04 2.00 2.29e-04 2.00 4.13e-02 1.0064 8.42e-05 2.05 1.49e-04 2.00 5.72e-05 2.00 2.06e-02 1.00 h (cid:107) p − p ∗ h (cid:107) (cid:107) q − q h (cid:107) (cid:107) r − r h (cid:107) error rate error rate error rate8 5.53e-03 – 9.51e-03 – 3.65e-03 –16 1.34e-03 2.04 2.39e-03 1.99 9.15e-04 2.0032 3.49e-04 1.94 5.98e-04 2.00 2.29e-04 2.0064 8.42e-05 2.05 1.49e-04 2.00 5.72e-05 2.00 Table 2.
Errors and convergence rates with V h ( K ) × Q h ( K ) × W h ( K ) = RTN ( K ) × P ( K ) × P ( K ; R d ) for the exact solutionin (4.1). h (cid:107) v − v h (cid:107) (cid:107) u − u h (cid:107) (cid:107) w − w h (cid:107) (cid:107) p − p h (cid:107) error rate error rate error rate error rate8 1.78e-01 – 7.50e-02 – 7.61e-02 – 1.65e-01 –16 8.86e-02 1.01 3.74e-02 1.00 3.81e-02 1.00 8.26e-02 1.0032 4.43e-02 1.00 1.87e-02 1.00 1.90e-02 1.00 4.13e-02 1.0064 2.21e-02 1.00 9.34e-03 1.00 9.52e-03 1.00 2.06e-02 1.00 h (cid:107) p − p ∗ h (cid:107) (cid:107) q − q h (cid:107) (cid:107) r − r h (cid:107) error rate error rate error rate8 1.78e-01 – 7.50e-02 – 7.61e-02 –16 8.86e-02 1.01 3.74e-02 1.00 3.81e-02 1.0032 4.43e-02 1.00 1.87e-02 1.00 1.90e-02 1.0064 2.21e-02 1.00 9.34e-03 1.00 9.52e-03 1.00 Table 3.
Errors and convergence rates with V h ( K ) × Q h ( K ) × W h ( K ) = BDM ( K ) × P ( K ) × P ( K ; R d ) for the exact solutionin (4.1). h (cid:107) v − v h (cid:107) (cid:107) u − u h (cid:107) (cid:107) w − w h (cid:107) (cid:107) p − p h (cid:107) error rate error rate error rate error rate8 1.18e-04 – 1.43e-04 – 5.91e-05 – 4.03e-03 –16 1.08e-05 3.44 9.90e-06 3.85 6.22e-06 3.25 1.01e-03 2.0032 1.38e-06 2.97 8.13e-07 3.61 7.37e-07 3.08 2.52e-04 2.0064 1.74e-07 2.98 8.33e-08 3.29 9.08e-08 3.02 6.30e-05 2.00 h (cid:107) p − p ∗ h (cid:107) (cid:107) q − q h (cid:107) (cid:107) r − r h (cid:107) error rate error rate error rate8 1.18e-04 – 1.43e-04 – 5.91e-05 –16 1.08e-05 3.44 9.90e-06 3.85 6.22e-06 3.2532 1.38e-06 2.97 8.13e-07 3.61 7.37e-07 3.0864 1.74e-07 2.98 8.33e-08 3.29 9.08e-08 3.02 We consider 4 different spatial discretizations such that the local finite elementspaces V h ( K ) × Q h ( K ) × W h ( K ) areBDM ( K ) × P ( K ) × P ( K ; R d ) , RTN ( K ) × P ( K ) × P ( K ; R d ) , (4.2) BDM ( K ) × P ( K ) × P ( K ; R d ) , RTN ( K ) × P ( K ) × P ( K ; R d ) . (4.3)For triangulation we use the structured meshes obtained by the bisection ofuniform N × N squares of Ω for N = 8 , , ,
64. Then the maximum mesh sizeis a multiple of h = 1 /N with a uniform constant independent of N . We usethe Crank–Nicolson scheme for time discretization with ∆ t = h for (4.2) and with∆ t = h for (4.3) in order to see optimal convergence rates of spatial discretization The datasets generated during and/or analyzed during the current study are available from thecorresponding author on reasonable request.
IXED METHOD FOR METAMATERIALS 19
Table 4.
Errors and convergence rates with V h ( K ) × Q h ( K ) × W h ( K ) = RTN ( K ) × P ( K ) × P ( K ; R d ) for the exact solutionin (4.1). h (cid:107) v − v h (cid:107) (cid:107) u − u h (cid:107) (cid:107) w − w h (cid:107) (cid:107) p − p h (cid:107) error rate error rate error rate error rate8 2.97e-03 – 2.23e-03 – 2.75e-03 – 4.10e-03 –16 6.93e-04 2.10 5.58e-04 2.00 6.88e-04 2.00 1.01e-03 2.0232 2.01e-04 1.79 1.40e-04 2.00 1.72e-04 2.00 2.55e-04 1.9964 4.56e-05 2.14 3.49e-05 2.00 4.30e-05 2.00 6.40e-05 1.99 h (cid:107) p − p ∗ h (cid:107) (cid:107) q − q h (cid:107) (cid:107) r − r h (cid:107) error rate error rate error rate8 2.97e-03 – 2.23e-03 – 2.75e-03 –16 6.93e-04 2.10 5.58e-04 2.00 6.88e-04 2.0032 2.01e-04 1.79 1.40e-04 2.00 1.72e-04 2.0064 4.56e-05 2.14 3.49e-05 2.00 4.30e-05 2.00 errors. For these four different discretization schemes the errors and convergencerates are presented in Tables 1–4.In these experiments the errors of (cid:107) u − u h (cid:107) , (cid:107) v − v h (cid:107) , (cid:107) w − w h (cid:107) , (cid:107) p − p h (cid:107) , (cid:107) q − q h (cid:107) , (cid:107) r − r h (cid:107) are computed at T = 0 .
25 whereas the error (cid:107) p − p ∗ h (cid:107) is computed at T − ∆ t . Although we used the weighted norms (cid:107)·(cid:107) ρ u , (cid:107)·(cid:107) ρ w , (cid:107)·(cid:107) ρ q , (cid:107)·(cid:107) ρ r in our error analysis for the errors of u , w , q , r , here we compute thestandard L norms for those errors. Since the L norms are the upper bounds ofthe weighted norms, the optimal convergence rates of the L norm errors are theresults stronger than the optimal convergence rates of the weighted norm errors.In all of these experiments we can see the convergence rates which are expected inour error analysis. Figure 1.
Wave propagation with p D in (4.4), µ f = 18, at t = 0 . , . Figure 2.
Wave propagation with p D in (4.4), µ f = 19, at t = 0 . , . Figure 3.
Wave propagation with p D in (4.4), µ f = 20, at t = 0 . , . , × [0 ,
2] with Ω = [3 / , / × [0 , ρ = Ω κ = (cid:40)
80 on Ω \ Ω , ω ρ = ω κ = 40 on Ω , so the medium is a metamaterial on Ω but is a conventional material on Ω \ Ω .We remark that the choices of these parameters are made without considerationof physical ranges of parameter values. We also set p D ( t, x, y ) = (cid:40)
10 sin( µ f π ( x + y − t ) if t > x + y and x < / , µ f is a constant. In the following experiments we impose the boundarycondition (2.2) with p D in the above and Γ D = ∂ Ω. Speaking more intuitively,
IXED METHOD FOR METAMATERIALS 21 this condition gives an incoming wave propagation from the bottom-left corner ofΩ with frequency 5 µ f . Figure 4.
Wave propagation with p D in (4.5) at t = 0 . , . , . , . µ f = 18 , ,
20. The finiteelements with V h ( K ) × Q h ( K ) × W h ( K ) = RTN ( K ) ×P ( K ) ×P ( K ; R d ) are usedfor spatial discretization and T h is the structured mesh obtained by bisecting 50 × t = 0 . µ f = 18, µ f = 19, and µ f = 20, respectively. We call the three regions theleft, the middle, and the right subdomains. In all of the figures in Figures 1–3 thewave propagation patterns look standard plane waves in the lower part of the leftsubdomain whereas they are more complicated due to the waves reflected by theinterface of the left and the middle subdomains. We can clearly see reversed wavepropagation patterns on the metamaterial layer Ω in the three figures at t = 0 . In addition, wave propagation patterns on the right subdomain in the figures at t = 0 .
4, are nearly plane waves with propagation directions similar to the patternson the left subdomain. One can see that the details of wave propagation patterns,particularly the shapes and directions of the reversed patterns on the metamateriallayer, depend on the frequency of waves.In the last experiment we set Ω = [0 , × [0 ,
2] with Ω = [3 / , × [0 , ρ = Ω κ = (cid:40)
80 on Ω \ Ω , ω ρ = ω κ = 80 on Ω . We also set p D ( t, x, y ) = (cid:40)
10 exp( − (1 + sin(20 π ( x + ( y − − t ))) if y − < . , . (4.5)We impose the boundary condition (2.2) with the above p D on the left-side { } × [0 ,
2] and with 0 on the other sides of Ω. The finite elements are V h ( K ) × Q h ( K ) × W h ( K ) = RTN ( K ) × P ( K ) × P ( K ; R d ) and T h is the structured mesh sameas the second set of experiments. ∆ t = 0 .
002 with the Crank–Nicolson scheme.The wave propagation patterns are presented in Figure 4. One can see that wavepropagation patterns are not conventional on the metamaterial layer Ω .5. Conclusion
In this paper we developed finite element methods for acoustic wave propaga-tion in the Drude-type metamaterials. We combined the mixed finite elements forthe Poisson equations and piecewise discontinuous finite element spaces for spatialdiscretization. For time discretization we use the Crank–Nicolson scheme. We car-ried out the a priori error analysis and proposed a local post-processing schemeto overcome low approximation property of the pressure for the BDM type finiteelements. The numerical experiments show the validity of our theoretical analysisas well as atypical wave propagation patterns in metamaterials.
References [1] Douglas N. Arnold, Richard S. Falk, and R. Winther,
Preconditioning in H (div) and appli-cations , Math. Comp. (1997), no. 219, 957–984. MR 1401938 (97i:65177)[2] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther, Multigrid in H (div) and H (curl),Numer. Math. (2000), no. 2, 197–217. MR 1754719 (2001d:65161)[3] C. Bellis and B. Lombard, Simulating transient wave phenomena in acoustic metamaterialsusing auxiliary fields , Wave Motion (2019), 175–194. MR 3904336[4] Daniele Boffi, Franco Brezzi, and Michel Fortin, Mixed finite element methods and applica-tions , Springer Series in Computational Mathematics, vol. 44, Springer, Heidelberg, 2013.MR 3097958[5] F. Brezzi and M. Fortin,
Mixed and hybrid finite element methods , Springer Series in com-putational Mathematics, vol. 15, Springer, 1992. MR MR2233925 (2008i:35211)[6] Franco Brezzi, Jr. Jim Douglas, and L. D. Marini,
Two families of mixed finite elementsfor second order elliptic problems , Numer. Math. (1985), no. 2, 217–235. MR 799685(87g:65133)[7] Lawrence C. Evans, Partial differential equations , Graduate Studies in Mathematics, vol. 19,American Mathematical Society, Providence, RI, 1998. MR 1625845 (99e:35001)[8] Tunc Geveci,
On the application of mixed finite element methods to the wave equations ,RAIRO Mod´el. Math. Anal. Num´er. (1988), no. 2, 243–250. MR 945124 (89i:65116) IXED METHOD FOR METAMATERIALS 23 [9] Yunqing Huang, Jichun Li, and Wei Yang,
Modeling backward wave propagation in metama-terials by the finite element time-domain method , SIAM J. Sci. Comput. (2013), no. 1,B248–B274. MR 3033069[10] Jichun Li, A literature survey of mathematical study of metamaterials , Int. J. Numer. Anal.Model. (2016), no. 2, 230–243. MR 3421776[11] Jichun Li and Yunqing Huang, Time-domain finite element methods for Maxwell’s equa-tions in metamaterials , Springer Series in Computational Mathematics, vol. 43, Springer,Heidelberg, 2013. MR 3013583[12] Jichun Li and Aihua Wood,
Finite element analysis for wave propagation in double negativemetamaterials , J. Sci. Comput. (2007), no. 2, 263–286. MR 2320572[13] Graeme W. Milton and Pierre Seppecher, Realizable response matrices of multi-terminalelectrical, acoustic and elastodynamic networks at a given frequency , Proc. R. Soc. Lond.Ser. A Math. Phys. Eng. Sci. (2008), no. 2092, 967–986. MR 2379501[14] Graeme W. Milton and John R. Willis,
On modifications of Newton’s second law and linearcontinuum elastodynamics , Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. (2007),no. 2079, 855–880. MR 2293080[15] J.-C. N´ed´elec,
Mixed finite elements in R , Numer. Math. (1980), no. 3, 315–341.MR 592160 (81k:65125)[16] , A new family of mixed finite elements in R , Numer. Math. (1986), no. 1, 57–81.MR 864305 (88e:65145)[17] A. N. Norris and A. L. Shuvalov, Elastic cloaking theory , Wave Motion (2011), no. 6,525–538. MR 2811881[18] P.-A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic prob-lems , Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delleRicerche (C.N.R.), Rome, 1975), Springer, Berlin, 1977, pp. 292–315. Lecture Notes in Math.,Vol. 606. MR 0483555 (58
Developing a time-domain finite element methodfor the Lorentz metamaterial model and applications , J. Sci. Comput. (2016), no. 2, 438–463. MR 3519188[20] Kosaku Yosida, Functional analysis , 6th ed., Springer Classics in Mathematics, Springer-Verlag, 1980.
Department of Mathematics, Baylor University, Waco, TX , USA
E-mail address ::