A Continuous Space-Time Finite Element Method for the Wave Equation
A Continuous Space-Time Finite Element Method for the Wave Equation Donald A. French
Department of Mathematical Sciences (ML 25), University of Cincinnati, Cincinnati, OH 45221
Todd E. Peterson
Department of Applied Mathematics, University of Virginia, Charlottesville, VA 22903
October 8, 1994 Abstract
We consider a nite element method for the nonhomogeneous second order wave equation which is formulated in terms of continuous approximation functions in both space and time, thereby giving a uni ed treatment of the spatial and temporal discretizations. Our analysis uses primarily energy arguments which are quite common for spatial discretizations but not for time. We present a priori estimates nodal (in time) superconvergence error estimates without any special time step restrictions. Our method is based on tensor-product spaces for the full discretization.
1 Introduction The continuous time Galerkin (CTG ) method is a nite element technique which provides time discretizations for evolution problems using approximation spaces of continuous functions. This approach is particularly appropriate for wave problems as it retains discrete versions of the important energy conservation properties provided by the initial/boundary value problem being approximated (see [FS]). Computations and analyses have shown this is especially useful in the approximations of solutions to nonlinear wave problems (see, for instance, [G], [GS], or [SV]). Recent work by DeFrutos and Sanz-Serna [DS] indicates that the constants in long-time estimates may be smaller for such methods. Another advantage of the CTG approach is that CTG methods of any desired order of accuracy are easily formulated.
Research supported in part by the University of Cincinnati through the University Research Council and Taft
Grants-in-aid as well as the Army Research Oce by grant 28535-MA.
1
The main purpose of this paper is to demonstrate new variational techniques to analyze these high order accurate space-time nite element methods. We will prove both global convergence and nodal in time superconvergence error estimates. The global error estimates we present have also been obtained by [BL] (see also [BL2]), however by non-variational arguments, and in earlier work of the authors [FP], but by di erent techniques which required a time step restriction. The approximation of the heat equation by CTG methods was studied by Aziz and Monk [AM]. Our report complements theirs; however, we note that the stability estimates for the wave equation are more complicated (see Section 3), and our proof of superconvergence is shorter and, we feel, more straightforward. The techniques we use would also apply to the heat equation. We remark that these CTG schemes in the homogeneous case (f = 0) are equivalent to GaussLegendre implicit Runge-Kutta methods (see [FS]). In this connection see also [BB]. For three other non-classical nite element treatments of the wave equation, see Babuska and Janik [BJ], Johnson [J], and Richter [R].
2 Preliminaries
Let be a bounded region in Rd with a smooth boundary @ , and let [0; T ] be a nite time interval. We consider the following initial/boundary value problem: nd U = U (x; t) such that Utt ? U = f in [0; T ] U = 0 on @ [0; T ] (1) U (; 0) = U0 and Ut (; 0) = V0 in
Our results easily generalize to the case where ? is replaced by any uniformly elliptic self-adjoint second order operator which is independent of t; the time-dependent case will be the subject of future work. For a domain S Rd , we will use the Lebesgue spaces L2 (S ) and L1 (S ), and the Sobolev spaces H s(S ) for s a positive integer, all de ned in the usual way. We will also use H01( ) and its dual H ?1( ). For H01( ), we take the norm to be jjv jjH01( ) = jjrv jjL2( ) . All of these spaces are Hilbert spaces except for L1 (S ). When S = , we will usually omit from our notation. For functions depending on both space and time variables, given a time interval [a; b] and H any of the above Hilbert spaces, we de ne the Hilbert space L2 ([a; b]; H ) by
jjvjjL2([a;b];H ) =
Zb
a
jjv(; t)jj2H dt
!1=2
:
There is an analogous de nition for L1 ([a; b]; H ). When [a; b] = [0; T ], these will be denoted simply by L2 (H ) and L1 (H ). We use C to denote a generic positive constant, not necessarily the same at di erent occurences, but always independent of all discretization parameters, solutions, and of T . 2
We reformulate (1) as a rst order system by introducing the function V = Ut . Letting ! ! " # U 0 0 ? I Y = V ; F = f ; and A = ? 0 then we have
Yt + AY = F in [0; T ] Y = 0 on @ [0; T ] (2) ! Y (; 0) = UV 0 in
0 where the domain of A is D(A) = (H 2 \ H01) H01. It will also be convenient to de ne a mapping T : H ?1 ! H01 by ?(Tv) = v in ; Tv = 0 on @ : Finally, we de ne
!
B = ?0I T0 ; with D(B ) = L2 H ?1 . Notice that BA and AB are identities on the appropriate domains. We next discuss the approximation spaces and their properties. Let Sph be a nite dimensional subspace of H01, depending on a discretization parameter h > 0. De ne the L2 projection x : L2 ! Sph by (x u; )L2 = (u; )L2 8 2 Sph ; and de ne the H01 projection Px : H01 ! Sph by (rPxu; r )L2 = (ru; r )L2 8 2 Sph : We assume for Sph the following properties:
jju ? xujjL2 Chr jjujjH
(3)
r
where u 2 H r \ H01 and 0 r p + 1; and
(4) jju ? PxujjH Chr?s jjujjH where u 2 H r \ H01, 0 s r p + 1 and s = 0; 1. De ne the discrete Laplace operator ?h : Sph ! Sph by (?h ; )L2 = (r; r )L2 8 2 Sph ; r
s
3
and de ne Th : H ?1 ! Sph by (rTh ; r )L2 = (; )L2 8 2 Sph : Note that Th restricted to Sph is the inverse of ?h . We also de ne the following two operators:
Ah =
0 ?I ?h 0
!
!
0 Th ?I 0 :
and Bh =
Notice that on Sph Sph these are inverses of each other. Let [0; T ] be partitioned by 0 = t0 < t1 < : : : < tN = T , and let
In = [tn?1 ; tn]; kn = tn ? tn?1 ; k = maxfkn : 1 n N g: For functions which depend continuously on time, we will often use the notation n = (tn ). The space of polynomials of degree q on an interval [a; b] is denoted by Pq ([a; b]). We de ne Sqk to be those continuous functions on [0; T ] whose restriction to any In belongs to Pq (In ). De ne an operator t : L2 (In ) ! Pq?1 (In ) by the equation (tu; )L2 (I ) = (u; )L2(I n
n
)
8 2 Pq?1 (In);
and also de ne Pt : H 1([0; T ]) ! Sqk by Pt u(0) = u(0) and (@t (Pt u); t)L2 ([0;T ]) = (ut ; t)L2 ([0;T ]) 8 2 Sqk : Note that
Ptu(tn ) = u(tn ) n = 0; 1; : : :; N and that there is no ambiguity if we talk of Pt : H 1 (In ) ! Pq (In ) (ie, Pt may be computed locally). Also note that t and Pt are projections into di erent spaces. We have the following approximation properties:
jju ? tujjL2(I ) Cknr jj@trujjL2(I ) n
where u 2 H r (In ) and 0 r q ; and
(5)
n
jj@ts(u ? Ptu)jjL2(I ) Cknr?sjj@trujjL2(I ) (6) where u 2 H r (In ) and 0 s r q + 1, s = 0; 1. Any function 2 Pq (In ) satis es the following n
inverse properties:
n
jj jjL1(I ) Ckn?1=2jj jjL2(I ) ; jj tjjL2(I ) Ckn?1jj jjL2(I ) ; n
n
n
n
4
(7) (8)
n
o
jj jjL2(I ) C kn1=2j (tn?1)j + jjt jjL2(I ) : (9) The space-time domains Q = [0; T ] and Sn = In will be used in this paper. Our approximate solutions will be de ned in the space Spqhk = Sph Sqk . The operators and estimates we n
n
have introduced for Sqk and Sph can be extended in obvious ways to the space Spqhk . We now introduce the approximation scheme. The method is based on the formulation (2), so U , Ut are approximated seperately by u; v 2 Spqhk . These approximates are de ned successively on each slab of space-time as follows: (ut ? v; )L2(S ) = 0 8 2 Sph Pq?1 (In ); n
(10)
(vt; )L2(S ) + (ru; r)L2(S ) = (f; )L2(S ) 8 2 Sph Pq?1 (In ): (11) Also take u(; 0) = u0 and v (; 0) = v0 where u0 ; v0 2 Sph are some suitable approximations of U0 ; V0. Note that the test functions ; are one degree lower (q ? 1) in time to account for the fact that u and v are xed a priori by continuity at t = tn?1 . Letting y = (u; v), we can reformulate this as follows: (yt + Ah y; )L2 (I ;H01L2 ) = (F; )L2(I ;H01 L2 ) 8 2 [Sph Pq?1 (In )]2: (12) One of the most appealing properties of this scheme is that it conserves energy in the same way as the continuous problem. Letting = vt in (10) and = ut in (11) we obtain n
n
n
n
n
En = En?1 + (f; ut)S where
n
En = 21 jjvnjj2L2 + 12 jjrunjj2L2 ;
and if f = 0 then the energy is conserved. We nish this section with an outline of the remainder of this paper. In Section 3 we analyse CTG methods for abstract problems in the form (2). In Section 4 we apply these results to the wave equation. Section 5 is devoted to the nal steps of our error estimate, and Section 6 contains some numerical results.
3 CTG approximation of an abstract IVP In this section we consider the discretization in time of an abstract initial value problem. Let H be a real Hilbert space, and let A be an operator de ned on a dense domain D(A) H , which generates a strongly continuous semi-group, which we will denote by etA . We assume that (AV; V ) 0 and that jjAV jjH C jjAV jjH for all V 2 D(A). Then jjetAV jjH jjV jjH for all V 2 H . In particular these assumptions are satis ed if A is skew-symmetric, which is the case for the wave equation; 5
however our analysis is more general, and would apply also for example to the heat equation. We consider the problem Yt + AY = F; Y (0) = Y0 : (13) Precise assumptions on Y0 and F will be stated below. In this section we will denote t and Pt simply by and P , respectively. The time-discrete CTG approximation to (13) is an element y of D(A) Sqk which satis es y(0) = Y0 and for 1 n N , (yt ; )L2(I ;H ) + (Ay; )L2(I ;H ) = (F; )L2 (I ;H ) 8 2 H Pq?1 (In ): We rst derive a basic stability estimate. Theorem 1 If y satis es (14), then n
n
(14)
n
n
o
a) jjytjjL2(H ) +jjAyjjL2(H ) C T 1=2jjAY0 jjH + T jjAF jjL2(H ) + jjF jjL2(H ) ; and for 0 t T ,
n
o
b) jjyt(t)jjH +jjAy(t)jjH C jjAY0 jjH + T 1=2jjAF jjL2(H ) + jjF jjL1 (H ) : proof: On each subinterval, (14) is equivalent to yt = ? Ay + F . Therefore, by (9), we have n
o
jjAyjjL2(I ;H ) C kn1=2jjAyn?1jjH + jjAyjjL2(I ;H ) n o C kn1=2jjAyn?1jjH + jjytjjL2(I ;H ) + jjF jjL2(I ;H ) : n
n
n
(15)
n
Taking = Ayt in (14) gives (yt ; Ayt )L2 (I ;H ) + (Ay; Ayt)L2 (I ;H ) = (F; Ayt )L2 (I ;H ); 1 jjAy jj2 ? 1 jjAy jj ?(AF; y ) t L2 (I ;H ): 2 n H 2 n?1 H Summing over n gives 1 jjAy jj2 1 jjAy jj2 ? (AF; y ) (16) t L2 ([0;t ];H ): 2 n H 2 0 H On In let w(t) = t ? tn?1 and write y = yn?1 + wy with y 2 H Pq?1 (In ). Choosing = y in (14) gives (yt ; y)L2 (I ;H ) + (Ay; y)L2 (I ;H ) = (F; y)L2 (I ;H ); (y ; y)L2 (I ;H ) + (wyt ; y)L2 (I ;H ) + (Ayn?1 ; y)L2 (I ;H ) + (wAy; y)L2 (I ;H ) = (F; y)L2 (I ;H ); jjyjj2L2(I ;H ) + (wyt; y)L2(I ;H ) ?(Ayn?1 ; y)L2(I ;H ) + (F; y)L2(I ;H ): n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
6
n
n
Integration by parts in time establishes that
(wyt; y)L2 (I ;H ) = ? 12 jjyjj2L2 (I ;H ) + 21 kn jjyn jj2H : n
Thus
n
1 jjyjj2 2 L2 (I ;H ) jjAyn?1 jjL2 (I ;H )jjyjjL2 (I ;H ) + jjF jjL2(I ;H )jjyjjL2(I ;H ) = kn1=2jjAyn?1 jjH jjyjjL2 (I ;H ) + jjF jjL2(I ;H )jjyjjL2(I ;H ) n
n
n
n
n
whence
n
n
n
jjyjjL2(I ;H ) 2 kn1=2jjAyn?1 jjH + jjF jjL2(I n
n
n
;H )
o
:
(17)
By the inverse estimate (8) and properties of w, we have jjytjjL2(I ;H ) = jjy + wytjjL2(I ;H ) jjyjjL2(I ;H ) + jjwytjjL2(I ;H ) C jjyjjL2(I ;H ): Equations (15)-(18) combine to give n
n
n
n
(18)
n
n
o
jjytjj2L2(I ;H ) + jjAyjj2L2(I ;H ) C kn jjAyn?1jj2H + jjF jj2L2(I ;H ) n o C kn jjAy0jj2H + knjjAF jjL2([0;t ];H )jjytjjL2([0;t ];H ) + jjF jj2L2(I ;H ) : (19) n
n
n
n
n
n
Summing over n yields jjytjj2L2([0;t ];H ) + jjAyjj2L2([0;t ];H ) n o C tn jjAy0jj2H + tnjjAF jjL2([0;t ];H )jjytjjL2([0;t ];H ) + jjF jj2L2([0;t ];H ) : (20) A simple kickback argument completes the proof of the rst result. To obtain the pointwise in time estimate, by (7) and (19), n
n
n
n
n
n
o
jjytjj2L1(I ;H ) + jjAyjj2L1(I ;H ) Ckn?1 jjytjj2L2(I ;H ) + jjAyjj2L2(I ;H ) n o C jjAy0jj2H + jjAF jjL2([0;t ];H )jjytjjL2([0;t ];H ) + kn?1 jjF jj2L2(I ;H ) n o C jjAy0jj2H + tnjjAF jj2L2([0;t ];H ) + t?n 1jjytjj2L2([0;t ];H ) + kn?1 jjF jj2L2(I ;H ) : n
n
n
n
n
n
n
n
Now by (20) we have
n
n
n
o
jjytjj2L1(I ;H ) + jjAyjj2L1(I ;H ) C jjAy0jj2H + tn jjAF jj2L2([0;t ];H ) + jjF jj2L1([0;t ];H ) : n
n
n
n
The desired result follows. //// If H is nite dimensional, then existence and uniqueness of the CTG approximation follow at once from the preceeding stability estimate. The next theorem shows that this holds true in general. 7
Theorem 2 Given Y0 2 D(A) and F which satis es jjAF jjL2(H ) + jjF jjL2(H ) < 1, there is a unique y 2 D(A) Sqk which satis es (14). proof: Let f n g be an orthonormal basis for H , with n 2 D(A), and set Hn = spanf 1 ; :::; ng. Let y0n be the orthogonal projection of Y0 into Hn . It suces to consider a generic time interval, such as I = [0; 1]. For each n there exists a unique solution to the ( nite-dimensional) problem: nd y n 2 Hn Pq (I ) such that y n (0) = y0n and (ytn + Ay n ; )L2 (I;H ) = (F; )L2 (I;H ) 8 2 Hn Pq?1 (I ): By the previous theorem and the inequality Z
jjyjjL2(I;H ) jjy(0)jjH + jjyt(t)jjH dt jjy(0)jjH + jjytjjL2(I;H ) I
it follows that fjjy n jjL2(I;H ) + jjAy n jjL2 (I;H )g is bounded. Because A is necessarily closed (being the generator of a strongly continuous semi-group), from this we can deduce that there is a subsequence, still denoted y n , such that y n converges weakly in L2 (I; H ) to some y 2 D(A), and further that Ayn converges weakly in L2(I; H ) to Ay. Since (ytn ; )L2 (I;H ) = ?(y n ; t )L2 (I;H ) ! ?(y; t )L2 (I;H ) = (yt ; )L2(I;H ) we also have that ytn converges weakly to yt . To show that y satis es (14), given 2 H Pq?1 (I ), let m be the orthogonal projection of into Hm Pq?1 (I ). Then for n m, (ytn + Ay n ; m )L2 (I;H ) = (F; m )L2 (I;H ) Fix m, and let n ! 1, and then let m ! 1. It only remains to show that the initial condition is satis ed. But this is trivial: by construction y n (0) converges in H to Y0 , and it also easy to deduce that y n (0) converges weakly in H to y (0); it follows that y (0) = Y0 . This proves existence. Uniqueness follows immediately from the stability estimate. //// The previous theorem guarantees that y (t) 2 D(A); standard arguments show that y (t) will have more regularity (i.e., lie in the domain of higher powers of A) under the appropriate assumptions on Y0 and F , and this fact will be tacitly used below. The stability estimate also allows us to prove the following error estimate: Theorem 3 Let Y be the solution of (13), and y the CTG approximation de ned by (14). Then for 0 t T , n
o
jjAy(t) ? AY (t)jjH Ckq+1 T 1=2jj@tq+1A2Y jjL2(H ) + jj@tq+1AY jjL1(H ) : 8
proof: Write y ? Y = (y ? PY ) ? (Y ? PY ) = + . Note that 2 H Sqk and (0) = 0. A short calculation establishes that satis es, for any n and any 2 H Pq?1 (In ),
(t + A; )L2 (I ;H ) = (A; )L2(I ;H ) : The stated estimate follows by applying Theorem 1b to , and estimating using (6). //// Our next goal is to derive a higher order estimate for the error at time nodes t = tn . For this we will need the following stability result: Lemma 1 The solution y of (14) satis es n
n
q X q q +1 1 = 2 2 q +1 2 q +1 jj@t A yjjL2(H ) CT jjA Y0jjH + CT jjA F jjL2(H ) + C jj@tq?j Aq+j F jjL2(H ): j =1
proof: Recall that on each subinterval In , y satisi es yt + Ay = F . Operating on this identity with @ti?1Aj gives
@tiAj y = ?@ti?1 Aj+1 y + @ti?1 (I ? )Aj+1 y + @ti?1 Aj F:
(21)
For the second term on the right hand side, we have by (5) and (8) that jj@ti?1(I ? )Aj+1yjjL2(I ;H ) Ckn?(i?1)jj(I ? )Aj+1yjjL2(I ;H ) n
Ck?(i?1) ki?1 jj@ i?1 Aj+1 yjj n
n
t
n
L2 (I ;H ); n
and thus taking norms in (21) gives
jj@tiAj yjjL2(I ;H ) C jj@ti?1Aj+1yjjL2(I ;H ) + jj@ti?1Aj F jjL2(I ;H ): n
n
n
(22)
By summing over n and repeated use of (22) we obtain
jj@tqAq+1 yjjL2(H ) C jj@tA2q yjjL2(H ) + C
qX ?1
jj@tq?j Aq+j F jjL2(H ):
j =1 The proof is now completed by applying Theorem 1 to A2q y , which is just the CTG solution to the problem with initial data A2q Y0 and nonhomogeneity A2q F . ////
Following is the nal result of this section. Theorem 4 Let Y be the solution of (13), and y the CTG approximation de ned by (14). Then, assuming Y0 and F have the indicated regularity, for 1 n N , 8
Department of Mathematical Sciences (ML 25), University of Cincinnati, Cincinnati, OH 45221
Todd E. Peterson
Department of Applied Mathematics, University of Virginia, Charlottesville, VA 22903
October 8, 1994 Abstract
We consider a nite element method for the nonhomogeneous second order wave equation which is formulated in terms of continuous approximation functions in both space and time, thereby giving a uni ed treatment of the spatial and temporal discretizations. Our analysis uses primarily energy arguments which are quite common for spatial discretizations but not for time. We present a priori estimates nodal (in time) superconvergence error estimates without any special time step restrictions. Our method is based on tensor-product spaces for the full discretization.
1 Introduction The continuous time Galerkin (CTG ) method is a nite element technique which provides time discretizations for evolution problems using approximation spaces of continuous functions. This approach is particularly appropriate for wave problems as it retains discrete versions of the important energy conservation properties provided by the initial/boundary value problem being approximated (see [FS]). Computations and analyses have shown this is especially useful in the approximations of solutions to nonlinear wave problems (see, for instance, [G], [GS], or [SV]). Recent work by DeFrutos and Sanz-Serna [DS] indicates that the constants in long-time estimates may be smaller for such methods. Another advantage of the CTG approach is that CTG methods of any desired order of accuracy are easily formulated.
Research supported in part by the University of Cincinnati through the University Research Council and Taft
Grants-in-aid as well as the Army Research Oce by grant 28535-MA.
1
The main purpose of this paper is to demonstrate new variational techniques to analyze these high order accurate space-time nite element methods. We will prove both global convergence and nodal in time superconvergence error estimates. The global error estimates we present have also been obtained by [BL] (see also [BL2]), however by non-variational arguments, and in earlier work of the authors [FP], but by di erent techniques which required a time step restriction. The approximation of the heat equation by CTG methods was studied by Aziz and Monk [AM]. Our report complements theirs; however, we note that the stability estimates for the wave equation are more complicated (see Section 3), and our proof of superconvergence is shorter and, we feel, more straightforward. The techniques we use would also apply to the heat equation. We remark that these CTG schemes in the homogeneous case (f = 0) are equivalent to GaussLegendre implicit Runge-Kutta methods (see [FS]). In this connection see also [BB]. For three other non-classical nite element treatments of the wave equation, see Babuska and Janik [BJ], Johnson [J], and Richter [R].
2 Preliminaries
Let be a bounded region in Rd with a smooth boundary @ , and let [0; T ] be a nite time interval. We consider the following initial/boundary value problem: nd U = U (x; t) such that Utt ? U = f in [0; T ] U = 0 on @ [0; T ] (1) U (; 0) = U0 and Ut (; 0) = V0 in
Our results easily generalize to the case where ? is replaced by any uniformly elliptic self-adjoint second order operator which is independent of t; the time-dependent case will be the subject of future work. For a domain S Rd , we will use the Lebesgue spaces L2 (S ) and L1 (S ), and the Sobolev spaces H s(S ) for s a positive integer, all de ned in the usual way. We will also use H01( ) and its dual H ?1( ). For H01( ), we take the norm to be jjv jjH01( ) = jjrv jjL2( ) . All of these spaces are Hilbert spaces except for L1 (S ). When S = , we will usually omit from our notation. For functions depending on both space and time variables, given a time interval [a; b] and H any of the above Hilbert spaces, we de ne the Hilbert space L2 ([a; b]; H ) by
jjvjjL2([a;b];H ) =
Zb
a
jjv(; t)jj2H dt
!1=2
:
There is an analogous de nition for L1 ([a; b]; H ). When [a; b] = [0; T ], these will be denoted simply by L2 (H ) and L1 (H ). We use C to denote a generic positive constant, not necessarily the same at di erent occurences, but always independent of all discretization parameters, solutions, and of T . 2
We reformulate (1) as a rst order system by introducing the function V = Ut . Letting ! ! " # U 0 0 ? I Y = V ; F = f ; and A = ? 0 then we have
Yt + AY = F in [0; T ] Y = 0 on @ [0; T ] (2) ! Y (; 0) = UV 0 in
0 where the domain of A is D(A) = (H 2 \ H01) H01. It will also be convenient to de ne a mapping T : H ?1 ! H01 by ?(Tv) = v in ; Tv = 0 on @ : Finally, we de ne
!
B = ?0I T0 ; with D(B ) = L2 H ?1 . Notice that BA and AB are identities on the appropriate domains. We next discuss the approximation spaces and their properties. Let Sph be a nite dimensional subspace of H01, depending on a discretization parameter h > 0. De ne the L2 projection x : L2 ! Sph by (x u; )L2 = (u; )L2 8 2 Sph ; and de ne the H01 projection Px : H01 ! Sph by (rPxu; r )L2 = (ru; r )L2 8 2 Sph : We assume for Sph the following properties:
jju ? xujjL2 Chr jjujjH
(3)
r
where u 2 H r \ H01 and 0 r p + 1; and
(4) jju ? PxujjH Chr?s jjujjH where u 2 H r \ H01, 0 s r p + 1 and s = 0; 1. De ne the discrete Laplace operator ?h : Sph ! Sph by (?h ; )L2 = (r; r )L2 8 2 Sph ; r
s
3
and de ne Th : H ?1 ! Sph by (rTh ; r )L2 = (; )L2 8 2 Sph : Note that Th restricted to Sph is the inverse of ?h . We also de ne the following two operators:
Ah =
0 ?I ?h 0
!
!
0 Th ?I 0 :
and Bh =
Notice that on Sph Sph these are inverses of each other. Let [0; T ] be partitioned by 0 = t0 < t1 < : : : < tN = T , and let
In = [tn?1 ; tn]; kn = tn ? tn?1 ; k = maxfkn : 1 n N g: For functions which depend continuously on time, we will often use the notation n = (tn ). The space of polynomials of degree q on an interval [a; b] is denoted by Pq ([a; b]). We de ne Sqk to be those continuous functions on [0; T ] whose restriction to any In belongs to Pq (In ). De ne an operator t : L2 (In ) ! Pq?1 (In ) by the equation (tu; )L2 (I ) = (u; )L2(I n
n
)
8 2 Pq?1 (In);
and also de ne Pt : H 1([0; T ]) ! Sqk by Pt u(0) = u(0) and (@t (Pt u); t)L2 ([0;T ]) = (ut ; t)L2 ([0;T ]) 8 2 Sqk : Note that
Ptu(tn ) = u(tn ) n = 0; 1; : : :; N and that there is no ambiguity if we talk of Pt : H 1 (In ) ! Pq (In ) (ie, Pt may be computed locally). Also note that t and Pt are projections into di erent spaces. We have the following approximation properties:
jju ? tujjL2(I ) Cknr jj@trujjL2(I ) n
where u 2 H r (In ) and 0 r q ; and
(5)
n
jj@ts(u ? Ptu)jjL2(I ) Cknr?sjj@trujjL2(I ) (6) where u 2 H r (In ) and 0 s r q + 1, s = 0; 1. Any function 2 Pq (In ) satis es the following n
inverse properties:
n
jj jjL1(I ) Ckn?1=2jj jjL2(I ) ; jj tjjL2(I ) Ckn?1jj jjL2(I ) ; n
n
n
n
4
(7) (8)
n
o
jj jjL2(I ) C kn1=2j (tn?1)j + jjt jjL2(I ) : (9) The space-time domains Q = [0; T ] and Sn = In will be used in this paper. Our approximate solutions will be de ned in the space Spqhk = Sph Sqk . The operators and estimates we n
n
have introduced for Sqk and Sph can be extended in obvious ways to the space Spqhk . We now introduce the approximation scheme. The method is based on the formulation (2), so U , Ut are approximated seperately by u; v 2 Spqhk . These approximates are de ned successively on each slab of space-time as follows: (ut ? v; )L2(S ) = 0 8 2 Sph Pq?1 (In ); n
(10)
(vt; )L2(S ) + (ru; r)L2(S ) = (f; )L2(S ) 8 2 Sph Pq?1 (In ): (11) Also take u(; 0) = u0 and v (; 0) = v0 where u0 ; v0 2 Sph are some suitable approximations of U0 ; V0. Note that the test functions ; are one degree lower (q ? 1) in time to account for the fact that u and v are xed a priori by continuity at t = tn?1 . Letting y = (u; v), we can reformulate this as follows: (yt + Ah y; )L2 (I ;H01L2 ) = (F; )L2(I ;H01 L2 ) 8 2 [Sph Pq?1 (In )]2: (12) One of the most appealing properties of this scheme is that it conserves energy in the same way as the continuous problem. Letting = vt in (10) and = ut in (11) we obtain n
n
n
n
n
En = En?1 + (f; ut)S where
n
En = 21 jjvnjj2L2 + 12 jjrunjj2L2 ;
and if f = 0 then the energy is conserved. We nish this section with an outline of the remainder of this paper. In Section 3 we analyse CTG methods for abstract problems in the form (2). In Section 4 we apply these results to the wave equation. Section 5 is devoted to the nal steps of our error estimate, and Section 6 contains some numerical results.
3 CTG approximation of an abstract IVP In this section we consider the discretization in time of an abstract initial value problem. Let H be a real Hilbert space, and let A be an operator de ned on a dense domain D(A) H , which generates a strongly continuous semi-group, which we will denote by etA . We assume that (AV; V ) 0 and that jjAV jjH C jjAV jjH for all V 2 D(A). Then jjetAV jjH jjV jjH for all V 2 H . In particular these assumptions are satis ed if A is skew-symmetric, which is the case for the wave equation; 5
however our analysis is more general, and would apply also for example to the heat equation. We consider the problem Yt + AY = F; Y (0) = Y0 : (13) Precise assumptions on Y0 and F will be stated below. In this section we will denote t and Pt simply by and P , respectively. The time-discrete CTG approximation to (13) is an element y of D(A) Sqk which satis es y(0) = Y0 and for 1 n N , (yt ; )L2(I ;H ) + (Ay; )L2(I ;H ) = (F; )L2 (I ;H ) 8 2 H Pq?1 (In ): We rst derive a basic stability estimate. Theorem 1 If y satis es (14), then n
n
(14)
n
n
o
a) jjytjjL2(H ) +jjAyjjL2(H ) C T 1=2jjAY0 jjH + T jjAF jjL2(H ) + jjF jjL2(H ) ; and for 0 t T ,
n
o
b) jjyt(t)jjH +jjAy(t)jjH C jjAY0 jjH + T 1=2jjAF jjL2(H ) + jjF jjL1 (H ) : proof: On each subinterval, (14) is equivalent to yt = ? Ay + F . Therefore, by (9), we have n
o
jjAyjjL2(I ;H ) C kn1=2jjAyn?1jjH + jjAyjjL2(I ;H ) n o C kn1=2jjAyn?1jjH + jjytjjL2(I ;H ) + jjF jjL2(I ;H ) : n
n
n
(15)
n
Taking = Ayt in (14) gives (yt ; Ayt )L2 (I ;H ) + (Ay; Ayt)L2 (I ;H ) = (F; Ayt )L2 (I ;H ); 1 jjAy jj2 ? 1 jjAy jj ?(AF; y ) t L2 (I ;H ): 2 n H 2 n?1 H Summing over n gives 1 jjAy jj2 1 jjAy jj2 ? (AF; y ) (16) t L2 ([0;t ];H ): 2 n H 2 0 H On In let w(t) = t ? tn?1 and write y = yn?1 + wy with y 2 H Pq?1 (In ). Choosing = y in (14) gives (yt ; y)L2 (I ;H ) + (Ay; y)L2 (I ;H ) = (F; y)L2 (I ;H ); (y ; y)L2 (I ;H ) + (wyt ; y)L2 (I ;H ) + (Ayn?1 ; y)L2 (I ;H ) + (wAy; y)L2 (I ;H ) = (F; y)L2 (I ;H ); jjyjj2L2(I ;H ) + (wyt; y)L2(I ;H ) ?(Ayn?1 ; y)L2(I ;H ) + (F; y)L2(I ;H ): n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
6
n
n
Integration by parts in time establishes that
(wyt; y)L2 (I ;H ) = ? 12 jjyjj2L2 (I ;H ) + 21 kn jjyn jj2H : n
Thus
n
1 jjyjj2 2 L2 (I ;H ) jjAyn?1 jjL2 (I ;H )jjyjjL2 (I ;H ) + jjF jjL2(I ;H )jjyjjL2(I ;H ) = kn1=2jjAyn?1 jjH jjyjjL2 (I ;H ) + jjF jjL2(I ;H )jjyjjL2(I ;H ) n
n
n
n
n
whence
n
n
n
jjyjjL2(I ;H ) 2 kn1=2jjAyn?1 jjH + jjF jjL2(I n
n
n
;H )
o
:
(17)
By the inverse estimate (8) and properties of w, we have jjytjjL2(I ;H ) = jjy + wytjjL2(I ;H ) jjyjjL2(I ;H ) + jjwytjjL2(I ;H ) C jjyjjL2(I ;H ): Equations (15)-(18) combine to give n
n
n
n
(18)
n
n
o
jjytjj2L2(I ;H ) + jjAyjj2L2(I ;H ) C kn jjAyn?1jj2H + jjF jj2L2(I ;H ) n o C kn jjAy0jj2H + knjjAF jjL2([0;t ];H )jjytjjL2([0;t ];H ) + jjF jj2L2(I ;H ) : (19) n
n
n
n
n
n
Summing over n yields jjytjj2L2([0;t ];H ) + jjAyjj2L2([0;t ];H ) n o C tn jjAy0jj2H + tnjjAF jjL2([0;t ];H )jjytjjL2([0;t ];H ) + jjF jj2L2([0;t ];H ) : (20) A simple kickback argument completes the proof of the rst result. To obtain the pointwise in time estimate, by (7) and (19), n
n
n
n
n
n
o
jjytjj2L1(I ;H ) + jjAyjj2L1(I ;H ) Ckn?1 jjytjj2L2(I ;H ) + jjAyjj2L2(I ;H ) n o C jjAy0jj2H + jjAF jjL2([0;t ];H )jjytjjL2([0;t ];H ) + kn?1 jjF jj2L2(I ;H ) n o C jjAy0jj2H + tnjjAF jj2L2([0;t ];H ) + t?n 1jjytjj2L2([0;t ];H ) + kn?1 jjF jj2L2(I ;H ) : n
n
n
n
n
n
n
n
Now by (20) we have
n
n
n
o
jjytjj2L1(I ;H ) + jjAyjj2L1(I ;H ) C jjAy0jj2H + tn jjAF jj2L2([0;t ];H ) + jjF jj2L1([0;t ];H ) : n
n
n
n
The desired result follows. //// If H is nite dimensional, then existence and uniqueness of the CTG approximation follow at once from the preceeding stability estimate. The next theorem shows that this holds true in general. 7
Theorem 2 Given Y0 2 D(A) and F which satis es jjAF jjL2(H ) + jjF jjL2(H ) < 1, there is a unique y 2 D(A) Sqk which satis es (14). proof: Let f n g be an orthonormal basis for H , with n 2 D(A), and set Hn = spanf 1 ; :::; ng. Let y0n be the orthogonal projection of Y0 into Hn . It suces to consider a generic time interval, such as I = [0; 1]. For each n there exists a unique solution to the ( nite-dimensional) problem: nd y n 2 Hn Pq (I ) such that y n (0) = y0n and (ytn + Ay n ; )L2 (I;H ) = (F; )L2 (I;H ) 8 2 Hn Pq?1 (I ): By the previous theorem and the inequality Z
jjyjjL2(I;H ) jjy(0)jjH + jjyt(t)jjH dt jjy(0)jjH + jjytjjL2(I;H ) I
it follows that fjjy n jjL2(I;H ) + jjAy n jjL2 (I;H )g is bounded. Because A is necessarily closed (being the generator of a strongly continuous semi-group), from this we can deduce that there is a subsequence, still denoted y n , such that y n converges weakly in L2 (I; H ) to some y 2 D(A), and further that Ayn converges weakly in L2(I; H ) to Ay. Since (ytn ; )L2 (I;H ) = ?(y n ; t )L2 (I;H ) ! ?(y; t )L2 (I;H ) = (yt ; )L2(I;H ) we also have that ytn converges weakly to yt . To show that y satis es (14), given 2 H Pq?1 (I ), let m be the orthogonal projection of into Hm Pq?1 (I ). Then for n m, (ytn + Ay n ; m )L2 (I;H ) = (F; m )L2 (I;H ) Fix m, and let n ! 1, and then let m ! 1. It only remains to show that the initial condition is satis ed. But this is trivial: by construction y n (0) converges in H to Y0 , and it also easy to deduce that y n (0) converges weakly in H to y (0); it follows that y (0) = Y0 . This proves existence. Uniqueness follows immediately from the stability estimate. //// The previous theorem guarantees that y (t) 2 D(A); standard arguments show that y (t) will have more regularity (i.e., lie in the domain of higher powers of A) under the appropriate assumptions on Y0 and F , and this fact will be tacitly used below. The stability estimate also allows us to prove the following error estimate: Theorem 3 Let Y be the solution of (13), and y the CTG approximation de ned by (14). Then for 0 t T , n
o
jjAy(t) ? AY (t)jjH Ckq+1 T 1=2jj@tq+1A2Y jjL2(H ) + jj@tq+1AY jjL1(H ) : 8
proof: Write y ? Y = (y ? PY ) ? (Y ? PY ) = + . Note that 2 H Sqk and (0) = 0. A short calculation establishes that satis es, for any n and any 2 H Pq?1 (In ),
(t + A; )L2 (I ;H ) = (A; )L2(I ;H ) : The stated estimate follows by applying Theorem 1b to , and estimating using (6). //// Our next goal is to derive a higher order estimate for the error at time nodes t = tn . For this we will need the following stability result: Lemma 1 The solution y of (14) satis es n
n
q X q q +1 1 = 2 2 q +1 2 q +1 jj@t A yjjL2(H ) CT jjA Y0jjH + CT jjA F jjL2(H ) + C jj@tq?j Aq+j F jjL2(H ): j =1
proof: Recall that on each subinterval In , y satisi es yt + Ay = F . Operating on this identity with @ti?1Aj gives
@tiAj y = ?@ti?1 Aj+1 y + @ti?1 (I ? )Aj+1 y + @ti?1 Aj F:
(21)
For the second term on the right hand side, we have by (5) and (8) that jj@ti?1(I ? )Aj+1yjjL2(I ;H ) Ckn?(i?1)jj(I ? )Aj+1yjjL2(I ;H ) n
Ck?(i?1) ki?1 jj@ i?1 Aj+1 yjj n
n
t
n
L2 (I ;H ); n
and thus taking norms in (21) gives
jj@tiAj yjjL2(I ;H ) C jj@ti?1Aj+1yjjL2(I ;H ) + jj@ti?1Aj F jjL2(I ;H ): n
n
n
(22)
By summing over n and repeated use of (22) we obtain
jj@tqAq+1 yjjL2(H ) C jj@tA2q yjjL2(H ) + C
qX ?1
jj@tq?j Aq+j F jjL2(H ):
j =1 The proof is now completed by applying Theorem 1 to A2q y , which is just the CTG solution to the problem with initial data A2q Y0 and nonhomogeneity A2q F . ////
Following is the nal result of this section. Theorem 4 Let Y be the solution of (13), and y the CTG approximation de ned by (14). Then, assuming Y0 and F have the indicated regularity, for 1 n N , 8
