A genuinely stable Lagrange-Galerkin scheme for convection-diffusion problems
aa r X i v : . [ m a t h . NA ] M a y A GENUINELY STABLE LAGRANGE–GALERKIN SCHEME FORCONVECTION-DIFFUSION PROBLEMS
MASAHISA TABATA
Department of Mathematics, Waseda University, 3-4-1, Ohkubo, Shinjuku, Tokyo169-8555, Japan tabata@ waseda. jp
SHINYA UCHIUMI
Research Fellow of Japan Society for the Promotion of ScienceGraduate School of Fundamental Science and Engineering, Waseda University,3-4-1, Ohkubo, Shinjuku, Tokyo 169-8555, Japan su48@ fuji. waseda. jp
Abstract.
We present a Lagrange–Galerkin scheme free from numerical quad-rature for convection-diffusion problems. Since the scheme can be implementedexactly as it is, theoretical stability result is assured. While conventionalLagrange–Galerkin schemes may encounter the instability caused by numer-ical quadrature error, the present scheme is genuinely stable. For the P k -element we prove error estimates of O (∆ t + h + h k +1 ) in ℓ ∞ ( L )-norm and of O (∆ t + h + h k ) in ℓ ∞ ( H )-norm. Numerical results reflect these estimates. Introduction
The Lagrange–Galerkin method, which is also called characteristics finite ele-ment method or Galerkin-characteristics method, is a powerful numerical methodfor flow problems such as the convection-diffusion equations and the Navier–Stokesequations. In this method the material derivative is discretized along the charac-teristic curve, which originates the robustness for convection-dominated problems.Although, as a result of the discretization along the characteristic curve, a compos-ite function term at the previous time appears, it is converted to the right-handside in the system of the linear equations. Thus, the coefficient matrix in the left-hand side is symmetric, which allows us to use efficient linear solvers for symmetricmatrices such as the conjugate gradient method and the minimal residual method[2, 15].
Date : April 18, 2015.2000
Mathematics Subject Classification.
Key words and phrases.
Lagrange–Galerkin scheme and Finite element method andConvection-diffusion problems and Exact integration . tability and error analysis of LG schemes has been done in [1, 3, 4, 6, 9, 10,11, 12, 13, 14, 16]; see also the bibliography therein. Pironneau [11] analyzedconvection-diffusion problems and the Navier–Stokes equations to obtain subopti-mal convergence results. Optimal convergence results were obtained by Douglas–Russell [6] for convection-diffusion problems and by S¨uli [16] for the Navier–Stokesequations. Optimal convergence results of second order in time were obtained byBoukir et al. [4] for the Navier–Stokes equations in multi-step method and byRui–Tabata [14] for convection-diffusion problems in single-step method. All thesetheoretical results are derived under the condition that the integration of the com-posite function term is computed exactly. Since, in real problems, it is difficultto get the exact integration value, numerical quadrature is usually employed. Itis, however, reported that instability may occur caused by numerical quadratureerror in [9, 17, 18]. That is, the theoretical stability results may collapse by theintroduction of numerical quadrature.Several methods have been studied to avoid the instability. The map of a particlefrom a time to the previous time along the trajectory, which is nothing but tosolve a system of ordinary differential equations (ODEs), is simplified in [3, 9, 13].Morton–Priestley–Suli [9] solved the ODEs only at the centroids of the elements,and Priestley [13] did only at the vertices of the elements. The map of the otherpoints is approximated by linear interpolation of those values. It becomes possibleto perform the exact integration of the composite function term with the simplifiedmap. Bermejo–Saavedra [3] used the same simplified map as [13] to employ anumerical quadrature of high accuracy to the composite function term. Tanaka–Suzuki–Tabata [19] approximated the map by a locally linearized velocity and thebackward Euler approximation for the solution of the ODEs in P -element. Theapproximate map makes possible the exact integration of the composite functionterm with the map. Pironneau–Tabata [12] used mass lumping in P -element todevelop a scheme free from quadrature for convection-diffusion problems.In this paper we prove the stability and convergence for the scheme with the sameapproximate map as [19] in P k -element for convection-diffusion problems. Since weneither solve the ODEs nor use numerical quadrature, our scheme can be preciselyimplemented to realize the theoretical results. It is, therefore, a genuinely stableLagrange–Galerkin scheme. Our convergence results are of O (∆ t + h + h k +1 ) in ℓ ∞ ( L )-norm and of O (∆ t + h + h k ) in ℓ ∞ ( H )-norm. They are best possible inboth norms for P -element and in ℓ ∞ ( H )-norm for P -elementThe contents of this paper are as follows. In the next section we describe theconvection-diffusion problem and some preparation. In section 3, after recalling theconventional Lagrange–Galerkin scheme, we present our genuinely stable Lagrange–Galerkin scheme. In section 4 we show stability and convergence results, which areproved in section 5. In section 6 we show some numerical results, which reflect thetheoretical convergence order. In section 7 we give conclusions.2. Preliminaries
We state the problem and prepare notation used throughout this paper.Let Ω be a polygonal or polyhedral domain of R d ( d = 2 ,
3) and
T > L p (Ω) with the norm k·k ,p , W s,p (Ω) and W s,p (Ω) withthe norm k·k s,p and the semi-norm |·| s,p for 1 ≤ p ≤ ∞ and a positive integer s .We will write H s (Ω) = W s, (Ω) and drop the subscript p = 2 in the corresponding orms. The L -norm k·k is simply denoted by k·k . The dual space of H (Ω) isdenoted by H − (Ω). For the vector-valued function w ∈ W , ∞ (Ω) d we define thesemi-norm | w | , ∞ by (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:26) d X i,j =1 (cid:16) ∂w i ∂x j (cid:17) (cid:27) / (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) , ∞ . The parenthesis ( · , · ) shows the L -inner product ( f, g ) ≡ R Ω f g dx . For a Sobolevspace X (Ω) we use abbreviations H m ( X ) = H m (0 , T ; X (Ω)) and C ( X ) = C ([0 , T ]; X (Ω)).We define a function space Z m ( t , t ) by Z m ( t , t ) ≡ { f ∈ H j ( t , t ; H m − j (Ω)); j = 0 , . . . , m, k f k Z m ( t ,t ) < ∞} , k f k Z m ( t ,t ) ≡ (cid:26) m X j =0 k f k H j ( t ,t ; H m − j ) (cid:27) / and denote Z m (0 , T ) by Z m .We consider the convection-diffusion problem: find φ : Ω × (0 , T ) → R such that ∂φ∂t + u · ∇ φ − ν ∆ φ = f, ( x, t ) ∈ Ω × (0 , T ) , (1a) φ = 0 , ( x, t ) ∈ ∂ Ω × (0 , T ) , (1b) φ = φ , x ∈ Ω , t = 0 , (1c)where ∂ Ω is the boundary of Ω and ν > ν . Functions u : Ω × (0 , T ) → R d , f ∈ C ( L ) and φ ∈ C ( ¯Ω)are given. Remark . As usual, in place of (1b), we can deal with the inhomogeneous boundarycondition φ = g by replacing the unknown function φ by e φ ≡ φ − e g if the function g defined on ∂ Ω × (0 , T ) can be extended to a function e g in Ω × (0 , T ) appropriately.Let ∆ t > N T ≡ ⌊ T / ∆ t ⌋ , t n ≡ n ∆ t and ψ n ≡ ψ ( · , t n ) fora function ψ defined in Ω × (0 , T ). For a set of functions ψ = { ψ n } N T n =0 , two norms k·k ℓ ∞ ( L ) and k·k ℓ ( n ,n ; L ) are defined by k ψ k ℓ ∞ ( L ) ≡ max {k ψ n k ; n = 0 , . . . , N T } , k ψ k ℓ ( n ,n ; L ) ≡ ∆ t n X n = n k ψ n k ! / and denote k ψ k ℓ (1 ,N T ; L ) by k ψ k ℓ ( L ) .Let u be smooth. The characteristic curve X ( t ; x, s ) is defined by the solutionof the system of the ordinary differential equations, dXdt ( t ; x, s ) = u ( X ( t ; x, s ) , t ) , t < s, (2a) X ( s ; x, s ) = x. (2b)Then, we can write the material derivative term ∂φ∂t + u · ∇ φ as (cid:18) ∂φ∂t + u · ∇ φ (cid:19) ( X ( t ) , t ) = ddt φ ( X ( t ) , t ) . or w : Ω → R d we define the mapping X ( w ) : Ω → R d by( X ( w ))( x ) ≡ x − w ( x )∆ t. (3) Remark . The image of x by X ( u ( · , t )) is nothing but the backward Euler ap-proximation of X ( t − ∆ t ; x, t ).The symbol ◦ stands for the composition of functions, e.g., ( g ◦ f )( x ) ≡ g ( f ( x )).Let T h ≡ { K } be a triangulation of ¯Ω and h ≡ max K ∈T h diam( K ) be the maxi-mum element size. Throughout this paper we consider a regular family of triangu-lations {T h } h ↓ . Let k be a fixed positive integer and V h ⊂ H (Ω) be the P k -finiteelement space, V h ≡ { v h ∈ C ( ¯Ω) ∩ H (Ω); v h | K ∈ P k ( K ) , ∀ K ∈ T h } , where P k ( K ) is the set of polynomials on K whose degrees are less than or equalto k . Let b φ h ∈ V h be the Poisson projection of φ ∈ H (Ω) defined by( ∇ ( b φ h − φ ) , ∇ ψ h ) = 0 , ∀ ψ h ∈ V h . (4)We use c to represent a generic positive constant independent of h , ∆ t , ν , f and φ which may take different values at different places. The notation c ( A ) means that c depends on a positive parameter A and that c increases monotonically when A in-creases. The constants c , c and c stand for c = c ( k u k C ( L ∞ ) ), c = c ( k u k C ( W , ∞ ) )and c = c ( k u k C ( W , ∞ ) ). We also use fixed positive constants α ∗ and δ ∗ defined inLemma 1 in the next section and in Lemma 5 in Section 5, respectively.3. A genuinely stable Lagrange–Galerkin scheme
The conventional Lagrange–Galerkin scheme, which we call Scheme LG, is de-scribed as follows.
Scheme LG.
Let φ h = b φ h . Find { φ nh } N T n =1 ⊂ V h such that for n = 1 , . . . , N T (cid:18) φ nh − φ n − h ◦ X n ∆ t , ψ h (cid:19) + ν ( ∇ φ nh , ∇ ψ h ) = ( f n , ψ h ) , ∀ ψ h ∈ V h , (5)where X n = X ( u n ).For this scheme error estimates k φ h − φ k ℓ ∞ ( L ) ≤ c ( h k + ∆ t ) , c (1 /ν )( h k +1 + ∆ t ) , k φ h − φ k ℓ ∞ ( H ) ≤ c (1 /ν )( h k + ∆ t ) (6)are proved in [6], where the composite function term ( φ n − h ◦ X n , ψ h ) is assumed tobe exactly integrated.Although the function φ n − h is a polynomial on each element K , the compositefunction φ n − h ◦ X n is not a polynomial on K in general since the image X n ( K )of an element K may spread over plural elements. Hence, it is hard to calculatethe composite function term ( φ n − h ◦ X n , ψ h ) exactly. In practice, the followingnumerical quadrature has been used. Let g : K → R be a continuous function. Anumerical quadrature I h [ g ; K ] of R K g dx is defined by I h [ g ; K ] ≡ meas( K ) N q X i =1 w i g ( a i ) , (7) here N q is the number of quadrature points and ( w i , a i ) ∈ R × K is a pair ofweight and point for i = 1 , . . . , N q . We call the practical scheme using numericalquadrature Scheme LG ′ . Scheme LG ′ . Let φ h = b φ h . Find { φ nh } N T n =1 ⊂ V h such that for n = 1 , . . . , N T t ( φ nh , ψ h ) − t X K ∈T h I h [( φ n − h ◦ X n ) ψ h ; K ] + ν ( ∇ φ nh , ∇ ψ h )= ( f n , ψ h ) , ∀ ψ h ∈ V h , (8)where X n = X ( u n ).It is reported that the results (6) do not hold for Scheme LG ′ [9, 17, 18, 19].We denote by Π (1) h the Lagrange interpolation operator to the P -finite elementspace. The following lemma is well-known [5]. Lemma 1. (i)
There exists a positive constant c Π such that for w ∈ W , ∞ (Ω) d k Π (1) h w − w k , ∞ ≤ c Π h | w | , ∞ . (ii) There exists a positive constant α ∗ ≥ such that for w ∈ W , ∞ (Ω) d | Π (1) h w | , ∞ ≤ α ∗ | w | , ∞ . We now present our genuinely stable scheme GSLG, which is free from quad-rature and exactly computable. We define a locally linearized velocity u h and amapping X n h by u h ≡ Π (1) h u, X n h ≡ X ( u nh ) . Scheme GSLG.
Let φ h = b φ h . Find { φ nh } N T n =1 ⊂ V h such that for n = 1 , . . . , N T (cid:18) φ nh − φ n − h ◦ X n h ∆ t , ψ h (cid:19) + ν ( ∇ φ nh , ∇ ψ h ) = ( f n , ψ h ) , ∀ ψ h ∈ V h . (9)We show that the integration ( φ n − h ◦ X n h , ψ h ) can be calculated exactly. At Figure 1.
Elements K , K and a polygon E first we prepare two lemmas. The next lemma on the mapping (3) is proved in [14]. Lemma 2 ([14, Proposition 1]) . Suppose w ∈ W , ∞ (Ω) d and ∆ t | w | , ∞ < . (10) Let F ≡ X ( w ) be the mapping defined in (3). Then, F : Ω → Ω is bijective. Lemma 3.
Let K , K ∈ T h and F : K → R d be linear and one-to-one. Let E ≡ K ∩ F − ( K ) and meas( E ) > . Then, the following hold. i) E is a polygon ( d = 2 ) or a polyhedron ( d = 3 ). (ii) φ h ◦ F | E ∈ P k ( E ) , ∀ φ h ∈ P k ( K ) .Proof. (i) Since both K and F − ( K ) are triangles ( d = 2) or tetrahedra ( d = 3),the intersection is a polygon or a polyhedron. See Fig. 1.(ii) F ∈ P ( K ) d implies that F ∈ P ( E ) d and it holds that F ( E ) ⊂ K .Hence, φ h ◦ F | E is well defined and φ h ◦ F | E ∈ P k ( E ). (cid:3) Proposition 1.
Let φ h , ψ h ∈ V h , w ∈ W , ∞ (Ω) and X h ≡ X (Π (1) h w ) , where X is the operator defined in (3). Suppose α ∗ ∆ t | w | , ∞ < . Then, R Ω ( φ h ◦ X h ) ψ h dx is exactly computable.Proof. It is sufficient to show that R K ( φ h ◦ X h ) ψ h dx can be computable exactlyfor any K ∈ T h . The mapping X h : Ω → Ω is bijective since we can apply Lemma2 thanks to ∆ t | Π (1) h w | , ∞ ≤ α ∗ ∆ t | w | , ∞ < . (11)Let Λ( K ) ≡ (cid:8) l ; K ∩ X − h ( K l ) = ∅ (cid:9) and E l ≡ K ∩ X − h ( K l ) for l ∈ Λ( K ). Notingthat [ l ∈ Λ( K ) E l = K ∩ [ l ∈ Λ( K ) X − h ( K l ) = K and that meas( E l ∩ E m ) = 0 for l = m , l, m ∈ Λ( K ), we can divide the integrationon K into the sum of those on E l for l ∈ Λ( K ), Z K ( φ h ◦ X h ) ψ h dx = X l ∈ Λ( K ) Z E l ( φ h ◦ X h ) ψ h dx. Since Lemma 3 with F = X h implies that both φ h ◦ X h and ψ h are polynomialson E l , we can execute the exact integration. (cid:3) Remark . In the case of d = 2, Priestley [13] approximated X ( t n − ; x, t n ) by e X h ( x ) = B λ ( x ) + B λ ( x ) + B λ ( x ) , x ∈ K on each K ∈ T h , where B i = X ( t n − ; A i , t n ), { A i } i =1 are vertices of K and { λ i } i =1 are the barycentric coordinates of K with respect to { A i } i =1 . Since e X h ( x )is linear in K , the decomposition Z K ( φ h ◦ e X h ) ψ h dx = X l ∈ Λ( K ) Z E l ( φ h ◦ e X h ) ψ h dx, Λ( K ) ≡ n l ; K ∩ e X − h ( K l ) = ∅ o , E l ≡ K ∩ e X − h ( K l )makes the exact integration possible. However, B i = X ( t n − ; A i , t n ) are the solu-tions of a system of ordinary differential equations and they cannot be solved ex-actly in general. In practice, some numerical method, e.g., Runge–Kutta method,is required, which introduces another error.4. Main results
We show the main results, the stability and convergence of Scheme GSLG.
Hypothesis 1. (i) u ∈ C (( W , ∞ ) d ), (ii) u ∈ C (( W , ∞ ∩ W , ∞ ) d ). Hypothesis 2. φ ∈ H ( H k +1 ) ∩ Z . ypothesis 3. The time increment ∆ t satisfies 0 < ∆ t ≤ ∆ t , where∆ t ≡ δ ∗ α ∗ | u | C ( W , ∞ ) , (12)and α ∗ and δ ∗ are the constants stated in Lemma 1 (Section 3) and Lemma 5(Section 5), respectively. Hypothesis 4.
There exists a positive constant c P such that, for ψ ∈ H k +1 (Ω) ∩ H (Ω), k b ψ h − ψ k ≤ c P h k +1 k ψ k k +1 , (13)where b ψ h is the Poisson projection defined in (4). Remark . (i) It is well-known that the H -estimate k b ψ h − ψ k ≤ c P h k k ψ k k +1 (14)holds without any specific condition. On the other hand, Hypothesis 4 holds, forexample, if Ω is convex, by Aubin–Nitsche lemma [5].(ii) Hypothesis 2 implies φ ∈ C ( H k +1 ) and φ ∈ H k +1 (Ω). Theorem 1.
Suppose Hypotheses 1-(i) and 3. Then, there exists a positive constant c ∗ independent of h, ∆ t, ν, φ and f such that k φ h k ℓ ∞ ( L ) + √ ν k∇ φ h k ℓ ( L ) ≤ c ∗ (cid:16)(cid:13)(cid:13) φ h (cid:13)(cid:13) + k f k ℓ ( L ) (cid:17) . Theorem 2.
Suppose Hypotheses 1-(ii), 2 and 3.(i) There exists a positive constant c ∗ independent of h, ∆ t, ν and φ such that k φ − φ h k ℓ ∞ ( L ) + √ ν k∇ ( φ − φ h ) k ℓ ( L ) ≤ c ∗ (cid:26) ∆ t k φ k Z + h k (cid:16)(cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( H k +1 ) + k φ k ℓ ∞ ( H k +1 ) + k φ k ℓ (0 ,N T ; H k +1 ) (cid:17) + h k∇ φ k ℓ (0 ,N T − L ) (cid:27) . (15) (ii) There exists a positive constant c ∗∗ independent of h, ∆ t, φ (but dependent on /ν ) such that k φ − φ h k ℓ ∞ ( H ) ≤ c ∗∗ (cid:26) ∆ t k φ k Z + h k (cid:18)(cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( H k +1 ) + k φ k ℓ ∞ ( H k +1 ) + k φ k ℓ (0 ,N T − H k +1 ) (cid:19) + h k∇ φ k ℓ (0 ,N T − L ) (cid:27) . (16) (iii) Moreover, suppose Hypothesis 4. Then, there exists a positive constant c ∗∗∗ independent of h, ∆ t, ν and φ such that k φ − φ h k ℓ ∞ ( L ) ≤ c ∗∗∗ (cid:26) ∆ t k φ k Z + h k +1 (cid:16)(cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( H k +1 ) + k φ k ℓ ∞ ( H k +1 ) + ν − / k φ k ℓ (0 ,N T − H k +1 ) (cid:17) + h k∇ φ k ℓ (0 ,N T − L ) (cid:27) . (17) emark . From Theorem 2, we have k φ − φ h k ℓ ∞ ( L ) ≤ c (∆ t + h + h k ) , c (cid:18) ∆ t + h + 1 √ ν h k +1 (cid:19) k φ − φ h k ℓ ∞ ( H ) ≤ c (cid:18) ν (cid:19) (∆ t + h + h k ) . In the case of P k -element, k = 1 ,
2, the estimate (15) shows the optimal L -convergence rate O (∆ t + h k ) independent of ν . The dependency on ν in (16)and (17) is also inevitable in Scheme LG.5. Proofs of main theorems
We recall some results used in proving main theorems. For their proofs we onlyshow outlines or refer to the bibliography.
Lemma 4 ([14, Lemma 1]) . Suppose w ∈ W , ∞ (Ω) d and ∆ t | w | , ∞ < . (18) Let F ≡ X ( w ) be the mapping defined in (3). Then, there exists a positive constant c ( | w | , ∞ ) such that for ψ ∈ L (Ω) k ψ ◦ F k ≤ (1 + c ∆ t ) k ψ k . The proof is given in [14].
Lemma 5.
There exists a constant δ ∗ ∈ (0 , such that, for w ∈ W , ∞ (Ω) d and ∆ t satisfying ∆ t | w | , ∞ ≤ δ ∗ , ≤ (cid:12)(cid:12)(cid:12)(cid:12) ∂X ( w ) ∂x (cid:12)(cid:12)(cid:12)(cid:12) ≤ , where | ∂X ( w ) /∂x | is the Jacobian of the mapping X ( w ) defined in (3). Lemma 5 is easily proved by the fact, (cid:18) ∂X ( w ) ∂x (cid:19) ij = δ ij − ∆ t ∂w i ∂x j . Lemma 6.
Let w i ∈ W , ∞ (Ω) d and F i ≡ X ( w i ) be the mapping defined in (3)for i = 1 , . Under the condition ∆ t | w i | , ∞ ≤ δ ∗ , i = 1 , , we have for ψ ∈ H (Ω) k ψ ◦ F − ψ ◦ F k ≤ √ t k w − w k , ∞ k∇ ψ k . Lemma 6 is a direct consequence of [1, Lemma 4.5] and Lemma 5.
Lemma 7.
Let w ∈ W , ∞ (Ω) d and F ≡ X ( w ) be the mapping defined in (3).Under the condition ∆ t | w | , ∞ ≤ δ ∗ , there exists a positive constant c ( k w k , ∞ ) such that for ψ ∈ L (Ω) k ψ − ψ ◦ F k H − (Ω) ≤ c ∆ t k ψ k . Lemma 7 is obtained from [6, Lemma 1] and Lemma 5. emma 8 (discrete Gronwall inequality) . Let a and a be non-negative numbers, ∆ t ∈ (0 , a ] be a real number, and { x n } n ≥ , { y n } n ≥ and { b n } n ≥ be non-negativesequences. Suppose x n − x n − ∆ t + y n ≤ a x n + a x n − + b n , ∀ n ≥ . Then, it holds that x n + ∆ t n X i =1 y i ≤ exp { (2 a + a ) n ∆ t } x + ∆ t n X i =1 b i ! , ∀ n ≥ . Lemma 8 is shown by using the inequalities11 − a ∆ t ≤ a ∆ t ≤ exp(2 a ∆ t ) . Outline of the proof of Theorem 1 . We substitute φ nh into ψ h in (9). We can applyLemma 4 with w = u nh and ψ = φ n − h by virtue of ∆ t | u h | C ( W , ∞ ) <
1. The rest ofthe proof is similar to [14, Theorem 1]. We, therefore, omit it.
Proof of Theorem 2.
We first show the estimate (15). Let e h ≡ φ h − b φ h , η ≡ φ − b φ h , (19)where b φ h is the Poisson projection defined in (4). From (1) and (9) we have (cid:18) e nh − e n − h ◦ X n h ∆ t , ψ h (cid:19) + ν ( ∇ e nh , ∇ ψ h ) = X i =1 ( R ni , ψ h ) (20)for ψ h ∈ V h , where R n ≡ ∂φ n ∂t + u n · ∇ φ n − φ n − φ n − ◦ X n ∆ t ,R n ≡ φ n − ◦ X n h − φ n − ◦ X n ∆ t ,R n ≡ η n − η n − ∆ t , R n ≡ η n − − η n − ◦ X n h ∆ t . (21)Substituting e nh into ψ h , applying Lemma 4 with F = X n h and ψ = e n − h , andevaluating the first term of the left-hand side as (cid:18) e nh − e n − h ◦ X n h ∆ t , e nh (cid:19) ≥ t ( k e nh k − (cid:13)(cid:13) e n − h ◦ X n h (cid:13)(cid:13) ) ≥ t ( k e nh k − (1 + c ∆ t ) (cid:13)(cid:13) e n − h (cid:13)(cid:13) )= 12∆ t ( k e nh k − (cid:13)(cid:13) e n − h (cid:13)(cid:13) ) − c c ∆ t ) (cid:13)(cid:13) e n − h (cid:13)(cid:13) , we have 12∆ t ( k e nh k − (cid:13)(cid:13) e n − h (cid:13)(cid:13) ) + ν k∇ e nh k ≤ c (cid:13)(cid:13) e n − h (cid:13)(cid:13) + X i =1 ε i k R ni k + X i =1 ε i ! k e nh k , (22)where { ε i } i =1 are positive constants satisfying ∆ t ≤ ε , ε ≡ P i =1 ε i . e evaluate R i , i = 1 , · · · ,
4. Setting y ( x, s ) = x + ( s − t u n ( x ) , t ( s ) = t n − + s ∆ t, we have φ n − φ n − ◦ X n ∆ t = 1∆ t (cid:2) φ ( y ( · , s ) , t ( s )) (cid:3) s =0 , which implies R n = ∂φ n ∂t + u n · ∇ φ n − Z (cid:26) u n ( · ) · ∇ φ + ∂φ∂t (cid:27) ( y ( · , s ) , t ( s )) ds = ∆ t Z ds Z s ((cid:18) u n ( · ) · ∇ + ∂∂t (cid:19) φ ) ( y ( · , s ) , t ( s )) ds = ∆ t Z s ((cid:18) u n ( · ) · ∇ + ∂∂t (cid:19) φ ) ( y ( · , s ) , t ( s )) ds . Hence, we have k R n k ≤ ∆ t Z s (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)((cid:18) u n ( · ) · ∇ + ∂∂t (cid:19) φ ) ( y ( · , s ) , t ( s )) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ds ≤ c √ ∆ t k φ k Z ( t n − ,t n ) , (23)where we have used the transformation of independent variables from x to y and s to t , and the estimate | ∂x/∂y | ≤ t | u | C ( W , ∞ ) , ∆ t | u h | C ( W , ∞ ) ≤ δ ∗ , and Lemmas 6 and 1 it holds that k R n k ≤ √ k∇ φ n − kk Π (1) h u n − u n k , ∞ ≤ c h (cid:13)(cid:13) ∇ φ n − (cid:13)(cid:13) . (24) R n is evaluated as k R n k = (cid:13)(cid:13)(cid:13)(cid:13)Z ∂η∂t ( · , t ( s )) ds (cid:13)(cid:13)(cid:13)(cid:13) ≤ c P h k √ ∆ t (cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; H k +1 ) , (25)where we have used (14).From ∆ t | u h | C ( W , ∞ ) ≤ δ ∗ and Lemma 6 it holds that k R n k ≤ √ k∇ η n − kk Π (1) h u n k , ∞ ≤ c h k (cid:13)(cid:13) φ n − (cid:13)(cid:13) k +1 . (26)Combining (22)–(26), we have12∆ t (cid:16) k e nh k − (cid:13)(cid:13) e n − h (cid:13)(cid:13) (cid:17) + ν k∇ e nh k ≤ ε k e nh k + c (cid:13)(cid:13) e n − h (cid:13)(cid:13) + c (cid:26) ∆ t k φ k Z ( t n − ,t n ) + h (cid:13)(cid:13) ∇ φ n − (cid:13)(cid:13) + h k ∆ t (cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; H k +1 ) + h k (cid:13)(cid:13) φ n − (cid:13)(cid:13) k +1 (cid:27) . rom Lemma 8 we obtain for n = 1 , . . . , N T k e nh k + 2 ν ∆ t N T X j =1 (cid:13)(cid:13)(cid:13) ∇ e jh (cid:13)(cid:13)(cid:13) ≤ c (cid:18)(cid:13)(cid:13) e h (cid:13)(cid:13) + ∆ t k φ k Z + h k (cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( H k +1 ) + h k ∆ t N T − X j =0 (cid:13)(cid:13) φ j (cid:13)(cid:13) k +1 + h ∆ t N T − X j =0 (cid:13)(cid:13) ∇ φ j (cid:13)(cid:13) (cid:19) , which implies (15) by virtue of e h = 0 and the triangle inequalities, k φ − φ h k ℓ ∞ ( L ) ≤ k e h k ℓ ∞ ( L ) + k η k ℓ ∞ ( L ) ≤ k e h k ℓ ∞ ( L ) + c P h k k φ k ℓ ∞ ( H k +1 ) , (27) k∇ ( φ − φ h ) k ℓ ( L ) ≤ k∇ e h k ℓ ( L ) + k∇ η k ℓ ( L ) ≤ k∇ e h k ℓ ( L ) + c P h k k φ k ℓ ( H k +1 ) . We show the estimate (16). The equation (20) can be rewritten as1∆ t ( e nh − e n − h , ψ h ) + ν ( ∇ e nh , ∇ ψ h ) = X i =1 ( R ni , ψ h ) , where R n ≡ t ( e n − h ◦ X n h − e n − h ) . From Lemma 6 it holds that k R n k ≤ √ k∇ e n − h kk Π (1) h u n k , ∞ ≤ c k∇ e n − h k . Substituting D ∆ t e nh ≡ t ( e nh − e n − h ) into ψ h , and using (23)–(26) for R , . . . , R ,we have (cid:13)(cid:13) D ∆ t e nh (cid:13)(cid:13) + 1∆ t (cid:16) ν k∇ e nh k − ν (cid:13)(cid:13) ∇ e n − h (cid:13)(cid:13) (cid:17) + ν t (cid:13)(cid:13) ∇ ( e nh − e n − h ) (cid:13)(cid:13) ≤ c (cid:26) ∆ t k φ k Z ( t n − ,t n ) + h (cid:13)(cid:13) ∇ φ n − (cid:13)(cid:13) + h k ∆ t (cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; H k +1 ) + h k (cid:13)(cid:13) φ n − (cid:13)(cid:13) k +1 (cid:27) + c ν (cid:16) ν (cid:13)(cid:13) ∇ e n − h (cid:13)(cid:13) (cid:17) + 12 (cid:13)(cid:13) D ∆ t e nh (cid:13)(cid:13) . From Lemma 8 we have for n = 1 , . . . , N T ∆ t N T X j =1 (cid:13)(cid:13) D ∆ t e nh (cid:13)(cid:13) + ν k∇ e nh k ≤ c exp (cid:18) c Tν (cid:19) (cid:18)(cid:13)(cid:13) ∇ e h (cid:13)(cid:13) + ∆ t k φ k Z + h k (cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( H k +1 ) + h k ∆ t N T − X j =0 (cid:13)(cid:13) φ j (cid:13)(cid:13) k +1 + h ∆ t N T − X j =0 (cid:13)(cid:13) ∇ φ j (cid:13)(cid:13) (cid:19) , which implies (16) by virtue of e h = 0, the triangle inequality, k∇ ( φ − φ h ) k ℓ ∞ ( L ) ≤ k∇ e h k ℓ ∞ ( L ) + k∇ η k ℓ ∞ ( L ) ≤ k∇ e h k ℓ ∞ ( L ) + c P h k k φ k ℓ ∞ ( H k +1 ) and the Poincar´e inequality, k v k ≤ c k∇ v k , ∀ v ∈ H (Ω) . (28) ow we show the estimate (17). We return to the error equation (20). Using(13) in place of (14) in the estimate of R n , we can evaluate (25) as k R n k ≤ c P h k +1 √ ∆ t (cid:13)(cid:13)(cid:13)(cid:13) ∂φ∂t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; H k +1 ) . From Lemma 7 we have k R n k H − (Ω) ≤ c (cid:13)(cid:13) η n − (cid:13)(cid:13) ≤ c h k +1 (cid:13)(cid:13) φ n − (cid:13)(cid:13) k +1 . Hence, it holds that( R n , e nh ) ≤ k R n k H − (Ω) k e nh k ≤ c ν h k +1) (cid:13)(cid:13) φ n − (cid:13)(cid:13) k +1 + ν k∇ e nh k , where we have used the Poincar´e inequality (28). Using this inequality instead of ε k R n k + ε k e nh k in (22) and replacing the last term of (27) by c P h k +1 k φ k ℓ ∞ ( H k +1 ) ,we obtain (17). (cid:3) Numerical results
We show numerical results in d = 2. We compare the conventional scheme(Scheme LG ′ ) with the present one (Scheme GSLG). We use FreeFem++ [8] forthe triangulation of the domain. Both P - and P -elements are used. For SchemeLG ′ we use the seven points quadrature formula of degree five [7]. A relative error E X is defined by E X ≡ k Π ( k ) h φ − φ h k ℓ ∞ ( X ) k Π ( k ) h φ k ℓ ∞ ( X ) , where Π ( k ) h is the Lagrange interpolation operator to the P k -finite element spaceand X = L (Ω) or H (Ω). Example 1 (The rotating Gaussian hill [14]) . In (1), Ω is an unit disk, and we set T = 2 π, ν = 10 − , u ( x, t ) ≡ ( − x , x ) , f ≡ , φ ≡ φ e ( · , , where φ e ( x, t ) ≡ σσ + 4 νt exp (cid:26) − (¯ x ( t ) − x ,c ) + (¯ x ( t ) − x ,c ) σ + 4 νt (cid:27) , (¯ x , ¯ x )( t ) ≡ ( x cos t + x sin t, − x sin t + x cos t ) , ( x ,c , x ,c ) ≡ (0 . , , σ ≡ . . In this problem the identity Π (1) h u = u holds. This problem does not satisfyour setting because Ω is not a polygon and u = 0 on ∂ Ω. The function φ e in Fig.2 (left) satisfies (1a) and (1c) but does not satisfy the boundary condition (1b).However, we may apply the schemes and treat φ e as the solution since the value of φ e on ∂ Ω is almost equal to zero, less than 10 − , and we may neglect the effect ofthe boundary value and the term R K ( φ n − h ◦ X n ) ψ h dx and R K ( φ n − h ◦ X n h ) ψ h dx onthe element K touching the boundary.Let N be the division number of the circle. We set h ≡ π/N, N = 32 , ,
128 and256. Figure 2 (right) shows the triangulation of ¯Ω for N = 64. The time increment∆ t is set to be c h and c h for P -element ( c = π ; . , c = π ; . igure 2. The function φ e ( · ,
0) (left) and the triangulation of ¯Ωfor N = 64 (right) in Example 1 c h and c h for P -element ( c = π ; . , c = π ; .
21) so that we canobserve the convergence behavior of O ( h k ) for E H , and O ( h k ) and O ( h k +1 ) for E L when P k -element is employed.In the following figures we use the symbols shown in Table 1. Figure 3 showsthe log-log graphs of E L and E H versus h . The left graph shows the results of P -element and Tables 2 shows the values of them. The convergence order of E L with ∆ t = O ( h ) is less than 1 in Scheme LG ′ ( ) and more than 1 in SchemeGSLG ( ). The orders of E L with ∆ t = O ( h ) are almost 2 for small h in bothschemes ( , ). The convergence of E H is not observed in Scheme LG ′ ( ) whilethe order is almost 1 in Scheme GSLG ( ). The right graph of Fig. 3 shows theresults of P -element and Tables 3 shows the values of them. The errors E L with∆ t = O ( h ) are too large at N = 128 and 256 to be plotted in the graph in SchemeLG ′ ( ) while the convergence order is almost 2 in Scheme GSLG ( ). The error E L with ∆ t = O ( h ) is large at N = 128, but it becomes small again at N = 256in Scheme LG ′ ( ). We will discuss the reason why such a behavior occurs in aforthcoming paper. The order is greater than 2.5 in Scheme GSLG ( ). The errors E H are too large at N = 128 and 256 to be plotted in the graph in Scheme LG ′ ( )while we can observe the convergence of E H but the order is less than 2 in SchemeGSLG ( ). The errors of Scheme GSLG are smaller than those of Scheme LG ′ inboth cases of P - and P -element. Table 1.
Symbols used in Figs. 3, 6 and 7 and Tables 2–7
X ℓ ∞ ( L ) ℓ ∞ ( L ) ℓ ∞ ( H )∆ t O ( h k ) O ( h k +1 ) O ( h k )Scheme LG ′ Scheme GSLGFigure 4 shows the solutions φ nh for h = 2 π/ ; . , ∆ t = 0 . n ∆ t ; π .In the case of P -element, the solution of Scheme LG ′ is oscillatory while that ofScheme GSLG is much better though a small ruggedness is observed. In the caseof P -element, the solution of Scheme LG ′ is quite oscillatory while that of SchemeGSLG is stable. .0250.05 0.1 0.20.0010.010.11. h E L , E H h E L , E H Figure 3.
Graphs of E L and E H versus h in Example 1 by P k -element. k = 1 (left) and k = 2 (right) - - - - - - - - - - Figure 4.
Solutions φ nh ( n ∆ t ; π ) in Example 1 by SchemeLG ′ (top left) and Scheme GSLG(top right) for P -element, and byScheme LG ′ (bottom left) and Scheme GSLG(bottom right) for P -element able 2. The values of errors and orders of the graph in Fig. 3by P -element N order order order32 7.58E-01 7.58E-01 9.52E-0164 5.65E-01 0.42 5.53E-01 0.45 1.04E+00 -0.13128 3.93E-01 0.52 1.87E-01 1.56 9.72E-01 0.10256 2.04E-01 0.95 4.15E-02 2.17 5.54E-01 0.81 N order order order32 7.26E-01 7.26E-01 9.01E-0164 4.28E-01 0.76 4.18E-01 0.80 7.34E-01 0.30128 1.36E-01 1.65 1.04E-01 2.01 3.07E-01 1.26256 5.62E-02 1.27 2.43E-02 2.10 1.45E-01 1.08 Table 3.
The values of errors and orders of the graph in Fig. 3by P -element N order order order32 6.86E-01 6.86E-01 9.22E-0164 4.06E-01 0.76 3.97E-01 0.79 1.31E+00 -0.51128 1.67E+02 -8.68 8.30E-01 -1.06 1.72E+03 -10.36256 1.42E+27 -82.81 3.05E-03 8.09 3.10E+28 -83.90 N order order order32 6.03E-01 6.03E-01 7.38E-0164 2.09E-01 1.53 1.17E-01 2.37 3.33E-01 1.15128 5.48E-02 1.93 1.59E-02 2.88 9.86E-02 1.76256 1.38E-02 1.99 2.66E-03 2.58 3.97E-02 1.31 Example 2.
In (1), Ω is the square (0 , × (0 , T = 1 , ν = 10 − and10 − , u ( x, t ) ≡ (sin πx sin πx , sin πx sin πx ) , f ≡ ∂φ e ∂t + u · ∇ φ e − ν ∆ φ e ,φ ≡ φ e ( · , , where φ e ( x, t ) ≡ cos(2 πt ) sin ( πx ) sin(2 πx ) . In this problem, Π (1) h u = u . Let N be the division number of each side of ¯Ω.We set h ≡ /N, N = 8 , ,
32 and 64. Figure 5 shows the triangulation of ¯Ωfor N = 16. The time increment ∆ t is set to be c h and c h for P -element( c = 0 . , c = 1), c h and c h for P -element ( c = 1 , c = 5 .
12) so that wecan observe the convergence behavior of O ( h k ) for E H , and O ( h k ) and O ( h k +1 )for E L when P k -element is employed.Figure 6 shows the log-log graphs of E L and E H versus h with ν = 10 − . Theleft graph shows the results of P -element and Table 4 shows the values of them.The convergence orders of E L with ∆ t = O ( h ) are almost 1 in both schemes ( ,). The orders of E L with ∆ t = O ( h ) are almost 2 in both schemes ( , ). Theorders of E H are almost 1 in both schemes ( , ). The right graph of Fig. 6 shows igure 5. The triangulation of ¯Ω for N = 16 in Example 2
164 132 116 18 h E L , E H
164 132 116 18 h E L , E H Figure 6.
Graphs of E L and E H versus h in Example 2 with ν = 10 − by P k -element. k = 1 (left) and k = 2 (right)
164 132 116 18 h E L , E H
164 132 116 18 h E L , E H Figure 7.
Graphs of E L and E H versus h in Example 2 with ν = 10 − by P k -element. k = 1 (left) and k = 2 (right) able 4. The values of errors and orders of the graph in Fig. 6by P -element N order order order8 8.14E-02 8.14E-02 1.10E-0116 3.64E-02 1.16 1.90E-02 2.10 4.36E-02 1.3432 1.70E-02 1.10 4.58E-03 2.05 1.87E-02 1.2264 8.53E-03 0.99 1.19E-03 1.94 9.18E-03 1.03 N order order order8 8.97E-02 8.97E-02 1.09E-0116 3.68E-02 1.29 2.13E-02 2.07 4.23E-02 1.3732 1.78E-02 1.05 4.83E-03 2.14 1.92E-02 1.1464 8.90E-03 1.00 1.29E-03 1.90 9.43E-03 1.03 Table 5.
The values of errors and orders of the graph in Fig. 6by P -element N order order order8 7.01E-02 4.50E-02 7.37E-0216 1.77E-02 1.99 5.72E-03 2.98 1.85E-02 1.9932 4.44E-03 2.00 7.11E-04 3.01 4.63E-03 2.0064 1.11E-03 2.00 8.86E-05 3.00 1.15E-03 2.01 N order order order8 6.31E-02 3.87E-02 6.60E-0216 1.58E-02 2.00 6.89E-03 2.49 1.64E-02 2.0132 3.98E-03 1.99 1.42E-03 2.28 4.12E-03 1.9964 9.90E-04 2.01 3.37E-04 2.08 1.02E-03 2.01 Table 6.
The values of errors and orders of the graph in Fig. 7by P -element N order order order8 9.53E-02 9.53E-02 2.90E-0116 3.94E-02 1.27 4.17E-02 1.19 2.00E-01 0.5432 1.82E-02 1.11 1.33E-02 1.65 1.09E-01 0.8864 9.07E-03 1.00 3.44E-03 1.95 5.12E-02 1.09 N order order order8 9.70E-02 9.70E-02 2.23E-0116 3.93E-02 1.30 2.86E-02 1.76 1.23E-01 0.8632 1.89E-02 1.06 7.42E-03 1.95 6.02E-02 1.0364 9.45E-03 1.00 1.91E-03 1.96 3.08E-02 0.97the results of P -element and Table 5 shows the values of them. The convergenceorders of E L with ∆ t = O ( h ) are almost 2 in both schemes ( , ). The orderof E L with ∆ t = O ( h ) is almost 3 in Scheme LG ′ ( ) and almost 2 in SchemeGSLG ( ). The orders of E H are almost 2 in both schemes ( , ). These results able 7. The values of errors and orders of the graph in Fig. 7by P -element N order order order8 7.45E-02 4.74E-02 1.88E-0116 1.81E-02 2.04 6.77E-03 2.81 1.25E-01 0.5932 3.94E+00 -7.77 1.16E-03 2.55 1.22E+02 -9.9364 1.10E+00 1.84 1.17E-04 3.31 8.17E+01 0.58 N order order order8 6.76E-02 4.17E-02 1.03E-0116 1.69E-02 2.00 8.80E-03 2.24 4.68E-02 1.1432 4.24E-03 1.99 2.04E-03 2.11 1.67E-02 1.4964 1.05E-03 2.01 4.43E-04 2.20 5.84E-03 1.52are consistent with the theoretical ones of Scheme GSLG, E L = O (∆ t + h + h k +1 )and E H = O (∆ t + h + h k ).Figure 7 shows the log-log graphs of E L and E H versus h with ν = 10 − . Theleft graph shows the results of P -element and Table 6 shows the values of them.The convergence orders of E L with ∆ t = O ( h ) are almost 1 in both schemes ( ,). The orders of E L with ∆ t = O ( h ) are almost 2 for small h in both schemes( , ). The orders of E H are almost 1 in both schemes ( , ). The right graphof Fig. 7 shows the results of P -element and Table 7 shows the values of them.The errors E L with ∆ t = O ( h ) are too large at N = 32 and 64 to be plotted inthe graph in Scheme LG ′ ( ) while the convergence order is almost 2 in SchemeGSLG ( ). The order E L with ∆ t = O ( h ) is almost 3 for small h in Scheme LG ′ ( ) and almost 2 in Scheme GSLG ( ). The errors E H are too large at N = 32and 64 to be plotted in the graph in Scheme LG ′ ( ) while we can observe theconvergence but the order is less than 2 in Scheme GSLG ( ). In order to obtainthe theoretical convergence order O ( h ), it seems that finer mesh will be necessary.7. Conclusions
We have presented a genuinely stable Lagrange–Galerkin scheme for convection-diffusion problems. In the scheme locally linearized velocities are used and theintegration is executed exactly without numerical quadrature. For the P k -elementwe have shown error estimates of O (∆ t + h + h k +1 ) in ℓ ∞ ( L )-norm and of O (∆ t + h + h k ) in ℓ ∞ ( H )-norm. We have also obtained error estimate, c (∆ t + h + h k )in ℓ ∞ ( L )-norm, where the coefficient c is dependent on the exact solution φ butindependent of the diffusion constant ν . Numerical results have reflected theseestimates. The extension to the Navier–Stokes equations will be discussed in aforthcoming paper. Acknowledgements.
The first author was supported by JSPS (Japan Societyfor the Promotion of Science) under Grants-in-Aid for Scientific Research (C)No.25400212 and (S)No. 24224004 and under the Japanese-German Graduate Ex-ternship (Mathematical Fluid Dynamics) and by Waseda University under Projectresearch, Spectral analysis and its application to the stability theory of the Navier-Stokes equations of Research Institute for Science and Engineering. The secondauthor was supported by JSPS under Grant-in-Aid for JSPS Fellows No. 26 · eferences [1] Achdou, Y., Guermond, J.: Convergence analysis of a finite element projection/Lagrange–Galerkin method for the incompressible Navier–Stokes equations. SIAM Journal on NumericalAnalysis (3), 799–826 (2000)[2] Barrett, R., Berry, M., Chan, T.F., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo,R., Romine, C., Van der Vorst, H.: Templates for the Solution of Linear Systems: BuildingBlocks for Iterative Methods, 2nd Edition. SIAM, Philadelphia, PA (1994)[3] Bermejo, R., Saavedra, L.: Modified Lagrange–Galerkin methods of first and second order intime for convection–diffusion problems. Numerische Mathematik (4), 601–638 (2012)[4] Boukir, K., Maday, Y., M´etivet, B., Razafindrakoto, E.: A high-order characteristics/finiteelement method for the incompressible Navier-Stokes equations. International Journal forNumerical Methods in Fluids (12), 1421–1454 (1997)[5] Ciarlet, P.G.: The Finite Element Method for Elliptic Problems. Classics in Applied Mathe-matics. SIAM (2002)[6] Douglas Jr., J., Russell, T.: Numerical methods for convection-dominated diffusion problemsbased on combining the method of characteristics with finite element or finite differenceprocedures. SIAM Journal on Numerical Analysis (5), 871–885 (1982)[7] Hammer, P.C., Marlowe, O.J., Stroud, A.H.: Numerical integration over simplexes and cones.Mathematics of Computation , 130–137 (1956)[8] Hecht, F.: New development in FreeFem++. Journal of Numerical Mathematics (3-4),251–265 (2012)[9] Morton, K.W., Priestley, A., Suli, E.: Stability of the Lagrange-Galerkin method with non-exact integration. Mod´elisation math´ematique et analyse num´erique (4), 625–653 (1988)[10] Notsu, H., Tabata, M.: Error estimates of a pressure-stabilized characteristics finite el-ement scheme for the Oseen equations. Journal of Scientific Computing (2015). DOI10.1007/s10915-015-9992-8. URL http://dx.doi.org/10.1007/s10915-015-9992-8 [11] Pironneau, O.: On the transport-diffusion algorithm and its applications to the Navier-Stokesequations. Numerische Mathematik , 309–332 (1982)[12] Pironneau, O., Tabata, M.: Stability and convergence of a Galerkin-characteristics finiteelement scheme of lumped mass type. International Journal for Numerical Methods in Fluids (10-12), 1240–1253 (2010)[13] Priestley, A.: Exact projections and the Lagrange-Galerkin method: a realistic alternative toquadrature. Journal of Computational Physics (2), 316–333 (1994)[14] Rui, H., Tabata, M.: A second order characteristic finite element scheme for convection-diffusion problems. Numerische Mathematik , 161–177 (2002)[15] Saad, Y.: Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia (2003)[16] S¨uli, E.: Convergence and nonlinear stability of the Lagrange-Galerkin method for the Navier-Stokes equations. Numerische Mathematik (4), 459–483 (1988)[17] Tabata, M.: Discrepancy between theory and real computation on the stability of some finiteelement schemes. Journal of computational and applied mathematics (2), 424–431 (2007)[18] Tabata, M., Fujima, S.: Robustness of a characteristic finite element scheme of second orderin time increment. In: Computational Fluid Dynamics 2004, pp. 177–182. Springer (2006)[19] Tanaka, K., Suzuki, A., Tabata, M.: A characteristic finite element method using the ex-act integration. Annual Report of Research Institute for Information Technology, KyushuUniversity , 11–18 (2002). (in Japanese), 11–18 (2002). (in Japanese)