Numer. Math. (2002) 92: 357–381 Digital Object Identifier (DOI) 10.1007/s002110100330
Numerische Mathematik
A numerical approach to degenerate parabolic equations Iuliu Sorin Pop1 , Wen-An Yong2 1 2
Eindhoven University of Technology, Faculty of Mathematics and Computer Science, P.O.Box 513, 5600 MB Eindhoven, The Netherlands; e-mail:
[email protected] Institut f¨ur Angewandte Mathematik, Universit¨at Heidelberg, Im Neuenheimer Feld 294, 69120 Heidelberg, Germany, e-mail:
[email protected]
Received January 20, 1999 / Revised version received February 28, 2000 / c Springer-Verlag 2001 Published online July 25, 2001 –
Summary. In this work we propose a numerical approach to solve some kind of degenerate parabolic equations. The underlying idea is based on the maximum principle. More precisely, we locally perturb the (initial and boundary) data instead of the nonlinear diffusion coefficients, so that the resulting problem is not degenerate. The efficiency of this method is shown analytically as well as numerically. The numerical experiments show that this new approach is comparable with the existing ones. Mathematics Subject Classification (1991): 65M20, 35K65
1. Introduction In this work we propose a numerical approach to a class of degenerate parabolic equations. This kind of equations typically describe the gas flow through porous media. The numerical difficulties in treating such problems arise from the nonlinear degeneracy of the equations. Let Ω be a bounded domain in Rd (d ≥ 1) with a Lipschitz continuous boundary and QT ≡ (0, T ) × Ω, with 0 < T < ∞, fixed. We deal here with the following nonlinear degenerate parabolic problem (1.1)
∂t u − ∇ · (∇β(u) + F (u)) = r(u), u(0, x) = u0 (x) ≥ 0, u = 0,
Correspondence to: I. Sorin Pop
in QT ≡ (0, T ) × Ω, in Ω, on ∂Ω,
358
I.S. Pop, W.-A. Yong
β : R → R being a strictly increasing smooth function. By degeneracy we mean a vanishing diffusion, namely β (u) = 0 for some u. Throughout this paper we assume, without loss of generality, that β (0) = 0. Such problems can be investigated through a parabolic regularization, which can be constructed in different ways. A widely used one is to modify β or β so that the resulting β does not vanish. This suggests various numerical schemes to approximate the above degenerate problems (see [9], [10], [15], [21] and references therein). On the other hand, it is well known that, if positive solutions are sought, the difficulties due to the degeneracy can be overcome by perturbing the data such that the corresponding solutions do admit the value 0—the degenerate point (see [17] for the classical porous medium equation, or, more recently [11] for more general scalar equations). To the authors’ knowledge, such a technique has not been employed for numerical purposes yet. The aim of this paper is to seek the possibility to use this technique numerically. More precisely, our approach is to solve the equations with the locally perturbed (initial and boundary) data, instead of the prescribed ones. The resulting problem has a solution taking values away from the degenerate point. This is guaranteed by the maximum principle obeyed by such equations. Thus, the initially degenerate problem is reduced to a regular parabolic one, which can be solved numerically by different methods. Based on the idea presented in [24], this paper extends the class of problems taken into consideration there. Let n be an integer and τ = T /n. The numerical scheme reads (1.2)
β −1 (θk ) − β −1 (θk−1 ) = τ ∇ · (∇θk + F (β −1 (θk )) +τ r(β −1 (θk−1 )), θk |∂Ω = β(ε),
for k = 1, n, where ε is a positive small perturbation parameter. The initial value is θ0 = β(u0 ) + β(ε)v0 , where v0 ∈ H 1 (Ω) is a perturbation function. Here θk is an approximation to β(u(tk )) (tk = kτ ), so β −1 (θk ) approximates u(tk ). The use of β −1 is justified by the assumptions below. The choice of v0 depends on the given initial data u0 (x), so that (1.3)
v0 |{u0 ≥ε} ≡ 0, v0 ≤ 1 and β(u0 ) + β(ε)v0 ≥ β(ε).
Here the inequalities hold almost everywhere. A typical example is β(u0 ) v0 := 1 − (1.4) , β(ε) + where [x]+ = x if x ≥ 0, otherwise [x]+ = 0. It is worth here to mention some aspects concerning the above scheme. First, this is based essentially on the Euler implicit method. A higher order method would be effective only if the solution would be more regular, which
A numerical approach to degenerate parabolic equations
359
- in general - is not the case. Next, since β(u) is more regular than u - as resulting from the weak formulation of the problem, the numerical scheme approximates β(u) instead of u itself. This idea has been suggested in other papers concerned with regularization procedures, like [15] or [9]. However, while in these cases the regularization step involves also a cut-off procedure for the derivative of β or its inverse in order to avoid the degeneracy, the perturbation of the data makes this step unnecessary in our approach. Finally, the scheme in (1.2) is nonlinear. For simplicity, the reaction term is discretized explicitly. Since the maximum principle plays a crucial role in our approach, a similar treatment of the transport is not possible. However, a semi-implicit linearization is possible, as given in [19]. About the problem in (1.1), we make the following assumptions. (A1) β is continuous and piece-wise differentiable, β(0) = 0, β (u) ≥ 0 and β may vanish at a single point at most (which we assume to be also 0). (A2) u0 ≥ 0, u0 ∈ L∞ (Ω) and β(u0 ) ∈ H01 (Ω). (A3) r : R → R and F : R → Rd are continuous in u and satisfy the condition |r(u) − r(v)|2 + |F (u) − F (v)|2 ≤ C(u − v)(β(u) − β(v)) for any u, v ∈ R, where C > 0 does not depend on x, t, u and v. Moreover, it is assumed here that r is positive for all positive arguments and the graphs of both functions contain the origin, hence r(0) = 0 and F (0) = ¯ 0. The solutions we are interested in here are positive. In (A2) the initial data are positive and therefore the same property holds for the solution of the problem in (1.1). If necessary, the constants appearing in the assumptions (A1) - (A3) can be given explicitly. The Lipschitz constant of β is denoted by Lβ , while the growth of F and r is controlled by CF , respectively Cr . Remark 1.1. In (A3) we have weakened the usual assumption on F and r, namely the Lipschitz continuity with respect to β(u) |r(u) − r(v)| + |F (u) − F (v)| ≤ C|β(u) − β(v)|. The weaker hypotheses in (A3) are enough for obtaining uniqueness of the solution (see [1], [23], or [18]). Non-homogeneous Dirichlet or natural boundary conditions may also be considered here, but only if they provide a positive solution. This restriction is fulfilled in many of the cases of practical interest.
360
I.S. Pop, W.-A. Yong
The assumption in (A1) implies that β has an inverse β −1 , which is continuously differentiable everywhere excluding 0. Then there are two positive constants such that 0 < C1 (ε) ≤ (β −1 ) (θ) ≤ C2 (ε) < ∞, for any real θ in [β(ε), β M eCT ], if this interval is included in the range of β. Here M is a large constant satisfying θ0 ≤ β(M ) − β(ε) almost everywhere and C is bigger than the constant in (A3). Moreover, C1 (ε) often depends on ε in a regular way. To see this, we refer to the typical porous medium case, where β(u) = um (m ≥ 1). This kind of dependence holds if β is Lipschitz continuous at least on bounded intervals. So we will use C1 to denote C1 (ε) and this is assumed in what follows. A usual definition of the solution for (1.1) is (1.5)
Definition 1.1. A function u is called a solution of the problem in (1.1) if u ∈ H 1 (0, T ; H −1 (Ω)), β(u) ∈ L2 (0, T ; H01 (Ω)), u(0) = u0 in H −1 and, for all ϕ ∈ L2 (0, T ; H01 (Ω)) the equation holds T T (∂t u(t), ϕ(t)) dt + 0 (∇β(u(t)) + F (u(t)), ∇ϕ(t))dt 0 (1.6) T = 0 (r(u(t)), ϕ(t))dt. Existence, uniqueness and boundedness of the solution for the above problem is studied in several papers (see, e.g., [1], [12], [18] and the references therein). Since the regularization approach considered here makes use of the maximum principle, we assume that the solution is uniformly bounded almost everywhere in the whole cylinder QT . This implies better regularity properties for u. Particularly, we get u ∈ C(0, T ; H −1 (Ω)), thus Problem (1.1) can be re-formulated in a stronger sense than it was considered in [1] and the initial condition holds in H −1 (Ω). In the rest of this paper we show the efficiency of the scheme in (1.2) from both theoretical and numerical point of view. Our numerical experiments presented later confirm that our approach is competitive. The main theoretical result is Theorem. Under the assumptions (A1), (A2) and (A3), if u is the weak solution of the problem in (1.1), then ε2 + (C2 (ε)β(ε))2 2 β(u) − θ∆ L2 (QT ) ≤ C τ + , τ where θ∆ (t) = θk for t ∈ (tk−1 , tk ] and k = 1, n. Remark 1.2. In the special case of a semigroup generated by a sub-gradient in a Hilbert space, based on the results in [20], the error estimates stated above can be improved.
A numerical approach to degenerate parabolic equations
361
Remark 1.3. If β satisfies ε · inf {β (x) : x ∈ [ε, M ]} ≥ Cβ(ε), which is true if β(u) = um with m ≥ 1, then the error estimate becomes β(u) − θ∆ 2L2 (QT ) ≤ C τ + ε2 + β 2 (ε) . Remark 1.4. As in (1.3), the initial data are perturbed only locally. This is essential for the method to be efficient. However, the above result remains valid also for a global perturbation (for example, v0 ≡ 1). We let (·, ·) stand for the inner product on L2 (Ω) or the duality pairing between H01 (Ω) and H −1 (Ω), · for the norm in L2 (Ω), · 1 and | · |1 for the respective norm and seminorm in H 1 (Ω). In addition, we often write u or u(t) instead of u(t, x) and use C to denote a generic positive constant independent of τ and ε. The paper is organized as follows. In the next section we study the nonlinear problems in (1.2). Sections 2 to 5 are devoted to the proof of the above theorem and in Sect. 6 some numerical experiments are presented. 2. Solving the elliptic problems In this section we study the nonlinear problems in (1.2), which, for each k, are elliptic and of Dirichlet type. First of all, the nonlinear problems are written in a certain variational form. As in the continuous case, it is necessary to give weak forms of the time discrete approximation schemes. Definition 2.1. Let k between 1 and n and θk−1 given. Find θk ∈ H01 (Ω) + β(ε) such that for all ϕ ∈ H01 (Ω), the equation (2.1)
(β −1 (θk ) − β −1 (θk−1 ), ϕ) + τ (∇θk + F (β −1 (θk )), ∇ϕ) = τ (r(β −1 (θk−1 )), ϕ)
holds. Remark 2.1. If F is assumed differentiable, the nonlinear convection ∇ · F (β −1 (θ)) can be rewritten as (β −1 ) (θ)f (β −1 (θ)) · ∇θ, which can be interpreted as the derivative of F with respect to θ. This term may generally go to infinity because of the derivative of β −1 . In connection with the maximum principle, this situation will be avoided due to the lower bound for θ. Taking into account the assumption made in (A3) on F , a simple calculation shows that for any θ above β(ε) the following holds |(β −1 ) (θ)f (β −1 (θ))| ≤ CF C2 (ε). If F is Lipschitz continuous with respect to β(u), C2 (ε) does not appear in the bound above.
362
I.S. Pop, W.-A. Yong
As stated before, the maximum principle plays a crucial role in the discretization method. For the initial problem written in variational form, this property has been studied in several papers ([17], see also [11] and the references therein). Our task here is to consider the elliptic problems defined above. To this end, we define, for each k less than or equal to n,
(2.2) Vk = ϕ ∈ β(ε) + H01 (Ω) : β(ε) ≤ ϕ ≤ β M eCr Lβ kτ , a.e. . It is clear that, for any k > 0, Vk includes any of the previous sets. The maximum principle is shown by mathematical induction. The assumption (A2) guarantees that, for M suitably chosen (e.g. M = u0 ∞ ), the initial data is in V0 . Now, assuming that θk−1 is given in Vk−1 , our aim is to show that if a solution θk of the problem given in Definition 2.1 exists, then it belongs to Vk . Fixing now k ≥ 1 and assuming θk−1 ∈ Vk−1 , the maximum principle is stated in the following lemma Lemma 2.1. Assume (A1) and (A3). Then, if a solution of the semi-discrete problem in (2.1) exists, it belongs to Vk . Proof. We start with the proof of the lower bound because this is essential in our approach. Since a solution θk ∈ H01 (Ω)
+ β(ε) for the problem given in Definition 2.1 is assumed to exist, ϕ = β(ε) − θk + ∈ H01 (Ω) can be taken as a test function and the following holds (I ) + (II) :=
(β −1 (θk ) − ε, β(ε) − θk + ) + τ (∇θk , ∇ β(ε) − θk + )
= (β −1 (θk−1 ) − ε, β(ε) − θk + ) + τ (r(β −1 (θk−1 )), β(ε) − θk + )
−τ (F (β −1 (θk )), ∇ β(ε) − θk + ) =: (III) + (IV ) + (V ). (2.3)
Let Ωε be the support of β(ε) − θk + and · Ωε the L2 norm on this subdomain. We have (I) = −(ε − β −1 (θk ), β(ε) − θk ) ≤ 0, +
because of the monotonicity of β. The same holds for (II), since (∇θk , ∇ β(ε) − θk ) = − ∇θk 2Ωε ≤ 0. +
Under the assumption (A3), for any real number δ, a vector valued function Gδ : R → Rd can be defined by θ (2.4) Fi (β −1 (s))ds, Gi (θ) = β(ε)
A numerical approach to degenerate parabolic equations
363
for all i = 1, d. Therefore we have Gi (β(ε)) = 0 and ∇x · G(θ) = F (β −1 (θ))∇x θ. Hence (III) vanishes due to the Gauß integration formula. Since θk−1 ≥ ε almost everywhere and r(u) ≥ 0 for all u ≥ 0, (IV ) and (V ) are also positive. Recalling the identity in (2.3) we get 0 ≥ − Ωε (β −1 (θk ) − ε)(θk − β(ε)) − τ ∇θ 2Ωε = Ωε (β −1 (θk−1 ) − ε + τ r(β −1 (θk−1 ))), β(ε) − θk ) ≥ 0.
Thus, Ωε (β −1 (θk )−ε)(θk −β(ε)) = 0, so either θk ≡ β(ε), or meas{Ωε } = 0, both leading to the desired lower bound. The upper bounds are obtained in a similar fashion. Take ϕ = k [θ − β(M eCr Lβ kτ )]+ in (2.1) and obtain as before 0 ≤ ΩM β −1 (θk ) − M eCr Lβ kτ θk − β(M eCr Lβ kτ ) + τ ∇θ 2ΩM = ΩM (β −1 (θk−1 ) − M eCr Lβ (k−1)τ , θk − β(M eCr Lβ kτ )) + ΩM (τ r(β −1 (θk−1 )) + M eCr Lβ (k−1)τ − M eCr Lβ kτ , θk − β(M eCr Lβ kτ )),
where ΩM is the support of θk − β(M eCr Lβ kτ ) + . Since r satisfies the assumption (A3), β is Lipschitz and θk−1 ≤ β(M eCr Lβ (k−1)τ ) almost ecerywhere in Ω, the last two integrals in the above identity are negative. Thus we get
β −1 (θk ) − M eCr Lβ kτ θk − β(M eCr Lβ kτ ) = 0, ΩM
and consequently the upper bound for θk . Remark 2.2. Following the ideas in [7] for proving of the weak maximum principle, the construction of a primitive function G can be avoided. The proof is more complicated, but can be useful if the convective term is not discretized implicitly (see also [19]). Remark 2.3. The problem in (1.1) considered here may be degenerate only at u = 0. The above maximum principle guarantees that if the (initial and boundary) data are away from 0, then the solution of any of the semi-discrete approximation defined before stays always away from 0 at all time steps. This is just the underlying idea of our approach to treat the degeneracy. Existence and uniqueness for the solution of the problems in (2.1) are a consequence of the theory of monotone operators, if τ is reasonably small. In this sense, let A : H01 (Ω) → H −1 be defined as Aθ = β −1 (θ + β(ε)) − τ ∇(∇θ + F (β −1 (θ + β(ε)))).
364
I.S. Pop, W.-A. Yong
If τ satisfies (2.5)
2 > CF τ,
A is strongly monotone. This results from (A θ − Aψ, θ − ψ) = (β −1 (θ + β(ε)) − β −1 (ψ + β(ε)), θ − ψ) + τ ∇(θ − ψ) 2 +τ (F (β −1 (θ + β(ε))) − F (β −1 (ψ + β(ε))), ∇(θ − ψ)) ≥ (β −1 (θ + β(ε)) − β −1 (ψ + β(ε)), θ − ψ) + τ ∇(θ − ψ) 2 −τ F (β −1 (θ + β(ε))) − F (β −1 (ψ + β(ε))) ∇(θ − ψ) ≥ (β −1 (θ + β(ε)) − β −1 (ψ + β(ε)), θ − ψ) + τ ∇(θ − ψ) 2 − τ2 CF (β −1 (θ + β(ε))) − β −1 (ψ + β(ε)), θ − ψ) − τ2 ∇(θ − ψ) 2 ≥ τ2 ∇(θ − ψ) 2 , where the inequality of Cauchy, (2.5) and (A3) have been used. Under the same assumptions, coercivity of A can be shown. Now, since β −1 is continuous, it is easy to see that A is hemicontinuous and therefore onto H −1 , which yields existence for the solution of problem Aθk = b, where the linear and bounded functional b is defined as b(ϕ) = (β −1 (θk−1 ), ϕ) + τ (r(β −1 (θk−1 )), ϕ). Moreover, uniqueness is a direct consequence of the strong monotonicity property of A. 3. Apriori estimates Convergence of the numerical approach is shown by proving the error estimates. To this aim, we are using some elementary identities, which are valid for all ak , bk ∈ Rq (q ≥ 1) (3.1)
2
m
2
2
ak (ak − ak−1 ) = |am | − |a0 | +
k=1 m k k=1 j=1
(3.3)
m
|ak − ak−1 |2 ,
k=1
2
(3.2)
m
m m 2 ak aj = ak + |ak |2 , k=1
ak (bk − bk−1 ) = am bm − a0 b0 −
k=1
The argument starts with a stability result.
k=1
m k=1
(ak − ak−1 )bk−1 .
A numerical approach to degenerate parabolic equations
365
Theorem 3.1. Assume (A1), (A2) and (A3). Then, for p ≤ n, if θk solves the problem given in Definition 2.1, we have τ
(3.4)
p
∇θk 2 ≤ C,
k=1
− β −1 (θk−1 ), θk − θk−1 ) + pk=1 θk − θk−1 2 +τ ∇θp 2 + τ pk=1 ∇(θk − θk−1 ) 2 ≤ Cτ C2 (ε).
p
k=1 (β
(3.5)
−1 (θ k )
Proof. By taking ϕ = θk − β(ε) ∈ H01 (Ω) in (2.1) and summing up for k = 0, p we get p p −1 k −1 k−1 ), θ k − β(ε)) + τ k 2 k=1 (β (θ ) − β (θ k=1 ∇θ (3.6)
= τ pk=1 −(F (β −1 (θk )), ∇θk ) + (r(β −1 (θk−1 )), θk − β(ε)) . Now, each of the resulting terms is estimated separately. To begin with, observe that b b C1 2 −1 −1 −1 b s(β ) (s)ds ≤ b(β (b) − β (a)), s(β −1 ) (s)ds ≥ 2 a 0 hold true for any real numbers a, b (due to the properties of β). Hence, for almost everywhere x ∈ Ω, it follows that k p p −1 k −1 k−1 ))θ k = k θ −1 k=1 (β (θ ) − β (θ k=1 θ θk−1 (β ) (s)ds θp θ0 p θk −1 −1 −1 ≥ k=1 θk−1 s(β ) (s)ds = 0 s(β ) (s)ds − 0 s(β ) (s)ds ≥ C(θp )2 − β −1 (θ0 )θ0 . Since u0 , θ0 ∈ L∞ , integrating over Ω the above inequality gives p −1 k −1 k−1 ), θ k ) ≥ C θ p 2 − (β −1 (θ 0 ), θ 0 ) k=1 (β (θ ) − β (θ ≥ C θp 2 − C ≥ −C. Using the maximum principle in Lemma 2.1, the remaining part of the first sum can be bounded as follows p (β −1 (θk ) − β −1 (θk−1 ), β(ε)) = (β −1 (θp ) − β −1 (θ0 ), β(ε)) k=1
≤ Cβ(ε). Next, the first term in the right hand side vanishes (as in (2.4)). For the last one, since θk is uniformly bounded (also independently on k) p −1 k−1 )), θ k − β(ε)) k=1 (r(β (θ ≤ τ pk=1 r(β −1 (θk−1 )) θk − β(ε) ≤ C.
366
I.S. Pop, W.-A. Yong
Since ε is small, we have β(ε) ≤ C. The above inequalities give the first part of the conclusion. For the second estimate, ϕ in (2.1) is replaced by θk − θk−1 ∈ H01 (Ω). Recalling the identity in (3.1) and summing again over k from 1 to p yields (I) + (II) := pk=1 (β −1 (θk ) − β −1 (θk−1 ), θk − θk−1 ) + τ2 ∇θp − ∇θ0 + pk=1 ∇(θk − θk−1 ) 2 (3.7) = −τ pk=1 (F (β −1 (θk )), ∇(θk − θk−1 )) +τ pk=1 (r(β −1 (θk−1 )), θk − θk−1 ) =: (III) + (IV ). For (I), because of the assumptions on β, the following inequality holds p p C1 k 1 −1 k k−1 2 (I) ≥ θ − θ + (β (θ ) − β −1 (θk−1 ), θk − θk−1 ). 2 2 k=1
k=1
Since ∇ · F (β −1 (θ)) = f (β −1 (θ))(β −1 ) (θ)∇θ, recalling the boundedness of f from Remark 2.1, for (III) we have |(III)| ≤ τ pk=1 |(f (β −1 (θk ))(β −1 ) (θk )∇θk , θk − θk−1 )| ≤ τ pk=1 CF C2 (ε) ∇θk θk − θk−1 . Applying here the inequality ab ≤ ηa2 + first estimates gives |(III)| ≤
1 2 4η b
for η =
2τ C1
and using the
p C1 k θ − θk−1 2 + Cτ C2 (ε). 8 k=1
Now, recalling the maximum principle again, the last term leads in a similar manner to p C1 k |(IV )| ≤ θ − θk−1 2 + Cτ. 8 k=1
Finally, since θ0 ∈ of the theorem.
H 1 (Ω),
the inequalities above show the remaining part
Remark 3.1. The large constant C2 (ε) appears in the last stability estimate due to the condition imposed in (A3) on F . As we will see below, this does not affect the error estimates. In the case F is Lipschitz continuous in β(u), since C2 (ε) disappears in (3.5), the apriori estimates become optimal (as obtained, e.g. in [16]).
A numerical approach to degenerate parabolic equations
367
Remark 3.2. By a different treatment of (III), apriori estimates without C2 (ε) can be obtained |(III)| ≤ τ pk=1 |(F (β −1 (θk )), ∇(θk − θk−1 ))| ≤ Cτ pk=1 ∇(θk − θk−1 ) ≤ C + τ4 pk=1 ∇(θk − θk−1 ) 2 , where the uniform bounds for θk have been used again. The rest of the proof follows as before and leads to p p −1 k −1 k−1 ), θ k − θ k−1 ) + k k−1 2 k=1 (β (θ ) − β (θ k=1 θ − θ (3.8) +τ ∇θp 2 + τ pk=1 ∇(θk − θk−1 ) 2 ≤ C. In this case the results are worse, since τ is lost. In addition we need the following lemma. Lemma 3.1. For u ≥ 0, define uε = β −1 (β(u) + β(ε)). If β satisfies (A1), then (3.9)
0 < uε − u ≤ ε + C2 (ε)β(ε),
where C2 (ε) is defined in (1.5). Proof. Because β −1 as well as β is increasing and ε > 0, we have u = β −1 (β(u)) < β −1 (β(u) + β(ε)) ≡ uε . On the other hand, it follows from Taylor’s theorem that uε ≡ β −1 (β(u) + β(ε)) = β −1 (β(max(u, ε))) + β(min(u, ε))(β −1 ) (θ), where θ is a convex combination of β(u) + β(ε) and β(max(u, ε)), which are both not less than β(ε) due to u ≥ 0. Thus, it is clear that β −1 (β(max(u, ε))) ≤ u + ε, β(min(u, ε)) ≤ β(ε), (β −1 ) (θ) ≤ C2 (ε) and hence
uε ≤ u + ε + C2 (ε)β(ε).
Remark 3.3. If β satisfies ε · inf β (x) : x ∈ [ε, M ] ≥ Cβ(ε), which is true if β(u) = um with m ≥ 1, then the above inequality becomes 0 < uε − u ≤ Cε. Similarly, in case that β is super-linear - satisfying β(u) + β(v) ≤ β(u + v) for any u, v - which holds again for the typical case mentioned above, Lemma 3.1 becomes 0 < uε − u ≤ ε. Both cases will simplify the error estimates.
368
I.S. Pop, W.-A. Yong
4. Error estimates For the error estimates let us set for any function f which is integrable in time and defined in QT 1 kτ k ¯ f (s, ·)ds, f := τ (k−1)τ if k ≥ 1 and f¯0 := f (0, ·). We will make use of the following notations k
ε − β −1 (θ k ), ¯k − β −1 (θk ), eε,k eku := u u := u k
k
k ε ekθ := β(u) − θk , eε,k θ := β(u ) − θ ,
where k ≥ 0. Using Lemma 3.1, because of the maximum principle, we have for any k ≥ 0 (4.1)
k eku ≤ eε,k u ≤ eu + ε + C2 (ε)β(ε), k ekθ ≤ eε,k θ ≤ eθ + β(ε).
Moreover, recalling the definition of θ0 in (1.2), Lemma 3.1 shows that initially these errors satisfy (4.2)
−ε + C2 (ε)β(ε) ≤ e0u ≤ 0 ≤ eε,0 u ≤ C2 (ε)β(ε), −β(ε) ≤ e0θ ≤ 0 ≤ eε,0 θ ≤ β(ε).
Because of the definition of uε , for k ≥ 0 ∇ekθ = ∇eε,k θ holds. The analysis below combines the approach in [15] with those in [6] and [16]. We have not considered the possibility given by the non-degeneracy property of the solution, since this kind of results are shown up to now only in some particular cases. In the sequel G : H −1 (Ω) → H01 (Ω) denotes the Green operator defined by (4.3) (∇Gψ, ∇ϕ) = (ψ, ϕ), for all ϕ ∈ H01 (Ω), where ψ is taken in H −1 (Ω). In the right hand side (·, ·) stands for the duality pairing between H −1 and H01 . Obviously, G is linear. Due to the Poincar´e-Friedrichs inequality, the expression ψ −1 =
|(ψ, ϕ)| ϕ∈H 1 , ϕ=0 ∇ϕ sup 0
A numerical approach to degenerate parabolic equations
369
is a norm in H −1 (Ω) equivalent to the usual one. Recalling the definition of the operator G we obtain ∇Gψ 2 = (∇Gψ, ∇Gψ) = (ψ, Gψ) ≤ ψ −1 ∇Gψ , hence ∇Gψ ≤ ψ −1 . The inequality in the inverse sense is a direct consequence of the same definitions. Moreover, if ψ ∈ L2 (Ω), the inequality ψ −1 ≤ C ψ is a direct consequence of the inequality mentioned above. Thus we have shown that ∇Gψ = ψ −1 , ψ −1 ≤ C ψ (4.4) (where the last inequality applies only if ψ ∈ L2 (Ω)). The error estimates are obtained in the following theorem. Theorem 4.1. Assume (A1), (A2) and (A3) and let u be the weak solution of the problem in Definition 1.1 and θk , k = 1, n solve the problems in Definition 2.1, then T 2 ε ε −1 sup eε,k u −1 + 0 (β(u (t)) − θ∆ (t), u (t) − β (θ∆ (t)))dt k=1,n ε2 + (C2 (ε)β(ε))2 2 , + β(u) − θ∆ L2 (Q) ≤ C τ + τ where θ∆ (t) = θk for t ∈ (tk−1 , tk ]. Proof. Let χ|I be the characteristic function of the time interval I ⊂ [0, T ]. Choosing ϕχ|[tj−1 ,tj ] for an arbitrary ϕ ∈ H01 (Ω) as a test function in (1.6) leads to the semi-discrete equations satisfied by the solution
tj tj (u(tj ) − u(tj−1 ), ϕ) + ∇ tj−1 β(u(t))dt + tj−1 F (u(t))dt, ∇ϕ
(4.5) tj r(u(t))dt, ϕ , = tj−1 for all ϕ ∈ H01 (Ω) and j > 0. Subtracting (2.1) from (4.5), taking ϕ = 1 Geε,j u ∈ H0 in the resulting difference and summing up for j = 1, k yields (I1 ) + (I2 ) := k −1 j −1 j−1 ), Geε,j ) u j=1 (u(tj ) − u(tj−1 ) − β (θ ) + β (θ k ε,j ε,j +τ j=1 (∇eθ , ∇Geu )
(4.6) tj F (u(t)) − F (β −1 (θj ))dt, ∇Geε,j = − kj=1 tj−1 u k tj ε,j −1 j−1 + j=1 tj−1 r(u(t)) − r(β (θ ))dt, Geu =: (I3 ) + (I4 ).
370
I.S. Pop, W.-A. Yong
Now we have to estimate each of the terms in (4.6). First we consider (I1 ) and decompose it as follows. k
(I1 ) =
j=1 (u(tj )
k
−u ¯j − u(tj−1 ) + u ¯j−1 , Geε,j u )
j j−1 uj − uε − u ¯j−1 − uε , Geε,j u ) j=1 (¯ k ε,j−1 + j=1 (eε,j , Geε,j u − eu u )
+
=: (I11 ) + (I12 ) + (I13 ). Recalling the elementary identity in (3.3), since u(t0 ) = u ¯0 , (I11 ) can be rewritten as (I11 ) = (u(tk ) − u ¯k , Geε,k u )−
k
j=2 (u(tj−1 )
ε,j−1 −u ¯j−1 , Geε,j ) u − Geu
=: (I111 ) + (I112 ). δb2 a2 + Because ∂t u ∈ L2 (0, T ; H −1 (Ω)), using the inequality ab ≤ δ 4 with δ > 0 we have tk ε,k (u(t ) − u(t), Ge ) u dt k tk−1 tk tk ε,k ≤ τ1 tk−1 t (∂s u(s), Geu )ds dt tk tk ε,k ≤ τ1 tk−1 t ∂s u −1 ∇Geu dsdt
|(I111 )| ≤ (4.7)
1 τ
≤ τ δ111 ∂t u 2L2 (tk−1 ,tk ;H −1 (Ω)) +
ε,k 2 1 4δ111 eu −1 ,
with δ111 > 0 given below. The estimation of (I112 ) follows similarly |(I112 )| ≤ (4.8)
1 τ
tj−1 ε,j ε,j−1 (u(t ) − u(t), Ge − Ge ) dt u u j−1 j=2 tj−2
k
≤ τ δ112 ∂t u 2L2 (0,tk−1 ;H −1 (Ω)) +
k 1 ε,j eu − eε,j−1 2−1 , u 4δ112 j=2
The bounds for (I12 ) are obtained in the same manner. k
0
uk − uε , Geε,k u0 − uε , Geε,0 (I12 ) = (¯ u ) − (¯ u ) k j−1 ε,j−1 uj−1 − uε , Geε,j ) − j=1 (¯ u − Geu =: (I121 ) + (I122 ) + (I123 ).
A numerical approach to degenerate parabolic equations
371
Employing the relations in (4.1), (4.2), Lemma 3.1 yields k
ε,k |(I121 )| ≤ ¯ uk − uε Geε,k u ≤ C(ε + C2 (ε)β(ε)) eu −1
≤ C(ε + C2 (ε)β(ε))2 δ121 +
ε,k 2 1 4δ121 eu −1 ,
(4.9) |(I122 )| ≤ C(ε + C2 (ε)β(ε))2 , j−1 ε,j−1 uj−1 − uε Geε,j |(I123 )| ≤ kj=1 ¯ u − Geu 2 ε,j−1 2 δ123 + 4δ1123 kj=1 eε,j −1 . ≤ C (ε+C2 (ε)β(ε)) u − eu τ Here we have used the boundedness of Ω and C is a generic constant. Using the identity in (3.1) and the properties of G, the estimate of (I13 ) is a simple matter, ε,j−1 (I13 ) = kj=1 (∇Geε,j , ∇Geε,j u − ∇Geu u )
(4.10) ε,0 2 ε,j ε,j−1 2 k 2 e − e = 12 eε,k u −1 − eu −1 + u u −1 . j=1 Now we proceed with the second term in (4.6), ε,j (I2 ) = τ kj=1 (eε,j θ , eu ) tj tj = kj=1 ( tj−1 (β(uε (t)) − θj )dt, τ1 tj−1 (uε (s) − β −1 (θj ))ds) tj = kj=1 tj−1 (β(uε (t)) − θj , uε (t) − β −1 (θj ))dt tj tj + kj=1 tj−1 (β(uε (t)) − θj , τ1 tj−1 (uε (s) − u(s))ds)dt tj ε j=1 tj−1 (β(u (t)) tj + kj=1 tj−1 (β(uε (t)) −
k
− θj , uε (t) − u(t))dt tj − θj , τ1 tj−1 (u(s) − u(t))ds)dt
=: (I21 ) + (I22 ) + (I23 ) + (I24 ). For (I21 ), recalling (1.5), we have tj (β(uε (t)) − θj , uε (t) − β −1 (θj ))dt (I21 ) ≥ 12 kj=1 tj−1 (4.11) tj + C21 kj=1 tj−1 β(uε (t)) − θj 2 dt. The estimates of (I22 ) and (I23 ) are identical and can be obtained as for (I123 ). |(I22 )|, |(I23 )| (4.12)
tj ≤ C(ε + C2 (ε)β(ε)) kj=1 tj−1 β(uε (t)) − θj dt tj ≤ δC2 kj=1 tj−1 β(uε (t)) − θj 2 dt + δ42 (ε + C2 (ε)β(ε))2 .
372
I.S. Pop, W.-A. Yong
In order to obtain upper bounds for (I24 ) we repeat the steps for (I112 ). Remembering the apriori estimates in Theorem 3.1 it follows that θj − β(ε) ∈ L2 (tj−1 , tj ; H01 (Ω)). Since ∂t u ∈ L2 (0, T ; H −1 (Ω)) and β(u) ∈ L2 (0, T ; H01 (Ω)), we have tj tj ε (t)) − θ j , s ∂ u(r)dr dsdt |I24 | ≤ τ1 kj=1 tj−1 β(u r t tj−1 k tj tj (4.13) ≤ j=1 tj−1 tj−1 ∇(β(uε (t)) − θj ) ∂r u(r) −1 drdt ≤ Cτ. Considering now the right hand side in (4.6), (I3 ) can be splitted into two terms. tj (F (u(t)) − F (uε (t)), ∇Geε,j (I3 ) = − kj=1 tj−1 u )dt k tj − j=1 tj−1 (F (uε (t)) − F (β −1 (θj )), ∇Geε,j u )dt =: (I31 ) + (I32 ). For the first term, because of the assumption (A3) on F , Lemma 3.1 can be used again in order to get tj F (u(t)) − F (uε (t)) ∇Geε,j |(I31 )| ≤ kj=1 tj−1 u dt 1 tj ε ε (4.14) ≤ C kj=1 eε,j u −1 tj−1 (u(t) − u (t), β(u(t)) − β(u (t))) 2 dt 2 ≤ δC31 β(ε)(ε + C2 (ε)β(ε)) + Cτ δ31 kj=1 eε,j u −1 . Analogous, for (I32 ) we can obtain |(I32 )| (4.15)
tj ε,j ε −1 j ε j 12 j=1 eu −1 tj−1 (u (t) − β (θ ), β(u (t)) − θ ) dt tj C k ε −1 j ε j j=1 tj−1 (u (t) − β (θ ), β(u (t)) − θ )dt 4δ32 2 +Cτ δ32 kj=1 eε,j u −1 .
≤C ≤
k
(I4 ) can be decomposed as
tj r(u(t)) − r(β −1 (θj ))dt, Geε,j (I4 ) = kj=1 tj−1 u +τ kj=1 (r(β −1 (θj−1 )) − r(β −1 (θj )), Geε,j u ) =: (I41 ) + (I42 ). Proceeding exactly in the same manner as before, the inequalities in (4.4) lead, for (I41 ), to estimates similar to the ones in (4.15). Moreover, the assumption on r in (A3) and the apriori estimates in (3.8) can be used for bounding (I42 ) −1 j −1 j−1 ), θ j − θ j−1 ) 12 |(I42 )| ≤ Cτ kj=1 eε,j u −1 (β (θ ) − β (θ (4.16) 2 ≤ Cτ + Cτ kj=1 eε,j u −1 .
A numerical approach to degenerate parabolic equations
373
Inserting the inequalities in (4.7) - (4.16) in (4.6), recalling (4.2) and choosing the δ’s properly, we get k tj 2 ε j ε −1 j eε,k u −1 + j=1 tj−1 (β(u (t)) − θ , u (t) − β (θ ))dt tj ε,j−1 2 + kj=1 eε,j −1 + C1 kj=1 tj−1 β(uε (t)) − θj 2 dt u − eu 2 2 + Cτ kj=1 eε,j ≤ Cτ ∂t u 2L2 (0,tk ;H −1 (Ω)) + Cτ + C (ε+C2 (ε)β(ε)) u −1 . τ Since u ∈ H 1 (0, T ; H −1 (Ω)), application of the discrete Gronwall inequality yields 2 eε,k u −1 +
tk
(β(uε ) − θ∆ , uε − β −1 (θ∆ )) (4.17) (ε + C2 (ε)β(ε))2 ε 2 + β(u ) − θ∆ L2 (0,tk ;L2 (Ω)) ≤ C τ + τ 0
for any k ≥ 0. Finally, we notice that β(uε ) − θj = β(u) − θj + β(ε) ≥ β(u) − θj − β(ε) ≥ β(u) − θj − Cβ(ε) and therefore (4.18)
2 β(uε ) − θj 2 ≥ β(u) − θj 2 − Cβ 2 (ε).
This, together with (4.17), gives the desired result.
Remark 4.1. An alternative proof is given in [24] which relies exclusively on the method proposed in [15] and can be applied only if the Problem (1.1) does not contain a convective term. Having now the result of Theorem 4.1 available, this procedure can be repeated in order to get error estimates involving H 1 (Ω) norms (see, e.g., [16]). Remark 4.2. In the situation described in Remark 3.3, the errors in Theorem 4.1 are bounded by C(τ + ε2 /τ ). Remark 4.3. The error estimates obtained above show that the numerical scheme considered here behaves at least as good as the algorithms in [15], [16], [9], [2] or [5]. In fact, the schemes in [15], [16] or [2] are proved to have a convergence order τ 1/2 , while the order τ 1/4 is demonstrated in [8], for the method proposed in [9]. For the nonlinear scheme in [5] the same error behaves asymptotically like τ 0.3 . However, under certain circumstances, optimal error estimates are obtained in the next section.
374
I.S. Pop, W.-A. Yong
Remark 4.4. The error estimates for u are obtained in the norm of L2 (0, T ; H −1 (Ω)). This can be improved if an inequality of the form (β −1 (θ) − β −1 (ψ))(θ − ψ) ≥ C(β −1 (θ) − β −1 (ψ))p holds for a positive constant C and an exponent p > 1. Then, using the estimate for the scalar product in (4.17), an error estimate in the better Lp (Q) norm can be obtained. For example, if β(u) = um , recalling the inequalities in (4.1), the estimate reads 2 2 ≤ C τ + ε + (C (ε)β(ε)) u − β −1 (θ∆ ) m+1 2 m+1 L (QT ) (see, e.g. [16], [5]).
5. Optimal estimates In this section we simplify the problem in (1.1), so that no convection or reaction is present. Then it is possible to improve the error estimates. We deal here with the problem ∂t u − ∆β(u) = 0, (5.1)
u(x, 0) = u0 (x) ≥ 0, u = 0,
in QT ≡ (0, T ) × Ω, in Ω, on ∂Ω,
the assumptions on β and u0 remaining the same. This problem can be formulated in terms of the semigroup theory As shown in [20], if the semigroup is generated by a sub-gradient in a Hilbert space, the backward Euler method has an O(τ ) order of convergence. Since β is maximal monotone, supposing its range is all of R, Rulla’s result applies also here. Recall now the Euler implicit discretization for the simplified problem in (5.1) (5.2)
uk − uk−1 = τ ∆β(uk ), β(uk )|∂Ω = 0
for k = 1, n with u0 from (5.1). In this case, uk approximates u(kτ ). Replacing uk with β −1 (θk ), the scheme in 1.2 is identical to the one above up to the initial and boundary data. In (1.2) these are subject to a small shift in order to get a regular parabolic problem, which is not the case for (5.1). The result in [20] establishing the optimal rate of convergence for the implicit scheme reads
A numerical approach to degenerate parabolic equations
375
Theorem 5.1. Let {uk , k = 0, n} be the solutions of the implicit scheme in (5.2) and uτ , respectively θτ the piece-wise constant functions interpolating {uk , k = 1, n} and {β(uk ), k = 1, n}. Then, for any k ≥ 0, T sup u(tk ) − uk 2−1 + 0 (u(t) − uτ (t), β(u(t)) − θτ (t))dt k=1,n
+τ β(u) − θτ 2L2 (0,T ;H 1 (Ω)) ≤ Cτ 2 , 0
where C depends on the
H1
norm of β(u0 ).
It is worth here to notice that the second term in the above inequality is positive due to the monotonicity of β. Now, if u0 is a positive L∞ (Ω) function, a semi-discrete maximum principle can be obtained easily, thus uk stays above 0 for any k > 1. Based on the previous estimates, we can get similar results for the numerical scheme in (1.2). Theorem 5.2. In the setting of Theorem 4.1, if θk , k > 0 solve the problems given in Definition 2.1 without convection or reaction, then β(u) − θ∆ 2L2 (Q) + τ β(u) − θ∆ 2L2 (0,T ;H 1 (Ω)) 0 2 2 2 ≤ C τ + ε + (C2 (ε)β(ε)) . Proof. First, the scheme in (5.2) is written in a weak form (uj − uj−1 , ϕ) + τ (∇β(uj ), ∇ϕ) = 0 (for all ϕ ∈ H01 (Ω), with j = 1, n). We want to compare the approximations generated by the two implicit numerical schemes. Subtracting (2.1) from the above equality, summing up for j = 1, k, taking ϕ = β(uε,k ) − θk ≡ β(uk ) + β(ε) − θk ∈ H01 in the resulting difference and summing up again for k = 1, n gives (I1 ) + (I2 ) := nk=1 (uk − β −1 (θk ), β(uε,k ) − θk ) (5.3) + τ nk=1 kj=1 (∇(β(uε,j ) − θj ), ∇(β(uε,k ) − θk )) = nk=1 (u0 − β −1 (θ0 ), β(uε,k ) − θk ) =: (I3 ). We have to estimate each of the terms above. First, (I1 ) can be decomposed into (I1 ) = nk=1 (uε,k − β −1 (θk ), β(uε,k ) − θk ) + nk=1 (uk − uε,k , β(uε,k ) − θk ) =: (I11 ) + (I12 ). Recalling the relation in (1.5), (I11 ) gives (I11 ) ≥ C1 nk=1 β(uε,k ) − θk 2 (5.4) ≥ C21 nk=1 β(uk ) − θk 2 −
C 2 τ β(ε) ,
376
I.S. Pop, W.-A. Yong
where the inequality in (4.18) has been applied here. As a consequence of Lemma 3.1 and of the maximum principle for the solution of the implicit scheme, (I11 ) can be bounded as |(I12 )| ≤ C(ε + C2 (ε)β(ε)) nk=1 β(uε,k ) − θk (5.5) ≤ Cτ (ε + C2 (ε)β(ε))2 + C81 nk=1 β(uk ) − θk 2 . For (I2 ) we make use of the identity in (3.2) in order to obtain n n τ τ ε,k k 2 (5.6) (I2 ) = ∇(β(u ) − θ ) + ∇ (β(uε,k ) − θk ) 2 . 2 2 k=1
k=1
Now, because of (4.2), (I3 ) can be bounded as (I12 ). Inserting all the inequalities in (5.4) - (5.6) in (5.3) and multiplying everything with τ we arrive at C1 n τ 2 n k k 2 ε,k k 2 k=1 τ β(u ) − θ + 2 k=1 ∇(β(u ) − θ ) 4 (5.7) ≤ C(ε + C2 (ε)β(ε))2 . The remaining part of proof can be completed using the above relation, Theorem 5.1 and the triangle inequality of norms ( a + b ≤ a + b ). Remark 5.1. This result shows that if the implicit Euler discretization has a linear order of convergence, the same holds also for the maximum principle based numerical scheme. Remark 5.2. Even though the above estimates can be extended to other cases (for example, monotone perturbations of a subgradient, [20]), they cannot cover the whole range of problems considered in this paper. Therefore we have maintained the analysis in the previous sections and complete it here with optimal results. 6. Numerical results In this section we give some numerical results to show the efficiency of the numerical approach in (1.2). The test problem is the well-known porous medium equation in the d-dimensional case, (6.1)
∂t u = ∆(um ),
with m = 1 + α. Homogeneous boundary data of Dirichlet type are considered in the domain Ω = [−L, L]d . For a particular choice of the initial data, an explicit solution was given by G. I. Barenblatt in [3], namely 1 1 2 m−1 1 1 2 − 2m x 1 − (6.2) u(t, x) = 2 d (md + 2 − d)(t + 1) md+2−d + (t + 1) (md+2−d)
A numerical approach to degenerate parabolic equations
377
Since the time-discrete problems are nonlinear, the following iteration scheme is used in the numerical examples for linearization purposes
(6.3)
θ¯i ∈ β(ε) + H01 (Ω), (σ(θ¯i−1 , θk−1 )(θ¯i − θk−1 ), ϕ) + τ (∇θ¯i , ∇ϕ) = τ ((β −1 ) (θ¯i−1 )f (β −1 (θ¯i−1 )) · ∇θ¯i , ϕ) + τ (r(β −1 (θk−1 )), ϕ), 1 σ(θ¯i , θk−1 ) = 0 β −1 (sθ¯i + (1 − s)θk−1 )ds
for all ϕ ∈ H01 (Ω) and i ≥ 1, where θ¯0 = θk−1 , σ(θ¯0 , θk−1 ) = β −1 (θk−1 ). This scheme is a linearization of the implicit time discretization in (1.2) and reproduces the one proposed in [9] and [10] for a particular choice of the parameters. However, in that case a cut-off approach is necessary in order to control the derivative of β −1 , which can become infinite. This situation is avoided here through the perturbation of the data, so that all the iterative solutions θ¯i , i ≥ 0, are bounded away from the degeneracy point, 0. This is due to the maximum principle, which holds also in this case. Since our approach relies on the maximum principle, an appropriate spatial discretization has to be considered. Here we use the lumped mass version of the finite element method (see [22]) on a regular decomposition of the domain Ω into d-simplices. If the decomposition is weakly acute, such a discretization preserves the maximum principle, which is essential in our approach. Although other discretization techniques may be combined with the time discretization in (1.2) (e.g. the adaptive mesh refinement near the free boundary), we do not discuss them here since the purpose of this paper is to provide the basic idea. However, if convection terms are also present, some upwinding techniques have to be considered ([19]). Convergence of the iteration in (6.3) is shown in [13], while the fully discrete case is studied in [19]. In our computations, at most 6 iterations were sufficient for approximating the solution of the nonlinear problem at each time step. The examples are implemented in U G ([4]) on an SGI O2 computer. The fully discrete linear problems are solved by multigrid procedures provided by the above mentioned package. The computations are done in two spatial dimensions (d = 2), we refer to [24] for one-dimensional examples. We have chosen L = 7, so the domain Ω includes the support of the solution. Two exponential nonlinearities are considered, m = 2 and m = 6. v0 is chosen as in (1.4). Since an exact solution is known in this case, the convergence order of the method can be evaluated. To this aim we approximate the L2 (QT ) errors
378
I.S. Pop, W.-A. Yong Table 1. L2 errors in u and θ, t = 1.0, m = 2 i
Eu
Eθ
αu
αθ
0 1 2
5.04 · 10−1 2.46 · 10−1 1.23 · 10−1
3.06 · 10−2 9.31 · 10−3 3.93 · 10−3
∗ 1.03 1.00
∗ 1.71 1.24
Table 2. L2 errors in u and θ, t = 1.0, m = 6. i
Eu
Eθ −1
0 1 2
7.47 · 10 5.38 · 10−1 3.10 · 10−1
−2
7.89 · 10 4.33 · 10−2 1.47 · 10−2
αu
αθ
∗ 0.47 0.80
∗ 0.86 1.56
in terms of u and β(u) by 1 N h
2 2 τ |uk,h (Ai ) − u(kτ, Ai )| meas(BAi ) , Eu := i=1
respectively 1 N h
2 2 τ |θk,h (Ai ) − β(u)(kτ, Ai )| meas(BAi ) , Eθ := i=1
where {Ai , i = 1, Nh } are the nodes of the triangularization and BA stands for the dual box centered in A, while meas(BA ) stands for the volume of BA . Assuming the above errors are of order τ α - more precisely Eu = Cu τ αu and Eθ = Cθ τ αθ for some exponents αu , αθ and constants Cu , Cθ - a reasonable evaluation is offered by the computation of approximations corresponding to different sets of discretization parameters τj , hj and εj . Now the convergence orders can be estimated from the relations αu =
ln(Eu1 /Eu2 ) ln(Eθ1 /Eθ2 ) and αθ = . ln(τ1 /τ2 ) ln(τ1 /τ2 )
We have avoided a comparison with other discretization methods because this always depends on the particular choices of the parameters. However, tests in one spatial dimension show that for the porous medium equation, the results obtained with the methods considered here are similar to the ones produced by the J¨ager-Kaˇcur scheme (more details can be found in [24]). All the results presented here are obtained at t = 1.0. The largest time step is τ = 0.04, the same as the perturbation parameter ε. We start with a uniform spatial grid which is refined two times. A finer computation halves τ
A numerical approach to degenerate parabolic equations
379
Fig. 1. Approximation of u (left) and β(u) (right), t = 1.0 and m = 2
and ε and refines the mesh once more. The first column of the tables contain the level of refinement, denoted by i. Table 1 and Table 2 present the errors Eu and Eθ for the maximum principle based approach together with the estimated order of convergence for the fully discrete scheme. Here αu and αθ stand for the order of the errors Eu , respectively Eθ . As we can see in the tables, the convergence rates for the more regular variable (β(u)) lie above the estimated one, 1. A better behaviour can be noticed if the nonlinearity is milder, namely m = 2. Then u itself is H 1 with respect to the spatial variable (see, e.g., [11]), implying better convergence orders also for Eu . The situation is different for m = 6, when only Eθ is of order O(τ ). The numerical approximation for the solutions u and β(u) are presented in Fig. 1 for m = 2 and Fig. 2 for m = 6. The pictures are zoomed 15 times in the z direction. We give here only the results on the finest mesh. Figure 3 displays the absolute error (zoomed 30 times in the z-direction) for the fully discrete implicit scheme after 1.0s, where the nonlinearity is β(u) = u2 .
Fig. 2. Approximation of u (left) and β(u) (right), t = 1.0 and m = 6
380
I.S. Pop, W.-A. Yong
Fig. 3. Errors in u (left) and θ (right), m = 2
Fig. 4. Errors in u (left) and θ (right), m = 6
These range until 2 · 10−2 , respectively 1 · 10−3 . The results for m = 6 are shown in Fig. 4. The errors are less than 0.13, respectively 2 · 10−2 . It becomes clear that in this case the errors are larger, because of the missing regularity of the solution u. It is worth here to notice that the errors appear mainly around the free boundary (namely inside a region covered by two, maximum three elements). Acknowledgements. We would like to express our thanks to Prof. W. J¨ager for suggesting us to consider this kind of problems, to Prof. J. Kaˇcur, Dr. N. Neuß and Mr. S. Schnadt for useful discussions. Particular thanks to Mr. Florin Radu for pointing us out a correction needed in Theorem 4.1. The suggestions from the referees helped us in improving the final form of the paper. This work was supported by the Deutsche Forschungsgemeinschaft through SFB 359 and by the Interdiscliplinary Center for Scientific Computing (IWR) at the University of Heidelberg.
References 1. H. W. Alt, S. Luckhaus: Quasilinear elliptic-parabolic differential equations. Math. Z. 183, 311–341 (1983) 2. T. Arbogast, M. F. Wheeler, N.-Y. Zhang: A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media. SIAM J. Numer. Anal. 33, 1669–1687 (1996)
A numerical approach to degenerate parabolic equations
381
3. G. I. Barenblatt: On some unsteady motion of a liquid or a gas in a porous medium. Prikl. Math. Meh. 16, 67–78 (1952) 4. P. Bastian, K. Birken, K. Johannsen, S. Lang, N. Neuss, H. Rentz-Reichert, C. Wieners: UG – a flexible Software toolbox for solving partial differential equations. Computing and Visualization in Science 1, 27–40 (1997) 5. C. Ebmeyer: Error estimates for a class of degenerate parabolic equations. SIAM J. Numer. Anal. 35, 1095–1113 (1998) 6. C. M. Elliott: Error analysis of the enthalpy method for the stefanproblem. IMA J. Numer. Anal. 7, 61–71 (1987) 7. D. Gilbarg, N. S. Trudinger: Elliptic Partial Differential Equations of Second Order. Grundlehren der Mathematischen Wissenschaften, 224, 2nd ed., Berlin; Springer 1983 8. A. Handloviˇcov´a: Error estimates of a linear approximation scheme for nonlinear diffusion problems. Acta Math. Univ. Comenianae 61, 27–39 (1992) 9. W. J¨ager, J. Kaˇcur: Solution of porous medium type systems by linear approximation schemes. Numer. Math. 60, 407–427 (1991) 10. W. J¨ager, J. Kaˇcur: Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes, M2 AN. Math. Model. Numer. Anal. 29, 605–627 (1995) 11. W. J¨ager, Y. G. Lu: H¨older continuity of solutions of degenerate parabolic equations, to appear 12. J. W. Jerome: Approximation of Nonlinear Evolution Systems. New York: Academic Press 1983 13. J. Kaˇcur: Solution to strongly nonlinear parabolic problems by a linear approximation scheme. IMA J. Numer. Anal. 19, 119–145 (1999) 14. J. Kaˇcur, A. Handloviˇcov´a, M. Kaˇcurov´a: Solution of nonlinear diffusion problems by linear approximation schemes. SIAM J. Numer. Anal. 30, 1703–1722 (1993) 15. E. Magenes, R. H. Nochetto, C. Verdi: Energy error estimates for a linear scheme to approximate nonlinear parabolic problems, M2 AN. Math. Model. Numer. Anal. 21, 655–678 (1987) 16. R. H. Nochetto, C. Verdi: Approximation of degenerate parabolic problems using numerical integration. SIAM J. Numer. Anal.25, 784–814 (1988) 17. O. A. Oleinik, A. S. Kalashnikov, Y.-L. Zhou: The Cauchy problem and boundary value problems for equations of the type of unsteady filtration. Izv. Akad. Nauk SSSR Ser. Mat. 22, 667–704 (1958) 18. F. Otto: L1 -contraction and uniqueness for quasilinear elliptic-parabolic equations. J. Differ. Equations 131, 20–38 (1996) 19. I. S. Pop: Regularization methods in the numerical analysis of some degenerate parabolic equations. SFB 359 Preprint 98-43, (1998), IWR, University of Heidelberg 20. J. Rulla: Error analysis for implicit approximations to solutions to Cauchy problems. SIAM J. Numer. Anal. 33, 68–87 (1996) 21. M. Slodiˇcka: Solution of nonlinear parabolic problems by linearization. Preprint M392, (1992), Comenius University Bratislava, Faculty of Mathematics and Physics 22. V. Thom´ee: Galerkin Finite Element Methods for Parabolic Problems. Lecture Notes in Math. 1054, Berlin: Springer 1984 23. M. Watanabe: An approach by difference to the porous medium equation with convection. Hiroshima Math. J. 25, 623–645 (1995) 24. W. A. Yong, I. S. Pop: A numerical approach to porous medium equations. SFB-359 Preprint 96-50, (1996), IWR, University of Heidelberg
