A Posteriori Error Estimates for Nonlinear ... - Semantic Scholar
A Posteriori Error Estimates for Nonlinear Problems. Lr (0, T ; Lρ (Ω))-Error Estimates for Finite Element Discretizations of Parabolic Equations ¨rth R. Verfu Fakult¨ at f¨ ur Mathematik, Ruhr-Universit¨ at Bochum, D-44780 Bochum, Germany
Summary: Using the abstract framework of [10] we analyze a residual a posteriori error estimator for space-time finite element discretizations of parabolic pdes. The estimator gives global upper and local lower bounds on the error of the numerical solution. The finite element discretizations in particular cover the so-called θ-scheme, which includes the implicit and explicit Euler methods and the Crank-Nicholson scheme. As particular examples we consider scalar quasilinear parabolic pdes of 2nd order and the time-dependent incompressible Navier-Stokes equations. Key words: A posteriori error estimates; nonlinear parabolic pdes; time-dependent Navier-Stokes equations; space-time finite elements; θ-scheme. AMS Subject Classification: 65N30, 65N15, 65J15, 76D05
1. Introduction We analyze a residual a posteriori error estimator for space-time finite element discretizations of parabolic pdes. Each space-time element K × J contributes the weighted sum of three terms: (1) The residual of the computed numerical solution with respect to the strong form of the differential operator evaluated on K × J, (2) the jump across ∂K × J of that trace operator which naturally connects the strong and the weak formulation of the differential equation, and (3) the jump of the numerical solution across K × ∂J. Here, K stands for an arbitrary element of the spatial mesh and J denotes an arbitrary interval of the time mesh. We could also extend our analysis to error estimators which are based on the solution of auxiliary local time-dependent problems. We do not follow this line here, in order not to overload the presentation. In order to construct our a posteriori error estimator and to prove that it yields upper and lower bounds on the error, we use the techniques introduced in [10] and consider in Section 2 abstract nonlinear problems of the form F (u) = 0 1
(1.1)
and corresponding discretizations of the form Fh (u0h ) = 0.
(1.2)
Here, F ∈ C 1 (X, Y ∗ ) and Fh ∈ C(Xh , Yh∗ ), Xh ⊂ X and Yh ⊂ Y are finite dimensional subspaces of the Banach spaces X and Y , and ∗ denotes the dual of a Banach space. If u0 ∈ X is a solution of problem (1.1) such that DF (u0 ) is an isomorphism of X onto Y ∗ and DF is Lipschitz continuous at u0 , we know from Proposition 2.1 in [10] that ckF (u)kY ∗ ≤ ku − u0 kX ≤ ckF (u)kY ∗ (1.3) holds for all u in a suitable neighbourhood of u0 . The constants c and c depend on DF (u0 ) and DF (u0 )−1 . They measure the sensitivity of the infinite dimensional problem (1.1) with respect to small perturbations. For a simple model problem we derive explicit bounds for c and c in Section 4. When applying estimate (1.3) to an approximate solution uh ∈ Xh of problem (1.2) one must evaluate the residual kF (uh )kY ∗ . This is as expensive as the solution of the original problem (1.1) since it amounts in the solution of an infinite dimensional maximization problem. In order to obtain an approximation of kF (uh )kY ∗ which is much easier to evaluate, we introduced in [10] an approximation F˜h : Xh → Y ∗ of F at uh , a finite dimensional subspace Y˜h ⊂ Y , and a restriction operator Rh : Y → Yh . For practical applications, Rh is a suitable local quasi-interpolation operator (cf. (3.10)), Feh (uh ) is obtained by locally projecting F (uh ) onto a suitable finite element space, and Yeh consists of judiciously chosen local cut-off functions (cf. Section 3, in particular Lemma 3.5). When applying these results to parabolic pdes, we obtain error estimates in a L (0, T ; W01,ρ (Ω))-norm (cf. Section 3 for a difinition of these spaces and norms). The ′ ′ space Y then consists of functions in Lp (0, T ; W01,π (Ω)) having their time derivative ′ ′ ′ in Lp (0, T ; W −1,π (Ω)). Due to the non-local nature of the space W −1,π (Ω) we get into troubles when deriving lower bounds on the error. This problem will be tackled in a subsequent paper. Here, we circumvent this difficulty by imposing a weaker ′ ′ ′ Lr (0, T ; Lρ (Ω))-norm on X and a stronger Lp (0, T ; W 2,π (Ω)∩W01,π (Ω))-norm on Y . The corresponding spaces will be denoted by X− and Y+ , respectively. In particular, ′ ′ the functions in Y+ now have time derivatives in Lp (0, T ; Lπ (Ω)). r
In Sections 4 and 5 we apply the general results of the previous sections to scalar quasilinear parabolic pdes of 2nd order and the time-dependent incompressible Navier-Stokes equations. We consider Petrov-Galerkin space-time finite elements. These discretizations are inspired by the concept of very weak solutions. The trial space Xh consists of functions which are discontinuous in time, piecewise polynomials of degree k, k ≥ 0. The test space Yh consists of functions which are continuous 2
in time, piecewise polynomials of degree k + 1 and which vanish at the final time T . This discretization corresponds to an implicit (k + 1)-stage Runge-Kutta-scheme. When applied to a linear problem its stability function is the (k + 1)-st diagonal Pad´e-approximation. For k = 0 we in particular obtain the Crank-Nicholson scheme. By slightly modifying the basis-functions of Yh we may also recover the popular θscheme for all θ ∈ [0, 1]. This in particular covers the explicit (θ = 0) and implicit (θ = 1) Euler schemes. We obtain global upper and local lower bounds for the error measured in an L (0, T ; Lρ (Ω))-norm. The upper and lower bounds differ by a factor 1+τ −1 h2 +τ h−2 . Here, h and τ are the local mesh-sizes in space and time, respectively. This factor reflects the fact that the differential operator is of 2nd order with respect to the space variables but only of 1st order with respect to the time variable. The local lower bounds may be combined to a global lower bound of the same type. When applied to the corresponding particular examples, our error estimates are similar to those obtained in [7, 8] using a different approach and a different discretization. r
2. Abstract Error Estimates Let X, Y be two Banach spaces with norms k.kX and k.kY . For any element u ∈ X and any real number R > 0 set BX (u, R) := {v ∈ X : ku − vkX < R}. We denote by L(X, Y ) and Isom(X, Y ) ⊂ L(X, Y ) the Banach space of continuous linear maps of X in Y equipped with the operator norm k.kL(X,Y ) and the open subset of linear homeomorphisms of X onto Y . By Y ∗ := L(Y, IR) and < ., . >Y we denote the dual space of Y and the corresponding duality pairing. Finally, A∗ ∈ L(Y ∗ , X ∗ ) denotes the adjoint of a given operator A ∈ L(X, Y ).
Let F ∈ C 1 (X, Y ∗ ) be a given continuously differentiable function. Given a solution u0 ∈ X of problem (1.1) and an arbitrary element u ∈ X ”close” to u0 , we may estimate the error ku − u0 kX by the residual kF (u)kY ∗ (cf. Proposition 2.1 in [10]). For parabolic pde’s we thus obtain control on the Lr (0, T ; W01,ρ (Ω))-norm of the error (cf. Section 3 for the definition of these spaces and norms). However, we are interested in controlling the Lr (0, T ; Lρ (Ω))-norm of the error. In order to achieve this within the present abstract framework, we must enlarge the space X and reduce the space Y . We therefore consider three additional Banach spaces X− , X+ , and Y+ such that X+ ⊂ X ⊂ X− and Y+ ⊂ Y with continuous and dense injections. Here, the +/− sign indicates a space with a stronger / weaker norm. We assume that X− is reflexive. 3
2.1 Proposition. Let u0 ∈ X be a solution of problem (1.1). Assume that u0 ∈ X+ , ∗ that DF (u0 )∗ ∈ Isom(Y+ , X− ), and that there are two numbers R0 > 0 and β > 0 such that k[DF (u0 ) − DF (u0 + tw)]wkY+∗ ≤ βtkwkX+ kwkX− (2.1) ∀w ∈ BX+ (0, R0 ), t ∈ [0, 1]. Set
−1
−1 R := min{R0 , β −1 kDF (u0 )∗ k−1 kDF (u0 )∗ kL(Y+ ,X−∗ ) }. L(X ∗ ,Y+ ) , 2β −
Then the following error estimate holds for all u ∈ BX+ (u0 , R) : 1 ∗−1 kDF (u0 )∗ k−1 kL(X−∗ ,Y+ ) kF (u)kY+∗ . ∗ ) kF (u)kY ∗ ≤ ku − u0 kX− ≤ 2kDF (u0 ) L(Y ,X + + − 2 (2.2) ∗ Proof. Let u ∈ BX+ (u0 , R). Consider an arbitrary element w ∈ X− and set ϕ := −1 ∗ DF (u0 ) w ∈ Y+ . We then have
< u − u0 , w >X− = < DF (u0 )(u − u0 ), ϕ >Y+
= < F (u), ϕ >Y+ Z 1 < [DF (u0 ) − DF (u0 + t(u − u0 ))](u − u0 ), ϕ >Y+ dt. + 0
−1
Inequality (2.1) and the continuity of DF (u0 )∗ imply that Z 1 < [DF (u0 ) − DF (u0 + t(u − u0 ))](u − u0 ), ϕ >Y+ dt| | 0
1 ≤ βku − u0 kX+ ku − u0 kX− kϕkY+ 2 −1 1 ≤ βkDF (u0 )∗ kL(X−∗ ,Y+ ) Rku − u0 kX− kwkX−∗ . 2 Combined with the above representation of < u − u0 , w >X− this yields −1 1 < u − u0 , w >X− ≤kF (u)kY+∗ kϕkY+ + βkDF (u0 )∗ kL(X−∗ ,Y+ ) Rku − u0 kX− kwkX−∗ 2 −1 1 ≤{kDF (u0 )∗ kL(X−∗ ,Y+ ) kF (u)kY+∗ + ku − u0 kX− }kwkX−∗ . 2 ∗ Since X− is reflexive and w ∈ X− was arbitrary, this implies the upper bound of estimate (2.2). In the same way, we obtain < F (u), ϕ >Y+ = < u − u0 , w >X− −
Z
0
1
< [DF (u0 ) − DF (u0 + t(u − u0 ))](u − u0 ), ϕ >Y+ dt
1 ≤ku − u0 kX− kwkX−∗ + βku − u0 kX+ ku − u0 kX− kϕkY+ 2 1 ≤{kDF (u0 )∗ kL(Y+ ,X−∗ ) ku − u0 kX− + βRku − u0 kX− }kϕkY+ 2 ∗ ≤2kDF (u0 ) kL(Y+ ,X−∗ ) ku − u0 kX− kϕkY+ . 4
Since ϕ ∈ Y+ is arbitrary, this proves the lower bound of estimate (2.2).
∗ The condition DF (u0 )∗ ∈ Isom(Y+ , X− ) of Proposition 2.1 is more restrictive ∗ than the assumption DF (u0 ) ∈ Isom(X, Y ) which is needed to bound ku − u0 kX by kF (u)kY ∗ (cf. Proposition 2.1 in [10]). For pdes, it is equivalent to an additional regularity condition. For linear problems, i.e. when DF is constant, one may extend F by continuity to a continuously differentiable map of X− to Y+∗ . Then the space X+ is not needed. For nonlinear problems, however, this extension often is impossible or the derivative of the extension is no longer Lipschitz continuous. This is the place where the space X+ comes into play.
Let Xh ⊂ X+ and Yh ⊂ Y be finite dimensional subspaces and Fh ∈ C(Xh , Yh∗ ) be an approximation of F . Given an approximate solution uh ∈ Xh of problem (1.2), Proposition 2.1 allows us to estimate the error ku0 −uh kX− by the residual kF (uh )kY+∗ . The evaluation of the latter, however, is a difficult task since it amounts in the solution of an infinite dimensional maximization problem. In order to obtain an approximation of the residual which is easier to compute, we introduce a restriction operator Rh ∈ L(Y, Yh ), a finite dimensional subspace Y˜h ⊂ Y+ , and an approximation F˜h : Xh → Y ∗ of F at uh . In the context of pdes F˜h is obtained by locally freezing the coefficients of the differential operator. We equip Yeh with the norm of Y+ .
2.2 Proposition. The following estimates hold:
kF (uh )kY+∗ ≤k(IdY+ − Rh )∗ F˜h (uh )kY+∗ + k(IdY+ − Rh )∗ [F (uh ) − F˜h (uh )]kY+∗ + kRh∗ [F (uh ) − Fh (uh )]kY+∗ + kRh∗ Fh (uh )kY+∗
(2.3)
and kF˜h (uh )kY˜ ∗ ≤ kF (uh )kY˜ ∗ + kF (uh ) − F˜h (uh )kY˜ ∗ . h
h
h
(2.4)
If there is a constant c0 , which does not depend on h, such that k(IdY+ − Rh )∗ F˜h (uh )kY+∗ ≤ c0 kF˜h (uh )kY˜ ∗ , h
(2.5)
then both k(IdY+ − Rh )∗ F˜h (uh )kY+∗ and kF˜h (uh )kY˜ ∗ yield upper and lower bounds for h the residual kF (uh )kY+∗ . Proof. Estimate (2.3) follows from the identity
< F (uh ), ϕ >Y+ = < F˜h (uh ), ϕ − Rh ϕ >Y+ + < F (uh ) − F˜h (uh ), ϕ − Rh ϕ >Y+ + < F (uh ) − Fh (uh ), Rh ϕ >Y+ + < Fh (uh ), Rh ϕ >Y+
which holds for all ϕ ∈ Y+ . Estimate (2.4) follows from the triangle inequality. The statement concerning upper and lower bounds for the residual is an obvious consequence of inequality (2.5). 5
A non-optimal estimate of the third and fourth term on the right-hand side of estimate (2.3) is given by kRh∗ [F (uh ) − Fh (uh )]kY+∗ + kRh∗ Fh (uh )kY+∗
≤kIdY+ kL(Y+ ,Y ) kRh kL(Y,Yh ) {kF (uh ) − Fh (uh )kYh∗ + kFh (uh )kYh∗ }. Here, Yh is equipped with the norm of Y . The second terms on the right-hand sides of estimates (2.3) and (2.4) measure the quality of the approximation F˜h (uh ) to F (uh ). Usually, they are higher order terms when compared with k(IdY+ − Rh )∗ F˜h (uh )kY+∗ and kF˜h (uh )k ˜ ∗ . The term kR∗ [F (uh ) − Fh (uh )]kY ∗ is the consistency error of the h
Yh
+
discretization. The term kRh∗ Fh (uh )kY+∗ measures the residual of the algebraic equation (1.2) and can be evaluated by standard methods. Inequality (2.5) is a non-trivial condition since it claims that a supremum with respect to an infinite dimensional space may be bounded from above by a supremum with respect to a finite dimensional space.
Propositions 2.1 and 2.2 will be used in the following way. Any upper bound for k(IdY+ − Rh )∗ F˜h (uh )kY+∗ yields a reliable error estimator. If in addition the error estimator yields a lower bound for kF˜h (uh )k ˜ ∗ , it is also efficient. Here, we adopt Yh
the standard convention that an error estimator is called reliable (resp. efficient) if it yields an upper (resp. lower) bound for the error. Moreover, upper and lower bounds may contain multiplicative constants which do not depend on the discretization parameter h. The lower bound of estimate (2.2) follows from the identity < F (u), ϕ >Y+ = < DF (u0 )(u − u0 ), ϕ >Y+ Z 1 < [DF (u0 ) − DF (u0 + t(u − u0 ))](u − u0 ), ϕ >Y+ dt − 0
with ϕ ∈ Y+ . In many applications Y+ will be a closed subspace of a suitable Sobolev space. When inserting in the above identity functions ϕ with a local support, one then obtains lower bounds for u − u0 restricted to the support of ϕ.
The factor kDF (u0 )∗ kL(Y+ ,X−∗ ) in the lower bound of estimate (2.2) is rather harmless since it corresponds to a differential operator which is local and the norm of which can easily be estimated in terms of its coefficients. On the other hand, −1 the factor kDF (u0 )∗ kL(X−∗ ,Y+ ) in the upper bound of estimate (2.2) is much more severe. In some applications, in particular when X− is a Hilbert space, a function ∗ w ∈ X− satisfying kwkX−∗ = 1
and hu − u0 , wiX− = ku − u0 kX− −1
can explicitly be given in terms of u and u0 . Then kDF (u0 )∗ kL(X−∗ ,Y+ ) can be −1
replaced by kDF (u0 )∗ wkY+ . The latter quantity may be estimated numerically by approximately solving a discrete analogue of the corresponding adjoint pde. 6
The Lipschitz condition of Proposition 2.1 may be replaced by weaker H¨older conditions of the form k[DF (u0 ) − DF (u0 + tw)]wkY+∗ ≤ βtα kwkX− kwkα X+ ∀w ∈ BX+ (0, R0 ), t ∈ [0, 1].
3. Auxiliary results Function spaces. Let Ω be a bounded, connected, open domain in IRn , n ≥ 2, with polyhedral boundary Γ. For any open subset ω of Ω with Lipschitz boundary γ, we denote by W k,p (ω), k ∈ IN, 1 ≤ p ≤ ∞, Lp (ω) := W 0,p (ω), and Lp (γ) the usual Sobolev and Lebesgue spaces equipped with the standard norms (cf. [1] and Vol. 3, Chap. IV in [6]). Let W01,p (Ω) := {u ∈ W 1,p (Ω) : u = 0 on Γ} and set for 1 < p < ∞
′
W −1,p (Ω) := W01,p (Ω)∗ .
Here, p′ denotes the dual exponent of p defined by p1′ + p1 = 1. In what follows, a prime will always denote the dual of a given Lebesgue exponent. Let V and W be two Banach spaces such that V ⊂ W with continuous and dense injection. Given two real numbers a and b with a < b, we denote by Lp (a, b; V ), 1 ≤ p < ∞, the space of measurable functions u defined on (a, b) with values in V such that the function t → ku(., t)kV is in Lp ((a, b)). Lp (a, b; V ) is a Banach space equipped with the norm Z b
kukLp (a,b;V ) := {
a
ku(., t)kpV dt}1/p
(cf. [6], Vol. 5, Chap. XVIII, §1). Slightly changing the notation of [6], we further consider the Banach space W p (a, b; V, W ) := {u ∈ Lp (a, b; V ) : ∂t u ∈ Lp (a, b; W )} equipped with the norm kukW p (a,b;V,W )
Z := {
a
b
ku(., t)kpV
dt +
Z
a
b
k∂t u(., t)kpW dt}1/p .
Here, the partial derivative ∂t u must be interpreted in the distributional sense (cf.[6]; loc. cit.). For all smooth functions ϕ ∈ D((a, b)) it satisfies the identity Z
a
b
∂t u(., t)ϕ(t)dt = − 7
Z
a
b
u(., t)ϕ′ (t)dt,
where the integrals are taken in W . From Proposition 9 in [6; loc. cit.] it follows that the traces u(., a) and u(., b) are defined as elements of W , provided p > 1. We therefore set for 1 < p < ∞ Wbp (a, b; V, W ) := {u ∈ W p (a, b; V, W ) : u(., b) = 0}. Given an interval I ⊂ IR, a suitable subset ω ⊂ Ω, and numbers 1 ≤ p, π < ∞, we will use the following abbreviation: Lpπ (ω × I) := Lp (I; Lπ (ω)),
Wπp (ω × I) := W p (I; W 2,π (ω), Lπ (ω)).
Finite element partition. Denote by T > 0 an arbitrary but fixed time. Let Iτ = {[tj , tj+1 ] : 1 ≤ j ≤ Nτ }, τ > 0, be a family of partitions of [0, T ] with 0 = t1 < t2 < ... < tNτ +1 = T. For 1 ≤ j ≤ Nτ set Jj := [tj , tj+1 ] , τj := tj+1 − tj . We assume that the family Iτ is shape regular, i.e., the ratios τj /τj+1 and τj+1 /τj are bounded from above independently of j and τ . With each 1 ≤ j ≤ Nτ , we associate a partition Tj of Ω into n-simplices. We denote by Ej the set of the interior faces of Tj . For K ∈ Tj and E ∈ Ej let hK , ρK , and hE be the diameter of K, the diameter of the largest ball inscribed into K, and the diameter of E. We assume that the partitions Tj satisfy the following two conditions: (1) admissibility: Any two simplices of Tj are either disjoint or share a complete smooth submanifold of their boundaries.
(2) shape regularity: The ratio hK /ρK is bounded from above independently of K ∈ Tj , j, and τ .
Condition (2) allows the use of locally refined meshes. It implies that the ratio hK /hE , for all K ∈ Tj , all faces E of K, and all j, is bounded from above and from below by constants which do not depend on K, E, j, and τ . Denote by Pτ := {K × Jj : 1 ≤ j ≤ Nτ , K ∈ Tj } the partition of the space-time cylinder Ω × [0, T ] into prisms which is induced by Iτ and the Tj ’s. For any E ∈ Ej , 1 ≤ j ≤ Nτ , and any piecewise continuous function u, we denote by [u]E the jump of u across E in an arbitrary, but fixed direction nE orthogonal to 8
E. Finally, we introduce the following neighbourhoods of elements and points U (tj ) := [tmax{1,j−1} , tmin{Nτ +1,j+1} ] , 1 ≤ j ≤ Nτ ,
U (Jj ) := U (tj ) ∪ U (tj+1 ) = [tmax{1,j−1} , tmin{Nτ +1,j+2} ] , 1 ≤ j ≤ Nτ , ′ U (K) := ∪K ′
, K ∈ Tj , 1 ≤ j ≤ Nτ ,
U (E) := ∪K′ ′
, E ∈ Ej , 1 ≤ j ≤ Nτ ,
K ∩K6=∅ K ′ ∈Tj
E∩K 6=∅ K ′ ∈Tj
U (P ) := ∪P′
′
P ∩P 6=∅ P ′ ∈Pτ
ωK :=
∪K ′
E⊂∂K K ′ ∈Tj
, P ∈ Pτ , , K ∈ Tj , 1 ≤ j ≤ Nτ ,
K ′ ∩K∈Ej K ′ ∈Tj
ωE := ∪K ′′
(3.1)
, E ∈ Ej , 1 ≤ j ≤ Nτ .
Here, K ∩K ′ ∈ Ej means that K and K ′ share a complete (n−1)-dimensional smooth submanifold of their boundaries. Finite element spaces. Denote by IPk , k ≥ 0, the space of polynomials (in x) of degree at most k. Given an admissible partition Th of Ω, we define the finite element spaces (in x) as usual: Shk,−1 := {u : Ω → IR : u|K ∈ IPk
Shk,0 k,0 Sh,0
:=
Shk,−1
:= {u ∈
∩ C(Ω)
Shk,0
, k ≥ 1,
∀K ∈ Th } , k ≥ 0,
(3.2)
: u = 0 on Γ} , k ≥ 1.
k,0 For any k ≥ 1 and any 1 ≤ p ≤ ∞ we obviously have Shk,0 ⊂ W 1,p (Ω) and Sh,0 ⊂ 1,p W0 (Ω).
Let V1 , ..., VNτ be finite element subspaces of C(Ω) associated with the partitions T1 , ..., TNτ introduced above. We then define finite element spaces in space and time by Sτk,−1 (Vh(τ ) ) := span{χJj (t)tµ vj,µ (x) : 0 ≤ µ ≤ k, 1 ≤ j ≤ Nτ , vj,µ ∈ Vj } , k ≥ 0. (3.3) Here, χJj denotes the characteristic function of the interval Jj . The elements of Sτk,−1 (Vh(τ ) ) are discontinuous at the intermediate points t2 , ..., tNτ . But the leftsided limits u(., tj − 0) := lim u(., tj − t) exist for all 2 ≤ j ≤ Nτ + 1; and the 0:=− < u0 , ϕ(., 0) >W 1,π Z T < u(., t), ∂t ϕ(., t) >W 1,π dt − 0
+
Z
T
0
Z
Ω
(4.2)
{a(x, u, ∇u)∇ϕ − b(x, u, ∇u)ϕ}dxdt.
We recall that a prime denotes the dual of a Lebesgue exponent, i.e. In order to ensure that F is well defined, we alsways assume that
1 p
+
1 p′
ρ ≥ π.
= 1 etc. (4.3)
Within the framework of Section 2, we further set ′
′
′
′
Y+ := WTp (0, T ; W 2,π (Ω) ∩ W01,π (Ω), Lπ (Ω)),
X− := Lr (0, T ; Lρ (Ω)),
(4.4)
X+ := Lr (0, T ; W01,σ (Ω)) where ρ ≤ σ < ∞.
In order to better understand the flavour of problem (1.1) and definition (4.2), we recall the notions of weak and very weak solutions of problem (4.1) (cf. [2], Sections 11 and 13). A function u ∈ W r (0, T ; W01,ρ (Ω), W −1,π (Ω)) is called a weak solution of problem (4.1) if u(., 0) = u0 in W −1,π (Ω) and Z
0
T
< ∂t u(., t), ϕ(., t) >W 1,π′ + ′
′
Z
0
T
Z
Ω
{a(x, u, ∇u)∇ϕ − b(x, u, ∇u)ϕ}dxdt = 0
for all ϕ ∈ Lp (0, T ; W 1,π (Ω)). It is a very weak solution if u ∈ X− and Z T − < u0 , ϕ(., 0) >Lρ − < u(., t), ∂t ϕ(., t) >Lρ dt +
Z
0
0
T
{< a(., u, ∇u), ∇ϕ >W 1,π′ − < b(x, u, ∇u), ϕ >W 1,π′ }dt = 0 19
for all ϕ ∈ Y+ . Obviously, every solution of problem (1.1) is a very weak solution of problem (4.1). Conversely, every very weak solution of problem (4.1) is a solution of problem (1.1) if it is contained in X. Using integration by parts with respect to the time variable, one sees that every weak solution of problem (4.1) is a solution of problem (1.1) and that the converse is true if the solution of problem (1.1) is contained in W r (0, T ; W01,ρ (Ω), W −1,π (Ω)). In this sense, a solution of problem (1.1) is weaker than a weak solution and stronger than a very weak solution of problem (4.1). ∗ One easily checks that DF (u)∗ ∈ Isom(Y+ , X− ) if and only if the adjoint linearized problem
−∂t w − ∇ · (A(x, u, ∇u)∇w) + ∂y a(x, u, ∇u) · ∇w
+∇ · (∇z b(x, u, ∇u)w) − ∂y b(x, u, ∇u)w = g in Ω × (0, T )
w = 0 on Γ × (0, T )
w(., T ) = 0 in Ω
∗ admits, for each g ∈ X− , a unique solution w ∈ Y+ such that kwkY+ ¹ kgkX−∗ .
Examples. We consider two particular examples: (1) a heat equation with nonlinear source term a(x, u, ∇u) = ∇u,
b(x, u, ∇u) = f (u),
f ∈ C 1 (IR, IR),
|f ′ (s)| ≤ γ
∀s ∈ IR,
p = r = π = ρ = 2. (2) a nonlinear convection-diffusion equation a(x, u, ∇u) := k(u)∇u,
b(x, u, ∇u) := f − c(x, u) · ∇u,
f ∈ L∞ (Ω), c ∈ C 1 (Ω × IR, IRn ), k ∈ C 2 (IR, IR), k(s) ≥ α > 0, |k (l) (s)| ≤ γ |∂s c(x, s)| ≤ γ
∀s ∈ IR, l ∈ {0, 1, 2},
∀x ∈ Ω, s ∈ IR, n , r ≥ 2p. ρ = π ∈ (n, 4), p > n−1
If in example (1) the constant γ is sufficiently small, we may use an energy estimate −1 and a perturbation argument to obtain explicit bounds on kDF (u)∗ kL(X−∗ ,Y+ ) in terms of the norm of the inverse Laplacian with homogeneous Dirichlet boundary conditions. More precisely, denote by c∆ := k(−∆)−1 kL(L2 (Ω),W 1,2 (Ω)∩W 2,2 (Ω)) 0
20
the norm of the inverse Laplacian with homogeneous Dirichlet boundary conditions. This quantity only depends on the geometry of Ω. Inserting a and b given above into the adjoint linearized problem, we immediadely see that DF (u)∗ = L − N where L is the operator associated with the time-reversed heat equation −∂t v − ∆v = g
in Ω × (0, T )
on Γ × (0, T )
v=0
v(., T ) = 0
(4.5)
in Ω
and where the operator N is given by < N v, ϕ >Y+ =
Z
T
0
Z
f ′ (u)vϕ.
(4.6)
Ω
From equations (4.5) and (4.6) we deduce that kLkL(Y+ ,X−∗ ) ≤ 2,
kN kL(Y+ ,X−∗ ) ≤ γ
and, hence, kDF (u)∗ kL(Y+ ,X−∗ ) ≤ 2 + γ. Multiplying the first equation of (4.5) with −∂t v, integrating over Ω × (0, T ), and using integration by parts with respect to the space variable, we conclude that k∂t vkL22 (Ω×(0,T )) ≤ kgkL22 (Ω×(0,T )) .
(4.7)
Writing the first equation of (4.5) in the form −∆v = g + ∂t v and using the estimate (4.7), we obtain on the other hand that √ kvkL2 (0,T ;W 2,2 (Ω)) ≤c∆ 2{kgkL22 (Ω×(0,T )) + k∂t vkL22 (Ω×(0,T )) } √ ≤c∆ 8kgkL22 (Ω×(0,T )) . Estimates (4.7) and (4.8) yield kL−1 kL(X−∗ ,Y+ ) ≤ 1 + Assume that γ(1 +
√
√
8c∆ ) < 1.
21
8c∆ .
(4.8)
A standard perturbation argument then gives
−1
kDF (u)∗ kL(X−∗ ,Y+ ) =k(L − N )−1 kL(X−∗ ,Y+ ) £ ¤−1 ≤kL−1 kL(X−∗ ,Y+ ) 1 − kL−1 kL(X−∗ ,Y+ ) kN kL(Y+ ,X−∗ ) √ 1 + 8c∆ √ . ≤ 1 − γ(1 + 8c∆ )
Finite element discretization. For the discretization of problem (4.1) we proceed as in Section 3. We choose a family Iτ of shape regular partitions of the interval [0, T ]. With each time tj , 1 ≤ j ≤ Nτ , we associate an admissible and shape regular partition Tj of Ω into n-simplices and a finite element space Vj ⊂ W01,σ (Ω) corresponding to Tj and consisting of affine equivalent finite elements in the sense of [4]. We choose an integer k and a parameter θ ∈ [0, 1] and set
Xh :=Sτk,−1 (Vh(τ ) ), Yh :=Sτθ;k+1,0 (Vh(τ ) ), < Fh (uh ), ϕh >Yh := < F (uh ), ϕh >Y
(4.9) ∀uh ∈ Xh , ϕh ∈ Yh .
For simplicity, we use in (4.9) the parameter h as an acronym for the meshsizes both in space and in time. We recall that the spaces on the right-hand side of (4.9) are defined (θ) in (3.3) and (3.4) with λj replaced by λj and that Xh ⊂ X+ and Yh ⊂ Y . Hence, the discretization (4.9) is conforming. It is also consistent, i.e kRh∗ [F (uh ) − Fh (uh )]kY+∗ = 0. It is a Petrov-Galerkin discretization since the test and trial spaces are different: The trial functions are discontinuous in time, piecewise polynomials of degree k, the test functions are continuous in time, piecewise polynomials of degree k + 1.
Relation to Runge-Kutta schemes. In order to better understand the flavour of problem (1.2) and definition (4.9) we rewrite Fh . Recalling that the functions in Yh are continuous at the intermediate times t2 , ..., tNτ and vanish at the final time T 22
and using integration by parts on each time interval, we get for uh ∈ Xh and ϕh ∈ Yh < F (uh ), ϕh >Y Z u0 (x)ϕh (x, 0)dx =− +
Ω N τ X j=1
{
Z
tj+1
tj
=
+ Z
Z
0
Ω
+
T
Ω
uh (x, tj + 0)ϕh (x, tj + 0)dx −
Z
Ω
uh (x, tj+1 − 0)ϕh (x, tj+1 − 0)dx}
{a(x, uh , ∇uh )∇ϕh − b(x, uh , ∇uh )ϕh }dxdt
[uh (x, 0 + 0) − u0 (x)]ϕh (x, 0)dx
Nτ Z X j=2
+
Ω
Z
∂t uh (x, t)ϕh (x, t)dxdt
Ω
Z
+
Z
Ω
Nτ Z X j=1
[uh (x, tj + 0) − uh (x, tj − 0)]ϕh (x, tj )dx
tj+1
tj
Z
Ω
[∂t uh ϕh + a(x, uh , ∇uh )∇ϕh − b(x, uh , ∇uh )ϕh ]dxdt.
Using the convention that uh (., 0 − 0) := u0 ,
(4.10)
we may write this in the compact form < F (uh ), ϕh >Y Nτ nZ X = [uh (x, tj + 0) − uh (x, tj − 0)]ϕh (x, tj )dx j=1
Ω
+
Z
tj+1
tj
Z
Ω
(4.11)
o [∂t uh ϕh + a(x, uh , ∇uh )∇ϕh − b(x, uh , ∇uh )ϕh ]dxdt .
We first consider the case k = 0 and set ujh := uh (., tj + 0) = uh (., tj+1 − 0) ∀1 ≤ j ≤ Nτ . (θ)
Observing that uh is piecewise constant on the time intervals, inserting ϕh = λj vj , 2 ≤ j ≤ Nτ , vj ∈ Vj , as a test function in (4.11) and using Lemma 3.1, we obtain < Fh (uh ), ϕh >Yh Z = {ujh − uj−1 h }vj dx Ω Z + θτj {a(x, ujh , ∇ujh )∇vj − b(x, ujh , ∇ujh )vj }dx Ω Z j−1 j−1 j−1 + (1 − θ)τj−1 {a(x, uj−1 h , ∇uh )∇vj − b(x, uh , ∇uh )vj }dx. Ω
23
(θ)
Inserting ϕh = λ1 v1 , v1 ∈ V1 , as a test function in (4.11) we similarly get Z < Fh (uh ), ϕh >Yh = {u1h − u0 }v1 dx Ω Z + θτ1 {a(x, u1h , ∇u1h )∇v1 − b(x, u1h , ∇u1h )v1 }dx. Ω
Hence, in the case k = 0, problem (1.2) yields the popular θ-scheme. In particular, the parameters θ = 0, θ = 1, and θ = 21 correspond to the explicit Euler scheme, the implicit Euler scheme, and the trapezoidal rule (Crank-Nicholson scheme). Thus the time discretization is of first order unless θ = 21 ; in this case it is of second order. Moreover, the time discretization is A-stable if θ ≥ 12 .
Next we consider the case k ≥ 1. Denote by pˆl , 0 ≤ l ≤ k−1, a set of orthonormal polynomials of degree l with respect to the weight function 4t(1 − t) on [0, 1]. Let Rt qˆl (t) := 0 pˆl (s)ds, 0 ≤ l ≤ k − 1, and qˆ−1 (t) = 1. For 1 ≤ j ≤ Nτ , set pl,j := pˆl ◦ FJ−1 , qm,j := qˆm ◦ FJ−1 , j j
0 ≤ l ≤ k − 1, −1 ≤ m ≤ k − 1.
Then every uh ∈ Xh and every ϕh ∈ Yh have unique representations of the form uh =
Nτ X j=1
ϕh =
Nτ X j=1
χJj (t){
k−1 X
qµ,j (t)vµ,j (x)},
µ=−1 (θ) λj (t)wj (x)
+
Nτ X
ψj (t)χJj (t){
k−1 X
pµ,j (t)wµ,j (x)}
µ=0
j=1
with vµ,j , wµ,j , wj ∈ Vj , 1 ≤ j ≤ Nτ , Here, χJj denotes the characteristic function of Jj and ψj (t) = τ42 (t − tj )(tj+1 − t). Consider a fixed j ∈ {2, ..., Nτ } and insert j
ϕh = ψj (t)χJj (t)pµ,j (t)wµ,j (x), 0 ≤ µ ≤ k − 1, as a test function in (4.11). We then get < Fh (uh ), ϕh >Yh Z tj+1 Z {∂t uh ϕh + a(x, uh , ∇uh )∇ϕh − b(x, uh , ∇uh )ϕh }dxdt = tj
=
Ω
Z
vµ,j (x)wµ,j (x)dx Z tj+1 Z + { [a(x, uh , ∇uh )∇wµ,j − b(x, uh , ∇uh )wµ,j ]dx}ψj (t)pµ,j (t)dt. Ω
tj
Ω
(θ)
Inserting ϕh = λj (t)wj (x) as a test function in (4.11), we obtain on the other hand < Fh (uh ), ϕh >Yh Z = {v−1,j (x) − uh (x, tj − 0)}wj (x)dx Ω Z tj+1 Z (θ) { [∂t uh wj + a(x, uh , ∇uh )∇wj − b(x, uh , ∇uh )wj ]dx}λj (t)dt. + tj−1
Ω
24
With the obvious modifications these expressions also hold for j = 1. Hence, the coefficients v0,j , ..., vk−1,j are functions of the coefficient v−1,j and the latter one is determined by the values of uh on the previous time interval. Thus problem (1.2) amounts in a (k + 1)-stage implicit Runge-Kutta scheme. A lengthy, but straightforward calculation shows that, for linear problems and k ∈ {1, 2}, θ = 12 , this scheme corresponds to the (k + 1)-st diagonal Pad´e-approximation. In particular, the time discretization then is of order 2k + 2 and A-stable. 4.1 Remark. When writing problem (1.2) in the form (4.11) it strongly resembles the discontinuous Galerkin method (cf. e.g. [7, 8]). In the discontinuous Galerkin method, however, the test and trial spaces are identical and consist both of discontinuous in time, piecewise polynomials of degree k. In particular, the case k = 0 corresponds to the implicit Euler scheme. Due to the discontinuities at the intermediate times t2 , ..., tNτ −1 the discontinuous Galerkin method is non-conforming both with respect to the standard weak formulation of problem (4.1) and to the formulation (4.2). This complicates its analysis within the framework of Secction 2. This difficulty will be overcome in a subsequent paper. A different analysis of the discontinuous Galerkin method is given in, e.g., [7, 8]. Definition of Rh , F˜h , and Y˜h . In order to put the discretization in the framework of Section 2, we assume that Yh contains the space θτ defined in (3.6). This is equivalent to assuming that the space discretization at least consists of linear el1,0 ements, i.e. Vj ⊃ Sj,0 , 1 ≤ j ≤ Nτ . As restriction operator Rh we use the operator Iτ defined in (3.10). For the construction of F˜h and Y˜h we define integers µ, ν and approximations ah of a and bh of b as follows: if a(x, vh , ∇vh ) a(x, uh , ∇uh ), µ,−1 ∈ Sτµ,−1 (Sh(τ ) ) ∀vh ∈ Xh X ah (x, uh , ∇uh ) := π1,Q a(x, uh , ∇uh ), µ := 1, otherwise Q∈Pτ
b(x, uh , ∇uh ), if b(x, vh , ∇vh ) ν,−1 ∈ Sτν,−1 (Sh(τ ) ) ∀vh ∈ Xh X bh (x, uh , ∇uh ) := π0,Q b(x, uh , ∇uh ), ν := 0, otherwise. Q∈Pτ
(4.12) Here, uh ∈ Xh is arbitrary and π0,Q and π1,Q denote the L (Q)-projections onto the spaces of polynomials of degree at most 0 and 1 (in the variables x and t), respectively. Now, F˜h is defined in the same way as F with a and b replaced by ah and bh , respectively, and 2
Y˜h := span{ψj ψK v, ψj ψE PE σ, ψj ψK Pj w : 1 ≤ j ≤ Nτ , K ∈ Tj , E ∈ Ej , ˜ m|K×{t } }. ˜ m|E×J , w ∈ IP ˜ m|K×J , σ ∈ IP v ∈ IP j
j
25
j
(4.13)
˜ m denotes the space of polynomials of degree at Here, m := max{µ − 1, ν} and IP most m in the variables x and t. The estimators. Given Q = K × Jj ∈ Pτ , we recall the abbreviation (3.17) and set εQ,π :=(h2K + τj )k∇ · [a(., uh , ∇uh ) − ah (., uh , ∇uh )] 1
π + hK
−1
+ [b(., uh , ∇uh ) − bh (., uh , ∇uh )]kLpπ (Q)
(h2K + τj )knE · [a(., uh , ∇uh ) − ah (., uh , ∇uh )]E kLpπ (∂QL ) ,
ηQ,π :=(h2K + τj )k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (Q) 1
−1
1
−1
π + hK
+ τjp
(4.14)
(h2K + τj )knE · [ah (., uh , ∇uh )]E kLpπ (∂QL )
(h2K + τj )kuh (., tj + 0) − uh (., tj − 0)kLπ (K) .
The quantity εQ,π obviously measures the quality of the approximation of a and b by ah and bh respectively, and can be estimated explicitly. Below, we will show that it yields upper bounds on the second terms on the right-hand sides of estimates (2.3) and (2.4). Note, that in our second example
εQ,π ≤ h2K kf − π0,Q f kLpπ (Q) + h3K k∇uh kLpπ (Q) , 1,0 if k = 0 and Vj = Sj,0 , 1 ≤ j ≤ Nτ .
Estimation of k(IdY+ −Rh )∗ F˜h (uh )kY+∗ and k(IdY+ −Rh )∗ [F (uh )−F˜h (uh )]kY+∗ . Next, we will derive upper bounds for the first and second terms on the right-hand side of inequality (2.3). Recalling equation (4.11) and using, for the space variables, integration by parts elementwise, we obtain for all ϕ ∈ Y < F (uh ), ϕ >Y Nτ n X Z X { [uh (x, tj + 0) − uh (x, tj − 0)]ϕ(x, tj )dx = j=1 K∈Tj
K
+ +
X Z
E∈Ej
Z
tj+1
tj tj+1
tj
Z
Z
E
K
[∂t uh − ∇ · a(x, uh , ∇uh ) − b(x, uh , ∇uh )]ϕdxdt}
nE · [a(x, uh , ∇uh )]E ϕdsdt 26
o
(4.15)
and < F˜h (uh ), ϕ >Y Nτ n X Z X { [uh (x, tj + 0) − uh (x, tj − 0)]ϕ(x, tj )dx = j=1 K∈Tj
K
+ X Z
+
E∈Ej
Z
tj+1
tj
tj+1
tj
Z
Z
K
E
[∂t uh − ∇ · ah (x, uh , ∇uh ) − bh (x, uh , ∇uh )]ϕdxdt}
o nE · [ah (x, uh , ∇uh )]E ϕdsdt .
(4.16) Let ϕ ∈ Y+ = × (0, T )) be arbitrary. From Lemma 3.4 we obtain for Q = K × Ij ∈ Pτ , 1 ≤ j ≤ Nτ , K ∈ Tj , and E ⊂ ∂K\Γ the estimates Z [uh (x, tj + 0) − uh (x, tj − 0)][ϕ(x, tj ) − Iτ ϕ(x, tj )]dx ′ Wπp′ (Ω
K
≤ kuh (., tj + 0) − uh (., tj − 0)kLπ (K) kϕ − Iτ ϕkLp′ (∂Q
B)
π′
1
¹ τjp Z tj+1 Z tj
−1
(h2K + τj )kuh (., tj + 0) − uh (., tj − 0)kLπ (K) kϕkW p′ (K×U (t
j ))
π′
K
,
[∂t uh − ∇ · ah (x, uh , ∇uh ) − bh (x, uh , ∇uh )][ϕ − Iτ ϕ]dxdt
≤ k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (Q) kϕ − Iτ ϕkLp′ (Q) π′
¹ (h2K + τj )k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (Q) kϕkW p′ (U (Q)) , π′ Z tj+1 Z nE · [ah (x, uh , ∇uh )]E [ϕ − Iτ ϕ]dsdt tj
E
≤ knE · [ah (., uh , ∇uh )]E kLpπ (∂QL ) kϕ − Iτ ϕkLp′ (∂Q π′
1
¹ hEπ
−1
L)
(h2K + τj )knE · [ah (., uh , ∇uh )]E kLpπ (∂QL ) kϕkW p′ (U (∂Q
L ))
π
.
Inserting these estimates in (4.16), using H¨older’s inequality for finite sums, and recalling definition (4.14), we conclude that < F˜h (uh ), ϕ − Iτ ϕ >
Nτ Nτ X X X X ′ π p/π 1/p ¹{ { ηQ,π } } { { kϕkπ p′ j=1 K∈Tj
j=1 K∈Tj
′
Wπ′ (U (K×Jj ))
′
′
}p /π }1/p .
Assume that p′ ≤ π ′ or, equivalently, p ≥ π. Then, Jensen’s inequality implies that Nτ X X ′ { { kϕkπ p′ j=1 K∈Tj
Wπ′ (U (K×Jj ))
′
′
′
}p /π }1/p ¹ kϕkW p′ (Ω×(0,T )) . π′
27
Using the abbreviation Nτ X X π η := { { ηQ,π }p/π }1/p
(4.17)
j=1 K∈Tj
we have thus shown that k(IdY+ − Rh )∗ F˜h (uh )kY+∗ ¹ η.
(4.18)
Replacing ah and bh by ah − a and bh − b, respectively, we conclude with the same arguments that k(IdY+ − Rh )∗ [F˜h (uh ) − F (uh )]kY+∗ ¹ ε (4.19) where
Nτ X X ε := { { επQ,π }p/π }1/p .
(4.20)
j=1 K∈Tj
Estimation of kF˜h (uh )kY˜ ∗ and kF˜h (uh ) − F (uh )kY˜ ∗ . Now, we will bound the h h terms in inequality (2.4). Given a subset Q of Ω × (0, T ), we set for abbreviation Y˜h|Q := {ϕ ∈ Y˜h : suppϕ ⊂ Q}. In order to bound the second term on the right-hand side of estimate (2.4), we conclude from the shape regularity of the partitions and a standard scaling argument that the estimate 1
π′ kϕkLp′ (K ′ ×J ′ ) + hK kϕkLp′ (∂K ′ ×J ′ ) ¹ h2K kϕkLp′ (J ′ ,W 2,π′ (K ′ )) π′
π′
holds for all ϕ ∈ Y˜h|U (Q) , Q = K × Jj ∈ Pτ , J ′ ⊂ U (Jj ), K ′ ⊂ U (K). Combining this with equations (4.14) – (4.16), and using H¨older’s inequality, we obtain the following estimate for all Q ∈ Pτ X kF (uh ) − F˜h (uh )kY˜ ∗ ¹ εQ′ ,π , (4.21) h|U (Q)
Q′ ⊂U (Q)
where Y˜h equipped with the norm of Y+ . In order to derive lower bounds for the left-hand side of estimate (2.4), consider ˜ m and equation (4.14) an arbitrary Q = K × Jj ∈ Pτ . From Lemma 3.5 with VQˆ = IP we then obtain k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (Q) Z −1 ¹ sup kvk p′ [∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )]ψj ψK v ˜m v∈IP
= sup ˜m v∈IP
Lπ′ (Q)
kvk−1p′ Lπ′ (Q)
≤kF˜h (uh )kY˜ ∗
h|Q
Q
< F˜h (uh ), ψj ψK v >Y+
sup kvk−1 kψj ψK vkY+ Lp (Q)
˜m v∈IP
π′
−1 ˜ ¹{h−2 K + τj }kFh (uh )kY˜ ∗ . h|Q
28
(4.22)
˜ ˆ ˆ yield for all Estimate (4.22), equation (4.16) and Lemma 3.5 with V∂ Qˆ L = IP m|E×J E ⊂ ∂K\Γ knE · [ah (., uh , ∇uh )]E kLpπ (∂QL ) Z −1 ¹ sup kσk p′ nE · [ah (., uh , ∇uh )]E ψj ψE PE σ Lπ′ (∂QL )
˜ m|∂Q σ∈IP L
≤
Lπ′ (∂QL )
˜ m|∂Q σ∈IP L
− ≤
kσk−1p′
sup
X Z
˜ m|∂Q σ∈IP L
K′
kσk−1p′
sup +
Jj
K ′ ⊂ωE
Z
K ′ ⊂ωE
{< F˜h (uh ), ψj ψE PE σ >Y+
[∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )]ψj ψE PE σdxdt}
Lπ′ (∂QL )
X
∂QL
{kF˜h (uh )kY˜h|ω
E ×Jj
kψj ψE PE σkY+
k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (K ′ ×Jj )
· kψj ψE PE σkLp′ (K ′ ×J ) } j
π′
1 π′
−1 ˜ ¹hE {h−2 E + τj }kFh (uh )kY˜ ∗
h|ωE ×Jj
.
From estimate (4.22), equation (4.16) and Lemma 3.5 with V∂ Qˆ B further conclude that
(4.23) ˜ = IPm|K×{0} we ˆ
kuh (., tj + 0) − uh (., tj − 0)kLπ (K) Z −1 kwkLπ′ (K) ¹ sup [uh (., tj + 0) − uh (., tj − 0)]ψj ψK Pj wdx ˜ m|K×{t } w∈IP j
=
sup ˜ m|K×{t } w∈IP j
− ≤
{< F˜h (uh ), ψj ψK Pj w >Y+ kwk−1 Lπ′ (K)
X Z Z
J⊂U (tj )
sup ˜ m|K×{t } w∈IP j
+
K
X
J⊂U (tj )
J
K
[∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )]ψj ψK Pj wdxdt}
kwk−1 {kF˜h (uh )kY˜ ∗ Lπ′ (K)
h|K×U (tj )
kψj ψK Pj wkY+
k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (K×J)
· kψj ψK Pj wkLp′ (K×J) } π′
1 p′
−1 ˜ ¹τj {h−2 K + τj }kFh (uh )kY˜ ∗
h|K×U (tj )
. (4.24)
Estimates (4.22) – (4.24) and definition (4.14) yield −1 2 ˜ ηQ,π ¹ {1 + τj h−2 K + τj hK }kFh (uh )kY˜ ∗
h|U (Q)
29
(4.25)
where Y˜h is endowed with the norm of Y+ . A posteriori error estimates. Combining estimates (4.18), (4.19), (4.21), and (4.25) with Propositions 2.1, 2.2 and recalling that kRh∗ [F (uh ) − Fh (uh )]kY+∗ = 0, we obtain the following result. 4.2 Proposition. Let u be a regular solution of problem (4.1) in the sense of Proposition 2.1 and definition (4.2) and let uh be an approximation of u in the sense of Proposition 2.2 and definition (4.5). Suppose that p ≥ π. Then the following a posteriori error estimates hold: ku − uh kLrρ (Ω×(0,T )) ¹ {η + ε + kRh∗ Fh (uh )kY+∗ } and
ηQ,π ¹{1 + τj−1 h2K + τj h−2 K }{ku − uh kLrρ (U (Q)) + ∀Q = K × Jj ∈ Pτ .
X
Q′ ⊂U (Q)
εQ′ ,π }
The quantities εQ,π , ηQ,π , η, and ε are given by equations (4.14), (4.17), and (4.20). 4.3 Remark. The local lower bounds for ku − uh kLrρ (Ω×(0,T )) can be combined in the standard way to the global lower bound η ¹ max max {1 + τj−1 h2K + τj h−2 K } · {ku − uh kLrρ (Ω×(0,T )) + ε}. 1≤j≤Nτ K∈Tj
The factor 1 + τj−1 h2K + τj h−2 K in this estimate and the second one of Proposition 4.2 reflects the fact that the differential operator is of 2nd order with respect to the space variables but only of 1st order with respect to the time variable. Remark. If p < π one may still obtain upper bounds on the error. Since, in case, Jensen’s inequality cannot be used in estimating < F˜h (uh ), ϕ − Iτ ϕ >Y+ < F˜h (uh ) − F (uh ), ϕ − Iτ ϕ >Y+ , one must now proceed as follows: Bound the space-integrals by using H¨ older’s inequality and Lemma 3.2. On each time-level sum up all contributions to that level and apply H¨older’s inequality for finite sums. 3. Bound the remaining time-integrals by using H¨older’s inequality and Lemma 3.3. 4. Sum up over all time-levels and use H¨older’s inequality for finite sums.
4.4 this and 1. 2.
5. Non-stationary, incompressible Navier-Stokes equations Variational setting. In this section we consider the two- and three dimensional 30
non-stationary, incompressible Navier-Stokes equations ∂t u − ν∆u + (u · ∇)u + ∇p = f
in Ω × (0, T )
in Ω × (0, T )
∇·u=0 u=0
u(., 0) = u0
on Γ × (0, T )
(5.1)
in Ω.
Here, ν > 0 is the constant viscosity of the fluid; u and p denote its velocity and pressure, respectively. In order to cast problem (5.1) into the abstract framework of Section 2, we introduce some auxiliary spaces Z 2 2 p = 0}, L0 (Ω) := {p ∈ L (Ω) : Ω ½ L∞ (0, T ; L2 (Ω)2 ) ∩ L2 (0, T ; W01,2 (Ω)2 ) , if n = 2, X1 := 1,2 ∞ 2 3 2 3 8 4 3 L (0, T ; L (Ω) ) ∩ L (0, T ; W0 (Ω) ) ∩ L (0, T ; L (Ω) ) , if n = 3, X2 := L2 (0, T ; L20 (Ω)),
Y1 := WT2 (0, T ; W01,2 (Ω)n , W −1,2 (Ω)n ), Y2 := L2 (0, T ; L20 (Ω)), H1 := L2 (0, T ; L2 (Ω)n ), H2 := L2 (0, T ; (W 1,2 (Ω) ∩ L20 (Ω))∗ ),
Z1 := WT2 (0, T ; W 2,2 (Ω)n ∩ W01,2 (Ω)n , L2 (Ω)n ) Z2 := L2 (0, T ; W 1,2 (Ω) ∩ L20 (Ω)),
and set X :=X1 × X2 , Y :=Y1 × Y2 ,
X− :=H1 × H2 ,
X+ :=L2 (0, T ; W01,3 (Ω)n ) × X2 ,
Y+ :=Z1 × Z2 Z Z u0 (x)v(x, 0)dx − < F ([u, p]), [v, q] >Y := − +
Ω Z T 0
+
Z
0
T
Z
0
T
Z
(5.2) u∂t vdxdt
Ω
Ω
{ν∇u · ∇v + [(u · ∇)u]v − p∇ · v − f v}dxdt
Ω
∇ · uqdxdt.
Z
One easily checks, that every solution [u, p] of problem (1.1) is a weak solution of problem (5.1) in the sense of [9; Problem III.3.2] if ∂t u ∈ L1 (0, T ; W −1,2 (Ω)n ). 31
Conversely, every weak solution [u, p] of problem (5.1) in the sense of [9] is a solution of problem (1.1) if u ∈ X. From this and Theorems III.3.1, III.3.2, and III.3.4 in [9] it follows that problem (1.1) with the definition (5.2) has exactly one solution, if n = 2, and at most one solution, if n = 3. Since Z1 ⊂ L∞ (0, T ; W01,2 (Ω)n ) ⊂ L∞ (0, T ; L6 (Ω)n ) with continuous imbeddings (cf. [6] Vol. 5, Chap. XVIII, §1, Theorem 1), one easily checks that condition (2.1) is fulfilled if the solution of problem (1.1) is contained in X+ . Moreover, DF ([u, p])∗ ∈ ∗ Isom(Y+ , X− ) if and only if the adjoint linearized Navier-Stokes equations −∂t v − ν∆v − (u · ∇)v + v∇u − ∇q = w
in Ω × (0, T )
v=0
on Γ × (0, T )
−∇ · v = r
u(., T ) = 0
in Ω × (0, T ) in Ω
∗ admit for each [w, r] ∈ X− a unique weak solution [v, q] ∈ Y+ satisfying
k[v, q]kY+ ¹ k[w, r]kX−∗ . Standard energy estimates show that this is true if u ∈ L2 (0, T ; W 2,2 (Ω)n ∩W01,2 (Ω)n ) and kukL2 (0,T ;W 2,2 (Ω)n ) is sufficiently small, e.g., if ν −2 kf kL22 (Ω×(0,T )) < 1.
Finite element discretization. For the discretization of problem (5.1) we proceed as in Sections 3 and 4. We choose a family Iτ of shape regular partitions of the interval [0, T ]. With each time tj , 1 ≤ j ≤ Nτ , we associate an admissible and shape regular partition Tj of Ω into n-simplices and two finite element spaces Mj ⊂ W01,2 (Ω)n and Qj ⊂ L20 (Ω) corresponding to Tj and consisting of affine equivalent elements in the sense of [4]. In what follows we always assume that there are two 1,0 n µ,0 n integers µ ≥ 1 and ν ≥ 0 such that, for all 1 ≤ j ≤ Nτ , [Sj,0 ] ⊂ Mj ⊂ [Sj,0 ] and 1,0 ν,0 0,−1 ν,−1 either Sj ⊂ Qj ⊂ Sj or Sj ⊂ Qj ⊂ Sj . We then set for k ≥ 0 and θ ∈ [0, 1] Xh := Sτk,−1 [Mh(τ ) ] × Sτk,−1 [Qh(τ ) ],
Yh := Sτθ;k+1,0 [Mh(τ ) ] × Sτk,−1 [Qh(τ ) ],
< Fh ([uh , ph ]), [v h , qh ] >Yh :=< F ([uh , ph ]), [v h , qh ] >Y
(5.3)
∀[uh , ph ] ∈ Xh , [v h , qh ] ∈ Yh . Note, that the above discretization is consistent, i.e. kRh∗ [F ([uh , ph ])−Fh ([uh , ph ])]kY+∗ = 0 and that for the pressure the test- and trial-space are identical. As in Section 4 one easily checks that this discretization corresponds to a (k + 1)-stage Runge-Kutta scheme. For k = 0, in particular, it yields the θ-scheme. Definition of Rh , F˜h , and Y˜h . In order to cast the discretization (5.3) in the framework of Section 2, we now define F˜h , Y˜h and Rh . Denote by π0,Q f the L2 projection of f onto the space of piecewise constant functions corresponding to the 32
partition Pτ of Ω × (0, T ). Then, F˜h is defined in the same way as F with f replaced by π0,Q f . Next, we set Y˜h :=
span{ψj ψK v, ψj ψE PE σ, ψj ψK Pj w : ˜ m |K×J , σ ∈ IP ˜ m |E×J , w ∈ IP ˜ m |K×{t } , v ∈ IP 1 j 2 j 3 j
1 ≤ j ≤ Nτ , K ∈ Tj , E ∈ Ej } ˜ m |K×J , 1 ≤ j ≤ Nτ , K ∈ Tj } ×span{χK×Jj q : q ∈ IP 4 j where
m1 := max{2k + 2µ + 1, k + ν − 1} , m2 := max{k + µ − 1, k + ν}, , m4 := k + µ − 1,
m3 := k + µ
˜m and χK×Jj denotes the characteristic function of the set K × Jj . Recall that IP denotes space of polynomials (in x and t) of degree at most m. In order to define Rh , we set, for 1 ≤ j ≤ Nτ , ( 0,−1 ∩ L20 (Ω) if Qj ⊂ 6 C(Ω), ˆ j := Sj Q 1,0 2 Sj ∩ L0 (Ω) if Qj ⊂ C(Ω), ˆ j either the L2 -projection of L2 (Ω) onto Q ˆ j , if Qj 6⊂ C(Ω), or the and denote by Π 0 Cl´ement interpolation operator Ij , if Qj ⊂ C(Ω). Let π ˆj , 1 ≤ j ≤ Nτ , be the L2 projection of L2 (Jj ; L20 (Ω)) onto span{tµ qµ : 0 ≤ µ ≤ k, qµ ∈ L20 (Ω)}. We then define ˆ h(τ ) ] by the operator Iˆτ : L2 (0, T ; L20 (Ω)) → Sτk,−1 [Q Iˆτ q :=
Nτ X
ˆ j q, χJj π ˆj Π
j=1
and set Rh := (Iτ , Iˆτ ), where Iτ is defined in (3.10). The proofs of Lemmas 3.2 and 3.3 can easily be modified ˆ j and π such that one obtains similar approximation results for the operators Π ˆj as for Ij and πj , respectively. The estimators. For Q = K × Jj ∈ Pτ we define εQ :=(h2K + τj )kf − π0,Q f kL22 (Q) , − 21
ηQ :=τj
(h2K + τj )kuh (., tj + 0) − uh (., tj − 0)kL22 (∂QB )
+ (h2K + τj )k∂t uh − ν∆uh + (uh · ∇)uh − ∇ph − π0,Q f kL22 (Q) −1
+ hK 2 (h2K + τj )k[νnE · ∇uh − ph nE ]E kL22 (∂QL ) + hK k∇ · uh kL22 (Q)
33
(5.4)
and ε := {
X
Q∈Pτ
ε2Q }1/2
,
η := {
X
Q∈Pτ
2 1/2 ηQ } .
(5.5)
As in Section 4, we will show that εQ yields upper bounds on the second terms on the right-hand sides of estimates (2.3) and (2.4). Estimation of k(IdY − Rh )∗ F˜h ([u , ph ])kY ∗ and k(IdY − Rh )∗ [F ([u , ph ]) − +
h
+
+
h
F˜h ([uh , ph ])]kY+∗ . Using integration by parts elementwise and the convention (4.10), we obtain, for all [v, q] ∈ Y , the following representation of the residual < F˜h ([uh , ph ]), [v, q] >Y Nτ n X Z X = { {uh (x, tj + 0) − uh (x, tj − 0)}v(x, tj )dx K
j=1 K∈Tj
Z
+ +
X Z
E∈Ej
+
tj
tj+1 Z
tj
X Z
K∈Tj
Z
tj+1
E
tj+1
tj
K
Z
{∂t uh − ν∆uh + (uh · ∇)uh + ∇ph − π0,Q f }vdxdt}
(5.6)
[νnE · ∇uh − ph nE ]E · vdsdt
K
o ∇ · uh qdxdt .
Let [v, q] ⊂ Y+ be arbitrary. Since, on each sub-interval Jj , 1 ≤ j ≤ Nτ , ∇ · uh ∈ span{tµ qµ : 0 ≤ µ ≤ k, qµ ∈ L20 (Ω)}, we conclude from the definition of Iˆτ that for all Q = K × Jj , 1 ≤ j ≤ Nτ , K ∈ Tj , Z ∇ · uh (q − Iˆτ q)dxdt Q Z ˆ j q))dxdt ˆj q + π ˆj (q − Π ∇ · uh (q − π = Q Z ˆ j q)dxdt ∇ · uh π = ˆj (q − Π Q
ˆ j q)kL2 (Q) ≤k∇ · uh kL22 (Q) kˆ πj (q − Π 2 ˆ j qkL2 (Q) ≤k∇ · uh kL22 (Q) kq − Π 2
¹hK k∇ · uh kL22 (Q) kqkL2 (Jj ;W 1,2 (K)) .
The terms in < F˜h ([uh , ph ]), [v, q] − Rh [v, q] >Y involving v can be estimated as in Section 4. This proves that k(IdY+ − Rh )∗ F˜h ([uh , ph ])kY+∗ ¹ η. 34
(5.7)
With the same notations and arguments, we obtain < F ([uh , ph ]) − F˜h ([uh , ph ]), [v, q] − Rh [v, q] >Y Nτ X X Z tj+1 Z = { [f − π0,Q f ][v − Iτ v]dxdt} j=1 K∈Tj
tj
K
Nτ X X { ¹ kf − π0,Q f kL22 (K×Jj ) (h2K + τj )kvkW22 (U (K×Jj )) } j=1 K∈Tj
¹
X
Q∈Pτ
εQ kvkW22 (U (Q))
¹εkvkW22 (Ω×(0,T )) . This proves that k(IdY+ − Rh )∗ [F ([uh , ph ]) − F˜h ([uh , ph ])]kY+∗ ¹ ε.
(5.8)
Estimation of kF˜h ([uh , ph ])kY˜ ∗ and kF ([uh , ph ])− F˜h ([uh , ph ])kY˜ ∗ . The second h h term on the right-hand side of estimate (2.4) can be bounded in the same way as in Section 4. In analogy to estimate (4.21), we now obtain for all Q ∈ Pτ kF ([uh , ph ]) − F˜h ([uh , ph ])kY˜ ∗
h|U (Q)
¹
X
εQ′ .
(5.9)
Q′ ⊂U (Q)
Thanks to the definition of Yh and Y˜h , we may insert [0, χQ ∇ · uh ], Q = K × Jj , 1 ≤ j ≤ Nτ , K ∈ Tj , as a test function in equation (5.6). From a standard scaling argument we then obtain k∇ · uh kL22 (Q)
=k∇ · uh k−1 < F˜h ([uh , ph ]), [0, χQ ∇ · uh ] >Y L2 (Q) 2
k∇ · uh kL2 (Jj ;W 1,2 (K)) ≤kF˜h ([uh , ph ])kY˜ ∗ k∇ · uh k−1 L2 (Q) h|Q
(5.10)
2
˜ ¹h−1 ∗ . K kFh ([uh , ph ])kY˜h|Q The remaining terms in ηQ can be bounded in the same way as in Section 4. We therefore obtain ˜ ηQ ¹ {1 + τj−1 h2K + τj h−2 K }kFh ([uh , ph ])kY˜ ∗
h|U (Q)
.
(5.11)
A posteriori error estimates. Combining estimates (5.7) – (5.9), and (5.11) with Propositions 2.1, 2.2 and recalling that the discretization (5.3) is consistent, we obtain the following result. 35
5.1 Proposition. Let [u, p] be a regular solution of problem (5.1) in the sense of Proposition 2.1 and definition (5.2) and let [uh , ph ] be an approximation of [u, p] in the sense of Proposition 2.2 and definition (5.3). Then the following a posteriori error estimates hold: ku − uh kL2 (0,T ;L2 (Ω)n ) + kp − ph kL2 (0,T ;W −1,2 (Ω))
¹{η + ε + kRh∗ Fh ([uh , ph ])kY+∗ },
ηQ ¹{1 + τj−1 h2K + τj h−2 K }{ku − uh kL2 (U (Jj );L2 (U (K))n ) X + kp − ph kL2 (U (Jj );W −1,2 (U (K))) + εQ′ } Q′ ⊂U (Q)
∀Q = K × Jj ∈ Pτ ,
and
η ¹ max max {1 + τj−1 h2K + τj h−2 K }{ku − uh kL2 (0,T ;L2 (Ω)n ) 1≤j≤Nτ K∈Tj
+ kp − ph kL2 (0,T ;W −1,2 (Ω)) + ε}. The quantities εQ , ηQ , ε, and η are given by equations (5.4) and (5.5). 5.2 Remark. Proposition 5.1 can be extended to the case of the slip boundary condition u · n = T (νu, p) − [n · T (νu, p) · n]n = 0, where
1 T (u, p) := ( (∂i uj + ∂j ui ) − pδij )1≤i,j≤n 2 denotes the stress tensor. One then has to replace ν∇uh − ph I in equation (5.4) by T (νuh , ph ), and Γ by the part of the boundary on which the no-slip condition u = 0 is imposed. Here, I := (δij )1≤i,j≤n denotes the unit tensor. Other boundary conditions may be treated similarly. 5.3 Remark. The results of Propositions 4.2 and 5.1 can be combined to handle non-Newtonian fluids for which the viscosity is a given function of the velocity and its gradient (cf. [3, 10] for the stationary case).
References [1] R. A. Adams: Sobolev spaces. Academic Press, New York, 1975 [2] H. Amann: Nonhomogeneous linear and quasilinear elliptic and parabolic boundary value problems. Preprint, Universit¨ at Z¨ urich, 1994 36
[3] J. Baranger and H. El Amri: Estimateurs a posteriori d’erreur pour le calcul adaptif d’´ecoulements quasi-Newtoniens. RAIRO M2 AN 25, 31 – 48 (1991) [4] P. G. Ciarlet: The finite element method for elliptic problems. North Holland, Amsterdam, 1978 [5] P. Cl´ement: Approximation by finite element functions using local regularization. RAIRO Anal. Num´er. 9 R-2, 77 – 84 (1975) [6] R. Dautray and J.-L. Lions: Mathematical analysis and numerical methods for science and technology. Springer, Berlin - Heidelberg - New York, 1992 [7] K. Eriksson and C. Johnson: Adaptive finite element methods for parabolic problems IV. Nonlinear problems. Chalmers University of G¨ oteborg, Preprint 1992:44 (1992) : Adaptive finite element methods for parabolic problems V. Long-time [8] integration. Chalmers University of G¨ oteborg, Preprint 1993:04 (1993) [9] R. Temam: Navier-Stokes Equations. 3rd edition, North Holland, Amsterdam, 1984 [10] R. Verf¨ urth: A posteriori error estimates for nonlinear problems. Finite element discretizations of elliptic equations. Math. Comput. 62 (206), 445 – 475 (1994)
37
Summary: Using the abstract framework of [10] we analyze a residual a posteriori error estimator for space-time finite element discretizations of parabolic pdes. The estimator gives global upper and local lower bounds on the error of the numerical solution. The finite element discretizations in particular cover the so-called θ-scheme, which includes the implicit and explicit Euler methods and the Crank-Nicholson scheme. As particular examples we consider scalar quasilinear parabolic pdes of 2nd order and the time-dependent incompressible Navier-Stokes equations. Key words: A posteriori error estimates; nonlinear parabolic pdes; time-dependent Navier-Stokes equations; space-time finite elements; θ-scheme. AMS Subject Classification: 65N30, 65N15, 65J15, 76D05
1. Introduction We analyze a residual a posteriori error estimator for space-time finite element discretizations of parabolic pdes. Each space-time element K × J contributes the weighted sum of three terms: (1) The residual of the computed numerical solution with respect to the strong form of the differential operator evaluated on K × J, (2) the jump across ∂K × J of that trace operator which naturally connects the strong and the weak formulation of the differential equation, and (3) the jump of the numerical solution across K × ∂J. Here, K stands for an arbitrary element of the spatial mesh and J denotes an arbitrary interval of the time mesh. We could also extend our analysis to error estimators which are based on the solution of auxiliary local time-dependent problems. We do not follow this line here, in order not to overload the presentation. In order to construct our a posteriori error estimator and to prove that it yields upper and lower bounds on the error, we use the techniques introduced in [10] and consider in Section 2 abstract nonlinear problems of the form F (u) = 0 1
(1.1)
and corresponding discretizations of the form Fh (u0h ) = 0.
(1.2)
Here, F ∈ C 1 (X, Y ∗ ) and Fh ∈ C(Xh , Yh∗ ), Xh ⊂ X and Yh ⊂ Y are finite dimensional subspaces of the Banach spaces X and Y , and ∗ denotes the dual of a Banach space. If u0 ∈ X is a solution of problem (1.1) such that DF (u0 ) is an isomorphism of X onto Y ∗ and DF is Lipschitz continuous at u0 , we know from Proposition 2.1 in [10] that ckF (u)kY ∗ ≤ ku − u0 kX ≤ ckF (u)kY ∗ (1.3) holds for all u in a suitable neighbourhood of u0 . The constants c and c depend on DF (u0 ) and DF (u0 )−1 . They measure the sensitivity of the infinite dimensional problem (1.1) with respect to small perturbations. For a simple model problem we derive explicit bounds for c and c in Section 4. When applying estimate (1.3) to an approximate solution uh ∈ Xh of problem (1.2) one must evaluate the residual kF (uh )kY ∗ . This is as expensive as the solution of the original problem (1.1) since it amounts in the solution of an infinite dimensional maximization problem. In order to obtain an approximation of kF (uh )kY ∗ which is much easier to evaluate, we introduced in [10] an approximation F˜h : Xh → Y ∗ of F at uh , a finite dimensional subspace Y˜h ⊂ Y , and a restriction operator Rh : Y → Yh . For practical applications, Rh is a suitable local quasi-interpolation operator (cf. (3.10)), Feh (uh ) is obtained by locally projecting F (uh ) onto a suitable finite element space, and Yeh consists of judiciously chosen local cut-off functions (cf. Section 3, in particular Lemma 3.5). When applying these results to parabolic pdes, we obtain error estimates in a L (0, T ; W01,ρ (Ω))-norm (cf. Section 3 for a difinition of these spaces and norms). The ′ ′ space Y then consists of functions in Lp (0, T ; W01,π (Ω)) having their time derivative ′ ′ ′ in Lp (0, T ; W −1,π (Ω)). Due to the non-local nature of the space W −1,π (Ω) we get into troubles when deriving lower bounds on the error. This problem will be tackled in a subsequent paper. Here, we circumvent this difficulty by imposing a weaker ′ ′ ′ Lr (0, T ; Lρ (Ω))-norm on X and a stronger Lp (0, T ; W 2,π (Ω)∩W01,π (Ω))-norm on Y . The corresponding spaces will be denoted by X− and Y+ , respectively. In particular, ′ ′ the functions in Y+ now have time derivatives in Lp (0, T ; Lπ (Ω)). r
In Sections 4 and 5 we apply the general results of the previous sections to scalar quasilinear parabolic pdes of 2nd order and the time-dependent incompressible Navier-Stokes equations. We consider Petrov-Galerkin space-time finite elements. These discretizations are inspired by the concept of very weak solutions. The trial space Xh consists of functions which are discontinuous in time, piecewise polynomials of degree k, k ≥ 0. The test space Yh consists of functions which are continuous 2
in time, piecewise polynomials of degree k + 1 and which vanish at the final time T . This discretization corresponds to an implicit (k + 1)-stage Runge-Kutta-scheme. When applied to a linear problem its stability function is the (k + 1)-st diagonal Pad´e-approximation. For k = 0 we in particular obtain the Crank-Nicholson scheme. By slightly modifying the basis-functions of Yh we may also recover the popular θscheme for all θ ∈ [0, 1]. This in particular covers the explicit (θ = 0) and implicit (θ = 1) Euler schemes. We obtain global upper and local lower bounds for the error measured in an L (0, T ; Lρ (Ω))-norm. The upper and lower bounds differ by a factor 1+τ −1 h2 +τ h−2 . Here, h and τ are the local mesh-sizes in space and time, respectively. This factor reflects the fact that the differential operator is of 2nd order with respect to the space variables but only of 1st order with respect to the time variable. The local lower bounds may be combined to a global lower bound of the same type. When applied to the corresponding particular examples, our error estimates are similar to those obtained in [7, 8] using a different approach and a different discretization. r
2. Abstract Error Estimates Let X, Y be two Banach spaces with norms k.kX and k.kY . For any element u ∈ X and any real number R > 0 set BX (u, R) := {v ∈ X : ku − vkX < R}. We denote by L(X, Y ) and Isom(X, Y ) ⊂ L(X, Y ) the Banach space of continuous linear maps of X in Y equipped with the operator norm k.kL(X,Y ) and the open subset of linear homeomorphisms of X onto Y . By Y ∗ := L(Y, IR) and < ., . >Y we denote the dual space of Y and the corresponding duality pairing. Finally, A∗ ∈ L(Y ∗ , X ∗ ) denotes the adjoint of a given operator A ∈ L(X, Y ).
Let F ∈ C 1 (X, Y ∗ ) be a given continuously differentiable function. Given a solution u0 ∈ X of problem (1.1) and an arbitrary element u ∈ X ”close” to u0 , we may estimate the error ku − u0 kX by the residual kF (u)kY ∗ (cf. Proposition 2.1 in [10]). For parabolic pde’s we thus obtain control on the Lr (0, T ; W01,ρ (Ω))-norm of the error (cf. Section 3 for the definition of these spaces and norms). However, we are interested in controlling the Lr (0, T ; Lρ (Ω))-norm of the error. In order to achieve this within the present abstract framework, we must enlarge the space X and reduce the space Y . We therefore consider three additional Banach spaces X− , X+ , and Y+ such that X+ ⊂ X ⊂ X− and Y+ ⊂ Y with continuous and dense injections. Here, the +/− sign indicates a space with a stronger / weaker norm. We assume that X− is reflexive. 3
2.1 Proposition. Let u0 ∈ X be a solution of problem (1.1). Assume that u0 ∈ X+ , ∗ that DF (u0 )∗ ∈ Isom(Y+ , X− ), and that there are two numbers R0 > 0 and β > 0 such that k[DF (u0 ) − DF (u0 + tw)]wkY+∗ ≤ βtkwkX+ kwkX− (2.1) ∀w ∈ BX+ (0, R0 ), t ∈ [0, 1]. Set
−1
−1 R := min{R0 , β −1 kDF (u0 )∗ k−1 kDF (u0 )∗ kL(Y+ ,X−∗ ) }. L(X ∗ ,Y+ ) , 2β −
Then the following error estimate holds for all u ∈ BX+ (u0 , R) : 1 ∗−1 kDF (u0 )∗ k−1 kL(X−∗ ,Y+ ) kF (u)kY+∗ . ∗ ) kF (u)kY ∗ ≤ ku − u0 kX− ≤ 2kDF (u0 ) L(Y ,X + + − 2 (2.2) ∗ Proof. Let u ∈ BX+ (u0 , R). Consider an arbitrary element w ∈ X− and set ϕ := −1 ∗ DF (u0 ) w ∈ Y+ . We then have
< u − u0 , w >X− = < DF (u0 )(u − u0 ), ϕ >Y+
= < F (u), ϕ >Y+ Z 1 < [DF (u0 ) − DF (u0 + t(u − u0 ))](u − u0 ), ϕ >Y+ dt. + 0
−1
Inequality (2.1) and the continuity of DF (u0 )∗ imply that Z 1 < [DF (u0 ) − DF (u0 + t(u − u0 ))](u − u0 ), ϕ >Y+ dt| | 0
1 ≤ βku − u0 kX+ ku − u0 kX− kϕkY+ 2 −1 1 ≤ βkDF (u0 )∗ kL(X−∗ ,Y+ ) Rku − u0 kX− kwkX−∗ . 2 Combined with the above representation of < u − u0 , w >X− this yields −1 1 < u − u0 , w >X− ≤kF (u)kY+∗ kϕkY+ + βkDF (u0 )∗ kL(X−∗ ,Y+ ) Rku − u0 kX− kwkX−∗ 2 −1 1 ≤{kDF (u0 )∗ kL(X−∗ ,Y+ ) kF (u)kY+∗ + ku − u0 kX− }kwkX−∗ . 2 ∗ Since X− is reflexive and w ∈ X− was arbitrary, this implies the upper bound of estimate (2.2). In the same way, we obtain < F (u), ϕ >Y+ = < u − u0 , w >X− −
Z
0
1
< [DF (u0 ) − DF (u0 + t(u − u0 ))](u − u0 ), ϕ >Y+ dt
1 ≤ku − u0 kX− kwkX−∗ + βku − u0 kX+ ku − u0 kX− kϕkY+ 2 1 ≤{kDF (u0 )∗ kL(Y+ ,X−∗ ) ku − u0 kX− + βRku − u0 kX− }kϕkY+ 2 ∗ ≤2kDF (u0 ) kL(Y+ ,X−∗ ) ku − u0 kX− kϕkY+ . 4
Since ϕ ∈ Y+ is arbitrary, this proves the lower bound of estimate (2.2).
∗ The condition DF (u0 )∗ ∈ Isom(Y+ , X− ) of Proposition 2.1 is more restrictive ∗ than the assumption DF (u0 ) ∈ Isom(X, Y ) which is needed to bound ku − u0 kX by kF (u)kY ∗ (cf. Proposition 2.1 in [10]). For pdes, it is equivalent to an additional regularity condition. For linear problems, i.e. when DF is constant, one may extend F by continuity to a continuously differentiable map of X− to Y+∗ . Then the space X+ is not needed. For nonlinear problems, however, this extension often is impossible or the derivative of the extension is no longer Lipschitz continuous. This is the place where the space X+ comes into play.
Let Xh ⊂ X+ and Yh ⊂ Y be finite dimensional subspaces and Fh ∈ C(Xh , Yh∗ ) be an approximation of F . Given an approximate solution uh ∈ Xh of problem (1.2), Proposition 2.1 allows us to estimate the error ku0 −uh kX− by the residual kF (uh )kY+∗ . The evaluation of the latter, however, is a difficult task since it amounts in the solution of an infinite dimensional maximization problem. In order to obtain an approximation of the residual which is easier to compute, we introduce a restriction operator Rh ∈ L(Y, Yh ), a finite dimensional subspace Y˜h ⊂ Y+ , and an approximation F˜h : Xh → Y ∗ of F at uh . In the context of pdes F˜h is obtained by locally freezing the coefficients of the differential operator. We equip Yeh with the norm of Y+ .
2.2 Proposition. The following estimates hold:
kF (uh )kY+∗ ≤k(IdY+ − Rh )∗ F˜h (uh )kY+∗ + k(IdY+ − Rh )∗ [F (uh ) − F˜h (uh )]kY+∗ + kRh∗ [F (uh ) − Fh (uh )]kY+∗ + kRh∗ Fh (uh )kY+∗
(2.3)
and kF˜h (uh )kY˜ ∗ ≤ kF (uh )kY˜ ∗ + kF (uh ) − F˜h (uh )kY˜ ∗ . h
h
h
(2.4)
If there is a constant c0 , which does not depend on h, such that k(IdY+ − Rh )∗ F˜h (uh )kY+∗ ≤ c0 kF˜h (uh )kY˜ ∗ , h
(2.5)
then both k(IdY+ − Rh )∗ F˜h (uh )kY+∗ and kF˜h (uh )kY˜ ∗ yield upper and lower bounds for h the residual kF (uh )kY+∗ . Proof. Estimate (2.3) follows from the identity
< F (uh ), ϕ >Y+ = < F˜h (uh ), ϕ − Rh ϕ >Y+ + < F (uh ) − F˜h (uh ), ϕ − Rh ϕ >Y+ + < F (uh ) − Fh (uh ), Rh ϕ >Y+ + < Fh (uh ), Rh ϕ >Y+
which holds for all ϕ ∈ Y+ . Estimate (2.4) follows from the triangle inequality. The statement concerning upper and lower bounds for the residual is an obvious consequence of inequality (2.5). 5
A non-optimal estimate of the third and fourth term on the right-hand side of estimate (2.3) is given by kRh∗ [F (uh ) − Fh (uh )]kY+∗ + kRh∗ Fh (uh )kY+∗
≤kIdY+ kL(Y+ ,Y ) kRh kL(Y,Yh ) {kF (uh ) − Fh (uh )kYh∗ + kFh (uh )kYh∗ }. Here, Yh is equipped with the norm of Y . The second terms on the right-hand sides of estimates (2.3) and (2.4) measure the quality of the approximation F˜h (uh ) to F (uh ). Usually, they are higher order terms when compared with k(IdY+ − Rh )∗ F˜h (uh )kY+∗ and kF˜h (uh )k ˜ ∗ . The term kR∗ [F (uh ) − Fh (uh )]kY ∗ is the consistency error of the h
Yh
+
discretization. The term kRh∗ Fh (uh )kY+∗ measures the residual of the algebraic equation (1.2) and can be evaluated by standard methods. Inequality (2.5) is a non-trivial condition since it claims that a supremum with respect to an infinite dimensional space may be bounded from above by a supremum with respect to a finite dimensional space.
Propositions 2.1 and 2.2 will be used in the following way. Any upper bound for k(IdY+ − Rh )∗ F˜h (uh )kY+∗ yields a reliable error estimator. If in addition the error estimator yields a lower bound for kF˜h (uh )k ˜ ∗ , it is also efficient. Here, we adopt Yh
the standard convention that an error estimator is called reliable (resp. efficient) if it yields an upper (resp. lower) bound for the error. Moreover, upper and lower bounds may contain multiplicative constants which do not depend on the discretization parameter h. The lower bound of estimate (2.2) follows from the identity < F (u), ϕ >Y+ = < DF (u0 )(u − u0 ), ϕ >Y+ Z 1 < [DF (u0 ) − DF (u0 + t(u − u0 ))](u − u0 ), ϕ >Y+ dt − 0
with ϕ ∈ Y+ . In many applications Y+ will be a closed subspace of a suitable Sobolev space. When inserting in the above identity functions ϕ with a local support, one then obtains lower bounds for u − u0 restricted to the support of ϕ.
The factor kDF (u0 )∗ kL(Y+ ,X−∗ ) in the lower bound of estimate (2.2) is rather harmless since it corresponds to a differential operator which is local and the norm of which can easily be estimated in terms of its coefficients. On the other hand, −1 the factor kDF (u0 )∗ kL(X−∗ ,Y+ ) in the upper bound of estimate (2.2) is much more severe. In some applications, in particular when X− is a Hilbert space, a function ∗ w ∈ X− satisfying kwkX−∗ = 1
and hu − u0 , wiX− = ku − u0 kX− −1
can explicitly be given in terms of u and u0 . Then kDF (u0 )∗ kL(X−∗ ,Y+ ) can be −1
replaced by kDF (u0 )∗ wkY+ . The latter quantity may be estimated numerically by approximately solving a discrete analogue of the corresponding adjoint pde. 6
The Lipschitz condition of Proposition 2.1 may be replaced by weaker H¨older conditions of the form k[DF (u0 ) − DF (u0 + tw)]wkY+∗ ≤ βtα kwkX− kwkα X+ ∀w ∈ BX+ (0, R0 ), t ∈ [0, 1].
3. Auxiliary results Function spaces. Let Ω be a bounded, connected, open domain in IRn , n ≥ 2, with polyhedral boundary Γ. For any open subset ω of Ω with Lipschitz boundary γ, we denote by W k,p (ω), k ∈ IN, 1 ≤ p ≤ ∞, Lp (ω) := W 0,p (ω), and Lp (γ) the usual Sobolev and Lebesgue spaces equipped with the standard norms (cf. [1] and Vol. 3, Chap. IV in [6]). Let W01,p (Ω) := {u ∈ W 1,p (Ω) : u = 0 on Γ} and set for 1 < p < ∞
′
W −1,p (Ω) := W01,p (Ω)∗ .
Here, p′ denotes the dual exponent of p defined by p1′ + p1 = 1. In what follows, a prime will always denote the dual of a given Lebesgue exponent. Let V and W be two Banach spaces such that V ⊂ W with continuous and dense injection. Given two real numbers a and b with a < b, we denote by Lp (a, b; V ), 1 ≤ p < ∞, the space of measurable functions u defined on (a, b) with values in V such that the function t → ku(., t)kV is in Lp ((a, b)). Lp (a, b; V ) is a Banach space equipped with the norm Z b
kukLp (a,b;V ) := {
a
ku(., t)kpV dt}1/p
(cf. [6], Vol. 5, Chap. XVIII, §1). Slightly changing the notation of [6], we further consider the Banach space W p (a, b; V, W ) := {u ∈ Lp (a, b; V ) : ∂t u ∈ Lp (a, b; W )} equipped with the norm kukW p (a,b;V,W )
Z := {
a
b
ku(., t)kpV
dt +
Z
a
b
k∂t u(., t)kpW dt}1/p .
Here, the partial derivative ∂t u must be interpreted in the distributional sense (cf.[6]; loc. cit.). For all smooth functions ϕ ∈ D((a, b)) it satisfies the identity Z
a
b
∂t u(., t)ϕ(t)dt = − 7
Z
a
b
u(., t)ϕ′ (t)dt,
where the integrals are taken in W . From Proposition 9 in [6; loc. cit.] it follows that the traces u(., a) and u(., b) are defined as elements of W , provided p > 1. We therefore set for 1 < p < ∞ Wbp (a, b; V, W ) := {u ∈ W p (a, b; V, W ) : u(., b) = 0}. Given an interval I ⊂ IR, a suitable subset ω ⊂ Ω, and numbers 1 ≤ p, π < ∞, we will use the following abbreviation: Lpπ (ω × I) := Lp (I; Lπ (ω)),
Wπp (ω × I) := W p (I; W 2,π (ω), Lπ (ω)).
Finite element partition. Denote by T > 0 an arbitrary but fixed time. Let Iτ = {[tj , tj+1 ] : 1 ≤ j ≤ Nτ }, τ > 0, be a family of partitions of [0, T ] with 0 = t1 < t2 < ... < tNτ +1 = T. For 1 ≤ j ≤ Nτ set Jj := [tj , tj+1 ] , τj := tj+1 − tj . We assume that the family Iτ is shape regular, i.e., the ratios τj /τj+1 and τj+1 /τj are bounded from above independently of j and τ . With each 1 ≤ j ≤ Nτ , we associate a partition Tj of Ω into n-simplices. We denote by Ej the set of the interior faces of Tj . For K ∈ Tj and E ∈ Ej let hK , ρK , and hE be the diameter of K, the diameter of the largest ball inscribed into K, and the diameter of E. We assume that the partitions Tj satisfy the following two conditions: (1) admissibility: Any two simplices of Tj are either disjoint or share a complete smooth submanifold of their boundaries.
(2) shape regularity: The ratio hK /ρK is bounded from above independently of K ∈ Tj , j, and τ .
Condition (2) allows the use of locally refined meshes. It implies that the ratio hK /hE , for all K ∈ Tj , all faces E of K, and all j, is bounded from above and from below by constants which do not depend on K, E, j, and τ . Denote by Pτ := {K × Jj : 1 ≤ j ≤ Nτ , K ∈ Tj } the partition of the space-time cylinder Ω × [0, T ] into prisms which is induced by Iτ and the Tj ’s. For any E ∈ Ej , 1 ≤ j ≤ Nτ , and any piecewise continuous function u, we denote by [u]E the jump of u across E in an arbitrary, but fixed direction nE orthogonal to 8
E. Finally, we introduce the following neighbourhoods of elements and points U (tj ) := [tmax{1,j−1} , tmin{Nτ +1,j+1} ] , 1 ≤ j ≤ Nτ ,
U (Jj ) := U (tj ) ∪ U (tj+1 ) = [tmax{1,j−1} , tmin{Nτ +1,j+2} ] , 1 ≤ j ≤ Nτ , ′ U (K) := ∪K ′
, K ∈ Tj , 1 ≤ j ≤ Nτ ,
U (E) := ∪K′ ′
, E ∈ Ej , 1 ≤ j ≤ Nτ ,
K ∩K6=∅ K ′ ∈Tj
E∩K 6=∅ K ′ ∈Tj
U (P ) := ∪P′
′
P ∩P 6=∅ P ′ ∈Pτ
ωK :=
∪K ′
E⊂∂K K ′ ∈Tj
, P ∈ Pτ , , K ∈ Tj , 1 ≤ j ≤ Nτ ,
K ′ ∩K∈Ej K ′ ∈Tj
ωE := ∪K ′′
(3.1)
, E ∈ Ej , 1 ≤ j ≤ Nτ .
Here, K ∩K ′ ∈ Ej means that K and K ′ share a complete (n−1)-dimensional smooth submanifold of their boundaries. Finite element spaces. Denote by IPk , k ≥ 0, the space of polynomials (in x) of degree at most k. Given an admissible partition Th of Ω, we define the finite element spaces (in x) as usual: Shk,−1 := {u : Ω → IR : u|K ∈ IPk
Shk,0 k,0 Sh,0
:=
Shk,−1
:= {u ∈
∩ C(Ω)
Shk,0
, k ≥ 1,
∀K ∈ Th } , k ≥ 0,
(3.2)
: u = 0 on Γ} , k ≥ 1.
k,0 For any k ≥ 1 and any 1 ≤ p ≤ ∞ we obviously have Shk,0 ⊂ W 1,p (Ω) and Sh,0 ⊂ 1,p W0 (Ω).
Let V1 , ..., VNτ be finite element subspaces of C(Ω) associated with the partitions T1 , ..., TNτ introduced above. We then define finite element spaces in space and time by Sτk,−1 (Vh(τ ) ) := span{χJj (t)tµ vj,µ (x) : 0 ≤ µ ≤ k, 1 ≤ j ≤ Nτ , vj,µ ∈ Vj } , k ≥ 0. (3.3) Here, χJj denotes the characteristic function of the interval Jj . The elements of Sτk,−1 (Vh(τ ) ) are discontinuous at the intermediate points t2 , ..., tNτ . But the leftsided limits u(., tj − 0) := lim u(., tj − t) exist for all 2 ≤ j ≤ Nτ + 1; and the 0:=− < u0 , ϕ(., 0) >W 1,π Z T < u(., t), ∂t ϕ(., t) >W 1,π dt − 0
+
Z
T
0
Z
Ω
(4.2)
{a(x, u, ∇u)∇ϕ − b(x, u, ∇u)ϕ}dxdt.
We recall that a prime denotes the dual of a Lebesgue exponent, i.e. In order to ensure that F is well defined, we alsways assume that
1 p
+
1 p′
ρ ≥ π.
= 1 etc. (4.3)
Within the framework of Section 2, we further set ′
′
′
′
Y+ := WTp (0, T ; W 2,π (Ω) ∩ W01,π (Ω), Lπ (Ω)),
X− := Lr (0, T ; Lρ (Ω)),
(4.4)
X+ := Lr (0, T ; W01,σ (Ω)) where ρ ≤ σ < ∞.
In order to better understand the flavour of problem (1.1) and definition (4.2), we recall the notions of weak and very weak solutions of problem (4.1) (cf. [2], Sections 11 and 13). A function u ∈ W r (0, T ; W01,ρ (Ω), W −1,π (Ω)) is called a weak solution of problem (4.1) if u(., 0) = u0 in W −1,π (Ω) and Z
0
T
< ∂t u(., t), ϕ(., t) >W 1,π′ + ′
′
Z
0
T
Z
Ω
{a(x, u, ∇u)∇ϕ − b(x, u, ∇u)ϕ}dxdt = 0
for all ϕ ∈ Lp (0, T ; W 1,π (Ω)). It is a very weak solution if u ∈ X− and Z T − < u0 , ϕ(., 0) >Lρ − < u(., t), ∂t ϕ(., t) >Lρ dt +
Z
0
0
T
{< a(., u, ∇u), ∇ϕ >W 1,π′ − < b(x, u, ∇u), ϕ >W 1,π′ }dt = 0 19
for all ϕ ∈ Y+ . Obviously, every solution of problem (1.1) is a very weak solution of problem (4.1). Conversely, every very weak solution of problem (4.1) is a solution of problem (1.1) if it is contained in X. Using integration by parts with respect to the time variable, one sees that every weak solution of problem (4.1) is a solution of problem (1.1) and that the converse is true if the solution of problem (1.1) is contained in W r (0, T ; W01,ρ (Ω), W −1,π (Ω)). In this sense, a solution of problem (1.1) is weaker than a weak solution and stronger than a very weak solution of problem (4.1). ∗ One easily checks that DF (u)∗ ∈ Isom(Y+ , X− ) if and only if the adjoint linearized problem
−∂t w − ∇ · (A(x, u, ∇u)∇w) + ∂y a(x, u, ∇u) · ∇w
+∇ · (∇z b(x, u, ∇u)w) − ∂y b(x, u, ∇u)w = g in Ω × (0, T )
w = 0 on Γ × (0, T )
w(., T ) = 0 in Ω
∗ admits, for each g ∈ X− , a unique solution w ∈ Y+ such that kwkY+ ¹ kgkX−∗ .
Examples. We consider two particular examples: (1) a heat equation with nonlinear source term a(x, u, ∇u) = ∇u,
b(x, u, ∇u) = f (u),
f ∈ C 1 (IR, IR),
|f ′ (s)| ≤ γ
∀s ∈ IR,
p = r = π = ρ = 2. (2) a nonlinear convection-diffusion equation a(x, u, ∇u) := k(u)∇u,
b(x, u, ∇u) := f − c(x, u) · ∇u,
f ∈ L∞ (Ω), c ∈ C 1 (Ω × IR, IRn ), k ∈ C 2 (IR, IR), k(s) ≥ α > 0, |k (l) (s)| ≤ γ |∂s c(x, s)| ≤ γ
∀s ∈ IR, l ∈ {0, 1, 2},
∀x ∈ Ω, s ∈ IR, n , r ≥ 2p. ρ = π ∈ (n, 4), p > n−1
If in example (1) the constant γ is sufficiently small, we may use an energy estimate −1 and a perturbation argument to obtain explicit bounds on kDF (u)∗ kL(X−∗ ,Y+ ) in terms of the norm of the inverse Laplacian with homogeneous Dirichlet boundary conditions. More precisely, denote by c∆ := k(−∆)−1 kL(L2 (Ω),W 1,2 (Ω)∩W 2,2 (Ω)) 0
20
the norm of the inverse Laplacian with homogeneous Dirichlet boundary conditions. This quantity only depends on the geometry of Ω. Inserting a and b given above into the adjoint linearized problem, we immediadely see that DF (u)∗ = L − N where L is the operator associated with the time-reversed heat equation −∂t v − ∆v = g
in Ω × (0, T )
on Γ × (0, T )
v=0
v(., T ) = 0
(4.5)
in Ω
and where the operator N is given by < N v, ϕ >Y+ =
Z
T
0
Z
f ′ (u)vϕ.
(4.6)
Ω
From equations (4.5) and (4.6) we deduce that kLkL(Y+ ,X−∗ ) ≤ 2,
kN kL(Y+ ,X−∗ ) ≤ γ
and, hence, kDF (u)∗ kL(Y+ ,X−∗ ) ≤ 2 + γ. Multiplying the first equation of (4.5) with −∂t v, integrating over Ω × (0, T ), and using integration by parts with respect to the space variable, we conclude that k∂t vkL22 (Ω×(0,T )) ≤ kgkL22 (Ω×(0,T )) .
(4.7)
Writing the first equation of (4.5) in the form −∆v = g + ∂t v and using the estimate (4.7), we obtain on the other hand that √ kvkL2 (0,T ;W 2,2 (Ω)) ≤c∆ 2{kgkL22 (Ω×(0,T )) + k∂t vkL22 (Ω×(0,T )) } √ ≤c∆ 8kgkL22 (Ω×(0,T )) . Estimates (4.7) and (4.8) yield kL−1 kL(X−∗ ,Y+ ) ≤ 1 + Assume that γ(1 +
√
√
8c∆ ) < 1.
21
8c∆ .
(4.8)
A standard perturbation argument then gives
−1
kDF (u)∗ kL(X−∗ ,Y+ ) =k(L − N )−1 kL(X−∗ ,Y+ ) £ ¤−1 ≤kL−1 kL(X−∗ ,Y+ ) 1 − kL−1 kL(X−∗ ,Y+ ) kN kL(Y+ ,X−∗ ) √ 1 + 8c∆ √ . ≤ 1 − γ(1 + 8c∆ )
Finite element discretization. For the discretization of problem (4.1) we proceed as in Section 3. We choose a family Iτ of shape regular partitions of the interval [0, T ]. With each time tj , 1 ≤ j ≤ Nτ , we associate an admissible and shape regular partition Tj of Ω into n-simplices and a finite element space Vj ⊂ W01,σ (Ω) corresponding to Tj and consisting of affine equivalent finite elements in the sense of [4]. We choose an integer k and a parameter θ ∈ [0, 1] and set
Xh :=Sτk,−1 (Vh(τ ) ), Yh :=Sτθ;k+1,0 (Vh(τ ) ), < Fh (uh ), ϕh >Yh := < F (uh ), ϕh >Y
(4.9) ∀uh ∈ Xh , ϕh ∈ Yh .
For simplicity, we use in (4.9) the parameter h as an acronym for the meshsizes both in space and in time. We recall that the spaces on the right-hand side of (4.9) are defined (θ) in (3.3) and (3.4) with λj replaced by λj and that Xh ⊂ X+ and Yh ⊂ Y . Hence, the discretization (4.9) is conforming. It is also consistent, i.e kRh∗ [F (uh ) − Fh (uh )]kY+∗ = 0. It is a Petrov-Galerkin discretization since the test and trial spaces are different: The trial functions are discontinuous in time, piecewise polynomials of degree k, the test functions are continuous in time, piecewise polynomials of degree k + 1.
Relation to Runge-Kutta schemes. In order to better understand the flavour of problem (1.2) and definition (4.9) we rewrite Fh . Recalling that the functions in Yh are continuous at the intermediate times t2 , ..., tNτ and vanish at the final time T 22
and using integration by parts on each time interval, we get for uh ∈ Xh and ϕh ∈ Yh < F (uh ), ϕh >Y Z u0 (x)ϕh (x, 0)dx =− +
Ω N τ X j=1
{
Z
tj+1
tj
=
+ Z
Z
0
Ω
+
T
Ω
uh (x, tj + 0)ϕh (x, tj + 0)dx −
Z
Ω
uh (x, tj+1 − 0)ϕh (x, tj+1 − 0)dx}
{a(x, uh , ∇uh )∇ϕh − b(x, uh , ∇uh )ϕh }dxdt
[uh (x, 0 + 0) − u0 (x)]ϕh (x, 0)dx
Nτ Z X j=2
+
Ω
Z
∂t uh (x, t)ϕh (x, t)dxdt
Ω
Z
+
Z
Ω
Nτ Z X j=1
[uh (x, tj + 0) − uh (x, tj − 0)]ϕh (x, tj )dx
tj+1
tj
Z
Ω
[∂t uh ϕh + a(x, uh , ∇uh )∇ϕh − b(x, uh , ∇uh )ϕh ]dxdt.
Using the convention that uh (., 0 − 0) := u0 ,
(4.10)
we may write this in the compact form < F (uh ), ϕh >Y Nτ nZ X = [uh (x, tj + 0) − uh (x, tj − 0)]ϕh (x, tj )dx j=1
Ω
+
Z
tj+1
tj
Z
Ω
(4.11)
o [∂t uh ϕh + a(x, uh , ∇uh )∇ϕh − b(x, uh , ∇uh )ϕh ]dxdt .
We first consider the case k = 0 and set ujh := uh (., tj + 0) = uh (., tj+1 − 0) ∀1 ≤ j ≤ Nτ . (θ)
Observing that uh is piecewise constant on the time intervals, inserting ϕh = λj vj , 2 ≤ j ≤ Nτ , vj ∈ Vj , as a test function in (4.11) and using Lemma 3.1, we obtain < Fh (uh ), ϕh >Yh Z = {ujh − uj−1 h }vj dx Ω Z + θτj {a(x, ujh , ∇ujh )∇vj − b(x, ujh , ∇ujh )vj }dx Ω Z j−1 j−1 j−1 + (1 − θ)τj−1 {a(x, uj−1 h , ∇uh )∇vj − b(x, uh , ∇uh )vj }dx. Ω
23
(θ)
Inserting ϕh = λ1 v1 , v1 ∈ V1 , as a test function in (4.11) we similarly get Z < Fh (uh ), ϕh >Yh = {u1h − u0 }v1 dx Ω Z + θτ1 {a(x, u1h , ∇u1h )∇v1 − b(x, u1h , ∇u1h )v1 }dx. Ω
Hence, in the case k = 0, problem (1.2) yields the popular θ-scheme. In particular, the parameters θ = 0, θ = 1, and θ = 21 correspond to the explicit Euler scheme, the implicit Euler scheme, and the trapezoidal rule (Crank-Nicholson scheme). Thus the time discretization is of first order unless θ = 21 ; in this case it is of second order. Moreover, the time discretization is A-stable if θ ≥ 12 .
Next we consider the case k ≥ 1. Denote by pˆl , 0 ≤ l ≤ k−1, a set of orthonormal polynomials of degree l with respect to the weight function 4t(1 − t) on [0, 1]. Let Rt qˆl (t) := 0 pˆl (s)ds, 0 ≤ l ≤ k − 1, and qˆ−1 (t) = 1. For 1 ≤ j ≤ Nτ , set pl,j := pˆl ◦ FJ−1 , qm,j := qˆm ◦ FJ−1 , j j
0 ≤ l ≤ k − 1, −1 ≤ m ≤ k − 1.
Then every uh ∈ Xh and every ϕh ∈ Yh have unique representations of the form uh =
Nτ X j=1
ϕh =
Nτ X j=1
χJj (t){
k−1 X
qµ,j (t)vµ,j (x)},
µ=−1 (θ) λj (t)wj (x)
+
Nτ X
ψj (t)χJj (t){
k−1 X
pµ,j (t)wµ,j (x)}
µ=0
j=1
with vµ,j , wµ,j , wj ∈ Vj , 1 ≤ j ≤ Nτ , Here, χJj denotes the characteristic function of Jj and ψj (t) = τ42 (t − tj )(tj+1 − t). Consider a fixed j ∈ {2, ..., Nτ } and insert j
ϕh = ψj (t)χJj (t)pµ,j (t)wµ,j (x), 0 ≤ µ ≤ k − 1, as a test function in (4.11). We then get < Fh (uh ), ϕh >Yh Z tj+1 Z {∂t uh ϕh + a(x, uh , ∇uh )∇ϕh − b(x, uh , ∇uh )ϕh }dxdt = tj
=
Ω
Z
vµ,j (x)wµ,j (x)dx Z tj+1 Z + { [a(x, uh , ∇uh )∇wµ,j − b(x, uh , ∇uh )wµ,j ]dx}ψj (t)pµ,j (t)dt. Ω
tj
Ω
(θ)
Inserting ϕh = λj (t)wj (x) as a test function in (4.11), we obtain on the other hand < Fh (uh ), ϕh >Yh Z = {v−1,j (x) − uh (x, tj − 0)}wj (x)dx Ω Z tj+1 Z (θ) { [∂t uh wj + a(x, uh , ∇uh )∇wj − b(x, uh , ∇uh )wj ]dx}λj (t)dt. + tj−1
Ω
24
With the obvious modifications these expressions also hold for j = 1. Hence, the coefficients v0,j , ..., vk−1,j are functions of the coefficient v−1,j and the latter one is determined by the values of uh on the previous time interval. Thus problem (1.2) amounts in a (k + 1)-stage implicit Runge-Kutta scheme. A lengthy, but straightforward calculation shows that, for linear problems and k ∈ {1, 2}, θ = 12 , this scheme corresponds to the (k + 1)-st diagonal Pad´e-approximation. In particular, the time discretization then is of order 2k + 2 and A-stable. 4.1 Remark. When writing problem (1.2) in the form (4.11) it strongly resembles the discontinuous Galerkin method (cf. e.g. [7, 8]). In the discontinuous Galerkin method, however, the test and trial spaces are identical and consist both of discontinuous in time, piecewise polynomials of degree k. In particular, the case k = 0 corresponds to the implicit Euler scheme. Due to the discontinuities at the intermediate times t2 , ..., tNτ −1 the discontinuous Galerkin method is non-conforming both with respect to the standard weak formulation of problem (4.1) and to the formulation (4.2). This complicates its analysis within the framework of Secction 2. This difficulty will be overcome in a subsequent paper. A different analysis of the discontinuous Galerkin method is given in, e.g., [7, 8]. Definition of Rh , F˜h , and Y˜h . In order to put the discretization in the framework of Section 2, we assume that Yh contains the space θτ defined in (3.6). This is equivalent to assuming that the space discretization at least consists of linear el1,0 ements, i.e. Vj ⊃ Sj,0 , 1 ≤ j ≤ Nτ . As restriction operator Rh we use the operator Iτ defined in (3.10). For the construction of F˜h and Y˜h we define integers µ, ν and approximations ah of a and bh of b as follows: if a(x, vh , ∇vh ) a(x, uh , ∇uh ), µ,−1 ∈ Sτµ,−1 (Sh(τ ) ) ∀vh ∈ Xh X ah (x, uh , ∇uh ) := π1,Q a(x, uh , ∇uh ), µ := 1, otherwise Q∈Pτ
b(x, uh , ∇uh ), if b(x, vh , ∇vh ) ν,−1 ∈ Sτν,−1 (Sh(τ ) ) ∀vh ∈ Xh X bh (x, uh , ∇uh ) := π0,Q b(x, uh , ∇uh ), ν := 0, otherwise. Q∈Pτ
(4.12) Here, uh ∈ Xh is arbitrary and π0,Q and π1,Q denote the L (Q)-projections onto the spaces of polynomials of degree at most 0 and 1 (in the variables x and t), respectively. Now, F˜h is defined in the same way as F with a and b replaced by ah and bh , respectively, and 2
Y˜h := span{ψj ψK v, ψj ψE PE σ, ψj ψK Pj w : 1 ≤ j ≤ Nτ , K ∈ Tj , E ∈ Ej , ˜ m|K×{t } }. ˜ m|E×J , w ∈ IP ˜ m|K×J , σ ∈ IP v ∈ IP j
j
25
j
(4.13)
˜ m denotes the space of polynomials of degree at Here, m := max{µ − 1, ν} and IP most m in the variables x and t. The estimators. Given Q = K × Jj ∈ Pτ , we recall the abbreviation (3.17) and set εQ,π :=(h2K + τj )k∇ · [a(., uh , ∇uh ) − ah (., uh , ∇uh )] 1
π + hK
−1
+ [b(., uh , ∇uh ) − bh (., uh , ∇uh )]kLpπ (Q)
(h2K + τj )knE · [a(., uh , ∇uh ) − ah (., uh , ∇uh )]E kLpπ (∂QL ) ,
ηQ,π :=(h2K + τj )k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (Q) 1
−1
1
−1
π + hK
+ τjp
(4.14)
(h2K + τj )knE · [ah (., uh , ∇uh )]E kLpπ (∂QL )
(h2K + τj )kuh (., tj + 0) − uh (., tj − 0)kLπ (K) .
The quantity εQ,π obviously measures the quality of the approximation of a and b by ah and bh respectively, and can be estimated explicitly. Below, we will show that it yields upper bounds on the second terms on the right-hand sides of estimates (2.3) and (2.4). Note, that in our second example
εQ,π ≤ h2K kf − π0,Q f kLpπ (Q) + h3K k∇uh kLpπ (Q) , 1,0 if k = 0 and Vj = Sj,0 , 1 ≤ j ≤ Nτ .
Estimation of k(IdY+ −Rh )∗ F˜h (uh )kY+∗ and k(IdY+ −Rh )∗ [F (uh )−F˜h (uh )]kY+∗ . Next, we will derive upper bounds for the first and second terms on the right-hand side of inequality (2.3). Recalling equation (4.11) and using, for the space variables, integration by parts elementwise, we obtain for all ϕ ∈ Y < F (uh ), ϕ >Y Nτ n X Z X { [uh (x, tj + 0) − uh (x, tj − 0)]ϕ(x, tj )dx = j=1 K∈Tj
K
+ +
X Z
E∈Ej
Z
tj+1
tj tj+1
tj
Z
Z
E
K
[∂t uh − ∇ · a(x, uh , ∇uh ) − b(x, uh , ∇uh )]ϕdxdt}
nE · [a(x, uh , ∇uh )]E ϕdsdt 26
o
(4.15)
and < F˜h (uh ), ϕ >Y Nτ n X Z X { [uh (x, tj + 0) − uh (x, tj − 0)]ϕ(x, tj )dx = j=1 K∈Tj
K
+ X Z
+
E∈Ej
Z
tj+1
tj
tj+1
tj
Z
Z
K
E
[∂t uh − ∇ · ah (x, uh , ∇uh ) − bh (x, uh , ∇uh )]ϕdxdt}
o nE · [ah (x, uh , ∇uh )]E ϕdsdt .
(4.16) Let ϕ ∈ Y+ = × (0, T )) be arbitrary. From Lemma 3.4 we obtain for Q = K × Ij ∈ Pτ , 1 ≤ j ≤ Nτ , K ∈ Tj , and E ⊂ ∂K\Γ the estimates Z [uh (x, tj + 0) − uh (x, tj − 0)][ϕ(x, tj ) − Iτ ϕ(x, tj )]dx ′ Wπp′ (Ω
K
≤ kuh (., tj + 0) − uh (., tj − 0)kLπ (K) kϕ − Iτ ϕkLp′ (∂Q
B)
π′
1
¹ τjp Z tj+1 Z tj
−1
(h2K + τj )kuh (., tj + 0) − uh (., tj − 0)kLπ (K) kϕkW p′ (K×U (t
j ))
π′
K
,
[∂t uh − ∇ · ah (x, uh , ∇uh ) − bh (x, uh , ∇uh )][ϕ − Iτ ϕ]dxdt
≤ k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (Q) kϕ − Iτ ϕkLp′ (Q) π′
¹ (h2K + τj )k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (Q) kϕkW p′ (U (Q)) , π′ Z tj+1 Z nE · [ah (x, uh , ∇uh )]E [ϕ − Iτ ϕ]dsdt tj
E
≤ knE · [ah (., uh , ∇uh )]E kLpπ (∂QL ) kϕ − Iτ ϕkLp′ (∂Q π′
1
¹ hEπ
−1
L)
(h2K + τj )knE · [ah (., uh , ∇uh )]E kLpπ (∂QL ) kϕkW p′ (U (∂Q
L ))
π
.
Inserting these estimates in (4.16), using H¨older’s inequality for finite sums, and recalling definition (4.14), we conclude that < F˜h (uh ), ϕ − Iτ ϕ >
Nτ Nτ X X X X ′ π p/π 1/p ¹{ { ηQ,π } } { { kϕkπ p′ j=1 K∈Tj
j=1 K∈Tj
′
Wπ′ (U (K×Jj ))
′
′
}p /π }1/p .
Assume that p′ ≤ π ′ or, equivalently, p ≥ π. Then, Jensen’s inequality implies that Nτ X X ′ { { kϕkπ p′ j=1 K∈Tj
Wπ′ (U (K×Jj ))
′
′
′
}p /π }1/p ¹ kϕkW p′ (Ω×(0,T )) . π′
27
Using the abbreviation Nτ X X π η := { { ηQ,π }p/π }1/p
(4.17)
j=1 K∈Tj
we have thus shown that k(IdY+ − Rh )∗ F˜h (uh )kY+∗ ¹ η.
(4.18)
Replacing ah and bh by ah − a and bh − b, respectively, we conclude with the same arguments that k(IdY+ − Rh )∗ [F˜h (uh ) − F (uh )]kY+∗ ¹ ε (4.19) where
Nτ X X ε := { { επQ,π }p/π }1/p .
(4.20)
j=1 K∈Tj
Estimation of kF˜h (uh )kY˜ ∗ and kF˜h (uh ) − F (uh )kY˜ ∗ . Now, we will bound the h h terms in inequality (2.4). Given a subset Q of Ω × (0, T ), we set for abbreviation Y˜h|Q := {ϕ ∈ Y˜h : suppϕ ⊂ Q}. In order to bound the second term on the right-hand side of estimate (2.4), we conclude from the shape regularity of the partitions and a standard scaling argument that the estimate 1
π′ kϕkLp′ (K ′ ×J ′ ) + hK kϕkLp′ (∂K ′ ×J ′ ) ¹ h2K kϕkLp′ (J ′ ,W 2,π′ (K ′ )) π′
π′
holds for all ϕ ∈ Y˜h|U (Q) , Q = K × Jj ∈ Pτ , J ′ ⊂ U (Jj ), K ′ ⊂ U (K). Combining this with equations (4.14) – (4.16), and using H¨older’s inequality, we obtain the following estimate for all Q ∈ Pτ X kF (uh ) − F˜h (uh )kY˜ ∗ ¹ εQ′ ,π , (4.21) h|U (Q)
Q′ ⊂U (Q)
where Y˜h equipped with the norm of Y+ . In order to derive lower bounds for the left-hand side of estimate (2.4), consider ˜ m and equation (4.14) an arbitrary Q = K × Jj ∈ Pτ . From Lemma 3.5 with VQˆ = IP we then obtain k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (Q) Z −1 ¹ sup kvk p′ [∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )]ψj ψK v ˜m v∈IP
= sup ˜m v∈IP
Lπ′ (Q)
kvk−1p′ Lπ′ (Q)
≤kF˜h (uh )kY˜ ∗
h|Q
Q
< F˜h (uh ), ψj ψK v >Y+
sup kvk−1 kψj ψK vkY+ Lp (Q)
˜m v∈IP
π′
−1 ˜ ¹{h−2 K + τj }kFh (uh )kY˜ ∗ . h|Q
28
(4.22)
˜ ˆ ˆ yield for all Estimate (4.22), equation (4.16) and Lemma 3.5 with V∂ Qˆ L = IP m|E×J E ⊂ ∂K\Γ knE · [ah (., uh , ∇uh )]E kLpπ (∂QL ) Z −1 ¹ sup kσk p′ nE · [ah (., uh , ∇uh )]E ψj ψE PE σ Lπ′ (∂QL )
˜ m|∂Q σ∈IP L
≤
Lπ′ (∂QL )
˜ m|∂Q σ∈IP L
− ≤
kσk−1p′
sup
X Z
˜ m|∂Q σ∈IP L
K′
kσk−1p′
sup +
Jj
K ′ ⊂ωE
Z
K ′ ⊂ωE
{< F˜h (uh ), ψj ψE PE σ >Y+
[∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )]ψj ψE PE σdxdt}
Lπ′ (∂QL )
X
∂QL
{kF˜h (uh )kY˜h|ω
E ×Jj
kψj ψE PE σkY+
k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (K ′ ×Jj )
· kψj ψE PE σkLp′ (K ′ ×J ) } j
π′
1 π′
−1 ˜ ¹hE {h−2 E + τj }kFh (uh )kY˜ ∗
h|ωE ×Jj
.
From estimate (4.22), equation (4.16) and Lemma 3.5 with V∂ Qˆ B further conclude that
(4.23) ˜ = IPm|K×{0} we ˆ
kuh (., tj + 0) − uh (., tj − 0)kLπ (K) Z −1 kwkLπ′ (K) ¹ sup [uh (., tj + 0) − uh (., tj − 0)]ψj ψK Pj wdx ˜ m|K×{t } w∈IP j
=
sup ˜ m|K×{t } w∈IP j
− ≤
{< F˜h (uh ), ψj ψK Pj w >Y+ kwk−1 Lπ′ (K)
X Z Z
J⊂U (tj )
sup ˜ m|K×{t } w∈IP j
+
K
X
J⊂U (tj )
J
K
[∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )]ψj ψK Pj wdxdt}
kwk−1 {kF˜h (uh )kY˜ ∗ Lπ′ (K)
h|K×U (tj )
kψj ψK Pj wkY+
k∂t uh − ∇ · ah (., uh , ∇uh ) − bh (., uh , ∇uh )kLpπ (K×J)
· kψj ψK Pj wkLp′ (K×J) } π′
1 p′
−1 ˜ ¹τj {h−2 K + τj }kFh (uh )kY˜ ∗
h|K×U (tj )
. (4.24)
Estimates (4.22) – (4.24) and definition (4.14) yield −1 2 ˜ ηQ,π ¹ {1 + τj h−2 K + τj hK }kFh (uh )kY˜ ∗
h|U (Q)
29
(4.25)
where Y˜h is endowed with the norm of Y+ . A posteriori error estimates. Combining estimates (4.18), (4.19), (4.21), and (4.25) with Propositions 2.1, 2.2 and recalling that kRh∗ [F (uh ) − Fh (uh )]kY+∗ = 0, we obtain the following result. 4.2 Proposition. Let u be a regular solution of problem (4.1) in the sense of Proposition 2.1 and definition (4.2) and let uh be an approximation of u in the sense of Proposition 2.2 and definition (4.5). Suppose that p ≥ π. Then the following a posteriori error estimates hold: ku − uh kLrρ (Ω×(0,T )) ¹ {η + ε + kRh∗ Fh (uh )kY+∗ } and
ηQ,π ¹{1 + τj−1 h2K + τj h−2 K }{ku − uh kLrρ (U (Q)) + ∀Q = K × Jj ∈ Pτ .
X
Q′ ⊂U (Q)
εQ′ ,π }
The quantities εQ,π , ηQ,π , η, and ε are given by equations (4.14), (4.17), and (4.20). 4.3 Remark. The local lower bounds for ku − uh kLrρ (Ω×(0,T )) can be combined in the standard way to the global lower bound η ¹ max max {1 + τj−1 h2K + τj h−2 K } · {ku − uh kLrρ (Ω×(0,T )) + ε}. 1≤j≤Nτ K∈Tj
The factor 1 + τj−1 h2K + τj h−2 K in this estimate and the second one of Proposition 4.2 reflects the fact that the differential operator is of 2nd order with respect to the space variables but only of 1st order with respect to the time variable. Remark. If p < π one may still obtain upper bounds on the error. Since, in case, Jensen’s inequality cannot be used in estimating < F˜h (uh ), ϕ − Iτ ϕ >Y+ < F˜h (uh ) − F (uh ), ϕ − Iτ ϕ >Y+ , one must now proceed as follows: Bound the space-integrals by using H¨ older’s inequality and Lemma 3.2. On each time-level sum up all contributions to that level and apply H¨older’s inequality for finite sums. 3. Bound the remaining time-integrals by using H¨older’s inequality and Lemma 3.3. 4. Sum up over all time-levels and use H¨older’s inequality for finite sums.
4.4 this and 1. 2.
5. Non-stationary, incompressible Navier-Stokes equations Variational setting. In this section we consider the two- and three dimensional 30
non-stationary, incompressible Navier-Stokes equations ∂t u − ν∆u + (u · ∇)u + ∇p = f
in Ω × (0, T )
in Ω × (0, T )
∇·u=0 u=0
u(., 0) = u0
on Γ × (0, T )
(5.1)
in Ω.
Here, ν > 0 is the constant viscosity of the fluid; u and p denote its velocity and pressure, respectively. In order to cast problem (5.1) into the abstract framework of Section 2, we introduce some auxiliary spaces Z 2 2 p = 0}, L0 (Ω) := {p ∈ L (Ω) : Ω ½ L∞ (0, T ; L2 (Ω)2 ) ∩ L2 (0, T ; W01,2 (Ω)2 ) , if n = 2, X1 := 1,2 ∞ 2 3 2 3 8 4 3 L (0, T ; L (Ω) ) ∩ L (0, T ; W0 (Ω) ) ∩ L (0, T ; L (Ω) ) , if n = 3, X2 := L2 (0, T ; L20 (Ω)),
Y1 := WT2 (0, T ; W01,2 (Ω)n , W −1,2 (Ω)n ), Y2 := L2 (0, T ; L20 (Ω)), H1 := L2 (0, T ; L2 (Ω)n ), H2 := L2 (0, T ; (W 1,2 (Ω) ∩ L20 (Ω))∗ ),
Z1 := WT2 (0, T ; W 2,2 (Ω)n ∩ W01,2 (Ω)n , L2 (Ω)n ) Z2 := L2 (0, T ; W 1,2 (Ω) ∩ L20 (Ω)),
and set X :=X1 × X2 , Y :=Y1 × Y2 ,
X− :=H1 × H2 ,
X+ :=L2 (0, T ; W01,3 (Ω)n ) × X2 ,
Y+ :=Z1 × Z2 Z Z u0 (x)v(x, 0)dx − < F ([u, p]), [v, q] >Y := − +
Ω Z T 0
+
Z
0
T
Z
0
T
Z
(5.2) u∂t vdxdt
Ω
Ω
{ν∇u · ∇v + [(u · ∇)u]v − p∇ · v − f v}dxdt
Ω
∇ · uqdxdt.
Z
One easily checks, that every solution [u, p] of problem (1.1) is a weak solution of problem (5.1) in the sense of [9; Problem III.3.2] if ∂t u ∈ L1 (0, T ; W −1,2 (Ω)n ). 31
Conversely, every weak solution [u, p] of problem (5.1) in the sense of [9] is a solution of problem (1.1) if u ∈ X. From this and Theorems III.3.1, III.3.2, and III.3.4 in [9] it follows that problem (1.1) with the definition (5.2) has exactly one solution, if n = 2, and at most one solution, if n = 3. Since Z1 ⊂ L∞ (0, T ; W01,2 (Ω)n ) ⊂ L∞ (0, T ; L6 (Ω)n ) with continuous imbeddings (cf. [6] Vol. 5, Chap. XVIII, §1, Theorem 1), one easily checks that condition (2.1) is fulfilled if the solution of problem (1.1) is contained in X+ . Moreover, DF ([u, p])∗ ∈ ∗ Isom(Y+ , X− ) if and only if the adjoint linearized Navier-Stokes equations −∂t v − ν∆v − (u · ∇)v + v∇u − ∇q = w
in Ω × (0, T )
v=0
on Γ × (0, T )
−∇ · v = r
u(., T ) = 0
in Ω × (0, T ) in Ω
∗ admit for each [w, r] ∈ X− a unique weak solution [v, q] ∈ Y+ satisfying
k[v, q]kY+ ¹ k[w, r]kX−∗ . Standard energy estimates show that this is true if u ∈ L2 (0, T ; W 2,2 (Ω)n ∩W01,2 (Ω)n ) and kukL2 (0,T ;W 2,2 (Ω)n ) is sufficiently small, e.g., if ν −2 kf kL22 (Ω×(0,T )) < 1.
Finite element discretization. For the discretization of problem (5.1) we proceed as in Sections 3 and 4. We choose a family Iτ of shape regular partitions of the interval [0, T ]. With each time tj , 1 ≤ j ≤ Nτ , we associate an admissible and shape regular partition Tj of Ω into n-simplices and two finite element spaces Mj ⊂ W01,2 (Ω)n and Qj ⊂ L20 (Ω) corresponding to Tj and consisting of affine equivalent elements in the sense of [4]. In what follows we always assume that there are two 1,0 n µ,0 n integers µ ≥ 1 and ν ≥ 0 such that, for all 1 ≤ j ≤ Nτ , [Sj,0 ] ⊂ Mj ⊂ [Sj,0 ] and 1,0 ν,0 0,−1 ν,−1 either Sj ⊂ Qj ⊂ Sj or Sj ⊂ Qj ⊂ Sj . We then set for k ≥ 0 and θ ∈ [0, 1] Xh := Sτk,−1 [Mh(τ ) ] × Sτk,−1 [Qh(τ ) ],
Yh := Sτθ;k+1,0 [Mh(τ ) ] × Sτk,−1 [Qh(τ ) ],
< Fh ([uh , ph ]), [v h , qh ] >Yh :=< F ([uh , ph ]), [v h , qh ] >Y
(5.3)
∀[uh , ph ] ∈ Xh , [v h , qh ] ∈ Yh . Note, that the above discretization is consistent, i.e. kRh∗ [F ([uh , ph ])−Fh ([uh , ph ])]kY+∗ = 0 and that for the pressure the test- and trial-space are identical. As in Section 4 one easily checks that this discretization corresponds to a (k + 1)-stage Runge-Kutta scheme. For k = 0, in particular, it yields the θ-scheme. Definition of Rh , F˜h , and Y˜h . In order to cast the discretization (5.3) in the framework of Section 2, we now define F˜h , Y˜h and Rh . Denote by π0,Q f the L2 projection of f onto the space of piecewise constant functions corresponding to the 32
partition Pτ of Ω × (0, T ). Then, F˜h is defined in the same way as F with f replaced by π0,Q f . Next, we set Y˜h :=
span{ψj ψK v, ψj ψE PE σ, ψj ψK Pj w : ˜ m |K×J , σ ∈ IP ˜ m |E×J , w ∈ IP ˜ m |K×{t } , v ∈ IP 1 j 2 j 3 j
1 ≤ j ≤ Nτ , K ∈ Tj , E ∈ Ej } ˜ m |K×J , 1 ≤ j ≤ Nτ , K ∈ Tj } ×span{χK×Jj q : q ∈ IP 4 j where
m1 := max{2k + 2µ + 1, k + ν − 1} , m2 := max{k + µ − 1, k + ν}, , m4 := k + µ − 1,
m3 := k + µ
˜m and χK×Jj denotes the characteristic function of the set K × Jj . Recall that IP denotes space of polynomials (in x and t) of degree at most m. In order to define Rh , we set, for 1 ≤ j ≤ Nτ , ( 0,−1 ∩ L20 (Ω) if Qj ⊂ 6 C(Ω), ˆ j := Sj Q 1,0 2 Sj ∩ L0 (Ω) if Qj ⊂ C(Ω), ˆ j either the L2 -projection of L2 (Ω) onto Q ˆ j , if Qj 6⊂ C(Ω), or the and denote by Π 0 Cl´ement interpolation operator Ij , if Qj ⊂ C(Ω). Let π ˆj , 1 ≤ j ≤ Nτ , be the L2 projection of L2 (Jj ; L20 (Ω)) onto span{tµ qµ : 0 ≤ µ ≤ k, qµ ∈ L20 (Ω)}. We then define ˆ h(τ ) ] by the operator Iˆτ : L2 (0, T ; L20 (Ω)) → Sτk,−1 [Q Iˆτ q :=
Nτ X
ˆ j q, χJj π ˆj Π
j=1
and set Rh := (Iτ , Iˆτ ), where Iτ is defined in (3.10). The proofs of Lemmas 3.2 and 3.3 can easily be modified ˆ j and π such that one obtains similar approximation results for the operators Π ˆj as for Ij and πj , respectively. The estimators. For Q = K × Jj ∈ Pτ we define εQ :=(h2K + τj )kf − π0,Q f kL22 (Q) , − 21
ηQ :=τj
(h2K + τj )kuh (., tj + 0) − uh (., tj − 0)kL22 (∂QB )
+ (h2K + τj )k∂t uh − ν∆uh + (uh · ∇)uh − ∇ph − π0,Q f kL22 (Q) −1
+ hK 2 (h2K + τj )k[νnE · ∇uh − ph nE ]E kL22 (∂QL ) + hK k∇ · uh kL22 (Q)
33
(5.4)
and ε := {
X
Q∈Pτ
ε2Q }1/2
,
η := {
X
Q∈Pτ
2 1/2 ηQ } .
(5.5)
As in Section 4, we will show that εQ yields upper bounds on the second terms on the right-hand sides of estimates (2.3) and (2.4). Estimation of k(IdY − Rh )∗ F˜h ([u , ph ])kY ∗ and k(IdY − Rh )∗ [F ([u , ph ]) − +
h
+
+
h
F˜h ([uh , ph ])]kY+∗ . Using integration by parts elementwise and the convention (4.10), we obtain, for all [v, q] ∈ Y , the following representation of the residual < F˜h ([uh , ph ]), [v, q] >Y Nτ n X Z X = { {uh (x, tj + 0) − uh (x, tj − 0)}v(x, tj )dx K
j=1 K∈Tj
Z
+ +
X Z
E∈Ej
+
tj
tj+1 Z
tj
X Z
K∈Tj
Z
tj+1
E
tj+1
tj
K
Z
{∂t uh − ν∆uh + (uh · ∇)uh + ∇ph − π0,Q f }vdxdt}
(5.6)
[νnE · ∇uh − ph nE ]E · vdsdt
K
o ∇ · uh qdxdt .
Let [v, q] ⊂ Y+ be arbitrary. Since, on each sub-interval Jj , 1 ≤ j ≤ Nτ , ∇ · uh ∈ span{tµ qµ : 0 ≤ µ ≤ k, qµ ∈ L20 (Ω)}, we conclude from the definition of Iˆτ that for all Q = K × Jj , 1 ≤ j ≤ Nτ , K ∈ Tj , Z ∇ · uh (q − Iˆτ q)dxdt Q Z ˆ j q))dxdt ˆj q + π ˆj (q − Π ∇ · uh (q − π = Q Z ˆ j q)dxdt ∇ · uh π = ˆj (q − Π Q
ˆ j q)kL2 (Q) ≤k∇ · uh kL22 (Q) kˆ πj (q − Π 2 ˆ j qkL2 (Q) ≤k∇ · uh kL22 (Q) kq − Π 2
¹hK k∇ · uh kL22 (Q) kqkL2 (Jj ;W 1,2 (K)) .
The terms in < F˜h ([uh , ph ]), [v, q] − Rh [v, q] >Y involving v can be estimated as in Section 4. This proves that k(IdY+ − Rh )∗ F˜h ([uh , ph ])kY+∗ ¹ η. 34
(5.7)
With the same notations and arguments, we obtain < F ([uh , ph ]) − F˜h ([uh , ph ]), [v, q] − Rh [v, q] >Y Nτ X X Z tj+1 Z = { [f − π0,Q f ][v − Iτ v]dxdt} j=1 K∈Tj
tj
K
Nτ X X { ¹ kf − π0,Q f kL22 (K×Jj ) (h2K + τj )kvkW22 (U (K×Jj )) } j=1 K∈Tj
¹
X
Q∈Pτ
εQ kvkW22 (U (Q))
¹εkvkW22 (Ω×(0,T )) . This proves that k(IdY+ − Rh )∗ [F ([uh , ph ]) − F˜h ([uh , ph ])]kY+∗ ¹ ε.
(5.8)
Estimation of kF˜h ([uh , ph ])kY˜ ∗ and kF ([uh , ph ])− F˜h ([uh , ph ])kY˜ ∗ . The second h h term on the right-hand side of estimate (2.4) can be bounded in the same way as in Section 4. In analogy to estimate (4.21), we now obtain for all Q ∈ Pτ kF ([uh , ph ]) − F˜h ([uh , ph ])kY˜ ∗
h|U (Q)
¹
X
εQ′ .
(5.9)
Q′ ⊂U (Q)
Thanks to the definition of Yh and Y˜h , we may insert [0, χQ ∇ · uh ], Q = K × Jj , 1 ≤ j ≤ Nτ , K ∈ Tj , as a test function in equation (5.6). From a standard scaling argument we then obtain k∇ · uh kL22 (Q)
=k∇ · uh k−1 < F˜h ([uh , ph ]), [0, χQ ∇ · uh ] >Y L2 (Q) 2
k∇ · uh kL2 (Jj ;W 1,2 (K)) ≤kF˜h ([uh , ph ])kY˜ ∗ k∇ · uh k−1 L2 (Q) h|Q
(5.10)
2
˜ ¹h−1 ∗ . K kFh ([uh , ph ])kY˜h|Q The remaining terms in ηQ can be bounded in the same way as in Section 4. We therefore obtain ˜ ηQ ¹ {1 + τj−1 h2K + τj h−2 K }kFh ([uh , ph ])kY˜ ∗
h|U (Q)
.
(5.11)
A posteriori error estimates. Combining estimates (5.7) – (5.9), and (5.11) with Propositions 2.1, 2.2 and recalling that the discretization (5.3) is consistent, we obtain the following result. 35
5.1 Proposition. Let [u, p] be a regular solution of problem (5.1) in the sense of Proposition 2.1 and definition (5.2) and let [uh , ph ] be an approximation of [u, p] in the sense of Proposition 2.2 and definition (5.3). Then the following a posteriori error estimates hold: ku − uh kL2 (0,T ;L2 (Ω)n ) + kp − ph kL2 (0,T ;W −1,2 (Ω))
¹{η + ε + kRh∗ Fh ([uh , ph ])kY+∗ },
ηQ ¹{1 + τj−1 h2K + τj h−2 K }{ku − uh kL2 (U (Jj );L2 (U (K))n ) X + kp − ph kL2 (U (Jj );W −1,2 (U (K))) + εQ′ } Q′ ⊂U (Q)
∀Q = K × Jj ∈ Pτ ,
and
η ¹ max max {1 + τj−1 h2K + τj h−2 K }{ku − uh kL2 (0,T ;L2 (Ω)n ) 1≤j≤Nτ K∈Tj
+ kp − ph kL2 (0,T ;W −1,2 (Ω)) + ε}. The quantities εQ , ηQ , ε, and η are given by equations (5.4) and (5.5). 5.2 Remark. Proposition 5.1 can be extended to the case of the slip boundary condition u · n = T (νu, p) − [n · T (νu, p) · n]n = 0, where
1 T (u, p) := ( (∂i uj + ∂j ui ) − pδij )1≤i,j≤n 2 denotes the stress tensor. One then has to replace ν∇uh − ph I in equation (5.4) by T (νuh , ph ), and Γ by the part of the boundary on which the no-slip condition u = 0 is imposed. Here, I := (δij )1≤i,j≤n denotes the unit tensor. Other boundary conditions may be treated similarly. 5.3 Remark. The results of Propositions 4.2 and 5.1 can be combined to handle non-Newtonian fluids for which the viscosity is a given function of the velocity and its gradient (cf. [3, 10] for the stationary case).
References [1] R. A. Adams: Sobolev spaces. Academic Press, New York, 1975 [2] H. Amann: Nonhomogeneous linear and quasilinear elliptic and parabolic boundary value problems. Preprint, Universit¨ at Z¨ urich, 1994 36
[3] J. Baranger and H. El Amri: Estimateurs a posteriori d’erreur pour le calcul adaptif d’´ecoulements quasi-Newtoniens. RAIRO M2 AN 25, 31 – 48 (1991) [4] P. G. Ciarlet: The finite element method for elliptic problems. North Holland, Amsterdam, 1978 [5] P. Cl´ement: Approximation by finite element functions using local regularization. RAIRO Anal. Num´er. 9 R-2, 77 – 84 (1975) [6] R. Dautray and J.-L. Lions: Mathematical analysis and numerical methods for science and technology. Springer, Berlin - Heidelberg - New York, 1992 [7] K. Eriksson and C. Johnson: Adaptive finite element methods for parabolic problems IV. Nonlinear problems. Chalmers University of G¨ oteborg, Preprint 1992:44 (1992) : Adaptive finite element methods for parabolic problems V. Long-time [8] integration. Chalmers University of G¨ oteborg, Preprint 1993:04 (1993) [9] R. Temam: Navier-Stokes Equations. 3rd edition, North Holland, Amsterdam, 1984 [10] R. Verf¨ urth: A posteriori error estimates for nonlinear problems. Finite element discretizations of elliptic equations. Math. Comput. 62 (206), 445 – 475 (1994)
37
