A Finite-Element Model for the Hasegawa-Mima Wave Equation
AA Finite-Element Model for the Hasegawa-Mima WaveEquation
Hagop Karakazian * Sophie Moufawad † Nabil Nassif ‡ January 26, 2021
Abstract
In a recent work [1], two of the authors have formulated the non-linear space-time Hasegawa-Mima plasma equationas a coupled system of two linear PDEs, a solution of which is a pair ( u, w ) , with w = ( I − ∆) u . The first equationis of hyperbolic type and the second of elliptic type. Variational frames for obtaining weak solutions to the initialvalue Hasegawa-Mima problem with periodic boundary conditions were also derived. Using the Fourier basis in thespace variables, existence of solutions were obtained. Implementation of algorithms based on Fourier series leads tosystems of dense matrices.In this paper, we use a finite element space-domain approach to semi-discretize the coupled variational Hasegawa-Mima model, obtaining global existence of solutions in H on any time interval [0 , T ] , ∀ T .In the sequel, full-discretization using an implicit time scheme on the semi-discretized system leads to a nonlinearfull space-time discrete system with a nonrestrictive condition on the time step.Tests on a semi-linear version of the implicit nonlinear full-discrete system are conducted for several initial data,assessing the efficiency of our approach. Acknowledgements : The authors would like to express their gratitude to Prof. Ghassan Antar, AUB Physics Depart-ment, for the several discussions held on the physics of the Hasegawa-Mima wave phenomena, providing us suitabletesting initial conditions and for his positive comments on the overall numerical results.
Keywords : Hasegawa-Mima; Periodic Sobolev Spaces; Petrov-Galerkin Approximations; Finite-Element Method;Semi-Discrete systems; Implicit Finite-Differences
AMS Subject Classification : 35A01; 35M33; 65H10; 65M60; 78M10
Magnetic plasma confinement is one of the most promising ways in future energy production. The Hasegawa-Mima(HM) model is a simplified two-dimensions turbulent system model which describes the time evolution of drift wavescaused during plasma confinement. To understand the phenomena, several mathematical models can be found inliterature[2, 3, 4, 5], of which the simplest and powerful two dimensions turbulent system model is the HM equationthat describes the time evolution of drift waves in a magnetically-confined plasma. It was derived by Akira Hasegawaand Kunioki Mima during late 70s[3, 4]. When normalized, it can[6, 7] be put as the following PDE that is third orderin space and first order in time: − ∆ u t + u t = { u, ∆ u } + { p, u } (1)where { u, v } = u x v y − u y v x is the Poisson bracket, u ( x, y, t ) describes the electrostatic potential, p = ln n ω ci is a function depending on the background particle density n and the ion cyclotron frequency ω ci , which in turndepends on the initial magnetic field. In this context, p = 0 refers to homogeneous plasma, and p (cid:54) = 0 refers to * Hagop Karakazian, Department of Mathematics, American University of Beirut, Lebanon ([email protected]) † Sophie Moufawad, Department of Mathematics, American University of Beirut, Lebanon ([email protected]) ‡ Nabil Nassif, Department of Mathematics, American University of Beirut, Lebanon ([email protected]) a r X i v : . [ m a t h . NA ] J a n on-homogeneous plasma. As a cultural note, equation (1) is also referred as the Charney-Hasegawa-Mima equationin geophysical context that models the time-evolution of Rossby waves in the atmosphere[6].In this paper, we deal with the numerical solution to Hasegawa-Mima equation on a rectangular domain with thesolution u , satisfying periodic boundary conditions (PBCs). For that purpose, we consider Ω = (0 , L ) × (0 , L ) ⊂ R and use the frame of periodic Sobolev spaces which are closed subspace of H m (Ω) , and therefore itself a Hilbertspace. Specifically: H P (Ω) = L (Ω) and H ∞ P (Ω) := ∩ m ≥ { H mP } ,H P (Ω) = { u ∈ H (Ω) | u ( x,
0) = u ( x, L ) , x ∈ (0 , L ) a.e. , u (0 , y ) = u ( L, y ) , y ∈ (0 , L ) a.e. } ,H P (Ω) = (cid:8) u ∈ H (Ω) | u, u x , u y ∈ H P (Ω) (cid:9) H P (Ω) = (cid:8) u ∈ H (Ω) | u, u x , u y , u xx , u yy , u xy , u yx ∈ H P (Ω) (cid:9) (2)In addition, we use for p > , the periodic Banach-Sobolev Spaces: W ,pP (Ω) = { u ∈ W ,p (Ω) | u ( x,
0) = u ( x, L ) a.e. x ∈ (0 , L ) , u (0 , y ) = u ( L, y ) a.e. y ∈ (0 , L ) } , Given an initial data u : Ω → R , we seek u : Ω × [0 , T ] → R such that: − ∆ u t + u t = { u, ∆ u } + p x u y − p y u x on Ω × (0 , T ] (1) PBCs on u, u x , and u y on ∂ Ω × (0 , T ] (2) u ( x, y,
0) = u ( x, y ) on Ω (3) (3)In [1], without loss of generality for the proof of existence, we assume that the background particle density n is afunction of x only, such that p x = ˆ k is a constant and p y = 0 , i.e. n = e Ax + B for A, B ∈ R . When dealing with(3.1), the major difficulty to circumvent, both theoretically and computationally, is the Poisson bracket { u, ∆ u } . Toovercome this issue, we have formulated it in [1], as a coupled system of linear hyperbolic-elliptic PDEs that will benaturally amenable to provide a Finite Element scheme for obtaining a numerical approximation/simulation. For thispurpose, a new variable w = − ∆ u + u is introduced, leading to the identity { u, ∆ u } = { u, u − w } = { u, u } + { u, − w } = −{ u, w } = { w, u } = w x u y − w y u x = − (cid:126)V ( u ) · ∇ w where (cid:126)V ( u ) = − u y (cid:126) i + u x (cid:126) j is a divergence-free vector field (div ( (cid:126)V ( u )) = 0 ). Then system (3) with p x = ˆ k and p y = 0 , becomes equivalent to the coupled hyperbolic-elliptic PDE system, w t + (cid:126)V ( u ) · ∇ w = ˆ ku y on Ω × (0 , T ] (1) − ∆ u + u = w on Ω × (0 , T ] (2) PBC’s on u, u x , u y , w on ∂ Ω × [0 , T ] (3) u (0) = u and w (0) = w on Ω . (4) (4)Using equation (5) which is obtained through Green’s formula and the imposition of periodic boundary conditions, (cid:68) (cid:126)V ( u ) · ∇ v, w (cid:69) L = − (cid:68) (cid:126)V ( u ) · ∇ w, v (cid:69) L (5)system (4) can be put in the following strong semi-variational form (on the space variables) whereby one seeks a pair { u, w } ∈ C ([0 , T ] , H (Ω) ∩ H P (Ω)) × (cid:2) C ([0 , T ] , L (Ω)) ∩ C ((0 , T ) , L (Ω)) (cid:3) such that (cid:104) w t , v (cid:105) L = (cid:68) (cid:126)V ( u ) · ∇ v, w (cid:69) L + (cid:68) ˆ ku y , v (cid:69) L , (1) (cid:104) u, v (cid:105) H = (cid:104) w, v (cid:105) L , (2) ∀ v ∈ W , ∞ P (Ω) ∩ H P (Ω) , ∀ t ∈ (0 , T ] (6)with u (0) = u ∈ H P (Ω) , w (0) = w = u − ∆ u ∈ L (Ω) .A similar formulation to (6) has been handled in [1] where using Fourier series, we prove the existence of a solution { u, w } ∈ C ([0 , T ] , H P (Ω)) ∩ L (0 , T ; H ) × (cid:2) C ([0 , T ] , L (Ω)) ∩ C ((0 , T ) , L (Ω)) (cid:3) . Note that the restriction v ∈ W ,pP (Ω) ensures that (cid:126)V ( u ) · ∇ v ∈ L (Ω) , which justifies this formulation. This allowsreducing the regularity of the initial conditions ( u , w ) from H P (Ω) × H P (Ω) in [1] to H P (Ω) × L (Ω) .The formulation used in this paper intends to avoid dealing with the time-derivative w t . For that purpose, we derive(6) a time integral semi-variational formulation. 2 ime Integral Formulation By integrating (6.1) over the temporal interval [ t, t + τ ] , with ≤ t ≤ T − τ , one reaches the following L IntegralFormulation: u ∈ L (0 , T ; H (Ω) ∩ H P (Ω)) , w ∈ L (0 , T ; L (Ω)) (cid:104) w ( t + τ ) − w ( t ) , v (cid:105) L = (cid:82) t + τt (cid:68) (cid:126)V ( u ( s )) · ∇ v, w ( s ) (cid:69) L + (cid:68) ˆ ku y ( s ) , v (cid:69) L ds (1) (cid:104) u ( s ) , v (cid:105) H = (cid:104) w ( s ) , v (cid:105) L , (2) ∀ v ∈ W , ∞ (Ω) ∩ H P (Ω) , ∀ t, τ, ≤ t < t + τ ≤ T, ∀ s ∈ [ t, t + τ ] (7)with u (0) = u ∈ H ∩ H P , and w (0) = w = u − ∆ u .Such formulation is well-suited for semi and full discretization of the original system (3) and (4).For semi-discretization, we start defining the finite element spaces as follows. Finite-Element Space Semi-Discretization
Systems (6) and (7) lead to equivalent P Finite-element space semi-discretization constructed as follows:Let P x = { x i | i = 1 , ..., n } be a partition of (0 , L ) : x < x < ... < x n = L in the x direction and similarly inthe y direction, P y = { y j | j = 1 , ..., n } . Let now: N = { P I ( x i , y j ) | I = 1 , , ..., N = n } = P x × P y , be a structured set of nodes covering Ω . Based on N , and as indicated in Figure 1, one obtains a conforming (Delaunay)Figure 1: A two-dimension meshing of Ω structured triangulation T of Ω , i.e., T = { E J | J = 1 , , ..., M } , Ω = ∪ J E J . The P finite element subspace X N of H (Ω) is given by: X N = { v ∈ C (Ω) | v restricted to E J ∈ P , J = 1 , .., M } ⊂ W ,pP , ≤ p ≤ ∞ with ∪ N ≥ { X N } approximating functions in H (Ω) . For that purpose, we let B N = { ϕ I | I = 1 , , ...N } be a finiteelement basis of functions with compact support in Ω , i.e.,: ∀ v N ∈ X N : v N ( x, y ) = N (cid:88) I =1 V I ϕ I ( x, y ) , V I = v N ( x I , y I ) . We state now useful estimates used in this paper. 3 pproximation properties of X N in H (Ω) These can be found in section 3.1 of Ciarlet [8], specifically: ∀ v ∈ W ,p we define π N ( v ) := N (cid:88) I =1 V I ϕ I ( x, y ) ∈ X N to be the interpolant of v in X N (8)One has the following estimates: ∀ v ∈ W ,p , | v − π N ( v ) | ,p ≤ C n | v | ,p , p ∈ (1 , ∞ ] (9)and ∀ v ∈ W , ∞ , | v − π N ( v ) | , ∞ ≤ C n | v | , ∞ , (10)Since the Hasegawa-Mima equation is set in H P (Ω) , we let X N,P = X N ∩ H P (Ω) . To discretize (7), we start with u N (0) = π N ( u ) , w N (0) = π N ( w ) , and w = u − ∆ u , then given ( u N ( t ) , w N ( t )) ∈ X N,P × X N,P where u N ( t ) = π N ( u ) = (cid:80) NJ =1 U J ( t ) ϕ J ( x, y ) , U J ( t ) = u N ( x J , y J , t ) , w N ( t ) = π N ( w ) = (cid:80) NI =1 W I ( t ) ϕ I ( x, y ) , and W I ( t ) = w N ( x I , y I , t ) , one seeks: u N ( t + τ ) ∈ X N,P , w N ( t + τ ) ∈ × X N,P (cid:104) w N ( t + τ ) − w N ( t ) , v (cid:105) L = (cid:82) t + τt (cid:68) (cid:126)V ( u N ( s )) · ∇ v, w N ( s ) (cid:69) L + (cid:68) ˆ ku N,y ( s ) , v (cid:69) L ds (1) (cid:104) u N ( s ) , v (cid:105) H = (cid:104) w N ( s ) , v (cid:105) L , (2) ∀ v ∈ X N,P ∀ t, τ, ≤ t < t + τ ≤ T, ∀ s ∈ { t, t + τ } (11)We now rewrite equation (11.1) by dividing it by τ giving τ (cid:104) w N ( t + τ ) − w N ( t ) , v (cid:105) L = 1 τ (cid:90) t + τt (cid:68) (cid:126)V ( u N ( s )) · ∇ v, w N ( s ) (cid:69) L + (cid:68) ˆ ku N,y ( s ) , v (cid:69) L ds Letting τ tend to , implies that every solution u N ( t ) = π N ( u ) , and w N ( t ) = π N ( w ) of (11) is a solution to thesemi-discretization of the H formulation (6) given by (cid:104) w N,t , v (cid:105) L − (cid:68) (cid:126)V ( u N ) · ∇ v, w N (cid:69) L = (cid:68) ˆ ku N,y , v (cid:69) L , ∀ v ∈ X N,P , ∀ t ∈ (0 , T ] (1) (cid:104) u N ( t ) , v (cid:105) H = (cid:104) w N ( t ) , v (cid:105) L , (2) ∀ v ∈ X N,P , ∀ t ∈ [0 , T ] (12)with u N (0) = π N ( u ) , w N (0) = π N ( w ) , w = u − ∆ u .Defining the associated vectors W ( t ) = { W I ( t ) } I and U ( t ) = { U J ( t ) } J , system (12) is equivalent in vector form to M dWdt + S ( U ) W = RU, ∀ t ∈ (0 , T ] (1) KU ( t ) = M W ( t ) ∀ t ∈ (0 , T ] (2) U (0) = U ; W (0) = W (3) (13)with M , K , S ( U ) and R , N × N matrices, defined as follows:• M = {(cid:104) ϕ I , ϕ J (cid:105) L | ≤ I, J ≤ N } , K = {(cid:104) ϕ I , ϕ J (cid:105) H | ≤ I, J ≤ N } , R = (cid:110)(cid:68) ˆ kϕ I,y , ϕ J (cid:69) L | ≤ I, J ≤ N (cid:111) • S ( U ) = (cid:110) − (cid:68) (cid:126)V ( u N ) · ∇ ϕ J , ϕ I (cid:69) L | ≤ I, J ≤ N (cid:111) = (cid:110)(cid:68) (cid:126)V ( u N ) · ∇ ϕ I , ϕ J (cid:69) L | ≤ I, J ≤ N (cid:111) ,where u N ( t ) = (cid:80) NK =1 U K ( t ) ϕ K ( x, y ) and ( (cid:126)V ( u n ) · ∇ ϕ I ) ϕ J = − ∂ϕ I ∂x ϕ J N (cid:88) K =1 (cid:18) U K ( t ) ∂ϕ K ∂y (cid:19) + ∂ϕ I ∂y ϕ J N (cid:88) K =1 (cid:18) U K ( t ) ∂ϕ K ∂x (cid:19) (14)When implementing system (13) one takes periodicity into account, reducing the degrees of freedom from N = n to N = ( n − . 4 tatement of Results Using a compactness technique, we prove in Section 2, the existence of a limit point ( u, w ) to the pair ( u N , w N ) andaccordingly, the existence of a solution to the Hasegawa-Mima coupled system which states as follows. Theorem 1.1.
Let u ∈ H (Ω) ∩ H P (Ω) and w = u − ∆ u ∈ L (Ω) . Then for all T > , there exists a uniquesolution ( u N ( t ) , w N ( t )) to (12). Furthermore the sequence { ( u N ( t ) , w N ( t )) , N > } admits a subsequence thatconverges to a solution pair ( u ( t ) , w ( t )) ∈ L (0 , T ; H (Ω) ∩ H P (Ω)) × L (0 , T ; L (Ω)) , such that: (cid:104) w ( t ) − w , v (cid:105) L = (cid:82) t (cid:68) (cid:126)V ( u ( s )) · ∇ v, w ( s ) (cid:69) L + ˆ k (cid:104) u y ( s ) , v (cid:105) L ds (1) (cid:104) u ( t ) , v (cid:105) H = (cid:104) w ( t ) , v (cid:105) L , (2) ∀ v ∈ H P (Ω) ∩ W , ∞ (Ω) , ∀ t ∈ [0 , T ] (15)In Section 3 we introduce the fully implicit nonlinear discrete scheme (58), and prove the existence and uniqueness ofits solution under a restriction on the time step. In Section 4, we complete the discretization cycle by presenting theresulting algorithm and its implementation using FreeFem++ software, generating in the sequel the system matrices M, K, S ( U ) , R . The obtained numerical results indicate the robustness of this new software, particularly in handlingthe complex discretization of the Poisson bracket and circumventing the difficulties encountered in [7]. Concludingremarks are provided in Section 5. At the core of the proof of this result, are1. Existence and uniqueness of solutions to (12) and equivalently to (13) proven in Section 2.1 which uses astandard existence theorem for systems of ODE’s of the form (cid:126)Y (cid:48) ( t ) = (cid:126)F ( (cid:126)Y ) , (cid:126)F Lipschitzian.2. A-priori estimates on these solutions shown in Section 2.2.3. In Section 2.3, a compactness argument would allow passing to the limit for test functions v ∈ W , ∞ (Ω) ∩ H P (Ω) in the formulation (15) of Theorem 1.1. Finally, a density argument of W , ∞ (Ω) ∩ H P (Ω) in H P (Ω) allows us to complete the proof.This thread of items to be proven requires skew-symmetry results obtained in the following section. Preliminary Result: Skew-Symmetry on X N,P
For that purpose, we start by obtaining a skew-symmetry result, stated in the following proposition.
Theorem 2.1.
For all { v, z, φ } ∈ X N,P × X N,P × X N,P , one has: (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L = − (cid:68) (cid:126)V ( v ) · ∇ φ, z (cid:69) L . (16)To prove Theorem (2.1), we start by breaking (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L into a sum of integrals over each triangle E J ∈ T , J = 1 , ..., M . Specifically: (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L = (cid:88) E J ∈T (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L ( E J ) (17) Lemma 2.2.
Let E J ∈ T with vertices J , J , J . Let also ν being the unit outer normal to ∂E J , defined piecewiseon each of the sides of the triangle E J and denoted respectively by ν , ν , ν on [1 J , J ] , [2 J , J ] and [3 J , J ] . Thenone has: (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L ( E J ) = − (cid:68) (cid:126)V ( v ) · ∇ φ, z (cid:69) L ( E J ) + γ ,J (cid:90) J J zφ + γ ,J (cid:90) J J zφ + γ ,J (cid:90) J J zφ, ith: γ ,J = [( v (2 J ) − v (1 J ))( ν ,x y J − y J − ν ,y x J − x J )] γ ,J = [( v (3 J ) − v (2 J ))( ν ,x y J − y J − ν ,y x J − x J )] γ ,J = [( v (1 J ) − v (3 J ))( ν ,x y J − y J − ν ,y x J − x J )] (18) Proof.
Using Green’s formula, where ν is the outer normal on ∂E J , we obtain: (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L ( E J ) = (cid:90) ∂E J zν · ( φ(cid:126)V ( v )) − (cid:90) E J z ∇ . ( φ(cid:126)V ( v )) Since (cid:126)V ( v ) is divergence free, then: (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L ( E J ) = (cid:90) ∂E J zν · ( φ(cid:126)V ( v )) − (cid:90) E J z ∇ . ( φ(cid:126)V ( v )) = (cid:90) ∂E J zν · ( φ(cid:126)V ( v )) − (cid:68) (cid:126)V ( v ) · ∇ φ, z (cid:69) L ( E J ) (19)Now the triangle boundary integral can be expressed as follows: (cid:90) ∂E J zν · ( φ(cid:126)V ( v )) = (cid:90) J J zφν · (cid:126)V ( v ) + (cid:90) J J zφν · (cid:126)V ( v ) + (cid:90) J J zφν · (cid:126)V ( v ) Handling as a sample one of these line integrals, one has for (example on [1 J , J ] ), ν · (cid:126)V ( v ) = ν ,x .v y − ν ,y .v x .Furthermore as v ∈ P , ν ,x .v y − ν ,y .v x is a constant on ( −−−→ J J ) and given by: ν ,x .v y − ν ,y .v x = ν ,x v (2 J ) − v (1 J ) y J − y J − ν ,y v (2 J ) − v (1 J ) x J − x J = ( v (2 J ) − v (1 J ))( ν ,x y J − y J − ν ,y x J − x J ) . Hence: (cid:90) J J zφν · (cid:126)V ( v ) = ( v (2 J ) − v (1 J ))( ν ,x y J − y J − ν ,y x J − x J ) (cid:90) J J zφ = γ ,J (cid:90) J J zφ, with similar identities obtained for (cid:82) J J zφν · (cid:126)V ( v ) and (cid:82) J J zφν · (cid:126)V ( v ) . Replacing these integrals by their expres-sions in (19), one obtains the result of this lemma.The next step is to consider the sum (cid:80) E J ∈T (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L ( E J ) and demonstrate the identity. Lemma 2.3. (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L = − (cid:68) (cid:126)V ( v ) · ∇ φ, z (cid:69) L + (cid:90) ∂ Ω zφ ν.(cid:126)V ( v ) . Proof.
On the basis of (17) and of lemma 2.2, one has: (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L = − (cid:88) E J ∈T (cid:68) (cid:126)V ( v ) · ∇ φ, z (cid:69) L ( E J ) + (cid:88) E J ∈T [ γ ,J (cid:90) J J zφ + γ ,J (cid:90) J J zφ + γ ,J (cid:90) J J zφ ] . (20)Given that any internal side [ AB ] of any triangle E J is also common to another triangle E K , then the correspondingline integral from E J is given by γ AB (cid:82) BA zφ and that coming from E K is γ BA (cid:82) AB zφ , with γ AB = − γ BA , leading toa zero sum for integrals on [ AB ] .Consequently, one is left on the right hand side of (20) with line integrals over ∂ Ω , i.e. the result of the lemma.Using the results of Lemmas (2.2) and (2.3), the we can complete the proof of Theorem 2.1 as follows.6 roof. The periodicity of v and z on ∂ Ω results in (cid:82) ∂ Ω zφ ν.(cid:126)V ( v ) = 0 , because of the periodicity of zφ on ∂ Ω and thefact that ν.(cid:126)V ( v ) is given by ± v x on the horizontal sides and by ± v y on the vertical sides.For example, handling the vertical sides gives (cid:90) (1 , , [ zφ ν.(cid:126)V ( v )] dx + (cid:90) (0 , , [ zφ ν.(cid:126)V ( v )] dx = (cid:90) (1 , , − zφv x dx + (cid:90) (0 , , zφv x (cid:48) dx (cid:48) = (cid:90) (1 , , − zφv x dx + (cid:90) (1 , , zφv − x d ( − x )= (cid:90) (1 , , − zφv x dx + (cid:90) (1 , , zφv x d ( x ) = 0 Then using Lemma (2.3) we complete the proof of Skew-symmetry.This obviously leads to the following corollary.
Corollary 2.4. ∀{ v, z } ∈ X N,P × X N,P , we have that (cid:68) (cid:126)V ( v ) · ∇ z, z (cid:69) L = 0 (12) (& (13) ) In the rest of the paper, we will use the following norms in R N for V ∈ R N : || V || M := V T M V (21) || V || K := V T K V (22)By associating the function v N ( x, y ) = (cid:80) NI =1 V I ϕ I ( x, y ) ∈ X N,P to V ∈ R N , then we have the isometries: || V || M = || v N || (23) || V || K = || v N || (24)We begin by obtaining relations between solutions to (12) and (13). For W ( t ) , U ( t ) ∈ R N with w N ( x, y, t ) = N (cid:88) I =1 W I ( t ) ϕ I ( x, y ) ∈ X N,P and u N ( x, y ) = N (cid:88) I =1 U I ( t ) ϕ I ( x, y ) ∈ X N,P one has the following:
Lemma 2.5.
Any pair ( U, W ) that solves (13.2) satisfies the following: || U ( t ) || M ≤ || U ( t ) || K ≤ || W ( t ) || M Proof.
From the equation KU ( t ) = M W ( t ) ⇐⇒ < u N , v > = < w N , v >, ∀ v ∈ X N,P . Letting v = u N ( t ) , yields || u N ( t ) || = < w N , u N > ≤ || w N ( t ) || . || u N ( t ) || , leading to: || u N ( t ) || ≤ || u N ( t ) || = < w N , u N > ≤ || w N ( t ) || . || u N ( t ) || = ⇒ || u N ( t ) || ≤ || w N ( t ) || , which translates to the result via the isometries (23) and (24) .The existence of a unique solution to the semi-discrete system can be obtained by reducing (12) (or (13)) to a systemof non-linear ordinary differential equations in W ( t ) . Specifically, in (13), we eliminate the variable U ( t ) , using thematrix A = K − M and obtain the system: (cid:26) M dWdt + S ( AW ) W = RAW, ∀ t ∈ (0 , T ] (1) W (0) = W (2) (25)7hich is equivalent to: (cid:26) dWdt = F ( W ( t )) , ∀ t ∈ (0 , T ] (1) W (0) = W (2) (26)where F ( W ) = M − [ RAW − S ( AW ) W ] . Now we show that F is locally Lipschitz on the spaces X N = { V ∈ R N | || V || M,C ([0 ,T ]; R N ) ≤ C T } where C T := e || ˆ k || ∞ T || w || is determined in Lemma 2.8. Lemma 2.6.
For W , W ∈ X N , there exists a positive constant L T independent from h such that || F ( W ) − F ( W ) || M ≤ L T h √ h || W − W || M Proof.
Let Z k = F ( W k ) with U k = AW k for k = 1 , .Let z k,N ( x, y ) = (cid:80) NI =1 Z k,I ( t ) ϕ I ( x, y ) , u k,N ( x, y ) = (cid:80) NI =1 U k,I ( t ) ϕ I ( x, y ) , w k,N ( x, y ) = (cid:80) NI =1 W k,I ( t ) ϕ I ( x, y ) .Then, M ( Z − Z ) = R ( U − U ) − ( S ( U ) W − S ( U ) W ) is equivalent to < z ,N − z ,N , φ > = − < (cid:126)V ( u ,N ) . ∇ w ,N − (cid:126)V ( u ,N ) . ∇ w ,N , φ > + < ˆ k ( u ,N − u ,N ) y , φ > for all φ ∈ X N,P .Let φ = z ,N − z ,N then by Cauchy-Schwartz we get that || z ,N − z ,N || = − < (cid:126)V ( u ,N ) . ∇ w ,N − (cid:126)V ( u ,N ) . ∇ w ,N , z ,N − z ,N > + < ˆ k ( u ,N − u ,N ) y , z ,N − z ,N > ≤ || (cid:126)V ( u ,N ) . ∇ w ,N − (cid:126)V ( u ,N ) . ∇ w ,N || || z ,N − z ,N || + || ˆ k ( u ,N − u ,N ) y || || z ,N − z ,N || Simplifying by || z ,N − z ,N || we get || z ,N − z ,N || ≤ || (cid:126)V ( u ,N ) . ∇ w ,N − (cid:126)V ( u ,N ) . ∇ w ,N || + || ˆ k ( u ,N − u ,N ) y || (27)Note that || ˆ k ( u ,N − u ,N ) y || ≤ || ˆ k || ∞ || ( u ,N − u ,N ) || H ≤ || ˆ k || ∞ || ( w ,N − w ,N ) || (28)where we have used the fact that K ( U − U ) = M ( W − W ) = ⇒ || u ,N − u ,N || H ≤ || w ,N − w ,N || .On the other hand, using the triangle inequality || (cid:126)V ( u ,N ) . ∇ w ,N − (cid:126)V ( u ,N ) . ∇ w ,N || ≤ || (cid:126)V ( u ,N − u ,N ) . ∇ w ,N || + || (cid:126)V ( u ,N ) . ∇ ( w ,N − w ,N ) || (29)Note that for any ( α N , β N ) ∈ X N,p × X N,p , one has || (cid:126)V ( α N ) . ∇ β N || ≤ max φ ∈ X N,P , || φ || =1 | < (cid:126)V ( α N ) . ∇ β N , φ > | (30)Using skew-symmetry, < (cid:126)V ( α N ) . ∇ β N , φ > = − < (cid:126)V ( α N ) . ∇ φ, β N > , then (30) becomes || (cid:126)V ( α N ) . ∇ β N || ≤ max φ ∈ X N,P , || φ || =1 | < (cid:126)V ( α N ) . ∇ φ, β N > | (31)We need now the following Lemma. 8 emma 2.7. For some constant C independent from h we have: || (cid:126)V ( α N ) . ∇ β N || ≤ Ch √ h || α N || || β N || Proof. | < (cid:126)V ( α N ) . ∇ φ, β N > | ≤ || (cid:126)V ( α N ) . ∇ φ || || β N || ≤ | φ | , ∞ || α N || || β N || ≤ Ch √ h || α N || || β N || , as using a result in Ciarlet ([8], Theorem 3.2.6), one has | φ | , ∞ ≤ Ch − / || φ || . We complete the proof, by applying this lemma twice to the right-hand side of (29) thus obtaining || (cid:126)V ( u ,N ) . ∇ w ,N − (cid:126)V ( u ,N ) . ∇ w ,N || ≤ || (cid:126)V ( u ,N − u ,N ) . ∇ w ,N || + || (cid:126)V ( u ,N ) . ∇ ( w ,N − w ,N ) ||≤ Ch √ h || u ,N − u ,N || || w ,N || + Ch √ h || u ,N || || w ,N − w ,N || (32)Plugging this inequality and (28) in (27) leads to: || z ,N − z ,N || ≤ || (cid:126)V ( u ,N ) . ∇ w ,N − (cid:126)V ( u ,N ) . ∇ w ,N || + || ˆ k || ∞ || ( w ,N − w ,N ) ||≤ Ch √ h || u ,N − u ,N || || w ,N || + Ch √ h || u ,N || || w ,N − w ,N || + || ˆ k || ∞ || ( w ,N − w ,N ) || Using the assumption that W , W ∈ X N , the last inequality leads to || z ,N − z ,N || ≤ CC T h √ h || w ,N − w ,N || + CC T h √ h || w ,N − w ,N || + || ˆ k || ∞ || w ,N − w ,N ||≤ (2 CC T h √ h + || ˆ k || ∞ ) || w ,N − w ,N || ≤ L T h √ h || w ,N − w ,N || (33)Now squaring and multiplying both sides by M T from the left, and using the equivalence of the L norm on X N,p andthe M norm on R N (23) completes the proof of Lemma 2.6.Thus, the semi-discrete system (26) has a unique solution W ( t ) or a solution pair { U ( t ) , W ( t ) } , for which we derivesome a priori estimates. (12) We may now state some a priori estimates.
Lemma 2.8.
Every unique solution { u N , w N } to (12), satisfies the following estimates: (cid:40) || w N || C ([0 ,T ]; L (Ω)) = max t ∈ [0 ,T ] || w N || ( t ) ≤ C T := e || ˆ k || ∞ T || w || (1) || u N || C ([0 ,T ]; H P (Ω)) = max t ∈ [0 ,T ] || u N || ( t ) ≤ C T := e || ˆ k || ∞ T || w || (2) (34) Proof.
In (12.1), let v = w N . Then one has: (cid:104) w N,t , w N (cid:105) L + (cid:68) (cid:126)V ( u N ) · ∇ w N , w N (cid:69) L = (cid:68) ˆ ku N,y , w N (cid:69) L ∀ t ∈ (0 , T ] , (35)and letting v = u N in (12.2), one obtains also: (cid:104) u N ( t ) , u N ( t ) (cid:105) H = (cid:104) w N ( t ) , u N ( t ) (cid:105) L , ∀ t ∈ (0 , T ] . (36)9sing Lemma 2.2 and the fact that (cid:104) w N,t , w N (cid:105) L = ddt || w N || ( t ) , then equations (35) and (36) lead to: ddt || w N || ( t ) ≤ || ˆ k || ∞ || w N || ( t ) || u N || ( t ) , ∀ t ∈ (0 , T ] || u N || ( t ) ≤ || w N || ( t ) || u N || ( t ) , ≤ || w N || ( t ) || u N || ( t ) , ∀ t ∈ [0 , T ] given the initial choice of u N (0) Hence from these two inequalities, one gets ∀ t ∈ [0 , T ] : || u N || ( t ) ≤ || w N || ( t ) , (37) ddt || w N || ( t ) ≤ || ˆ k || ∞ || w N || ( t ) . (38)Integration of the differential inequality (38) gives: || w N || ( t ) ≤ e || ˆ k || ∞ t || w N || (0) , ∀ t ∈ [0 , T ] ≤ e || ˆ k || ∞ T || w N || (0) , ∀ t ∈ [0 , T ] ∴ || w N || ( t ) ≤ e || ˆ k || ∞ T || w N || (0) , ∀ t ∈ [0 , T ] (39)Thus inequalities (37) and (39) give the results of the lemma, where || w N || (0) = || π N ( w ) || ≤ || w || . At this point we introduce the sequence { z N } , defined by: z N ( t ) ∈ H P (Ω) ∩ H (Ω) : − ∆ z N ( t ) + z N ( t ) = w N ( t ) with || z N || ( t ) ≤ C || w N || ( t ) (40)Note that in this case, the finite element approximation to z n in X N,P is u N ( t ) . Using the well-known Cea’s estimatefor elliptic problems ([8], Theorem 3.2.2) in addition to (40), one has ∀ t ∈ [0 , T ] : || z N ( t ) − u N ( t ) || ≤ || z N ( t ) − π N ( u N )( t ) || ≤ Cn | z N ( t ) | , ≤ Cn || w N ( t ) || (41)with C a generic constant independent of N . Thus, instead of studying the convergence of the pair ( u N , w N ) , westudy that of ( z N , w N ) .Following (34.1) and (41.1), one concludes that: || z N ( t ) || C ([0 ,T ]; H ) ≤ Ce || ˆ k || ∞ T/ || π N ( w ) || . Lemma 2.9.
There exists an element u ∈ L (0 , T ; H P ) and a subsequence { z N i } ⊂ { z N } , such that: lim N i →∞ || u − z N i || L (0 ,T ; H ) = 0 . Proof.
This result follows from the Rellich-Kondrachov theorem that stipulates the compact injection of H (Ω) in H (Ω) ([9], p.285). .Let us now denote ( z N i , w N i ) by ( z N , w N ) and seek first a limit point to the sequence { w N } . Specifically, we havethe following result. Lemma 2.10.
There exists w ∈ L (0 , T ; L (Ω)) and a subsequence { w N j } of { w N } such that: w N j ( t ) (cid:42)w ( t ) in L (0 , T ; L (Ω)) (42) w N j ( t ) (cid:42)w ( t ) in L (Ω) for all t ∈ [0 , T ] (43) Furthermore || w || L (0 ,T ; L (Ω)) ≤ e || ˆ k || ∞ T || w || . roof. Observe via (34.1) that the sequence { w N } is uniformly bounded in the reflexive space L (0 , T ; L (Ω)) , andso it has a subsequence { w N j } which converges weakly, say to w ∈ L (0 , T ; L (Ω)) . i.e. (cid:90) T (cid:104) w N ( t ) , v (cid:105) dt −→ (cid:90) T (cid:104) w ( t ) , v (cid:105) dt for all v ∈ L (0 , T ; L (Ω)) Now fix v ∈ L (Ω) and consider the sequence F j ( t ) := (cid:82) t (cid:10) w N j ( s ) , v (cid:11) ds of functions on [0 , T ] . Observe that:1. F j ( t ) ∈ C (0 , T ) and F (cid:48) j ( t ) = (cid:10) w N j ( t ) , v (cid:11) for all j by the Fundamental Theorem of Calculus.2. F j ( t ) −→ F ( t ) := (cid:82) t (cid:104) w ( s ) , v (cid:105) ds pointwise on [0 , T ] .3. (cid:107) F j ( t ) (cid:107) ’s are uniformly bounded on [0 , T ] by (34.1).4. { F j ( t ) } are uniformly equicontinuous on [0 , T ] as | F j ( t ) − F j ( s ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) ts (cid:10) w N j ( τ ) , v (cid:11) dτ (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:90) ts (cid:12)(cid:12)(cid:10) w N j ( τ ) , v (cid:11)(cid:12)(cid:12) dτ ≤ (cid:90) ts (cid:13)(cid:13) w N j ( τ ) (cid:13)(cid:13) (cid:107) v (cid:107) dτ ≤ C T (cid:107) v (cid:107) | t − s | so that by Arzel`a-Ascoli theorem, F j ( t ) has a subsequence F j k ( t ) that converges to F ( t ) uniformly on [0 , T ] , and so (cid:68) w N jk ( t ) , v (cid:69) = F (cid:48) j k ( t ) −→ F (cid:48) ( t ) = (cid:104) w ( t ) , v (cid:105) for every t ∈ [0 , T ] , which gives (43) after relabelling.Finally, weakly lower-semicontinuity of norms implies that || w || L (0 ,T ; L (Ω)) ≤ lim inf N →∞ || w N || L (0 ,T ; L (Ω)) ≤ e || ˆ k || ∞ T || w || To complete the proof of Theorem 1.1, we denote the pair ( u N j , w N j ) by ( u N , w N ) and aim at proving that the limitpair ( u, w ) satisfies (cid:40) (cid:104) w ( t ) − w , v (cid:105) L = (cid:82) t (cid:68) (cid:126)V ( u ( s )) · ∇ v, w ( s ) (cid:69) L + ˆ k (cid:104) u y ( s ) , v (cid:105) L ds, (1) (cid:104) u, v (cid:105) H = (cid:104) w, v (cid:105) L , ∀ v ∈ H P (Ω) , (2) (44) ∀ v ∈ H P (Ω) ∩ W , ∞ (Ω) , ∀ t ∈ (0 , T ] with w (0) = π N ( w ) .For that purpose, we replace in (11) t by and t + τ by t and simultaneously use the skew-symmetry property inTheorem 2.1, getting consequently: (cid:104) w N ( t ) − w N (0) , v N (cid:105) L − (cid:82) t (cid:68) (cid:126)V ( u N ( s )) · ∇ v N , w N ( s ) (cid:69) L ds = (cid:82) t (cid:68) ˆ ku N,y ( s ) , v N (cid:69) L ds, (1) (cid:104) u N ( t ) , v N (cid:105) H = (cid:104) w N ( t ) , v N (cid:105) L , (2) ∀ v N ∈ X N,P , ∀ t ∈ (0 , T ] (45)with w N (0) = π N ( u − ∆ u ) , which is the direct semi-discretization of (15).To consider limit points when N → ∞ of each of the terms in (45.1) and (45.2), the following sequence of lemmas isneeded in which C ( T ) is a generic constant of the form ( aT + b ) e d || ˆ k || ∞ T independent from n where a, b, d ∈ N . Lemma 2.11.
For all v ∈ H P (Ω) ∩ W , ∞ (Ω) , and for all t ∈ [0 , T ] , one has: (cid:104) w N ( t ) , v (cid:105) L = (cid:104) w N ( t ) , π N ( v ) (cid:105) L + (cid:15) N, ( t ) with | (cid:15) N, ( t ) | ≤ C ( T ) 1 n || w || | v | (46) where (cid:15) N, ( t ) = (cid:104) w N ( t ) , π N ( v ) − v (cid:105) L roof. Given the identity: (cid:104) w N ( t ) , π N ( v ) (cid:105) L = (cid:104) w N ( t ) , v (cid:105) L + (cid:104) w N ( t ) , π N ( v ) − v (cid:105) L , and letting: (cid:15) N, ( t ) = (cid:104) w N ( t ) , v N − v (cid:105) L , one has using (9) and (34): | (cid:15) N, ( t ) | ≤ || w N ( t ) || . || v − π N ( v ) || , (47) ≤ e || ˆ k || ∞ T || w || C n | v | ≤ C ( T ) 1 n || w || | v | (48)where C ( T ) = C e || ˆ k || ∞ T .Similarly, one has: Lemma 2.12.
For all v ∈ H P (Ω) ∩ W , ∞ (Ω) , one has: (cid:90) t (cid:68) ˆ ku N,y ( s ) , v (cid:69) L ds = (cid:90) t (cid:68) ˆ ku N,y ( s ) , π N ( v ) (cid:69) L ds + (cid:15) N, ( t ) with | (cid:15) N, ( t ) | ≤ C ( T ) 1 n || w || | v | (49) Proof.
Given the identity: (cid:90) t (cid:68) ˆ ku N,y ( s ) , v (cid:69) L ds = (cid:90) t (cid:68) ˆ ku N,y ( s ) , π N ( v ) (cid:69) L ds + (cid:90) t (cid:68) ˆ ku N,y ( s ) , v − π N ( v ) (cid:69) L , and letting: (cid:15) N, ( t ) = (cid:82) t (cid:68) ˆ ku N,y ( s ) , v − π N ( v ) (cid:69) L , one has using (9) and (34),: | (cid:15) N, ( t ) | ≤ || v − π N ( v ) || . (cid:90) t || u N ( s ) || ds, (50) ≤ C n | v | T e || ˆ k || ∞ T || w || ≤ C ( T ) 1 n || w || | v | (51)where C ( T ) = T C e || ˆ k || ∞ T We turn now to the key term in (45) and prove the following.
Lemma 2.13.
For all v ∈ H P (Ω) ∩ W , ∞ (Ω) , one has: (cid:90) t (cid:68) (cid:126)V ( z N ( s )) · ∇ v, w N ( s ) (cid:69) L ds = (cid:90) t (cid:68) (cid:126)V ( u N ( s )) · ∇ π N ( v ) , w N ( s ) (cid:69) L ds + (cid:15) N, ( t ) (52) where (cid:15) N, ( t ) = (cid:90) t (cid:68) (cid:126)V ( z N ( s ) − u N ( s )) · ∇ v, w N ( s ) (cid:69) L ds + (cid:90) t (cid:68) (cid:126)V ( u N ( s )) · ∇ ( v − π N ( v )) , w N ( s ) (cid:69) L ds, with | (cid:15) N, ( t ) | ≤ C ( T ) 1 n || w || max {| v | , ∞ , | v | , ∞ } . Proof.
Adding up the following two identities and using the definition of (cid:15) N, ( t ) : (cid:90) t (cid:68) (cid:126)V ( z N ( s )) · ∇ v, w N ( s ) (cid:69) L ds = (cid:90) t (cid:68) (cid:126)V ( z N ( s ) − u N ( s )) · ∇ v, w N ( s ) (cid:69) L ds + (cid:90) t (cid:68) (cid:126)V ( u N ( s )) · ∇ v, w N ( s ) (cid:69) L ds (cid:90) t (cid:68) (cid:126)V ( u N ( s )) · ∇ v, w N ( s ) (cid:69) L ds = (cid:90) t (cid:68) (cid:126)V ( u N ( s )) · ∇ π N ( v ) , w N ( s ) (cid:69) L ds + (cid:90) t (cid:68) (cid:126)V ( u N ( s )) · ∇ ( v − π N ( v )) , w N ( s ) (cid:69) L ds. yields: (cid:90) t (cid:68) (cid:126)V ( z N ( s )) · ∇ v, w N ( s ) (cid:69) L ds = (cid:90) t (cid:68) (cid:126)V ( u N ( s )) · ∇ π N ( v ) , w N ( s ) (cid:69) L ds + (cid:15) N, ( t ) . | (cid:15) N, ( t ) | ≤ | v − π N ( v ) | , ∞ . (cid:90) t | u N | ( s ) . || w N ( s ) || ds + | v | , ∞ . (cid:90) t | u N − z N | ( s ) . || w N ( s ) || ds ≤ Cn | v | , ∞ (cid:90) t | u N | ( s ) . || w N ( s ) || ds + | v | , ∞ . (cid:90) t | u N − z N | ( s ) . || w N ( s ) || ds ≤ | v | , ∞ Cn T e || ˆ k || ∞ T || w || + | v | , ∞ Cn T e || ˆ k || ∞ T || w || ≤ C ( T ) 1 n || w || max {| v | , ∞ , | v | , ∞ } where C ( T ) = 2 T C e || ˆ k || ∞ T .In a similar way to Lemma 2.11, one can prove the following lemma: Lemma 2.14.
For all v ∈ H P (Ω) ∩ W , ∞ (Ω) , and for all t ∈ [0 , T ] one has: (cid:104) u N ( t ) , v (cid:105) H = (cid:104) u N ( t ) , π N ( v ) (cid:105) H + (cid:15) N, ( t ) with | (cid:15) N, ( t ) | ≤ C ( T ) 1 n || w || | v | (53) where (cid:15) N, ( t ) = (cid:104) u N ( t ) , π N ( v ) − v (cid:105) L . Synthesis: Completion of Proof of Existence (Theorem 1.1)
Using (46), (49), and (52) one gets: (cid:104) w N ( t ) − w N (0) , v (cid:105) L = (cid:104) w N ( t ) − w N (0) , π N ( v ) (cid:105) L + (cid:15) N, ( t ) − (cid:15) N, (0) (54) − (cid:90) t (cid:68) ˆ ku N,y ( s ) , v (cid:69) L ds = − (cid:90) t (cid:68) ˆ ku N,y ( s ) , π N ( v ) (cid:69) L ds − (cid:15) N, ( t ) (55) − (cid:90) t (cid:68) (cid:126)V ( z N ( s )) · ∇ v, w N ( s ) (cid:69) L ds = − (cid:90) t (cid:68) (cid:126)V ( u N ( s )) · ∇ π N ( v ) , w N ( s ) (cid:69) L ds − (cid:15) N, ( t ) (56)Then, by summing up (54), (56), (55), and using (45.1) , we get ∀ t ∈ (0 , T ) : (cid:104) w N ( t ) − w N (0) , v (cid:105) L − (cid:90) t (cid:68) (cid:126)V ( z N ( s )) · ∇ v, w N ( s ) (cid:69) L + (cid:68) ˆ ku N,y ( s ) , v (cid:69) L ds = (cid:15) N, ( t ) − (cid:15) N, (0) − (cid:88) k =2 (cid:15) N,k ( t ) (57)We may now let N → ∞ . Using the previous lemmas in this section, we have subsequently:1. For the right hand side of (57), using Lemmas 2.11, 2.12 and 2.13: lim N →∞ (cid:15) N, ( t ) − (cid:15) N, (0) − (cid:88) k =2 (cid:15) N,k ( t ) = 0 , ∀ v ∈ W , ∞ (Ω) ∩ H P (Ω)
2. For the left-hand side of (57):• lim N →∞ (cid:104) w N ( t ) − w N (0) , v (cid:105) L = (cid:104) w ( t ) − w (0) , v (cid:105) L , using Lemma 2.11.• lim N →∞ (cid:82) t (cid:68) ˆ ku N,y ( s ) , v (cid:69) L ds = (cid:82) t (cid:68) ˆ ku y ( s ) , v (cid:69) L ds , using Lemma 2.12.• For the term (cid:82) t (cid:68) (cid:126)V ( z N ( s )) · ∇ v, w N ( s ) (cid:69) L ds , note that for all s ∈ [0 , t ] and for all v ∈ W , ∞ (Ω) , (cid:126)V ( z N ( s )) . ∇ v converges to (cid:126)V ( u ( s )) . ∇ v strongly in L (0 , T ; L (Ω)) , as (cid:90) t (cid:68) (cid:126)V ( z N ( s ) − u ( s )) ∇ v, ϕ (cid:69) ds ≤ (cid:90) t (cid:107) z N ( s ) − u ( s ) (cid:107) (cid:107) v (cid:107) W , ∞ (cid:107) ϕ (cid:107) ds → as N → Combining this with the weak convergence of w N to w in L (0 , T ; L (Ω)) , we obtain lim N →∞ (cid:90) t (cid:68) (cid:126)V ( z N ( s )) · ∇ v, w N ( s ) (cid:69) L ds = (cid:90) t (cid:68) (cid:126)V ( u ( s )) · ∇ v, w ( s ) (cid:69) L ds
13. By applying similar techniques with respect to (cid:104) u N ( t ) , v N (cid:105) H = (cid:104) w N ( t ) , v N (cid:105) L , ∀ v N ∈ X N,P , ∀ t ∈ [0 , T ] in (45.2), using Lemmas 2.11 and 2.14 one obtains as N → ∞ : (cid:104) u ( t ) , v (cid:105) H = (cid:104) w ( t ) , v (cid:105) L , ∀ v ∈ H P (Ω) , ∀ t ∈ [0 , T ] . Thus these last 3 consecutive points prove formulation (15) for test functions v ∈ W , ∞ (Ω) ∩ H P (Ω) .Finally, the density of W , ∞ (Ω) ∩ H P (Ω) in W , ∞ (Ω) ∩ H P (Ω) completes the proof of Theorem 1.1. As to fully discretizing the Hasegawa-Mima system, starting with (7) and to avoid any constraint of the CFL type onthe choice of τ , the term (cid:82) t + τt (cid:68) (cid:126)V ( u ( s )) · ∇ v, w ( s ) (cid:69) L is first discretized using an implicit right rectangular rule: (cid:90) t + τt (cid:68) (cid:126)V ( u ( s )) · ∇ v, w ( s ) (cid:69) L = τ (cid:68) (cid:126)V ( u ( t + τ )) · ∇ v, w ( t + τ ) (cid:69) L + (cid:15) τ , leading to the following fully implicit Computational Model. Given ( u N ( t ) , w N ( t )) ∈ X N,P × X N,P , one seeks ( u N ( t + τ ) , w N ( t + τ )) ∈ X N,P × X N,P , such that: (cid:104) w N ( t + τ ) − w N ( t ) , v (cid:105) L = τ (cid:68) (cid:126)V ( u N ( t + τ )) · ∇ v, w N ( t + τ ) (cid:69) L + τ (cid:68) ˆ ku N,y ( t + τ ) , v (cid:69) L , (1) (cid:104) u N ( s ) , v (cid:105) H = (cid:104) w N ( s ) , v (cid:105) L , (2) ∀ v ∈ X N,P , ∀ s ∈ { t, t + τ } (58)In matrix notations and using the expressions: w N ( t ) = N (cid:88) I =1 W I ( t ) ϕ I ( x, y ) , and u N ( x, y, t ) = N (cid:88) J =1 U J ( t ) ϕ J ( x, y ) , where W I ( t ) = w N ( x I , y I , t ) , and U J ( t ) = w N ( x J , y J , t ) , then (58) can be rewritten as follows:Given ( U ( t ) , W ( t )) ∈ R N × R N , seek ( U ( t + τ ) , W ( t + τ )) ∈ R N × R N , such that: (cid:26) ( M + τ S ( U ( t + τ )) W ( t + τ ) − τ R U ( t + τ ) = M W ( t ) (1) KU ( s ) = M W ( s ) , ∀ s ∈ { t, t + τ } (2) (59)In Section 3.1, using a fixed point approach we start by showing the existence of solution to (59), i.e. (58). Then weprove uniqueness of this solution in Section 3.2. To prove the existence of a solution to the Fully Discrete System (59), we start by transforming it into the fixed pointproblem (65). Then, we prove the existence of a solution to (65) using the Leray-Schauder fixed point Theorem [10].Using equation (59.2) for s = t + τ , one gets U ( t + τ ) = K − M W ( t + τ ) . Substituting U ( t + τ ) in equation (59.1),system (59) can be rewritten as follows (cid:2) M + τ S ( K − M W ( t + τ )) − τ RK − M (cid:3) W ( t + τ ) = M W ( t ) , (60)which is equivalent to: (cid:2) M − τ RK − M (cid:3) W ( t + τ ) = M W ( t ) − τ S ( K − M W ( t + τ )) W ( t + τ ) , (61)Let B = M − τ RK − M , Z = W ( t ) and Y = W ( t + τ ) , then we get the first fixed point BY = M Z − τ S ( K − M Y ) Y, (62)14 heorem 3.1. Define the bilinear form on H (Ω) × H (Ω) : a τ ( w, v ) = < w, v > − τ < ˆ kw y , v > . Then, for τ ≤ || ˆ k || ∞ this bilinear form is coercive in the sense that | a τ ( w, w ) | ≥ || w || , ∀ w ∈ H (Ω) Proof.
This simply results from: | a τ ( w, w ) | = | < w, w > − τ < ˆ kw y , w > | ≥ (cid:16) − τ || ˆ k || ∞ (cid:17) || w || ≥ || w || , f or τ ≤ || ˆ k || ∞ As a consequence, one gets the following corollary.
Corollary 3.2.
The matrix B is invertible for τ ≤ || ˆ k || ∞ .Proof. Bα = β ⇐⇒ ( I − τ RK − ) M α = β ⇐⇒ ( I − τ RK − ) α (1) = β (63)where α (1) = M α . Let α (2) = K − α (1) , i.e. α (1) = Kα (2) , and β (1) = M − β . Hence, ( K − τ R ) α (2) = M β (1)
In variational form, this is equivalent to < α (2) N , v > − τ < ˆ kα (2) N,y , v > = < β (1) N , v > (64)for all v ∈ X N,P , with β (1) N = N (cid:80) i =1 β (1) i ϕ i , α (2) N = N (cid:80) i =1 α (2) i ϕ i , and α (2) N,y = N (cid:80) i =1 α (2) i,y ϕ i .Since from theorem 3.1, the bilinear form a τ ( w, v ) = < w, v > − τ < ˆ kw y , v > is coercive, then by applyingLax-Milgram on (64), one completes the proof of the corollary.Given that B is invertible for τ ≤ || ˆ k || ∞ , the fixed point format (62) becomes Y = G ( Y ) = B − (cid:2) M Z − τ S ( K − M Y ) Y (cid:3) . (65)Using the Leray-Schauder fixed point theorem in R N , we prove the existence of a solution to the fixed point problem(65). Theorem 3.3. (Leray-Schauder Theorem [10] ) Let T be a continuous and compact mapping of a Banach space X into itself, such that the set { x ∈ X : x = λ T x for some ≤ λ ≤ } is bounded. Then T has a fixed point. Theorem 3.4.
Let ≤ s ≤ and consider the fixed point problem BY s = s (cid:2) M Z − τ S ( K − M Y s ) Y s (cid:3) (66) then for τ ≤ || ˆ k || ∞ , the fixed point problem Y = G ( Y ) admits a solution. roof. Multiplying equation (66) by Y Ts we get: Y Ts BY s = sY Ts M Z − τ sY Ts S ( K − M Y s ) Y s Using skew-symmetry of the matrix S ( K − M Y s ) we get: Y Ts BY s = sY Ts M Z ≤ s || Y s || M || Z || M (67) Y Ts ( M − τ RK − M ) Y s = || Y s || M − τ Y Ts RK − M Y s ≤ s || Y s || M || Z || M (68) || Y s || M ≤ s || Y s || M || Z || M + τ Y Ts RK − M Y s = s || Y s || M || Z || M + τ Y Ts RW s (69)where KW s = M Y s . Let w s,N = N (cid:80) i =1 W s,i ϕ i and y s,N = N (cid:80) i =1 Y s,i ϕ i . Then, || w s,N || = || W s || M and || y s,N || = || Y s || M . Also, from KW s = M Y s one gets: < w s,N , v > = < y s,N , v > = ⇒ || w s,N || ≤ || y s,N || (70)and Y Ts RW s = < ˆ k ( w s,N ) y , y s,N > . Hence, using the last two equations, (69) becomes || Y s || M ≤ s || Y s || M || Z || M + τ Y Ts RW s ≤ || Y s || M || Z || M + τ || ˆ k || ∞ || y s,N || (71) ≤ || Y s || M || Z || M + τ || ˆ k || ∞ || Y s || M (72) = ⇒ (1 − τ || ˆ k || ∞ ) || Y s || M ≤ || Y s || M || Z || M (73)Then for τ ≤ || ˆ k || ∞ we get || Y s || M ≤ || Z || M . Thus all solutions Y s are uniformly bounded, leading to theexistence of a fixed point for all ≤ s ≤ , using the Leray-Schauder fixed point theorem. Now we turn to uniqueness and consider the ball B N = { Y ∈ R N | || Y || M ≤ C Z = 2 || Z || M } (74)we then prove the following theorem. Theorem 3.5.
For Y ( k ) ∈ B N , k = 1 , and τ ≤
18 min (cid:40) || ˆ k || ∞ , h cC Z (cid:41) , one has: || G ( Y (1) ) − G ( Y (2) ) || M ≤ || Y (1) − Y (2) || and hence under such restriction on τ and h , the fixed problem Y = G ( Y ) admits a unique solution.Proof. For k = 1 , , let Λ ( k ) = G ( Y ( k ) ) (75) K Σ ( k ) = M Y ( k ) (76) K ∆ ( k ) = M Λ ( k ) (77)Then by applying this equation for k = 1 and , followed by a subtraction, we get ( M − τ RK − M )(Λ (1) − Λ (2) ) = τ S ( K − M Y (2) ) Y (2) − τ S ( K − M Y (1) ) Y (1) = τ S (Σ ) Y (2) − τ S (Σ ) Y (1) ( K − τ R )(∆ (1) − ∆ (2) ) = τ S (Σ (2) ) Y (2) − τ S (Σ (1) ) Y (1) (78)Let Y ( k ) N = N (cid:80) i =1 Y ( k ) i ϕ i and (cid:15) N = Y (1) N − Y (2) N , then || (cid:15) N || = || Y (1) − Y (2) || M .Let Λ ( k ) N = N (cid:80) i =1 Λ ( k ) ϕ i and λ N = Λ (1) N − Λ (2) N , then || λ N || = || G ( Y (1) ) − G ( Y (2) ) || M .In variational form, equations , (76), (77) and (78) lead to: < ∆ ( k ) N , v > − τ < ˆ k (∆ ( k ) N ) y , v > = < Z N , v > − τ < (cid:126)V (Σ ( k ) N ) . ∇ Y ( k ) N , v > (79) < ∆ (1) N − ∆ (2) N , v > − τ < ˆ k (∆ (1) N − ∆ (2) N ) y , v > = − τ < (cid:126)V (Σ (1) N ) . ∇ Y (1) N − (cid:126)V (Σ (2) N ) . ∇ Y (2) N , v > (80) < Σ ( k ) N , v > = < Y ( k ) N , v > (81) < ∆ ( k ) N , v > = < Λ ( k ) N , v > (82)For all v ∈ X N,P , where Z N = N (cid:80) i =1 Z i ϕ i , Σ ( k ) N = N (cid:80) i =1 Σ ( k ) i ϕ i , ∆ ( k ) N = N (cid:80) i =1 ∆ ( k ) i ϕ i , and ∆ ( k ) N,y = N (cid:80) i =1 ∆ ( k ) i,y ϕ i , for k = 1 , .Let us define δ N = ∆ (1) N − ∆ (2) N , σ N = Σ (1) N − Σ (2) N , and . Using again the bilinear form a τ ( ., . ) defined in theorem3.1, a τ ( δ N , v ) = − τ < (cid:126)V ( σ N ) . ∇ Y (1) N − (cid:126)V (Σ (2) N ) . ∇ (cid:15) N , v > (83)Then, (80), (81), and (82) lead to: a τ ( δ N , v ) = − τ < (cid:126)V ( σ N ) . ∇ Y (1) N − (cid:126)V (Σ (2) N ) . ∇ (cid:15) N , v > (84) < σ N , v > = < (cid:15) N , v > (85) < δ N , v > = < λ N , v > (86)for all v ∈ X N,P .Now, let v = δ N in (84), then for τ ≤ || ˆ k || ∞ , we have: || δ N || ≤ τ | < (cid:126)V ( σ N ) . ∇ Y (1) N , δ N > | + τ | < (cid:126)V (Σ (2) N ) . ∇ (cid:15) N , δ N > | (87)Since, theorem (2.1) asserts that ∀{ v, z, φ } ∈ X N,P × X N,P × X N,P : (cid:68) (cid:126)V ( v ) · ∇ z, φ (cid:69) L = − (cid:68) (cid:126)V ( v ) · ∇ φ, z (cid:69) L , then (87) can be rewritten as: || δ N || ≤ τ | < (cid:126)V ( σ N ) . ∇ δ N , Y (1) N > | + τ | < (cid:126)V (Σ (2) N ) . ∇ δ N , (cid:15) N > | (88) ≤ τ || (cid:126)V ( σ N ) . ∇ δ N || || Y (1) N || + τ || (cid:126)V (Σ (2) N ) . ∇ δ N || || (cid:15) N || (89)Using Lemma 2.7, this last inequality leads to: || δ N || ≤ τh || δ N || (cid:104) C Z || σ N || + || Σ (2) N || . || (cid:15) N || (cid:105) . Therefore, given that || δ N || ≤ || δ N || , it results that: || δ N || ≤ τh (cid:104) C Z || σ N || + || Σ (2) N || . || (cid:15) N || (cid:105) . (90)Note that if in (85), v = σ N , one proves that: || σ N || ≤ || (cid:15) N || , k = 2 and v = Σ (2) N , we obtain: || Σ (2) N || ≤ || Y (2) N || ≤ C Z . Combining the last two inequalities with inequality (90), we conclude the inequality: || δ N || ≤ τ C Z h || (cid:15) N || . (91)Finally, if in (86), we let v = λ N , then: || λ N || = | < δ N , λ N > | ≤ || δ N || || λ N || Using a result from Ciarlet ([8], Theorem 3.2.6), one has: || λ N || ≤ ch || λ N || , then || λ N || ≤ ch || δ N || (92)When combining (91) and (92) we reach the final result: || G ( Y ) − G ( Y ) || M = || Λ (1) N − Λ (2) N || = || λ N || ≤ ch || δ N || ≤ τ cC Z h || (cid:15) N || = 4 τ cC Z h || Y − Y || M . On the other hand, considering (79), one writes: a τ (∆ ( k ) N , v ) = < Z N , v > − τ < (cid:126)V (Σ ( k ) N ) . ∇ Y ( k ) N , v >, k = 1 , Consequently, our Theorem is proved
Corollary 3.6.
Under the conditions of theorem 3.5, namely τ ≤
18 min (cid:40) || ˆ k || ∞ , h cC Z (cid:41) , the iteration Y ( k +1) = G ( Y ( k ) ) with Y (0) = Z in the ball B N converges to the unique solution of Y = G ( Y ) .Proof. Starting with || Y (1) − Y || M = || G ( Z ) − G ( Y ) || M ≤ || Z − Y || then by induction we get that || Y ( k +1) − Y || M = || G ( Y ( k ) ) − G ( Y ) || M ≤ || Y ( k ) − Y || ≤ k +1 || Y (0) − Y || where ≤ lim k →∞ || Y ( k +1) − Y || M ≤ lim k →∞ k +1 || Y (0) − Y || = 0 . Recall from (59) that in matrix notations and using the expressions: w N ( t ) = N (cid:88) I =1 W I ( t ) ϕ I ( x, y ) , and u N ( x, y, t ) = N (cid:88) J =1 U J ( t ) ϕ J ( x, y ) , where W I ( t ) = w N ( x I , y I , t ) , and U J ( t ) = w N ( x J , y J , t ) , then (58) can be rewritten as follows:Given ( U ( t ) , W ( t )) ∈ R N × R N , seek ( U ( t + τ ) , W ( t + τ )) ∈ R N × R N , such that: (cid:26) ( M + τ S ( U ( t + τ )) W ( t + τ ) − τ R U ( t + τ ) = M W ( t ) (1) KU ( s ) = M W ( s ) , ∀ s ∈ { t, t + τ } (2)
18o solve (59) we use in this paper a semi-linearized approach: (cid:40) (cid:104) w N ( t + τ ) , v (cid:105) L − τ (cid:68) (cid:126)V ( u N ( t )) · ∇ v, w N ( t + τ ) (cid:69) L = (cid:104) w N ( t ) , v (cid:105) L + τ (cid:68) ˆ ku N,y ( t ) , v (cid:69) L , (1) (cid:104) u N ( s ) , v (cid:105) H = (cid:104) w N ( s ) , v (cid:105) L , s ∈ { t, t + τ } (2) (93) ∀ v ∈ X N,P , which in matrix form is given by: (cid:26) ( M + τ S ( U ( t )) W ( t + τ ) = M W ( t ) + τ R U ( t ) (1) KU ( t + τ ) = M W ( t + τ ) , (2) (94)where M, K , S ( U ) and R are N × N matrices defined in (13). However, by taking periodicity into account, thedegrees of freedom are reduced from N = n to N = ( n − . Note that M , is the well-known Mass matrix forperiodic boundary conditions and K = M + A where A is the stiffness matrices for periodic boundary conditions.The nonlinearity of the problem originates from S ( U ) , which we derive its corresponding local matrix, in addition tothat of R , over each triangle for equally spaced nodes in Section 4.1. Then, deduce the block sparsity pattern of theglobal matrix, which is the same for M and K .Thus, to solve (94) these matrices should be generated for a given meshing. In Section 4.2, we implement Algorithm(1) using Freefem++ [11], a programming language and software focused on solving partial differential equationsusing the finite element method.Using different initial conditions, we test in Section 4.3 our algorithm for the case when p = ln n ω ci is a function of x ,such that ˆ k = p x is a constant and p y = 0 , and for the case when p, p x , and p y are functions of ( x, y ) . S ( U ) and R To compute the matrices S ( U ) and R , the square domain Ω is partitioned into n equally-spaced nodes in each of the x and y direction, leading to a set of n nodes. N = { P I ( x i , y j ) | I = 1 , , ..., N = n } = P x × P y The indexing of these nodes starts from left to right, and bottom to top as shown in Figure 2 for n = 5 . Moreover, theset of M = 2( n − triangles covering Ω are also indexed from left to right, and bottom to top T = { T J | J = 1 , , ..., M } , Ω = ∪ J T J The global indexing of the vertices of triangles T j − of type a are { j + c, j + n + 1 + c, j + n + c } , whereas that oftriangles T j of type b are { j + c, j + 1 + c, j + n + 1 + c } , for j = 1 , , .., M n − and c = (cid:24) jn − (cid:25) − .Figure 2: A two-dimension meshing of Ω with the corresponding nodes and triangle indexing S ( U ) and R are first computed locally on each triangle T J and then assembled globally. Triangle T J has three nodes with global indexing { α, β, γ } depending on its type, and local indexing { , , } . On triangle T J , theonly non-zero basis functions are ψ α , ψ β , and ψ γ are locally denoted by ψ , ψ , and ψ . Thus, the local S ( U ) and R matrices on triangle T J have at most nonzero entries (in rows and columns α, β, γ ) that can be computed using a × matrix, denoted by S ( J ) ( U ( J ) ) and R ( J ) , where U ( J ) = [ U α ( t ) , U β ( t ) , U γ ( t )] is a vector of length .By assembling the global matrices and imposing periodic boundary conditions, the degrees of freedom are reducedfrom N = n to N = ( n − . Letting { i , i , · · · , i N } ⊂ { , , · · · , N } , one defines the extracted vectors ˜ z ∈ R N from z ∈ R N , and extracted matrices ˜ E ∈ R N × N from E ∈ R N × N as shown in the appendices A.2 and B.2 for thematrices S ( U ) and R . Note that in all following sections we drop the tilde notation, and the matrices M, K , S ( U ) and R are assumed to be of size N × N , and the vectors U ( t ) ≡ U, W ( t ) ≡ W are of size N . Moreover, based on thisextraction of the minimum number of degrees of freedom, we consider w N ( x, y, t ) = N (cid:88) I =1 W I ( t ) ψ I ( x, y ) , and u N ( x, y, t ) = N (cid:88) J =1 U J ( t ) ψ J ( x, y ) where { ψ J | J = i , ..., i N } is the modified basis extracted from { ϕ I | I = 1 , ..., N } .In what follows, we compute the local matrices S ( J ) ( U ( J ) ) and R ( J ) and state the sparsity patterns of the resultingglobal matrices S ( U ) and R . Local Matrix S ( J ) ( U ( J ) ) The 9 entries S ( J ) i,j of the local matrix S ( J ) ( U ( J ) ) are defined as follows for i, j = 1 , , S ( J ) ( U ( J ) ) = S ( J )1 , S ( J )1 , S ( J )1 , S ( J )2 , S ( J )2 , S ( J )2 , S ( J )3 , S ( J )3 , S ( J )3 , S ( J ) i,j = (cid:82) T J (cid:126)V ( u n ) · ∇ ψ i ψ j dA = − (cid:82) T J ∂ψ i ∂x ψ j (cid:80) k =1 (cid:18) U ( J ) k ∂ψ k ∂y (cid:19) dA + (cid:82) T J ∂ψ i ∂y ψ j (cid:80) k =1 (cid:18) U ( J ) k ∂ψ k ∂x (cid:19) dA where ψ i ( x, y ) = a i + b i x + c i y with a = x y − x y Area ( T J ) , b = y − y Area ( T J ) , c = x − x Area ( T J ) a = x y − x y Area ( T J ) , b = y − y Area ( T J ) , c = x − x Area ( T J ) a = x y − x y Area ( T J ) , b = y − y Area ( T J ) , c = x − x Area ( T J ) and Area ( T J ) = 0 . x ( y − y ) + 0 . x ( y − y ) + 0 . x ( y − y ) . Moreover, ∂ψ i ∂y = c i , and ∂ψ i ∂x = b i . Thus, S ( J ) i,j = − (cid:90) T J b i ψ j (cid:88) k =1 (cid:16) c k U ( J ) k (cid:17) dA + (cid:90) T J c i ψ j (cid:88) k =1 (cid:16) b k U ( J ) k (cid:17) dA (95) = − b i (cid:32) (cid:88) k =1 c k U ( J ) k (cid:33) (cid:90) T J ψ j dA + c i (cid:32) (cid:88) k =1 b k U ( J ) k (cid:33) (cid:90) T J ψ j dA (96)Let b ( J ) = [ b , b , b ] and c ( J ) = [ c , c , c ] , then (cid:80) k =1 b k U ( J ) k = b ( J ) · U ( J ) and (cid:80) k =1 c k U ( J ) k = c ( J ) · U ( J ) areconstants per triangle. Let η i = (cid:82) T J ψ i dA , and η ( J ) = [ η , η , η ] then,20 ( J ) ( U ( J ) ) = − ( c ( J ) · U ( J ) ) b η b η b η b η b η b η b η b η b η + ( b ( J ) · U ( J ) ) c η c η v η c η c η c η c η c η c η . (97)Note that η i = (cid:82) T J ψ i dA = 13 Area ( T J ) . Thus, b η i = 16 ( y − y ) , c η i = 16 ( x − x ) , and S ( J ) ( U ( J ) ) = −
16 ( c ( J ) · U ( J ) ) y − y y − y y − y (cid:2) (cid:3) + 16 ( b ( J ) · U ( J ) ) x − x x − x x − x (cid:2) (cid:3) (98) = d J (cid:16)(cid:98) c ( J ) · U ( J ) (cid:98) b ( J ) (cid:2) (cid:3) − (cid:98) b ( J ) · U ( J ) (cid:98) c ( J ) (cid:2) (cid:3)(cid:17) (99)where (cid:98) b ( j ) = y − y y − y y − y , (cid:98) c ( j ) = x − x x − x x − x , and d J = 112 Area ( T J ) After computing S ( J ) ( U ( J ) ) its rows and columns are mapped from local indexing { , , } to the global indexing { α, β, γ } and added to the global matrix S ( U ) . Thus, for computing S ( U ) two difference matrices B, C of size M × have to be computed once and stored, where the J th row of B is (cid:98) b ( J ) and the J th row of C is (cid:98) c ( J ) ; in addition to an M × vector of triangle areas. Note that the matrix S ( U ) has to be assembled at each time iteration.Assuming that the set of nodes on Ω are equally spaced, i.e. x i +1 − x i = y i +1 − y i = h, ∀ i = 0 , , .., n − , then S ( J ) ( U ( J ) ) can be further simplified. In this case, there are 2 types of triangles with the local nodes numbering asshown in Figure 3.Figure 3: The 2 types of triangles with their local nodes’ indexing, assuming that the set of nodes on Ω are equally spaced in thex and y directions. The triangle on the left is denoted by type a, whereas that on the right by type b. For triangles of type a, (cid:98) b ( J ) = h (cid:2) − (cid:3) T and (cid:98) c ( J ) = h (cid:2) − (cid:3) T . Whereas, for triangles of type b, (cid:98) b ( J ) = h (cid:2) − (cid:3) T and (cid:98) c ( J ) = h (cid:2) − (cid:3) T . In both cases, Area ( T J ) = h and S ( J ) ( U ( J ) ) = 16 U ( J )3 − U ( J )2 U ( J )1 − U ( J )3 U ( J )2 − U ( J )1 (cid:2) (cid:3) (100)Thus, the computation of S ( U ) in the case of equally spaced nodes reduces to taking differences of the U vectorentries, without the need to store any values. In addition, S ( U ) is a block tridiagonal matrix with 2 additional blocksin the upper right and lower left corner. Moreover, it is a skew-symmetric matrix ( S ( U ) T = − S ( U ) ) that is linear in U , with 6 nonzero entries per row, 6 nonzero entries per column, and zeros on the diagonal assuming the meshing of21 shown in Figure 2. S ( U ) = 16 S , S , · · · S ,k S , S , S , · · · . . . . . . . . . . . . ...... . . . S j,l S j,j S j,i · · · S i,j S i,i S i,k A k, · · · S k,i S k,k with S i,j ≡ S i,j ( U ) where i = n − , j = n − , k = n − , l = n − , and the n − nonzero block matrices S i,j are of size ( n − × ( n − with n − nonzero entries each, and the following sparsity patterns:• S i,i for i = 1 , .., n − are tridiagonal matrices with zero diagonal entries, and nonzero S i,i (1 , n − , and S i,i ( n − , .• S ,n − and S i +1 ,i for i = 1 , , , .., n − are lower bidiagonal matrices, with nonzero entry in first row andcolumn n − .• S n − , and S i,i +1 for i = 1 , , .., n − are upper bidiagonal matrices with nonzero entry in first column androw n − .Thus, S ( U ) has a total of n − n −
1) = 6 N nonzero entries. As for the explicit expressions/values of theentries, refer to appendix (A.2). Local Matrix R ( J ) The 9 entries R ( J ) i,j of the local matrix R ( J ) are defined as follows for i, j = 1 , , R ( J ) = R ( J )1 , R ( J )1 , R ( J )1 , R ( J )2 , R ( J )2 , R ( J )2 , R ( J )3 , R ( J )3 , R ( J )3 , R ( J ) i,j = (cid:82) T J ˆ k ψ i,y ψ j dA = (cid:82) T J ˆ k ∂ψ i ∂y ψ j dA where ψ i ( x, y ) = a i + b i x + c i y, and ∂ψ i ∂y = c i with c = x − x Area ( T J ) c = x − x Area ( T J ) c = x − x Area ( T J ) and Area ( T J ) = 0 . x ( y − y ) + 0 . x ( y − y ) + 0 . x ( y − y ) . Thus, assuming ˆ k is constant, then R ( J ) i,j = ˆ k c i (cid:90) T J ψ j dA = ˆ k c i Area ( T J ) (101) R ( J ) = 16 ˆ k x − x x − x x − x x − x x − x x − x x − x x − x x − x = 16 ˆ k x − x x − x x − x (cid:2) (cid:3) = 16 ˆ k (cid:98) c ( J ) (cid:2) (cid:3) (102)After computing R ( J ) , its rows and columns are mapped from local indexing { , , } to the global indexing { α, β, γ } and added to the global matrix R . Thus, for computing R one difference matrix C of size M × has to be computedonce and stored, where the J th row of C is (cid:98) c ( J ) . Note that the matrix R is computed once.Assuming that the set of nodes on Ω are equally spaced, i.e. x i +1 − x i = h, ∀ i = 0 , , .., n − , then R ( J ) can befurther simplified. 22or triangles of type a shown in Figure 3, (cid:98) c ( J ) = h (cid:2) − (cid:3) T and R ( J ) a = h k − − −
10 0 01 1 1 .Whereas, for triangles of type b, (cid:98) c ( J ) = h (cid:2) − (cid:3) T and R ( J ) b = h k − − −
11 1 1 .Given that the matrix R is independent of U , its computation is straightforward. In addition, R has the same blocksparsity pattern as S ( U ) . For example, assuming that ˆ k is constant, then R is a skew-symmetric matrix ( R T = − R )with zeros on the diagonal, and 6 nonzero entries per row and 6 nonzero entries per column, of the form h kα where α = − , − , , or , with the sum of entries per row or column being zero. (refer to appendix (B.2)). R = h k R , R , · · · R ,n − R , R , R , · · · . . . . . . . . . . . . ...... . . . R j,m R j,j R j,i · · · R i,j R i,i R i,l A l, · · · R l,i R l,l where i = n − , j = n − , l = n − , m = n − , and the n − nonzero block matrices R i,j are of size ( n − × ( n − with n − nonzero entries each, and the following sparsity patterns:• R i,i for i = 1 , , .., n − are such that: R i,i ( j, j + 1) = 1 , R i,i ( j + 1 , j ) = − , for j = 1 , , ..., n − , R i,i (1 , n −
1) = − , and R i,i ( n − ,
1) = 1 .• R ,n − = R i +1 ,i for i = 1 , , .., n − are lower bidiagonal matrices ( R i +1 ,i ( j, j ) = 2 , R i +1 ,i ( j + 1 , j ) = 1 ),with R i +1 ,i (1 , n −
1) = 1 .• R n − , = R i,i +1 for i = 1 , .., n − are upper bidiagonal matrices ( R i,i +1 ( j, j ) = − , R i,i +1 ( j, j + 1) = − ),with R i,i +1 ( n − ,
1) = − .Thus, R has a total of N nonzero entries. (94) The semi linear scheme (94) can be solved at each time step as shown in algorithm (1), where the matrices need to begenerated as described previously, based on the meshing of the space domain given in Figure 2. For that purpose, weimplemented Algorithm 1 using FreeFem++.
Algorithm 1
Numerical Hasegawa-Mima semi-linearized Finite Element Scheme
Input: M : mass matrix ; K = M + A , A : stiffness matrix; S ( U ) : algorithm that builds S ( U ) as in Appendix A; R : matrix defined in Section 4.1; U , W : the discrete initial condition vectors; τ : time step; T : end time; U = U ; W = W ; for t = 0 : τ : T − τ do G = M ∗ W ; Solve for W : ( M + τ S ( U )) W = τ R U + G ; Solve for U : K U = M W ; end for
23e consider a square domain [ x , x n ] × [ y , y n ] with a uniform mesh in the x and y direction ( x i − x i − = x n − x n = y i − y i − = y n − y n for i = 1 , , .., n and n intervals in the x and y directions respectively) and the finite element P space with periodic boundary conditions, using appropriate Freefem++ functions. The function p of the initialHasegawa-Mima PDE is defined, where in most simulations it is assumed that p y = 0 , and ˆ k = p x is a constant,unless stated otherwise. The initial conditions u is given as input. As for the initial condition w = u − ∆ u itcould be given as input if u is a simple function. However, for any function u , we compute the vector W = W (0) by solving the linear system M ∗ W = K ∗ U where the vectors U = U (0) , W , the mass matrix M , and the matrix K = M + A are defined in (13), with A beingthe stiffness matrix. Note that the matrices M , A , R and S ( U j ) are generated in Freefem++ using the correspondingvariational formulations: a ( u, v ) = (cid:90) T h ( u x ∗ v x + u y ∗ v y ); where the matrix A = a ( V h, V h ); b ( u, v ) = (cid:90) T h ( u ∗ v ); where the matrix M = b ( V h, V h ); c ( u, v ) = (cid:90) T h ( p x ∗ u y ∗ v − p y ∗ u x ∗ v ); where the matrix R = c ( V h, V h ); d ( w, v ) = (cid:90) T h ( u jx ∗ w y ∗ v − u jy ∗ w x ∗ v ) where the matrix S ( ˜ U j ) = d ( V h, V h ) . (103)Another alternative in Freefem++ is to just define the variational formulation of (94) as problems that are solved ateach time iteration: hypo ( w j +1 , v ) = (cid:90) T h ( w j +1 ∗ v/τ − w j ∗ v/τ ) + (cid:90) T h ( u jx ∗ w j +1 y ∗ v − u jy ∗ w j +1 x ∗ v ) + (cid:90) T h ( p y ∗ u jx ∗ v − p x ∗ u jy ∗ v ) ellip ( u j +1 , v ) = (cid:90) T h ( u j +1 x ∗ v x + u j +1 y ∗ v y ) + (cid:90) T h ( u j +1 ∗ v ) − (cid:90) T h ( w j +1 ∗ v ) where u j +1 , w j +1 refer to the sought solutions at the current time iteration, i.e. u ( t + τ ) , w ( t + τ ) , and u j , w j refer to u ( t ) , w ( t ) the solutions at the previous time iteration. At each time iteration, Freefem++ will generate thecorresponding matrices and vectors and solve the linear systems accordingly. However, this implies that the fixedmatrices M, K, and R will be regenerated redundantly at each iteration. Thus, it is preferable timewise to generatethe fixed matrices once in the algorithm and solve the matrix form of the problem (Algorithm refalg:HMC-Newton).Table 1 validates this claim, where the execution time is reduced at least by half.Note that the simulation is stopped once the maximum value of u ( t ) at one of the mesh nodes is . , which correspondsto the maximum value attained physically. τ T FreeFem++ FreeFem++Variational Form Matrix Form1/8 34.2500 51.8857 25.61411/10 40.7000 78.3870 35.82611/16 61.3125 188.9680 85.79791/32 119.0000 676.0540 318.10101/64 235.0000 3245.3400 1272.4000Table 1:
Execution time of the semi linear scheme (94) in FreeFem++ using the Variational Form or the Matrix Form for n = 64 , u = 10 − sin (3 x ) , Ω = [0 , π ] × [0 , π ] , ˆ k = p x = 12 , p y = 0 , and a given time step τ . T denotes the end time at which thesimulation was stopped as u ( t ) reached the maximum value of . at one of the mesh nodes. .3 Testing We start by testing Algorithm 1 for n ω ci = e Ax + B , i.e. p = Ax + B where ˆ k = p x = A and p y = 0 . As explainedin [7], the solution is expected to be a traveling wave in the y-direction for a nonzero A . The speed of the motion andits direction depend on the magnitude and sign of A respectively. We consider different initial conditions for the sameexponential density profile n . We consider the following cases:1. Domain Ω = [0 , × [0 , with the number of intervals in the x and y direction n = 64 ( h = 1 / ≈ . ), A = 12 , the initial condition u ( x, y ) = 10 − sin (10 πy ) and τ = 0 . . Figure 4 shows the time evolution of thesolution of the Semi-Linear Scheme (94) at time t = 0 , , , , . It is clear that the sin function is movingin the y-direction as time proceeds without any perturbation in its initial form up till t = 200 . Afterwards, thesolution grows with time to reach || u || ∞ = 0 . at t = 260 . when the algorithm is stopped.Note that even though τ > || ˆ k || ∞ = 124 does not satisfy the sufficient condition for proving the existence anduniqueness of a solution to (59), however the algorithm still converges and produces the expected behavior.Figure 4: Time evolution of solution u for u = 10 − sin (10 πy ) , τ = 0 . , and a × grid on Ω = [0 , × [0 , .
25. Domain
Ω = [0 , π ] × [0 , π ] with A = 12 , the initial condition u ( x, y ) = 10 − sin (3 y ) , and τ = 0 . .In Figure 5 we consider the number of intervals in the x and y direction n = 32 ( h = π/ ≈ . ), and inFigure 6 n = 64 ( h = π/ ≈ . ). We notice the same behavior where the solution moves in y-directionand grows with time at a faster rate to reach || u || ∞ = 0 . at t = 9 . . The larger h value ( h ≈ . , n = 32 ) doesnot affect the solution with respect to that of h ≈ . , n = 64 , apart from the smoothness of the 3D surface.From this perspective, this shows the robustness of the algorithm for reasonable h < values.Figure 5: Time evolution of solution u for u = 10 − sin (3 y ) , τ = 0 . , and a × grid on Ω = [0 , π ] × [0 , π ] . Figure 6:
Time evolution of solution u for u = 10 − sin (3 y ) , τ = 0 . , and a × grid on Ω = [0 , π ] × [0 , π ] .
26. To test the fact that the solution will always converge to a sine function moving in the y direction when n ω ci = e Ax + B , i.e. p = Ax + B , we start with u ( x, y ) = 10 − sin (3 x ) for τ = 0 . , n = 32 , with A = 12 . Thesolution remains unchanged up till t = 26 , and after the transitional time where the solution shifts from a sinefunction in the x direction to a sine function in the y-direction, the same behavior is observed (Figure 7).Figure 7: Time evolution of solution u for u = 10 − sin (3 x ) , τ = 0 . , and a × grid on Ω = [0 , π ] × [0 , π ] .
27. It should be noted that even if we start with a random initial condition u ( x, y ) , after some time the solutionconverges to the same sine function moving in the y direction. Figure 8 shows the time evolution of the solutionfor u ( x, y ) = 10 − xy ( x − sin ( x ) , τ = 0 . , Ω = [0 , π ] × [0 , π ] , n = 32 ( h = π/ ≈ . ) with A = 12 .Figure 8: Time evolution of solution u for u = 10 − xy ( x − sin ( x ) , τ = 0 . , and a × grid on Ω = [0 , π ] × [0 , π ] . p y = 0 , however the algorithm works for any inputfunction p . Note that if we set p x = 0 and p y = 12 , then the solution will be moving in the x-direction in a similar man-ner as shown in the previous examples. We test the algorithm for the case where n = 10 e − ( x − / − ( y − / , w ci = 10 , and Ω = [0 , × [0 , . Since ∇ p = [ − ( x − / , − ( y − / , the solution is expected to have acircular motion around the center of the domain (10 , as shown in Figure 9.Figure 9: Time evolution of solution u for u = − − ( x − e − . x − − . y − , τ = 0 . , and a × grid on Ω = [0 , × [0 , . ∇ p = [( x − / , ( y − / for the same initial conditions, then the solution will be moving in theopposite direction. To sum up the results of this paper, we have proven an existence theorem for the Hasegawa-Mima wave equation in C (0 , T ; H ∩ H P ) × L ∞ (0 , T ; L ) . Uniqueness of the solution would require more regularity on the initial condi-tions as proven in [1]. On the other hand, we have also considered a full discretization scheme based on coupling theFinite-Element in space and a non-linear discretization of the time variable. The implementation of that scheme usesa semi-linear approach that provides a robust algorithm as revealed by early experiments.Overall, future avenues of research include the following:1. Proof of convergence of the solution to the nonlinear (58), (59) schemes and the semi-linear (93) one as τ and h go to zero.2. Testing other alternatives to solve (59). These include:(a) Using the fixed point approach discussed in Section 3.1, where the system can be written as follows:Since KU j = M W j , and if we let Z = M W ( t ) , Y = W ( t + τ ) , and Y = W ( t ) , then Y k +1 = G ( Y k ) ⇐⇒ ( M − τ RK − M ) Y k +1 = Z − τ S ( U k ) Y k which leads to the following predictor-corrector scheme to move from time t to t + τ U = U ( t ) , Y = W ( t ) , Z = M W ( t ) error = 1 , k = 1 While error > tol do: ( M − τ R K − M ) Y k +1 = Z − τ S ( U k ) Y k KU k +1 = M Y k +1 error = max || W k +1 − W k |||| W k +1 || , || U k +1 − U k |||| U k +1 || k = k + 1 end do U ( t + τ ) = U k and W ( t + τ ) = Y k (104)Following the results obtained in Section 3.1 , specifically corollary (3.6), this algorithm is convergent. Afirst look at this approach indicates the necessity to deal with dense linear systems, which matrix is ( M − τ R K − M ) . However, this difficulty can be lifted since ( M − τ R K − M ) Y = r is equivalentto ( K − τ R ) K − M Y = r where the solution is Y = M − K ( K − τ R ) − r can be obtained by solvingtwo time-independent sparse systems ( K − τ R ) ˜ Y = r followed by M Y = K ˜ Y .(b) A second approach to handle (59) would be based on Newton’s method.3. Another interesting problemhas to deal with the Modon Traveling Waves Solutions to (6). These solutions areobtained by considering the pair of variables ( ξ, η ) given by ξ = x η = y − ct , one looks for solutions to(4) in the form u ( x, y, t ) = φ ( ξ, η ) = φ ( x, y − ct ) and w ( x, y, t ) = ψ ( ξ, η ) = ψ ( x, y − ct ) . By defining ∀ t ∈ (0 , T ) : Ω t = { ξ, η | < ξ < L, − ct < η < L − ct } , then in terms of φ and ψ , the system (4) reduces tobe solved on Ω = Ω . Thus, with ∇ = ∇ ξ,η , one seeks { φ, ψ } : Ω → R , such that: − cψ η + (cid:126)V ( φ ) · ∇ ψ = kφ η on Ω (1) − ∆ φ + φ = ψ on Ω (2)
PBC’s on φ, φ ξ , φ η , ψ on ∂ Ω (3) (105)Undergoing research is also being carried out on this problem.30 eferences [1] H. Karakazian and N. Nassif, “Local existence of an h p solution to the hasegawa-mima plasma equation,” Sub-mitted. arXiv:1712.05524 , 2019.[2] ITER Physics Expert Groups on Confinement and Transport and Confinement Modeling and Database, “Plasmaconfinement and transport,”
Nuclear Fusion , vol. 39, no. 12, pp. 2175–2249, 1999.[3] A. Hasegawa and K. Mima, “Stationary spectrum of strong turbulence in magnetized nonuniform plasma,”
Physics of Fluids , vol. 39, pp. 205–208, Jul 1977.[4] A. Hasegawa and K. Mima, “Pseudo-three-dimensional turbulence in magnetized nonuniform plasma,”
Physicsof Fluids , vol. 21, pp. 87–92, Jan 1978.[5] A. Hasegawa and M. Wakatani, “Plasma edge turbulence,”
Phys. Rev. Lett. , vol. 50, pp. 682–686, Feb 1983.[6] B. Shivamoggi, “Charney-hasegawa-mima equation: A general class of exact solutions,”
Physics Letters A ,vol. 138, pp. 37–42, Jun 1989.[7] F. A. Hariri, “Simulating bi-dimensional turbulence in fusion plasma with linear geometry,” Master’s thesis,American University of Beirut, 2010.[8] P. G. Ciarlet,
The Finite Element Method for Elliptic Problems . SIAM, 1979.[9] H. Brezis,
Functional Analysis, Sobolev Spaces and Partial Differential Equations . Springer, 2011.[10] J. Mawhin, “Leray-schauder degree: a half century of extensions and applications,”
Topol. Methods NonlinearAnal. , vol. 14, no. 2, pp. 195–228, 1999.[11] F. Hecht, “New development in freefem++,”
J. Numer. Math. , vol. 20, no. 3-4, pp. 251–265, 2012.31 ppendices
A Assembling the Global matrix S ( U ) First, we derive the sparsity pattern of the global matrix S ( U ) without any assumptions related to boundary conditions,i.e. define the general matrix S ( U ) , in Section A.1. Then, we derive the sparsity pattern of the matrix S ( U ) assumingperiodic boundary conditions in Section A.2. A.1 The General Matrix S ( U ) Given the meshing of Figure 2 and without any assumptions on the boundary nodes, S ( U ) is an n × n sparse matrixwith at most 6 nonzero entries per row, as discussed below. To get those nonzero entries, note that each node has atmost 6 edges connecting it with its neighboring nodes, i.e. belongs to at most 6 triangles. We consider the 4 types ofnodes shown in different colors in Figure 2: the ( n − black internal nodes that belong to 6 triangles, the n − red boundary nodes that belong to 4 triangles, the 2 green corner nodes that belong to 2 triangles, and the 2 blue cornernodes that belong to 1 triangle.
1- Black Internal Nodes:
Each of the ( n − black internal nodes with index v = kn + i for k = 1 , , ..., n − and i = 2 , , ..., n − belongs to triangles T j +1 , T j +2 , T j +3 , T j + n ) , T j + n )+1 , T j + n )+2 , where j = ( n − k −
1) + ( i − .Thus, we first define the nonzero entries in row v per local triangle and then add them up.• Triangle T j +1 = T j +1) − has global vertices ( j + 1) + ( k −
1) = nk − n + i − j + 1) + ( n + 1) + ( k −
1) = nk + i = v ( j + 1) + n + ( k −
1) = nk + i − where node v is the second local vertex of triangle T j +1 . Thus row v of the matrix S ( U ) has the entry U nk − n + i − − U nk + i − in columns nk − n + i − , nk + i − , and v = nk + i .• Triangle T j +2 = T j +1) has global vertices ( j + 1) + ( k −
1) = nk − n + i − j + 1) + 1 + ( k −
1) = nk − n + i ( j + 1) + ( n + 1) + ( k −
1) = nk + i = v where node v is the third local vertex of triangle T j +2 . Thus row v of the matrix S ( U ) has the entry U nk − n + i − U nk − n + i − in columns nk − n + i − , nk − n + i , and v = nk + i .• Triangle T j +3 = T j +2) − has global vertices ( j + 2) + ( k −
1) = nk − n + i ( j + 2) + ( n + 1) + ( k −
1) = nk + i + 1( j + 2) + n + ( k −
1) = nk + i = v where node v is the third local vertex of triangle T j +3 . Thus row v of the matrix S ( U ) has the entry U nk + i +1 − U nk − n + i in columns nk − n + i, v = nk + i , and nk + i + 1 .• Triangle T j + n ) has global vertices ( j + n ) + ( k ) = nk + i − j + n ) + 1 + ( k ) = nk + i = v ( j + n ) + ( n + 1) + ( k ) = nk + i + n where node v is the second local vertex of triangle T j + n ) . Thus row v of the matrix S ( U ) has the entry U nk + i − − U nk + i + n in columns nk + i − , v = nk + i , and nk + i + n .• Triangle T j + n )+1 = T j + n +1) − has global vertices ( j + n + 1) + ( k ) = nk + i = v ( j + n + 1) + ( n + 1) + ( k ) = nk + i + n + 1( j + n + 1) + n + ( k ) = nk + i + n v is the first local vertex of triangle T j + n )+1 . Thus row v of the matrix S ( U ) has the entry U nk + i + n − U nk + i + n +1 in columns v = nk + i, nk + i + n , and nk + i + n + 1 .• Triangle T j + n )+2 = T j + n +1) has global vertices ( j + n + 1) + ( k ) = nk + i = v ( j + n + 1) + 1 + ( k ) = nk + i + 1( j + n + 1) + ( n + 1) + ( k ) = nk + i + n + 1 where node v is the first local vertex of triangle T j + n )+2 . Thus row v of the matrix S ( U ) has the entry U nk + i + n +1 − U nk + i +1 in columns v = nk + i, nk + i + 1 , and nk + i + n + 1 .Thus the nonzero entries in row v = nk + i of S ( U ) for k = 1 , , ..., n − and i = 2 , , ..., n − are in columns nk − n + i − U nk − n + i − − U nk + i − + 16 U nk − n + i − U nk − n + i − = 16 U nk − n + i − U nk + i − nk − n + i : 16 U nk − n + i − U nk − n + i − + 16 U nk + i +1 − U nk − n + i = 16 U nk + i +1 − U nk − n + i − nk + i − U nk − n + i − − U nk + i − + 16 U nk + i − − U nk + i + n = 16 U nk − n + i − − U nk + i + n v = nk + i : 16 U nk − n + i − − U nk + i − + 16 U nk − n + i − U nk − n + i − + 16 U nk + i +1 − U nk − n + i + 16 U nk + i − − U nk + i + n + 16 U nk + i + n − U nk + i + n +1 + 16 U nk + i + n +1 − U nk + i +1 = 0 nk + i + 1 : 16 U nk + i +1 − U nk − n + i + 16 U nk + i + n +1 − U nk + i +1 = 16 U nk + i + n +1 − U nk − n + i nk + i + n : 16 U nk + i − − U nk + i + n + 16 U nk + i + n − U nk + i + n +1 = 16 U nk + i − − U nk + i + n +1 nk + i + n + 1 : 16 U nk + i + n − U nk + i + n +1 + 16 U nk + i + n +1 − U nk + i +1 = 16 U nk + i + n − U nk + i +1
2- Red Boundary Nodes:
The red boundary nodes are of 4 types: a) the left boundary with index v = kn + 1 and k = 1 , , ..., n − , that belong to triangles T k − n − , T k ( n − , and T k ( n − .Thus, we first define the nonzero entries in row v = kn + 1 per local triangle and then add them up.• Triangle T k − n − = T nk − n − k +2) − has global vertices nk − n − k + 2 + ( k −
1) = nk − n + 1 nk − n − k + 2 + ( n + 1) + ( k −
1) = nk + 2 nk − n − k + 2 + n + ( k −
1) = nk + 1 = v where node v is the third local vertex of triangle T k − n − . Thus row v of the matrix S ( U ) has theentry U nk +2 − U nk − n +1 in columns nk − n + 1 , v = nk + 1 , and nk + 2 .• Triangle T k ( n − = T (cid:0) k ( n − (cid:1) − has global vertices k ( n −
1) + 1 + ( k ) = nk + 1 = vk ( n −
1) + 1 + ( n + 1) + ( k ) = nk + n + 2 k ( n −
1) + 1 + n + ( k ) = nk + n + 1 where node v is the first local vertex of triangle T k ( n − . Thus row v of the matrix S ( U ) has the entry U nk + n +1 − U nk + n +2 in columns v = nk + 1 , nk + n + 1 , , and nk + n + 2 .33 Triangle T k ( n − = T k ( n − has global vertices k ( n −
1) + 1 + ( k ) = nk + 1 = vk ( n −
1) + 1 + 1 + ( k ) = nk + 2 k ( n −
1) + 1 + ( n + 1) + ( k ) = nk + n + 2 where node v is the first local vertex of triangle T k ( n − . Thus row v of the matrix S ( U ) has the entry U nk + n +2 − U nk +2 in columns v = nk + 1 , nk + 2 , and nk + n + 2 .Thus the nonzero entries in row v = nk + 1 of S ( U ) for k = 1 , , ..., n − are in columns nk − n + 1 : 16 U nk +2 − U nk − n +1 = 16 U nk +2 − U nk − n +1 v = nk + 1 : 16 U nk +2 − U nk − n +1 + 16 U nk + n +1 − U nk + n +2 + 16 U nk + n +2 − U nk +2 = 16 U nk + n +1 − U nk − n +1 nk + 2 : 16 U nk +2 − U nk − n +1 + 16 U nk + n +2 − U nk +2 = 16 U nk + n +2 − U nk − n +1 nk + n + 1 : 16 U nk + n +1 − U nk + n +2 = 16 U nk + n +1 − U nk + n +2 nk + n + 2 : 16 U nk + n +1 − U nk + n +2 + 16 U nk + n +2 − U nk +2 = 16 U nk + n +1 − U nk +2 b) the right boundary with index v = kn and k = 2 , , ..., n − , that belong to triangles T k − n − , T k − n − − , and T k ( n − .Thus, we first define the nonzero entries in row v = kn per local triangle and then add them up.• Triangle T k − n − − has global vertices ( k − n −
1) + ( k −
2) = nk − n − k − n −
1) + ( n + 1) + ( k −
2) = nk = v ( k − n −
1) + n + ( k −
2) = nk − where node v is the second local vertex of triangle T k − n − − . Thus row v of the matrix S ( U ) hasthe entry U nk − n − − U nk − in columns nk − n − , nk − , and v = nk .• Triangle T k − n − has global vertices ( k − n −
1) + ( k −
2) = nk − n − k − n −
1) + 1 + ( k −
2) = nk − n ( k − n −
1) + ( n + 1) + ( k −
2) = nk = v where node v is the third local vertex of triangle T k − n − . Thus row v of the matrix S ( U ) has theentry U nk − n − U nk − n − in columns nk − n − , nk − n , and v = nk .• Triangle T k ( n − has global vertices k ( n −
1) + ( k −
1) = nk − k ( n −
1) + 1 + ( k −
1) = nk = vk ( n −
1) + ( n + 1) + ( k −
1) = nk + n where node v is the second local vertex of triangle T k ( n − . Thus row v of the matrix S ( U ) has the entry U nk − − U nk + n in columns nk − , v = nk , and nk + n .34hus the nonzero entries in row v = nk of S ( U ) for k = 2 , , ..., n − are in columns nk − n − U nk − n − − U nk − + 16 U nk − n − U nk − n − = 16 U nk − n − U nk − nk − n : 16 U nk − n − U nk − n − = 16 U nk − n − U nk − n − nk − U nk − n − − U nk − + 16 U nk − − U nk + n = 16 U nk − n − − U nk + n v = nk : 16 U nk − n − − U nk − + 16 U nk − n − U nk − n − + 16 U nk − − U nk + n = 16 U nk − n − U nk + n nk + n : 16 U nk − − U nk + n = 16 U nk − − U nk + n c) the lower boundary with index v = i and i = 2 , , ..., n − , that belong to triangles T i − , T i − , and T i − .Thus, we first define the nonzero entries in row v = i per local triangle and then add them up.• Triangle T i − has global vertices i − i − i −
1) + 1 = i = v ( i −
1) + ( n + 1) = i + n where node v = i is the second local vertex of triangle T i − . Thus row v = i of the matrix S ( U ) hasthe entry U i − − U i + n in columns i − , v = i , and i + n .• Triangle T i − = T i − has global vertices i = vi + ( n + 1) = i + n + 1 i + n = i + n where node v = i is the first local vertex of triangle T i − . Thus row v = i of the matrix S ( U ) has theentry U i + n − U i + n +1 in columns v = i, i + n , and i + n + 1 .• Triangle T i − = T i has global vertices i = vi + 1 = i + 1 i + ( n + 1) = i + n + 1 where node v = i is the first local vertex of triangle T i − . Thus row v = i of the matrix S ( U ) has theentry U i + n +1 − U i +1 in columns v = i, i + 1 , and i + n + 1 .Thus the nonzero entries in row v = i of S ( U ) for i = 2 , , ..., n − are in columns i − U i − − U i + n = 16 U i − − U i + n v = i : 16 U i − − U i + n + 16 U i + n − U i + n +1 + 16 U i + n +1 − U i +1 = 16 U i − − U i +1 i + 1 : 16 U i + n +1 − U i +1 = 16 U i + n +1 − U i +1 i + n : 16 U i − − U i + n + 16 U i + n − U i + n +1 = 16 U i − − U i + n +1 i + n + 1 : 16 U i + n − U i + n +1 + 16 U i + n +1 − U i +1 = 16 U i + n − U i +1 ) the upper boundary with index v = n ( n −
1) + i and i = 2 , , ..., n − , that belong to triangles T n − n − i − − , T n − n − i − , and T n − n − i − .Thus, we first define the nonzero entries in row v = n ( n −
1) + i per local triangle and then add them up.• Triangle T n − n − i − − has global vertices ( n − n −
2) + i − n −
2) = n ( n −
2) + i − n − n −
2) + i − n + 1) + ( n −
2) = n ( n −
1) + i = v ( n − n −
2) + i − n + ( n −
2) = n ( n −
1) + i − where node v = n ( n −
1) + i is the second local vertex of triangle T n − n − i − − . Thus row v = n ( n −
1) + i of the matrix S ( U ) has the entry U n ( n − i − − U n ( n − i − in columns n ( n −
2) + i − , n ( n −
1) + i − , and v = n ( n −
1) + i .• Triangle T n − n − i − has global vertices ( n − n −
2) + i − n −
2) = n ( n −
2) + i − n − n −
2) + ( i −
1) + 1 + ( n −
2) = n ( n −
2) + i ( n − n −
2) + ( i −
1) + ( n + 1) + ( n −
2) = n ( n −
1) + i = v where node v = n − n + i is the third local vertex of triangle T i − . Thus row v = i of the matrix S ( U ) has the entry U n ( n − i − U n ( n − i − in columns n ( n − i − , n ( n − i , and v = n ( n − i .• Triangle T n − n − i − = T n − n − i − has global vertices ( n − n −
2) + i + ( n −
2) = n ( n −
2) + i ( n − n −
2) + i + ( n + 1) + ( n −
2) = n ( n −
1) + i + 1( n − n −
2) + i + n + ( n −
2) = n ( n −
1) + i = v where node v = n ( n −
1) + i is the third local vertex of triangle T n − n − i − . Thus row v = n ( n − i of the matrix S ( U ) has the entry U n ( n − i +1 − U n ( n − i in columns n ( n − i, v = n ( n −
1) + i , and n ( n −
1) + i + 1 .Thus the nonzero entries in row v = n ( n −
1) + i of S ( U ) for i = 2 , , ..., n − , are in columns n ( n −
2) + i − U n ( n − i − − U n ( n − i − + 16 U n ( n − i − U n ( n − i − = 16 U n ( n − i − U n ( n − i − n ( n −
2) + i : 16 U n ( n − i − U n ( n − i − + 16 U n ( n − i +1 − U n ( n − i = 16 U n ( n − i +1 − U n ( n − i − n ( n −
1) + i − U n ( n − i − − U n ( n − i − = 16 U n ( n − i − − U n ( n − i − v = n ( n −
1) + i : 16 U n ( n − i − − U n ( n − i − + 16 U n ( n − i − U n ( n − i − + 16 U n ( n − i +1 − U n ( n − i = 16 U n ( n − i +1 − U n ( n − i − n ( n −
1) + i + 1 : 16 U n ( n − i +1 − U n ( n − i = 16 U n ( n − i +1 − U n ( n − i - Green Corner Nodes: There are 2 green corners with index : a) v = 1 that belongs to triangles T , and T . First, we define the nonzero entries in row v = 1 per local triangle.• Triangle T has global vertices v = 1 , n + 2 , n + 1 where node v = 1 is the first local vertex of triangle T .Thus row v = 1 of the matrix S ( U ) has the entry U n +1 − U n +2 in columns v = 1 , n + 1 , and n + 2 .• Triangle T has global vertices v = 1 , , n + 2 where node v = 1 is the first local vertex of triangle T .Thus row v = 1 of the matrix S ( U ) has the entry U n +2 − U in columns v = 1 , , and n + 2 .Thus the nonzero entries in row v = 1 of S ( U ) are in columns v = 1 : 16 U n +1 − U n +2 + 16 U n +2 − U = 16 U n +1 − U U n +2 − U = 16 U n +2 − U n + 1 : 16 U n +1 − U n +2 = 16 U n +1 − U n +2 n + 2 : 16 U n +1 − U n +2 + 16 U n +2 − U = 16 U n +1 − U b) v = n that belongs to triangles T n − n − − , and T n − n − .First, we define the nonzero entries in row v = n per local triangle and then add them up.• Triangle T n − n − − has global vertices n ( n − − , v = n , n − where node v = n is thesecond local vertex of triangle T n − n − − . Thus row v = n of the matrix S ( U ) has the entry U n ( n − − − U n − in columns n ( n − − , n − , and v = n .• Triangle T n − n − − has global vertices n ( n − − , n ( n − , v = n , where node v = n is the thirdlocal vertex of triangle T . Thus row v = n of the matrix S ( U ) has the entry U n ( n − − U n ( n − − in columns n ( n − − , n ( n − , and v = n .Thus the nonzero entries in row v = n of S ( U ) are in columns n ( n − − U n ( n − − − U n − + 16 U n ( n − − U n ( n − − = 16 U n ( n − − U n − n ( n −
1) : 16 U n ( n − − U n ( n − − = 16 U n ( n − − U n ( n − − n − U n ( n − − − U n − = 16 U n ( n − − − U n − v = n : 16 U n ( n − − − U n − + 16 U n ( n − − U n ( n − − = 16 U n ( n − − U n −
4- Blue Corner Nodes:
There are 2 blue corners with index : a) v = n that belongs to triangle T n − with global vertices n − , v = n, n , where node v = n is the secondlocal vertex of triangle T n − . Thus row v = n of the matrix S ( U ) has the entry U n − − U n in columns n − , v = n , and n . 37 ) v = n ( n −
1) + 1 that belongs to triangle T n − n − with global vertices n ( n −
2) + 1 , n ( n −
1) + 2 ,v = n ( n −
1) + 1 where node v = n ( n −
1) + 1 is the third local vertex of triangle T n − n − − . Thus row v = n ( n −
1) + 1 of the matrix S ( U ) has the entry U n ( n − − U n ( n − in columns n ( n −
2) + 1 ,v = n ( n −
1) + 1 , and n ( n −
1) + 2 .For n = 5 , the matrix S ( U ) corresponding to the mesh in Figure 2, has the following sparsity pattern: S ( U ) = ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
00 0 0 0 0 0 0 0 0 0 0 0 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
In general, the S ( U ) matrix is an n × n block tridiagonal matrix with at most 6 nonzero entries per row and 6 percolumn. S ( U ) = 16 A , A , · · · A , A , A , · · · . . . . . . . . . . . . ...... . . . A j,k A j,j A j,i · · · A i,j A i,i A i,n · · · A n,i A n,n where i = n − , j = n − , k = n − , and the block matrices are of size n × n with the following sparsity patterns:• A , and A n,n are tridiagonal matrices, each with n − nonzero entries.• A i,i for i = 2 , , .., n − are tridiagonal matrices with zero diagonal entries except for the first and the last.• A i +1 ,i for i = 1 , , , .., n − are lower bidiagonal matrices, with n − nonzero entries.• A i,i +1 for i = 1 , , , .., n − are upper bidiagonal matrices, with n − nonzero entries.Thus, S ( U ) has a total of n − n − n )+(2 n − n −
1) = 6 n − n − n +4 n − n +2 = 6 n − n − nonzero entries. As for the nonzero entries in S ( U ) , each is of order || U || ∞ since each nonzero entry is of the form
16 ( U i − U j ) ≤ | U i − U j | ≤
16 ( | U i | + | U j | ) ≤ || U || ∞ .2 The Matrix S ( U ) assuming Periodic Boundary Conditions on U Given the meshing of Figure 2 and assuming periodic boundary conditions, where the values of U are equal at theupper and lower red vertices, the left and right red boundary vertices, and the corner vertices i.e. U kn +1 = U kn + n , f or k = 1 , , .., n − (106) U i = U i + n ( n − , f or i = 2 , , .., n − (107) U = U n = U n ( n − = U n (108)and assuming that our domain is torus shaped and that the vector U is of size ( n − , U = [ U , .., U n − , U n +1 , ..., U n − , U n +1 , ..., U n − , ......, U n ( n − − ] then S ( U ) is an ( n − × ( n − sparse matrix with at most 6 nonzero entries per row, as discussed below. This“periodic” S ( U ) matrix, can be obtained from the general one described in the previous section by merging/addingthe rows corresponding to equal U entries and also the columns, and using the periodicity of U .Note that since rows/columns n, n, n, ..., ( n − n of the general S ( U ) matrix are merged with other rows/columns,then the indices have to be reindexed to get the corresponding rows/columns of the “periodic” matrix S ( U ) as such: [1 , · · · , n − , n + 1 , · · · , n − , n + 1 , · · · , n − , · · · · · · , ( n − n − ↓ (109) [1 , · · · , n − , n, · · · , n − , n − , · · · , n − , · · · · · · , ( n − (cid:3) But the indices of the U vector are not reindexed in what follows. The entries in the matrix S ( U ) that will be modifiedare the ones corresponding to the green left corner, lower and left red boundary nodes, and the upper and left boundaryblack nodes in Figure 2.
1- Green Left Node:
Assuming that the vertices , n , n ( n − and n coincide, then the first row of the “periodic” S ( U ) will be the sum of the entries in rows , n , n ( n −
1) + 1 and n of the general S ( U ) : a) Row v = n of the general matrix S ( U ) has the entry U n − − U n = U n − − U n +1 in columns n − , v = n = 1 , and n = n + 1 . b) Row v = n ( n − of the general matrix S ( U ) has the entry U n ( n − − U n ( n − = U − U n ( n − in columns n ( n −
2) + 1 , v = n ( n −
1) + 1 = 1 , and n ( n −
1) + 2 = 2 . c) Row v = 1 of the general matrix S ( U ) has nonzero entries in columns v = 1 : 16 U n +1 − U n +2 + 16 U n +2 − U = 16 U n +1 − U U n +2 − U = 16 U n +2 − U n + 1 : 16 U n +1 − U n +2 = 16 U n +1 − U n +2 n + 2 : 16 U n +1 − U n +2 + 16 U n +2 − U = 16 U n +1 − U d) Row v = n of the general matrix S ( U ) has nonzero entries in columns n ( n − − U n ( n − − U n − = 16 U n ( n − − U n − n ( n −
1) = n ( n −
2) + 1 : 16 U n ( n − − U n ( n − − = 16 U n ( n − − U n ( n − − n − n − U n ( n − − − U n − = 16 U n ( n − − − U n − v = n = 1 : 16 U n ( n − − U n − = 16 U n ( n − − U n − of the “periodic” S ( U ) with the following nonzero entries: U n − − U n +1 + 16 U − U n ( n − + 16 U n +1 − U + 16 U n ( n − − U n − = 02 : 16 U − U n ( n − + 16 U n +2 − U = 16 U n +2 − U n ( n − n − U n − − U n +1 + 16 U n ( n − − − U n − = 16 U n ( n − − − U n +1 n + 1 : 16 U n − − U n +1 + 16 U n +1 − U n +2 = 16 U n − − U n +2 n + 2 : 16 U n +1 − U n ( n −
2) + 1 : 16 U − U n ( n − + 16 U n ( n − − U n ( n − − = 16 U − U n ( n − − n ( n − − U n ( n − − U n − Recall that these column indices have to be reindex by (109).
2- Lower Red Boundary Nodes:
Assuming that the vertices i and i + n ( n − coincide for i = 2 , , .., n − , thenthe rows i will be the sum of the entries in rows i and i + n ( n − of the general S ( U ) : a) Row i of the general matrix S ( U ) for i = 2 , , ..., n − has nonzero entries in columns i − U i − − U i + n ) / v = i : ( U i − − U i +1 ) / i + 1 : ( U i + n +1 − U i +1 ) / i + n : ( U i − − U i + n +1 ) / i + n + 1 : 16 U i + n − U i +1 b) Row v = n ( n −
1) + i of the general matrix S ( U ) for i = 2 , , ..., n − , has nonzero entries in columns n ( n −
2) + i − U n ( n − i − U n ( n − i − = 16 U n ( n − i − U i − n ( n −
2) + i : 16 U n ( n − i +1 − U n ( n − i − = 16 U i +1 − U n ( n − i − n ( n −
1) + i − i − U n ( n − i − − U n ( n − i − = 16 U n ( n − i − − U i − v = n ( n −
1) + i = i : 16 U n ( n − i +1 − U n ( n − i − = 16 U i +1 − U i − n ( n −
1) + i + 1 = i + 1 : 16 U n ( n − i +1 − U n ( n − i = 16 U i +1 − U n ( n − i Adding up these 2 rows, we get row i of the “periodic” S ( U ) with the following nonzero entries for i = 2 , , ..., n − : i − U i − − U i + n + 16 U n ( n − i − − U i − = 16 U n ( n − i − − U i + n i : 16 U i − − U i +1 + 16 U i +1 − U i − = 0 i + 1 : 16 U i + n +1 − U i +1 + 16 U i +1 − U n ( n − i = 16 U i + n +1 − U n ( n − i i + n : 16 U i − − U i + n +1 i + n + 1 : 16 U i + n − U i +1 ( n −
2) + i − U n ( n − i − U i − n ( n −
2) + i : 16 U i +1 − U n ( n − i − Recall that these column indices have to be reindex by (109). Note that for i = n − we get the following nonzeroentries: n − U n ( n − n − − U n − n = 1 : 16 U n − U n ( n − n − = 16 U n +1 − U n ( n − n − n − U n − − U n = 16 U n − − U n +1 n = n + 1 : 16 U n − − U n = 16 U n − − U n ( n −
2) + n − U n ( n − n − − U n − n ( n −
2) + n − U n − U n ( n − n − = 16 U − U n ( n − n −
3- Left Red Boundary Nodes:
Assuming that the vertices kn + 1 and kn + n coincide for k = 1 , , .., n − , thenthe corresponding rows kn + 1 − k of the “periodic” S ( U ) will be the sum of the entries in rows kn + 1 and kn + n of the general S ( U ) : a) Row kn + 1 of the general matrix S ( U ) for k = 1 , , .., n − , has nonzero entries in columns nk − n + 1 : ( U nk +2 − U nk − n +1 ) / nk + 1 : 16 U nk + n +1 − U nk − n +1 nk + 2 : 16 U nk + n +2 − U nk − n +1 nk + n + 1 : 16 U nk + n +1 − U nk + n +2 nk + n + 2 : 16 U nk + n +1 − U nk +2 b) Row ( k + 1) n of the general matrix S ( U ) for k = 1 , , ..., n − , has nonzero entries in columns n ( k + 1) − n − nk − U n ( k +1) − n − U n ( k +1) − = 16 U nk − n +1 − U nk + n − n ( k + 1) − n = nk = nk − n + 1 : 16 U n ( k +1) − n − U n ( k +1) − n − = 16 U nk − n +1 − U nk − n ( k + 1) − nk + n − U n ( k +1) − n − − U n ( k +1)+ n = 16 U nk − − U nk + n +1 n ( k + 1) = nk + 1 : 16 U n ( k +1) − n − U n ( k +1)+ n = 16 U nk − n +1 − U nk + n +1 n ( k + 1) + n = nk + n + 1 : 16 U n ( k +1) − − U n ( k +1)+ n = 16 U nk + n − − U nk + n +1 Adding up these 2 rows, we get row nk + 1 − k of the “periodic” S ( U ) with the following nonzero entries for k = 1 , , ..., n − : nk − n + 1 : 16 U nk +2 − U nk − n +1 + 16 U nk − n +1 − U nk − = 16 U nk +2 − U nk − nk − U nk − n +1 − U nk + n − k + 1 : 16 U nk + n +1 − U nk − n +1 + 16 U nk − n +1 − U nk + n +1 = 0 nk + 2 : 16 U nk + n +2 − U nk − n +1 nk + n − U nk − − U nk + n +1 nk + n + 1 : 16 U nk + n +1 − U nk + n +2 + 16 U nk + n − − U nk + n +1 = 16 U nk + n − − U nk + n +2 nk + n + 2 : 16 U nk + n +1 − U nk +2 Recall that these column indices have to be reindex by (109). Note that for k = n − we get the following nonzeroentries: n ( n −
3) + 1 : 16 U n ( n − − U n ( n − − n ( n − − U n ( n − − U n ( n − − n ( n −
2) + 2 : 16 U n ( n − − U n ( n − = 16 U − U n ( n − n ( n − − U n ( n − − − U n ( n − = 16 U n ( n − − − U n ( n −
1) + 1 = 1 : 16 U n ( n − − − U n ( n − = 16 U n ( n − − − U n ( n −
1) + 2 = 2 : 16 U n ( n − − U n ( n − = 16 U − U n ( n −
4- Right Black Boundary Nodes:
Rows ( k +1) n − k − of the “periodic” S ( U ) matrix correspond to rows ( k +1) n − of the general S ( U ) for k = 1 , .., n − with nonzero entries in columns: nk − n + n − − nk − U nk − n + n − − U nk + n − − = 16 U nk − − U nk + n − nk − n + n − nk − U nk + n − − U nk − n + n − − = 16 U nk +1 − U nk − nk + n − − nk + n − U nk − n + n − − − U nk + n − n = 16 U nk − − U nk +2 n − nk + n − nk + n = nk + 1 : 16 U nk + n − n +1 − U nk − n + n − = 16 U nk + n +1 − U nk − nk + n − n = nk + 2 n − U nk + n − − − U nk + n − n +1 = 16 U nk + n − − U nk + n +1 nk + n − n + 1 = nk + 2 n = nk + n + 1 : 16 U nk + n − n − U nk + n − = 16 U nk +2 n − − U nk +1 Recall that these column indices have to be reindex by (109).
5- Upper Black Boundary Nodes:
Rows ( n − n −
2) + i of the “periodic” S ( U ) matrix correspond to rows n ( n −
2) + i of the general S ( U ) for i = 2 , , .., n − with nonzero entries in columns: n ( n − − n + i − n ( n −
3) + i − U n ( n − − n + i − U n ( n − i − = 16 U n ( n − i − U n ( n − i − n ( n − − n + i = n ( n −
3) + i : 16 U n ( n − i +1 − U n ( n − − n + i − = 16 U n ( n − i +1 − U n ( n − i − ( n −
2) + i − U n ( n − − n + i − − U n ( n − i + n = 16 U n ( n − i − − U i n ( n −
2) + i + 1 : 16 U n ( n − i + n +1 − U n ( n − − n + i = 16 U i +1 − U n ( n − i n ( n −
2) + i + n = n ( n −
1) + i = i : 16 U n ( n − i − − U n ( n − i + n +1 = 16 U n ( n − i − − U i +1 n ( n −
2) + i + n + 1 = n ( n −
1) + i + 1 = i + 1 : 16 U n ( n − i + n − U n ( n − i +1 = 16 U i − U n ( n − i +1 Recall that these column indices have to be reindex by (109). Note that for i = n − we get the following nonzeroentries: n ( n −
3) + n − − n ( n − − U n ( n − n − − U n ( n − n − − = 16 U n ( n − − − U n ( n − − n ( n −
3) + n − n ( n − − U n ( n − n − − U n ( n − n − − = 16 U n ( n − − U n ( n − − n ( n −
2) + n − − n ( n − − U n ( n − n − − − U n − = 16 U n ( n − − − U n − n ( n −
2) + n − n ( n −
1) = n ( n −
2) + 1 : 16 U n − − U n ( n − n − = 16 U − U n ( n − − n − U n ( n − n − − − U n − = 16 U n ( n − − − U n − n = 1 : 16 U n − − U n ( n − n − = 16 U n − − U n ( n − For n = 5 , the matrix S ( U ) corresponding to the mesh in Figure 2, has the following sparsity pattern: S ( U ) = ∗ ∗ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
00 0 0 0 0 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗
00 0 ∗ ∗ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗ with U = [ U , U , U , U , U , U , U , U , U , U , U , U , U , U , U , U ] In general, the S ( U ) matrix is an ( n − × ( n − block tridiagonal matrix with 2 additional blocks in the upperright and lower left corner. Moreover, it is a skew-symmetric matrix ( A T = − A ) with 6 nonzero entries per row, 6nonzero entries per column, and zeros on the diagonal. S ( U ) = 16 A , A , · · · A ,n − A , A , A , · · · . . . . . . . . . . . . ...... . . . A j,h A j,j A j,i · · · A i,j A i,i A i,k A k, · · · A k,i A k,k i = n − , j = n − , k = n − , h = n − , and the n − nonzero block matrices are of size ( n − × ( n − with n − nonzero entries each, and the following sparsity patterns:• A i,i for i = 1 , , .., n − are tridiagonal matrices with zero diagonal entries, and nonzero A i,i (1 , n − , A i,i ( n − , .• A ,n − and A i +1 ,i for i = 1 , , , .., n − are lower bidiagonal matrices, with nonzero entry in first row, column n − .• A n − , and A i,i +1 for i = 1 , , .., n − are upper bidiagonal matrices with nonzero entry in first column, row n − .Thus, S ( U ) has a total of n − n −
1) = 6( n − nonzero entries.As for the nonzero entries in S ( U ) , each is of order || U || ∞ since each nonzero entry is of the form
16 ( U i − U j ) ≤ | U i − U j | ≤
16 ( | U i | + | U j | ) ≤ || U || ∞ B Assembling the Global matrix R First, we derive the sparsity pattern of the global matrix R without any assumptions related to boundary conditions,i.e. define the general matrix R , in Section B.1. Then, we derive the sparsity pattern of the matrix R assuming periodicboundary conditions in Section B.2. B.1 The General Matrix R Given the meshing of Figure 2 and without any assumptions on the boundary nodes, R is an n × n sparse matrixwith at most 6 nonzero entries per row, as discussed below. To get those nonzero entries, note that each node has atmost 6 edges connecting it with its neighboring nodes, i.e. belongs to at most 6 triangles. We consider the 4 types ofnodes shown in different colors in Figure 2: the ( n − black internal nodes that belong to 6 triangles, the n − red boundary nodes that belong to 4 triangles, the 2 green corner nodes that belong to 2 triangles, and the 2 blue cornernodes that belong to 1 triangle.
1- Black Internal Nodes:
Each of the ( n − black internal nodes with index v = ln + i for l = 1 , , ..., n − and i = 2 , , ..., n − belongs to triangles T j +1 , T j +2 , T j +3 , T j + n ) , T j + n )+1 , T j + n )+2 , where j = ( n − l −
1) + ( i − .Thus, we first define the nonzero entries in row v per local triangle and then add them up.• Triangle T j +1 of type a, has node v as the second local vertex. Thus row v of the matrix R has the entry incolumns nl − n + i − , nl + i − , and v = nl + i .• Triangle T j +2 of type b, has node v as the third local vertex. Thus row v of the matrix R has the entry h k incolumns nl − n + i − , nl − n + i , and v = nl + i .• Triangle T j +3 of type a, has node v as the third local vertex. Thus row v of the matrix R has the entry h k incolumns nl − n + i, v = nl + i , and nl + i + 1 .• Triangle T j + n ) of type b, has node v as the second local vertex . Thus row v of the matrix R has the entry − h k in columns nl + i − , v = nl + i , and nl + i + n .• Triangle T j + n )+1 of type a, has node v as the first local vertex. Thus row v of the matrix R has the entry − h k in columns v = nl + i, nl + i + n , and nl + i + n + 1 .• Triangle T j + n )+2 of type b, has node v as the first local vertex. Thus row v of the matrix R has the entry incolumns v = nl + i, nl + i + 1 , and nl + i + n + 1 .44hus the nonzero entries in row v = nl + i of R for l = 1 , , ..., n − and i = 2 , , ..., n − are in columns nl − n + i − h k = h knl − n + i : h k + h k = h knl + i − − h k = − h kv = nl + i : h k + h k − h k − h k = 0 nl + i + 1 : h k = h knl + i + n : − h k − h k = − h knl + i + n + 1 : − h k = − h k
2- Red Boundary Nodes:
The red boundary nodes are of 4 types: a) the left boundary with index v = ln + 1 and l = 1 , , ..., n − , that belong to triangles T l − n − , T l ( n − , and T l ( n − .Thus, we first define the nonzero entries in row v = ln + 1 per local triangle and then add them up.• Triangle T l − n − of type a, has node v as the third local vertex. Thus row v of the matrix R has theentry h k in columns nl − n + 1 , v = nl + 1 , and nl + 2 .• Triangle T l ( n − of type a, has node v as the first local vertex. Thus row v of the matrix R has the entry − h k in columns v = nl + 1 , nl + n + 1 , , and nl + n + 2 .• Triangle T l ( n − has node v is the first local vertex. Thus row v of the matrix R has the entry incolumns v = nl + 1 , nl + 2 , and nl + n + 2 .Thus the nonzero entries in row v = nl + 1 of R for l = 1 , , ..., n − are in columns nl − n + 1 : h k = h kv = nl + 1 : h k − h k = 0 nl + 2 : h k = h knl + n + 1 : − h k = − h knl + n + 2 : − h k = − h k b) the right boundary with index v = ln and l = 2 , , ..., n − , that belong to triangles T l − n − , T l − n − − , and T l ( n − .Thus, we first define the nonzero entries in row v = ln per local triangle and then add them up.• Triangle T l − n − − of type a, has node v is the second local vertex. Thus row v of the matrix R hasthe entry in columns nl − n − , nl − , and v = nl .• Triangle T l − n − of type b, has node v as the third local vertex. Thus row v of the matrix R has theentry h k in columns nl − n − , nl − n , and v = nl .45 Triangle T l ( n − of type b, has node v is the second local vertex . Thus row v of the matrix R has theentry − h k in columns nl − , v = nl , and nl + n .Thus the nonzero entries in row v = nl of R for l = 2 , , ..., n − are in columns nl − n − h k = h knl − n : h k = h knl − − h k = − h kv = nl : h k − h k = 0 nl + n : − h k = − h k c) the lower boundary with index v = i and i = 2 , , ..., n − , that belong to triangles T i − , T i − , and T i − .Thus, we first define the nonzero entries in row v = i per local triangle and then add them up.• Triangle T i − of type b, has node v = i as the second local vertex. Thus row v = i of the matrix R hasthe entry − h k in columns i − , v = i , and i + n .• Triangle T i − of type a, has node v = i as the first local vertex. Thus row v = i of the matrix R hasthe entry − h k in columns v = i, i + n , and i + n + 1 .• Triangle T i − of type b, has node v = i as the first local vertex. Thus row v = i of the matrix R hasthe entry in columns v = i, i + 1 , and i + n + 1 .Thus the nonzero entries in row v = i of R for i = 2 , , ..., n − are in columns i − − h k = − h kv = i : − h k − h k = − h ki + n : − h k − h k = − h ki + n + 1 : − h k = − h k d) the upper boundary with index v = n ( n −
1) + i and i = 2 , , ..., n − , that belong to triangles T n − n − i − − , T n − n − i − , and T n − n − i − .Thus, we first define the nonzero entries in row v = n ( n −
1) + i per local triangle and then add them up.• Triangle T n − n − i − − of type a, has node v = n ( n −
1) + i as the second local vertex. Thusrow v = n ( n −
1) + i of the matrix R has the entry in columns n ( n −
2) + i − , n ( n −
1) + i − , and v = n ( n −
1) + i .• Triangle T n − n − i − of type b, has node v = n − n + i as the third local vertex. Thus row v = i of the matrix R has the entry h k in columns n ( n −
2) + i − , n ( n −
2) + i , and v = n ( n −
1) + i .• Triangle T n − n − i − of type a, has node v = n ( n −
1) + i as the third local vertex. Thus row v = n ( n −
1) + i of the matrix R has the entry h k in columns n ( n −
2) + i, v = n ( n −
1) + i , and n ( n −
1) + i + 1 . 46hus the nonzero entries in row v = n ( n −
1) + i of R for i = 2 , , ..., n − , are in columns n ( n −
2) + i − h k = h kn ( n −
2) + i : h k + h k = h kv = n ( n −
1) + i : h k + h k = h kn ( n −
1) + i + 1 : h k = h k
3- Green Corner Nodes:
There are 2 green corners with index : a) v = 1 that belongs to triangles T , and T . First, we define the nonzero entries in row v = 1 per local triangle.• Triangle T of type a, has node v = 1 is the first local vertex. Thus row v = 1 of the matrix R has theentry − h k in columns v = 1 , n + 1 , and n + 2 .• Triangle T of type b, has node v = 1 as the first local vertex of triangle T . Thus row v = 1 of the matrix R has the entry in columns v = 1 , , and n + 2 .Thus the nonzero entries in row v = 1 of R are in columns v = 1 : − h kn + 1 : − h kn + 2 : − h k b) v = n that belongs to triangles T n − n − − , and T n − n − .First, we define the nonzero entries in row v = n per local triangle and then add them up.• Triangle T n − n − − of type a, has node v = n as the second local vertex of triangle T n − n − − .Thus row v = n of the matrix R has the entry in columns n ( n − − , n − , and v = n .• Triangle T n − n − of type n, has node v = n as the third local vertex. Thus row v = n of the matrix R has the entry h k in columns n ( n − − , n ( n − , and v = n .Thus the nonzero entries in row v = n of R are in columns n ( n − − h kn ( n −
1) : h kv = n : h k
4- Blue Corner Nodes:
There are 2 blue corners with index : a) v = n that belongs to triangle T n − of type b where node v = n is the second local vertex. Thus row v = n of the matrix R has the entry − h k in columns n − , v = n , and n .47 ) v = n ( n −
1) + 1 that belongs to triangle T n − n − of type a, where node v = n ( n −
1) + 1 is the thirdlocal vertex. Thus row v = n ( n −
1) + 1 of the matrix R has the entry h k in columns n ( n −
2) + 1 ,v = n ( n −
1) + 1 , and n ( n −
1) + 2 .For n = 5 , the matrix R corresponding to the mesh in Figure 2, is given by h kA where A = − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − −
10 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 − −
10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 2 1 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 2 1 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 2 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1
In general, the R matrix is an n × n block tridiagonal matrix with at most 6 nonzero entries per row and 6 percolumn. R = h k A , A , · · · A , A , A , · · · . . . . . . . . . . . . ...... . . . A j,k A j,j A j,i · · · A i,j A i,i A i,n · · · A n,i A n,n where i = n − , j = n − , k = n − , and the block matrices are of size n × n with the following sparsity patterns:• A , is a lower bidiagonal matrix with n − nonzero entries.• A n,n = − A T , .• A i,i for i = 2 , , .., n − are tridiagonal matrices with zero diagonal entries with A ( i, i ) T = − A ( i, i ) .• A i +1 ,i = − A , for i = 1 , , , .., n − • A i,i +1 = A T , for i = 1 , , , .., n − Thus, R has a total of ( n − n −
2) + (2 n )(2 n −
1) = 2 n − n + 4 + 4 n − n = 6 n − n + 4 nonzero entries.As for the nonzero entries in R , their absolute values are less than or equal to h | ˆ k | .48 .2 The Matrix R assuming Periodic Boundary Conditions Given the meshing of Figure 2 and assuming periodic boundary conditions, where the values of U are equal at theupper and lower red vertices, the left and right red boundary vertices, and the corner vertices i.e. U kn +1 = U kn + n , f or k = 1 , , .., n − (110) U i = U i + n ( n − , f or i = 2 , , .., n − (111) U = U n = U n ( n − = U n (112)and assuming that our domain is torus shaped and that the vector U is of size ( n − , U = [ U , .., U n − , U n +1 , ..., U n − , U n +1 , ..., U n − , ......, U n ( n − − ] then R is an ( n − × ( n − sparse matrix with at most 6 nonzero entries per row, as discussed below. This“periodic” R matrix, can be obtained from the general one described in the previous section by merging/adding therows corresponding to equal U entries and also the columns.Note that since rows/columns n, n, n, ..., ( n − n of the general R matrix are merged with other rows/columns,then the indices have to be reindexed to get the corresponding rows/columns of the “periodic” matrix R as such: [1 , · · · , n − , n + 1 , · · · , n − , n + 1 , · · · , n − , · · · · · · , ( n − n − ↓ (113) [1 , · · · , n − , n, · · · , n − , n − , · · · , n − , · · · · · · , ( n − (cid:3) The entries in the matrix R that will be modified are the ones corresponding to the green left corner, lower and left redboundary nodes, and the upper and left boundary black nodes in Figure 2.
1- Green Left Node:
Assuming that the vertices , n , n ( n − and n coincide, then the first row of the “periodic” R will be the sum of the entries in rows , n , n ( n −
1) + 1 and n of the general R : a) Row v = 1 of the general matrix R has the entry − h k in columns , n + 1 , n + 2 . b) Row v = n of the general matrix R has the entry − h k in columns n − , v = n = 1 , and n = n + 1 . c) Row v = n ( n −
1) + 1 of the general matrix R has the entry has the entry h k in columns n ( n −
2) + 1 , v = n ( n −
1) + 1 = 1 , and n ( n −
1) + 2 = 2 . d) Row v = n of the general matrix R has the entry h k in columns n ( n − − , n ( n −
1) = n ( n −
2) + 1 , v = n = 1 .Adding up these rows, we get row of the “periodic” R with the following nonzero entries: − h k − h k + h k + h k = 02 : h kn − − h k n + 1 : − h k − h k = − h kn + 2 : − h k n ( n −
2) + 1 : h k + h k = h kn ( n − − h k Recall that these column indices have to be reindex by (113).
2- Lower Red Boundary Nodes:
Assuming that the vertices i and i + n ( n − coincide for i = 2 , , .., n − , thenthe rows i will be the sum of the entries in rows i and i + n ( n − of the general S ( U ) : a) Row i of the general matrix R for i = 2 , , ..., n − has nonzero entries in columns49 − − h ki + n : − h k v = i : − h ki + n + 1 : − h k b) Row v = n ( n −
1) + i of the general matrix R for i = 2 , , ..., n − , has nonzero entries in columns n ( n −
2) + i − h kn ( n −
2) + i : h k v = n ( n −
1) + i = i : h kn ( n −
1) + i + 1 = i + 1 : h k Adding up these 2 rows, we get row i of the “periodic” R with the following nonzero entries for i = 2 , , ..., n − : i − − h ki + 1 : h k i + n : − h ki + n + 1 : − h k n ( n −
2) + i − h kn ( n −
2) + i : h k i : h k − h k = 0 Recall that these column indices have to be reindex by (113). Note that for i = n − we get the following nonzeroentries: n − − h kn = 1 : h k n − − h k n = n + 1 : − h k n ( n −
2) + n − h kn ( n −
2) + n − h k
3- Left Red Boundary Nodes:
Assuming that the vertices ln + 1 and ln + n coincide for l = 1 , , .., n − , then thecorresponding rows ln + 1 − l of the “periodic” R will be the sum of the entries in rows ln + 1 and ln + n of thegeneral R : a) Row ln + 1 of the general matrix R for l = 1 , , .., n − , has nonzero entries in columns nl − n + 1 : h knl + 2 : h k nl + n + 1 : − h knl + n + 2 : − h k b) Row ( l + 1) n of the general matrix R for l = 1 , , ..., n − , has nonzero entries in columns nl − h knl = nl − n + 1 : h k nl + n − − h knl + n + n = nl + n + 1 : − h k Adding up these 2 rows, we get row nl + 1 − l of the “periodic” R with the following nonzero entries for l =1 , , ..., n − : nl − n + 1 : h k + h k = h knl + n + 1 : − h k − h k = − h k nl − h knl + 2 : h k nl + n − − h knl + n + 2 : − h k Recall that these column indices have to be reindex by (113). Note that for l = n − we get the following nonzeroentries: n ( n − − n + 1 : h kn ( n −
2) + n + 1 = 1 : − h k n ( n − − h kn ( n −
2) + 2 : h k n ( n −
2) + n − − h kn ( n −
2) + n + 2 = 2 : − h k - Right Black Boundary Nodes: Rows ( l + 1) n − l − of the “periodic” R matrix correspond to rows ( l + 1) n − of the general R for l = 1 , .., n − with nonzero entries in columns: nl − h knl − h k nl + n − − h knl + 2 n − − h k nl + n = nl + 1 : h knl + 2 n = nl + n + 1 : − h k Recall that these column indices have to be reindex by (113).
5- Upper Black Boundary Nodes:
Rows ( n − n − i of the “periodic” R matrix correspond to rows n ( n − i of the general R for i = 2 , , .., n − with nonzero entries in columns: n ( n −
3) + i − h kn ( n −
3) + i : h k n ( n −
2) + i − − h kn ( n −
2) + i + 1 : h k n ( n −
1) + i = i : − h kn ( n −
1) + i + 1 = i + 1 : − h k Recall that these column indices have to be reindex by (113). Note that for i = n − we get the following nonzeroentries: n ( n − − h kn ( n − − h k n ( n − − − h kn − − h k n ( n −
1) = n ( n −
2) + 1 : h kn − n = 1 : − h k For n = 5 , the matrix R corresponding to the mesh in Figure 2 has the following sparsity pattern R = h k − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − −
10 0 0 0 0 0 1 2 1 0 − − − − − − − − − − − − − − − In general, the R matrix is an ( n − × ( n − block tridiagonal matrix with 2 additional blocks in the upper rightand lower left corner. Moreover, it is a skew-symmetric matrix ( R T = − R ) with zeros on the diagonal, and 6 nonzeroentries per row of the form h kα , 6 nonzero entries per column, of the form h kα where α = − , − , − , , , . R = h k A , A , · · · A ,n − A , A , A , · · · . . . . . . . . . . . . ...... . . . A j,m A j,j A j,i · · · A i,j A i,i A i,l A l, · · · A l,i A l,l where i = n − , j = n − , l = n − , m = n − , and the n − nonzero block matrices are of size ( n − × ( n − with n − nonzero entries each, and the following sparsity patterns:51 A i,i for i = 1 , , .., n − are such that, A i,i ( j, j + 1) = 1 , A i,i ( j + 1 , j ) = − , for j = 1 , , ..., n − A i,i (1 , n −
1) = − , and A i,i ( n − ,
1) = 1 .• A ,n − = A i +1 ,i for i = 1 , , .., n − are lower bidiagonal matrices ( A i +1 ,i ( j, j ) = 2 , A i +1 ,i ( j + 1 , j ) = 1 ),with A i +1 ,i (1 , n −
1) = 1 ,• A n − , = A i,i +1 for i = 1 , .., n − are upper bidiagonal matrices ( A i,i +1 ( j, j ) = − , A i,i +1 ( j, j + 1) = − ),with A i,i +1 ( n − ,
1) = − ,Thus, R has a total of n − n −