A modified L-scheme to solve nonlinear diffusion problems

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A modified L-scheme to solve nonlinear diffusion problems ∗

K. Mitra a,b , , I.S. Pop b,c a b c

Eindhoven University of Technology, Department of Mathematics and Computer Science, Netherlands Hasselt University, Faculty of Science, Belgium University of Bergen, Department of Mathematics, Norway

article

info

Article history: Available online xxxx Keywords: Nonlinear diffusion problem Linearization Newton, Picard, L-scheme Unconditional convergence Stability Richards equation

a b s t r a c t In this work, we propose a linearization technique for solving nonlinear elliptic partial differential equations that are obtained from the time-discretization of a wide variety of nonlinear parabolic problems. The scheme is inspired by the L-scheme, which gives unconditional convergence of the linear iterations. Here we take advantage of the fact that at a particular time step, the initial guess for the iterations can be taken as the solution of the previous time step. First it is shown for quasilinear equations that have linear diffusivity that the scheme always converges, irrespective of the time step size, the spatial discretization and the degeneracy of the associated functions. Moreover, it is shown that the convergence is linear with convergence rate proportional to the time step size. Next, for the general case it is shown that the scheme converges linearly if the time step size is smaller than a certain threshold which does not depend on the mesh size, and the convergence rate is proportional to the square root of the time step size. Finally numerical results are presented that show that the scheme is at least as fast as the modified Picard scheme, faster than the L-scheme and is more stable than the Newton or the Picard scheme. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction In this paper, a linearization technique is considered for the generalized nonlinear advection diffusion equations of the type

[ ] ∂t b(u) + ∇ · F⃗ (⃗x, u) = ∇ · D(⃗x, u)∇ u + f (⃗x, t , u);

(1.1)

completed with suitable boundary and initial conditions. Eq. (1.1) appears as mathematical model for many real world applications, like flow through porous media or reactive transport. For the discretization in time, the backward Euler method is often used due to its stability. This changes (1.1) into a sequence of nonlinear elliptic equations. For solving these, linear iterative schemes are required. A commonly used linearization technique is the Newton scheme (NS). Being quadratically convergent, it is widely used for solving nonlinear equations [1,2]. However, this quadratic convergence is featured under certain restrictions. In particular, degenerate problems do not fulfil these restrictions [3]. Another drawback of the NS is that it is only locally convergent, meaning that the initial guess for the iterations has to be close enough to the actual solution for the scheme to be convergent [1,3]. In many cases, this requires extremely small time step sizes limiting the applicability of the NS. For this reason, pre-conditioners, line-search methods and different parametrizations are often used to enhance the robustness of the NS [4,5]. ∗ Corresponding author at: Eindhoven University of Technology, Department of Mathematics and Computer Science, Netherlands. E-mail address: [email protected] (K. Mitra). https://doi.org/10.1016/j.camwa.2018.09.042 0898-1221/© 2018 Elsevier Ltd. All rights reserved.

Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

2

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

An alternative to the NS is a modified Picard scheme (PS), proposed in [6]. In [2,7] it is shown that this scheme is quite fast despite having linear convergence. Another linearly converging scheme is the Jäger–Kačur scheme (see [8–10]). A sufficient condition for convergence was derived in [3] for all the schemes mentioned above. When applied to the Euler implicit discretization of (1.1) one needs to take

τ < Cmrb hd ,

(1.2)

to guarantee the convergence of these schemes. Here τ is the time step size, h is the mesh size, d is the spatial dimension of the problem, mb ≥ 0 is the lower bound of b′ and C , r > 0 are constants that depend on the nonlinear functions. For d ≥ 2 this imposes a severe restriction on the time step size, which can increase the computation time to a great extent. Moreover, in the degenerate case mb = 0, either a regularized version of the function b has to be used [8,9,11] or the initial/boundary data has to be shifted [12] to ensure convergence. For porous media applications, where all the associated functions are nonlinear and the problem might become degenerate, stability is an important issue. Also, extremely large timescales are involved for such processes and so condition (1.2) cannot always be satisfied. To address this, a fixed point iteration scheme, termed L-scheme or simply LS in the context of this discussion, was proposed in [13–15]. The LS is linearly convergent but it has the interesting property of unconditional convergence, meaning that it converges to the time-discrete solution irrespective of the choice of initial guess. This is due to the fact that in the LS, the stabilization terms are estimated globally as opposed to the local estimations used in the NS, the PS or the Jäger–Kačur scheme. However, as shown in [7,16], this increases the convergence rate of the LS, making it slower when compared to the NS or the PS. This has motivated authors to either use the LS to provide initial guesses for the NS [7] or to apply it in a domain decomposition type approach that boosts the speed of the LS [16]. All of the schemes mentioned above are mostly designed to solve nonlinear elliptic problems and do not use the fact that these are the outcomes of the time discretization of a nonlinear parabolic problem. In this context, the solution from the previous time step can be used as initial guess in the iterative process. For solutions that have a good regularity in time, as being the case for parabolic problems, the changes in the solution from one time step to the next one are limited. However, the standard schemes do not use this fact and thus their implementation for parabolic and elliptic problems is more or less the same. The main idea in this work is to exploit the fact that the nonlinear elliptic problems are the result of the time discretization of (1.1). The proposed scheme is a combination of the PS and the LS and it uses local estimations to improve the convergence behaviour of the LS without affecting its stability. After introducing the problem in Section 2 in Section 3 we propose the iterative scheme for simple quasilinear equations, and analyse its convergence. It can be observed that, apart from being unconditionally convergent, the scheme has a convergence rate proportional to the time step size for sufficiently small time steps. This is unlike the LS, where a lower time step increases the convergence rate. Section 4 generalizes these ideas for (1.1). The scheme converges also in the degenerate case, however for small enough time step sizes. For the non-degenerate case, the convergence is linear and the rate is proportional to the square root of the time step size. Finally in Section 5 some numerical experiments are presented. These show that the scheme is more stable than the NS or the PS and converges at least as fast as the PS and sometimes comparable to the NS. 2. The general problem and linearization schemes Let Ω be a bounded domain in Rd which has a Lipschitz boundary ∂ Ω and define Q = Ω × [0, T ) for some T > 0. For the rest of our discussion (·, ·) and ∥·∥ will represent L2 (Ω ) inner product and norm. Any other norm will be presented as ∥·∥V with V being the corresponding space. The Sobolev space W k,p (Ω )∑ is the set of functions u defined on Ω such that Dk u ∈ Lp (Ω ) for the multi-index k, equipped with the norm ∥u∥W k,p (Ω ) = ( q≤k ∫Ω |Dq u|p )1/p for 1 ≤ p < ∞ [17]. Further, H k (Ω ) = W k,2 (Ω ) and H0k (Ω ) represents the set of elements of H k (Ω ) which have 0 trace at the boundary ∂ Ω [17]. The Ho¨ lder space C ℓ,δ (Ω ) refers to the space of δ -Ho¨ lder continuous functions up to the ℓth space derivative for the metric ⃗) = |⃗x − y⃗|, ⃗x, y⃗ ∈ Rd . The associated norms of these spaces are defined in detail in [18]. With these basic definitions dist(⃗ x, y stated, we introduce the problem. 2.1. Time-discrete formulation and properties of functions We consider (1.1) in the space–time domain Q . For simplicity homogeneous Dirichlet boundary condition is assumed on

∂ Ω × [0, T ]. Further, let u0 (·) ∈ H01 (Ω ) be the initial condition. We refer to [19] and [20] for results on the existence and uniqueness of weak solutions to (1.1) for this case. For the time discretization of (1.1) we let τ = T /N, N ∈ N, be the time step size and n ∈ {1, . . . , N } represent the time step. Define tn = nτ . The Euler implicit discretization of (1.1) leads to the sequence of elliptic problems (n ∈ {1, . . . , N })

{ (P 1)

b(un ) − b(un−1 ) + τ ∇ · F⃗ (⃗ x, un ) = τ ∇ · [D(⃗ x, un )∇ un ] + τ f (⃗ x, tn , un ) in Ω ;

(a)

un = 0 on ∂ Ω .

(b)

(2.1)

Below we consider weak solutions to the time discrete problem (P 1), defined as: Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

3

Definition 0.1. Let n ∈ {1, . . . , N } and un−1 ∈ L2 (Ω ) be given. A weak solution to Problem (P 1) is a function un ∈ H01 (Ω ) satisfying for any test function φ ∈ H01 (Ω ) (b(un ) − b(un−1 ), φ ) + τ (D(⃗ x, un )∇ un , ∇φ ) = τ (F⃗ (⃗ x, un ), ∇φ ) + τ (f (⃗ x, tn , un ), φ ).

(2.2)

The existence of un ∈ H01 (Ω ) for a given un−1 ∈ L2 (Ω ) is shown in [14]. Note that the sequence {un }Nn=1 can be used to approximate the solution to the original parabolic problem (1.1) (the Rothe method, [21]). Below we assume the following: (P. 1) b : R → R is a C 2 (R) function such that b′ ≥ mb and |b′ |, |b′′ | ≤ Mb for some mb , Mb ≥ 0. (P. 2) f : Q × R → R is twice differentiable with respect to u and satisfies T ∂u f (⃗ x, t , u) ≤ b′ (u) and |f |, |∂u f |, |∂uu f | ≤ Mf (Mf ≥ 0) a.e. in (⃗ x, t) ∈ Q and u ∈ R. (P. 3) D : Ω × R → R+ is a C 2 (Ω × R) function. Moreover, there exists Dm , DM > 0 such that Dm ≤ D and q D, |∂u D|, |∂xj D|, |∂u ∂xj D| ≤ DM for j ∈ {1, . . . , d}, q ∈ {1, 2}. q

(P. 4) F⃗ : Ω × R → Rd , with xj -component denoted by Fj , admits partial derivatives that satisfy |∂u Fj |, |∂xj Fj |, |∂u ∂xj Fj | ≤ MF for j ∈ {1, . . . , d}, q ∈ {1, 2} and MF ≥ 0. (P. 5) u0 ∈ H01 (Ω ). In some cases we also use u0 ∈ W 2,∞ (Ω ).

In (P. 2) instead of using ∂u f ≤ 0 we have used the less restrictive condition b′ ≥ T ∂u f . At this stage we define the following function which will be used later, z : Ω × R → R,

z(⃗ x, u) = b(u) − τ f (⃗ x, tn , u).

(2.3)

Observe that actually z depends on the time step size τ and n as well. However the focus here is on constructing an iterative 2 scheme for a fixed time step size ( τ at )a fixed (time tn). Hence, we disregard this dependence. From (P. 1) and (P. 2), z ∈ C and satisfies ∂u z = b′ − τ ∂u f ≥ b′ 1 − Tτ ≥ mb 1 − Tτ ≥ 0. So by defining m = mb (1 − Tτ ) we get the inequality

∂u z ≥ m ≥ 0 a.e. for ⃗x ∈ Ω and u ∈ R.

(2.4)

Remark 2.1 (On the properties of the functions). The situation mb = 0 in (P. 1) gives rise to degeneracy, which will be discussed in detail later. The boundedness of b′′ is assumed in our analysis. This is a more severe condition compared to the Lipschitz continuity assumed in [7,13–15,22] and the Ho¨ lder continuity assumed in [23]. However, it is true for many porous media flow models, where the form of b is as in (5.6). D can also be a matrix D˜ , the only constraint being that the condition (P. 3) is satisfied for all the components D˜ ij and that D˜ is positive definite. The bounds assumed in (P. 4) guarantee that F⃗ , which corresponds to the ‘flux’ in physics, is bounded and varies smoothly within Ω . Remark 2.2 (Boundary conditions). For the ease of presentation a zero Dirichlet boundary condition is assumed at the boundary, but the results can be proved for more general boundary conditions, including non-homogeneous Dirichlet, Neumann, Robin or mixed type ones. 2.2. Standard linearization techniques For resolving the nonlinear terms in problem (P 1), Newton scheme (NS) uses iterations where the values of the nonlinear terms in the current iteration are approximated by Taylor series expansion. For i ∈ N let uin stand for the ith iterate at the nth time step and let δ uin = uin − uin−1 for i > 0. If φni is one of the functions b, f , z , D, Fj for u = uin and t = tn , then one takes φni ≈ φni−1 + ∂u φni−1 δ uin , ∂u φni being the partial derivative of φ with respect to u at u = uin . With this substitution and given uin−1 , NS resumes to find δ uin that solves

[ ] ⎧ ′ i−1 i−1 i i−1 i i−1 i−1 i−1 i ⃗ ⎪ (b (u ) − τ ∂ f ) δ u − τ ∇ · D ∇δ u + ( ∂ D ∇ u − ∂ F ) δ u ⎪ u u u n n n n n n n n n ⎪ ⎨ [ ] n,i (PN ) = −(b(ui−1 ) − b(un−1 )) + τ ∇ · (Di−1 ∇ ui−1 ) − ∇ · F⃗ i−1 + f i−1 in Ω , (a) n n n n n ⎪ ⎪ ⎪ ⎩ i δ un = 0 on ∂ Ω . (b)

(2.5)

n ,i

For Problem (PN ), one of the natural initial guesses is u0n = un−1 . The modified Picard scheme (PS) can be interpreted as a simplified version of the NS [6]. Here the Taylor expansion is used only for the function z = b − τ f . For i ∈ N and with given uni−1 one seeks uin solving

[ i−1 i ] ⎧ i−1 i ⎪ ⎨∂u zn un − τ ∇ · Dn ∇ un n,i (PP ) = ∂u zni−1 uin−1 − (b(uin−1 ) − b(un−1 )) + τ [−∇ · F⃗ni−1 + fni−1 ] in Ω , (a) ⎪ ⎩ i un = 0 on ∂ Ω . (b)

(2.6)

Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

4

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

The computational cost per iteration of the PS is lower than that for the NS, because the number of gradients to be calculated is less. However, the scheme is linearly convergent and also only conditionally stable. A convergence proof can be found in [3]. The PS takes advantage of the fact that the function z is increasing with respect to u. This insight is taken a step further in the L-scheme (LS). For a L ≥ maxu∈R,(⃗x,t)∈Q {∂u z } = maxu∈R,(⃗x,t)∈Q {|b′ (u) − τ ∂u f (⃗ x, t , u)|} and a given uin−1 the scheme reads (see [13–15])

[ i−1 i ] ⎧ i ⎪ ⎨Lun − τ ∇ · Dn ∇ un n ,i (PL ) = Luin−1 − (b(uin−1 ) − b(un−1 )) + τ [−∇ · F⃗ni−1 + fni−1 ] in Ω , (a) ⎪ ⎩ i un = 0 on ∂ Ω . (b)

(2.7)

The difference between the LS and the PS is that, instead of using the factor ∂u zni−1 , the former uses an upper bound of it, L. This can affect the convergence order [7,16], but transforms the scheme to a H 1 -contraction which converges irrespective of the spatial discretization or the initial guess. Convergence can also be achieved for smaller values of L, e.g. for L = 21 maxu∈R,(⃗x,t)∈Q {∂u z } convergence is guaranteed (see [7,16]) although the contraction property cannot be ensured in this case. In the next section we show how to improve the speed of the LS while preserving this stability. n,i n,i n ,i Observe that each of the problems, (PN ), (PP ) and (PL ), has unique solution in H01 (Ω ) as they are linear and coercive. i Also, for the schemes mentioned above, if un → un in strong sense in H01 (Ω ) then un is a weak solution of (P1 ) in the sense of Definition 0.1. For the convergence analysis we let ein denote the error of the iterate at the nth time step, ein = uin − un .

(2.8)

The convergence rate (linear) of the schemes is defined as

 i e  1 n H (Ω ) α = lim sup  i−1  . en  1 i→∞

(2.9)

H (Ω )

Notice that α < 1 implies convergence. Also worth pointing out is that the H 1 -norm has been used to define the convergence rate in (2.9) but the convergence rate will roughly be the same if defined with respect to the L2 norm [7,13]. Hence, in Section 5 the L2 -norm has been used to estimate α . Remark 2.3. The convergence rate of the LS has the form α = (L − m)/(L + C τ ) < 1, for some C > 0 and L > m (see [7,13,14,22]). However, for large L, or small τ or m, α approaches 1. This leads to extremely slow convergence of the L-scheme. This issue has been reported for a variety of problems in literature [7,16,24]. The proofs given below make use of the following notations. [·]+ = max(·, 0) and [·]− = min(·, 0) are the positive and negative cut functions respectively. For any a, b ∈ R the interval I (a, b) is defined as I (a, b) = {x : min(a, b) < x < max(a, b)}.

We will further use Young’s inequality: for a, b ∈ R and ρ > 0 one has ab ≤

1 2ρ

a2 +

ρ 2

b2 .

(2.10)

Finally, CΩ is the constant appearing in the Poincare´ inequality. 3. The modified L-scheme In what follows we discuss a modified form of the L-scheme that preserves its stability property while having an improved convergence rate. The idea is to replace the constant L with a function Lin : Ω → R+ defined at each iteration. This leads to the Modified L-scheme (MS). Let n ∈ {1, . . . , N } and i ∈ N be fixed and assume that un−1 , uin−1 ∈ H01 (Ω ) are given. Find uin ∈ H01 (Ω ) satisfying (Lin uin , φ ) + τ (D(⃗ x, uin−1 )∇ uin , ∇φ )

= (Lin uin−1 , φ ) − (b(uin−1 ) − b(un−1 ), φ ) + τ (F⃗ (⃗x, uin−1 ), ∇φ ) + τ (f (⃗x, tn , uni−1 ), φ ), for all φ ∈

H01 (Ω ),

Lin (x)

where

⃗ = max{[b

′

Lin

(3.1)

: Ω → R is defined as ⃗ − τ ∂u f (⃗x, tn , uin−1 (⃗x)) + Mτ ], 2Mτ }.

(uin−1 (x))

(3.2)

u0n

Here = un−1 , and M is a positive constant that will be specified later. Observe that M = 0 corresponds to the PS, and disregarding the b′ − τ ∂u f term in the definition of Lin leads to the LS. In this sense the scheme is a combination of the PS and the LS. Below we show that it inherits the qualities of both schemes. The convergence results and estimates are first obtained for a simple version of (P 1) where ∂u D ≡ 0 and ∂u Fj ≡ 0 and then for the general problem. Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

5

3.1. Quasilinear equations with linear diffusivity and flux To analyse the scheme we first look at the simpler quasilinear parabolic equation,

∂t b(u) + ∇ · F⃗ (⃗x) = ∇ · (D(⃗x)∇ u) + f

(3.3)

which corresponds to ∂u D = 0 and ∂u Fj = 0 in (1.1). Problems where ∂u D ̸ = 0 can also be reduced to this ∫case if D is u separable in the variables, i.e. D(⃗ x, u) = D1 (⃗ x)D2 (u). By (P. 3), D2 > 0 and by using the Kirchhoff transform U = D2 (ϱ)dϱ (see [13,19,25]) one obtains ∇ · [D(⃗ x, u)∇ u] = ∇ · (D1 (⃗ x)∇ U). Using U as the primary unknown one gets an equation similar to (3.3). Well-known examples, such as the porous medium equation or the Richards equation, can be reduced to this form. Recalling the homogeneous boundary conditions, a weak solution to the time discrete version of (3.3) satisfies (b(un ) − b(un−1 ), φ ) + τ (D∇ un , ∇φ ) = τ (fn , φ ) + τ (F⃗ , ∇φ ),

(3.4)

for all φ ∈ H01 (Ω ). For the subsequent analysis we assume the following: (A. 1) There exists a Λ > 0 such that ∥un − un−1 ∥L∞ (Ω ) ≤ Λτ for all n ∈ N. Note that since un is, in fact, the time discrete approximation of the solution to the parabolic problem (3.3), Assumption (A. 1) is similar to saying that ∂t u ∈ L∞ (Q ). Sufficient conditions for boundedness of ∂t u in the L∞ (Ω )-norm can be found in [18, Chapter 5]. Below we present a result that shows the validity of Assumption (A. 1), in this sense. Proposition 3.1. For a fixed n ∈ {1, . . . , N } let un solve (3.4) and assume that m > 0 and ∇ · (D∇ un−1 ) ∈ L∞ (Ω ). Then ∇ · (D∇ un ) ∈ L∞ (Ω ) and a Λ ≥ 0, independent of τ , exists such that ∥un − un−1 ∥L∞ (Ω ) ≤ Λτ . The proof is given in Appendix. Observe that Proposition 3.1 is valid if u0 ∈ W 2,∞ (Ω ), which extends then to all time steps. Based on Assumption (A. 1) we have Lemma 3.1. Assume (A. 1) and that Lin satisfies a.e. in Ω Lin ≥ sup{|b′ (ζ ) − τ ∂u f (⃗ x, tn , ζ )| : ζ ∈ I (uin−1 , un )}. Then ∥uin − un ∥L∞ (Ω ) ≤ Λτ for all i ∈ N, uin being the solution of (3.1). Before giving the proof we observe that the result is quite general with respect to the choice of the function Lin . However, it will be used for Lin either constant or as in (3.2), which corresponds to the LS and the MS. Proof. The proof is by induction. Assumption (A. 1) guarantees that the assertion holds for i = 0. Let the assertion hold for uin−1 , i.e. ∥uin−1 − un ∥L∞ (Ω ) = ∥ein−1 ∥L∞ (Ω ) ≤ Λτ . We prove that this implies ∥uin − un ∥L∞ (Ω ) = ∥ein ∥L∞ (Ω ) ≤ Λτ . As ∂u D, ∂u Fj = 0, subtracting (3.4) from (3.1), a ζ ∈ I (uin−1 , un ) exists such that (Lin ein , φ ) + τ (D∇ ein , ∇φ ) = (Lin ein−1 , φ ) − ((b(uin−1 ) − b(un )), φ ) + τ (fni−1 − fn , φ )

= ((Lin − ∂u z(ζ ))ein−1 , φ ).

(3.5)

Here we used the definition of z from (2.3). With φ = [

− Λτ ]+ one gets  2 − Λτ ]+ ) + τ Dm ∇[ein − Λτ ]+  ein

(Lin

≤

− Λτ ], [

ein (Lin ein

[

∫ ≤ Ω

ein

,[

− Λτ ]+ ) + Λτ

(Lin

,[

− Λτ ]+ ) + τ (D∇ , ∇[

ein ein

− Λτ ]+ ) = ((Lin − ∂u z(ζ ))eni−1 , [ein − Λτ ]+ ) ∫ |Lin − ∂u z(ζ )||ein−1 |[ein − Λτ ]+ ≤ (Lin − m)Λτ [ein − Λτ ]+ ≤ Λτ (Lin , [ein − Λτ ]+ ). ein

ein

Ω

In the last estimates we used the inequalities |ein−1 | ≤ Λτ a.e. and 0 ≤ (Lin − ∂u z(ζ )) ≤ Lin − m for ζ ∈ I (uni−1 , un ). Cancelling the common terms in both sides gives

2

(Lin , [ein − Λτ ]2+ ) + τ Dm ∇[ein − Λτ ]+  ≤ 0,



which implies that ein < Λτ a.e. in Ω . Similarly taking φ = [ein + Λτ ]− one gets that ein > −Λτ a.e. which concludes the proof. □ For the LS, the condition stated in Lemma 3.1 is satisfied as Lin = L ≥ sup{|∂u z(⃗ x, ζ )| : ζ ∈ R}. However this leads to an overestimation of L at most points. Below we show that this estimation is improved significantly for the MS, resulting in much better convergence rates. Theorem 3.1. Under Assumptions (A. 1) and (P. 1)–(P. 5) and for M0 = Λ maxu∈R {|b′′ | + τ |∂uu f |} ≥ 0, the MS for Eq. (3.4) converges linearly in H01 (Ω ) for all M ≥ M0 and τ > 0. More precisely, for M ≥ M0 and τ > 0 the limit un = limi→∞ uin exists and is a solution to (3.4), whereas for the convergence rate α in (2.9) one has α < 1. Moreover if m > 0 (the non-degenerate case) then α = O(τ ) for τ small enough. Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

6

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

Proof. By (P. 1) and (P. 2), M0 is well defined. Observe that, for any ζ ∈ I (uin−1 , un ) there exists a ζ1 ∈ I (uin−1 , ζ ) such that

|(b′ (uin−1 ) − τ ∂u f (⃗x, tn , uin−1 )) − (b′ (ζ ) − τ ∂u f (⃗x, tn , ζ ))| = |∂u z(uin−1 ) − ∂u z(ζ )| = |∂uu z(ζ1 )||(uin−1 − ζ )| ≤ max{|b′′ | + τ |∂uu f |}Λτ = M0 τ . u∈R

This implies that if M ≥ M0 and Lin = ∂u z(uin−1 ) + Mτ , then Lin − ∂u z(ζ ) ≥ 0. Moreover, if Lin = 2Mτ then ∂u z(uin−1 ) ≤ Mτ , which means that ∂u z(ζ ) ≤ ∂u z(uin−1 ) + M0 τ ≤ 2Mτ , giving Lin − ∂u z(ζ ) ≥ 0. Hence, for M ≥ M0 one has x, tn , ζ ) ≥ 0, Lin − b′ (ζ ) + τ ∂u f (⃗ which by Lemma 3.1 implies that ∥ein ∥L∞ (Ω ) < Λτ for all i ∈ N. Using similar arguments, if Lin = ∂u zni−1 + Mτ and M ≥ M0 then one gets Lin − b′ (ζ ) + τ ∂u f (⃗ x, tn , ζ ) = ∂uu z(ζ1 )(uin−1 − ζ ) + Mτ ≤ Λτ |∂uu z(ζ1 )| + Mτ ≤ 2Mτ . If Lin = 2Mτ then Lin − ∂u z(ζ ) ≤ Lin − m ≤ 2Mτ . Combining both gives x, tn , ζ ) ≤ 2Mτ , 0 ≤ Lin − b′ (ζ ) + τ ∂u f (⃗

(3.6)

or, more strongly, 0 ≤ (M − M0 )τ ≤ Lin − b′ (ζ ) + τ ∂u f (⃗ x, tn , ζ ) ≤ 2Mτ .

(3.7)

These two inequalities will be used repeatedly in the proofs that follow. First we prove the convergence of the scheme in H01 (Ω ). Taking φ = ein in (3.5) yields

 2

2

2Mτ ein  + τ Dm ∇ ein  ≤ (Lin ein , ein ) + τ (D∇ ein , ∇ ein )



 2  2 = ((Lin − ∂u z(ζ ))ein−1 , ein ) ≤ 2Mτ (|ein−1 |, |ein |) ≤ Mτ ein−1  + Mτ ein  .

(3.8)

Cancelling common terms on both sides and applying the Poincare´ inequality one gets

(

Dm

M+

2

) CΩ

 i 2 Dm  i 2   e  + ∇ e  ≤ Mei−1 2 , n n n 2

which gives

 i 2 e  + n

(

As ∥u∥2 +

Dm

(2M + CΩ Dm )

Dm (2M+CΩ Dm )

∥∇ u∥2

 i 2 ∇ e  ≤ n

)1/2

(

2M

 i−1 2 e  + n

(2M + CΩ Dm )

Dm

(2M + CΩ Dm )

)  i−1 2 ∇ e  . n

is equivalent to the H 1 (Ω )-norm it follows that {uin }i∈N converges in H 1 (Ω ). The conver-

gence rate in this equivalent norm is

α=

√

2M/(2M + CΩ Dm ) < 1.

(3.9)

Observe that the convergence does not depend on the spatial discretization and also holds in the degenerate case when ∂u z may vanish. In the non-degenerate case, when ∂u z ≥ m > 0, substituting φ = ein in (3.5) one gets

 2

2

2

(m + Mτ )ein  + τ Dm ∇ ein  ≤ (Lin ein , ein ) + Dm τ ∇ ein 





= ((L − ∂u z(ζ ))ein−1 , ein ) ≤ 2Mτ (|ein−1 |, |ein |) ≤

  2M2 τ 2  i−1 2 e  + m ei 2 . n m 2 n

Cancelling the common term in both sides we obtain

 i 2 e  + n

2τ Dm (m + 2Mτ )

 i 2 ∇ e  ≤ n

4M2 τ 2 m(m + 2Mτ )

(  i−1 2 e  + n

2τ Dm (m + 2Mτ )

)  i−1 2 ∇ e  . n

(3.10)

With τ small enough, one obtains as before, the convergence of uin in H 1 (Ω ). The convergence rate is

α= √

2Mτ

m(m + 2Mτ )

which is less than 1 for τ < with respect to a different

( α = min

√

<

2Mτ m

.

(3.11)

τ |

 2  2 |, |ein |) < Mτ (ein−1  +ein  ) to prove contraction

m . One can also use the inequality 2M ( ein−1 2M Mτ H 1 (Ω )-norm with . Hence, the actual m

2Mτ

m(m + 2Mτ )

α=

√ ,

Mτ m

√ ,

√

2M 2M + CΩ Dm

convergence rate is

) . □

(3.12)

Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

7

We conclude this section with a result that shows that for the non-degenerate case even the L∞ (Ω ) errors, i.e. ∥ein ∥L∞ (Ω ) , decrease linearly for τ sufficiently small. Proposition 3.2. For a fixed n ∈ {1, . . . , N }, let {uin }i∈N be the sequence of functions resulting from the modified L-scheme for (3.4). Assume (A. 1) and (P. 1)–(P. 5). If m > 0, then for small enough τ and M ≥ M0 ,

 i    u − un  ∞ < β ui−1 − un  ∞ , n n L (Ω ) L (Ω ) for i ∈ N and a constant β ∈ (0, 1). Moreover, β = O(τ ). The proof is given in Appendix. Remark 3.1 (Robustness of the MS). Theorem 3.1 shows that the MS converges irrespective of the spatial discretization and the time step. It converges in the degenerate case m = 0 too. Unlike the LS, the convergence rate α is independent of Mb , and scales with τ for small τ . This robustness is extremely helpful for computations, as it is difficult to satisfy condition (1.2), which guarantees the convergence of schemes like NS, PS and Jäger–Kačur, for higher dimensional computations, i.e. d > 1. Moreover, the α = O(τ ) property makes the scheme faster as the time step size is made smaller. 4. General nonlinear diffusion equation For the general problem, when F⃗ and D are functions of u, convergence can be shown for both the LS and the MS when τ is small enough. For both schemes,√ the convergence rate does not depend on the spatial discretization. Moreover, for the MS the convergence rate scales with τ for small τ values. For proving the main theorem of this section, we assume (A. 2) ∥∇ un ∥L∞ (Ω ) ≤ Λ1 for all n ∈ {1, . . . , N } and some Λ1 > 0. Due to (P. 3), this is equivalent to assuming that the flux is bounded. Similar to Section 3, the MS works if for all i ∈ N one has

∥uin − un ∥L∞ (Ω ) ≤ Λτ .

(4.1)

Though non-trivial, this condition is fulfilled under certain assumptions, as follows from Lemma 4.1. Let n ∈ {1, . . . , N } be fixed and Ω be a C 2 domain. Assume (A. 1), m > 0 and let Λ2 > 0 be such that ∥un ∥W 2,2q (Ω ) ≤ Λ2 for some q > 2d , q ∈ N. Further, assume that a Λ3 > 0 exists such that

∥un − un−1 ∥W 1,2q (Ω ) ≤ Λ3 τ .

(4.2)

Let {uin }i∈N be the sequence generated by the MS. Then there exists a τ˜ > 0 such that for all τ ≤ τ˜ and i ∈ N, ∥uin − un ∥L∞ (Ω ) ≤ Λτ , ∥uin − un ∥W 1,2q (Ω ) ≤ Λ3 τ and uin ∈ W 2,2q (Ω ). Proof. Similar to Lemma 3.1, we give a proof by induction. Assume that ∥ekn ∥L∞ (Ω ) ≤ Λτ , ∥∇ ekn ∥L2q (Ω ) ≤ Λ3 τ and ukn ∈ W 2,2q (Ω ) for k < i. This is true for k = 0 because of Assumptions (A. 1) and (4.2). We show that this implies ∥ein ∥L∞ (Ω ) ≤ Λτ , ∥∇ ein ∥L2q (Ω ) ≤ Λ3 τ and uin ∈ W 2,2q (Ω ) for small τ . More specifically, we show that there exists a C > 0 such that

∥ein ∥ ≤ C τ 2 ,

∥ein ∥L∞ (Ω ) ≤ C τ

1 1+ 2q

,

1

∥ein ∥W 1,p (Ω ) ≤ C τ 1+ p ,

∥ein ∥W 2,2q (Ω ) ≤ C τ ,

(4.3)

for all p ≥ 2 and τ > 0. This proves the lemma for τ small enough. First, observe that inequality (3.6) holds for Lin defined in (3.2). Moreover, the assumption un ∈ W 2,2q (Ω ) implies (A. 2), i.e. there exists a Λ1 > 0 such that ∥∇ un ∥L∞ (Ω ) ≤ Λ1 . This is a direct consequence of Morrey’s inequality [17, Chapter 5] and 2q > d. By the same argument uin−1 ∈ W 1,∞ (Ω ) ∈ H 1 (Ω ). Hence, by Theorem 5 of [17, Chapter 6] it follows that uin is a classical solution to Lin uin − τ ∇ · (Dni−1 ∇ uin ) = Lin uin−1 − (bin−1 − bn−1 ) + τ ∇ · F⃗ni−1 + τ fni−1 .

(4.4)

H01

uin

as well. In particular, is essentially bounded and lies in Subtracting (2.1) from (4.4) and rearranging the terms leads to Lin ein − τ ∇ · (Dni−1 ∇ ein )

= Lin ein−1 − ∂u z(ζ )ein−1 + τ ∇ · ((Dni−1 − Dn )∇ un ) − τ ∇ · (F⃗ni−1 − F⃗n ), where ζ ∈

I (uin−1

(4.5)

, un ). Denoting the terms on the right by I1 , I2 , I3 we have

|I1 | = |(Lin − ∂u z(ζ ))ein−1 | ≤ 2ΛMτ 2 ; Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

8

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

|I2 | = τ |(Dni−1 − Dn )∆un + (∂u Dni−1 ∇ uin−1 − ∂u Dn ∇ un ) · ∇ un +

)

–

d ∑

(∂xj Dni−1 − ∂xj Dn )∂xj un |

j=1

≤ ΛDM τ |∆un | + τ |∇ un ∥∂ 2

i−1 u Dn

∇

ein−1

+ (∂

i−1 u Dn

− ∂u Dn )∇ un | + τ |∇ un ∥

∑

(∂xj Dni−1 − ∂xj Dn )|

≤ ΛDM τ 2 |∆un | + τ DM Λ1 |∇ ein−1 | + DM ΛΛ21 τ 2 + dΛ1 DM Λτ 2 ; |I3 | = τ |∇ · (F⃗ni−1 − F⃗n )| = τ |(∂u F⃗ni−1 · ∇ uin−1 − ∂u F⃗n · ∇ un ) +

d ∑

(∂xj Fji,−n 1 − ∂xj Fj,n )|

j=1

≤ τ |∂

i−1 u Fn

≤ τ MF |∇

·∇

ein−1

ein−1

| + τ |(∂u F⃗ni−1 − ∂u F⃗n ) · ∇ un | + τ |

∑

(∂

i−1 xj Fj,n

− ∂xj Fj,n )|

| + ΛMF Λ1 τ + dΛMF τ . 2

2

Define Si := I1 + I2 + I3 . As ∥ein−1 ∥W 1,2q (Ω ) ≤ Λ3 τ and ∥∆un ∥L2q (Ω ) ≤ Λ2 it follows that a C1 > 0 exists such that

∥Si ∥L2q (Ω ) ≤ C1 τ 2 .

(4.6)

Now we test (4.5) with the test function φ = (ein )2q−1 ∈ H01 (Ω ). This gives 2q

m∥ein ∥L2q (Ω ) + τ Dm (2q − 1)

≤

1 2qm2q−1

∫

|ein |

Ω

2(q−1)

(2q − 1)m

∥Si ∥2q + L2q (Ω )

2q

2

|∇ ein | ≤ (|Si |, |ein |

2q−1

)

∥ein ∥2q . L2q (Ω )

The last inequality follows from repeated application of Young’s inequality. This implies

∥ein ∥L2q (Ω ) ≤

1 m

∥Si ∥L2q (Ω ) ≤

C1 m

τ 2.

Rewriting (4.5) as

( ) ∇ · Dni−1 ∇ ein = Lin

(

ein

) −

τ

1

τ

Si ,

(4.7)

we see that the L2q (Ω )-norm of the right hand side is bounded by some constant times τ . As |∇ Dni−1 | ≤ |∂u Dni−1 ∇ uni−1 | + ∑ i−1 i−1 i−1 2q j |∂xj Dn | ≤ DM (|∇ en | + Λ1 + d), |∇ Dn | ∈ L (Ω ). Hence, we apply Theorem 15.1 of [26, Chapter 2] to get that d i 2,2q 1,γ en ∈ W (Ω ) ∩ C (Ω ) for γ = 1 − 2q . Moreover, there exists a C > 0 such that

∥ein ∥W 2,2q (Ω ) ≤ C τ and ∥∇ ein ∥L∞ (Ω ) ≤ ∥ein ∥C1,γ (Ω ) ≤ C τ .

(4.8)

This proves the last statement of (4.3). Next, we test (4.5) with φ = ein to get m∥ein ∥2 + τ Dm ∥∇ ein ∥2 ≤ (Si , ein ) ≤

1 2m

∥Si ∥2 +

m 2

∥ein ∥2 .

(4.9)

As Si ∈ L2q (Ω ) ⊆ L2 (Ω ) we get using (4.6) that for some constant C2 > 0

∥∇ ein ∥2 ≤ C2 τ 3 .

(4.10)

This implies that for p ≥ 2,

∥∇ ein ∥Lp (Ω ) =

(∫ Ω

|∇ ein |

p

)1/p

( )1/p ∫ 1 −2 i 2 ≤ ∥∇ ein ∥pL∞ |∇ e | ≤ (C p−2 C2 )1/p τ 1+ p . n (Ω )

(4.11)

Ω

Finally, as C 0,γ (Ω ) ⊆ W 1,2q (Ω ), see Morrey’s inequality (Theorem 5 of [17, Chapter 5]), we get that

∥ein ∥L∞ (Ω ) ≤ ∥ein ∥C0,γ (Ω ) ≤ C τ

1 1+ 2q

.

From this we get ∥ein ∥L∞ (Ω ) ≤ Λτ if τ ≤ (Λ/C )2q . Similarly one obtains ∥∇ ein ∥L2q (Ω ) ≤ Λ3 τ for small τ .

(4.12) □

Before giving the main theorem of this section we give a context that can lead to improved convergence behaviour of the MS. Specifically, we assume

|∂u D(⃗x, u)| +

d ∑

|∂u Fj (⃗x, u)| ≤ Υ ∂u z(⃗x, u) a.e. in ⃗x ∈ Ω and u ∈ R.

(4.13)

j=1

Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

9

Remark 4.1. The bound (4.13) is true, e.g. for the Richards equation (see (5.1), Section 5). In this case the diffusivity and the flux terms are D = k(b(u)) and F⃗ = k(b(u))eˆ g (eˆ g is a constant unit vector) with b : R → [0, 1] and k ∈ C 1 ([0, 1]) giving |∂u D|, |∂u Fj | ≤ sups∈[0,1] {k′ (s)}b′ (u). Also the quasilinear system discussed in Section 3 is just a special case when Υ = 0. Theorem 4.1. For a fixed n ∈ {1, . . . , N } let {uin }, i ∈ N be the sequence provided by the MS. Assume (4.1) holds for i ∈ N. If Assumptions (A. 1)–(A. 2) and (P. 1)–(P. 5) are satisfied then the following hold: (a) If inequality (4.13) is satisfied then there exist a M1 > 0 and τ ∗ > 0 independent of m such that if M ≥ M1 and τ ≤ τ ∗ then uin → un in the strong sense in H 1 (Ω ). (b) If m > √ 0 then there exists a τ¯ > 0 such that for τ < τ¯ and M ≥ M0 , uin converges linearly to un in H 1 (Ω ). Moreover α = O( τ ) for this case. Proof ((a)). We follow the line of arguments presented in [7] for proving this part. Subtracting (2.2) from (3.1) and taking

φ = ein gives

(Lin (ein − ein−1 ), ein ) + τ (Dni−1 ∇ uin − Dn ∇ un , ∇ ein )

= −(z(uin−1 ) − z(un ), ein ) + τ (F⃗ni−1 − F⃗n , ∇ ein ).

(4.14)

This can be rearranged into the form T1 + τ T2 + T3 = T4 + τ T5 + τ T6 where the terms Tj are estimated as: T1 := (Lin (ein − ein−1 ), ein ) =

∫ T2 :=

2

Ω

1

∫

2

1

2

Lin |ein | −

Ω

∫

2

2

Ω

Lin |eni−1 | +

1 2

∫

2

Ω

Lin |ein − eni−1 | ;

2



Dni−1 |∇ ein | ≥ Dm ∇ ein  ;

(z(uin−1 )

− z(un ),

ein−1 )

∫

1

2

|z(uin−1 ) − z(un )| ; ζ0 ) Ω ∂u z(∫ ∫ 1 1 1 2 2 i−1 i−1 i−1 i |z(un ) − z(un )| + Lin |eni−1 − ein | ; T4 := (z(un ) − z(un ), en − en ) ≤ i T3 :=

T5 := −

((Dni−1

− Dn )∇ un , ∇

=

ein )

≤

2 1

Dm

Ω

Ln

(Dni−1

∥

2

2

− Dn )∇ un ∥ +

Dm

4

Ω

∥∇ ein ∥2

 2 2    Λ21  ∂u D(ζ1 ) i−1   + Dm ∥∇ ei ∥2 ≤ (Υ Λ1 ) z(ui−1 ) − z(un )2 + Dm ∥∇ ei ∥2 ; ≤ (z(u ) − z(u )) n n n n n  Dm  ∂u z(ζ1 ) 4 Dm 4 T6 := (F⃗ni−1 − F⃗n , ∇ ein ) ≤

≤

1 Dm

Dm ∥F⃗ni−1 − F⃗n ∥2 + ∥∇ ein ∥2

4

d 1 ∑  ∂u Fj (ζ2 )

Dm

j=1

 2 2    i−1   + Dm ∥∇ ei ∥2 ≤ Υ z(ui−1 ) − z(un )2 + Dm ∥∇ ei ∥2 . (z(u ) − z(u )) n n n n n  ∂u z(ζ2 )  4 Dm 4

For the functions ζ1 , ζ2 : Ω → R one has ζ1,2 ∈ I (uni−1 , un ). For T3 , if z(uni−1 ) = z(un ) then the equality holds for any ζ0 ∈ R for which ∂u z(ζ0 ) ̸= 0. If z(uin−1 ) ̸= z(un ) then by the Mean Value Theorem there exists a ζ0 ∈ I (uin−1 , un ) such that ∂u z(ζ0 ) ̸= 0. Hence, there exists a ζ0 : Ω → R such that equality for T3 is satisfied and ∂u z(ζ0 ) ̸= 0 a.e. In the following, T3 will only induce an upper bound for τ which will depend on the lower bound of ∂ z(1ζ ) . Therefore we claim that the choice u 0 of ζ0 , discussed above, has no influence on the result. Putting everything together we get the following inequality

) ∫ (  i 2 2 1 2τ Υ 2 2 2   | | + τ Dm ∇ en + − i − (1 + Λ1 ) |z(uni−1 ) − z(un )| ∂ z( ζ ) L D m Ω Ω n ∫ ∫ ∫u 0 i i−1 2 i−1 i−1 2 i i−1 i−1 2 ≤ Ln |en | ≤ Ln |en | + |Ln − Ln ||en | .

∫

Lin

2 ein

Ω

Ω

(4.15)

Ω

This inequality is useful if all the terms on the left are positive. This is achieved if will add some restrictions on τ . To see this we define

M1 = max{4M0 , 2M0 + G},

2

∂u z(ζ0 )

−

1 Lin

−

2τ Υ 2 (1 Dm

+ Λ21 ) ≥ 0, which (4.16)

where M0 = Λ maxu∈R {|b | + τ |∂uu f |} is defined just as in the proof of Theorem 3.1 and the value of G will be clarified ′′

later. Observe that ( ∂2z − u

1 ) Lin

=

1

∂u z

+

Lin −∂u z Lin ∂u z

. From (3.7) we get Lin − ∂u z(ζ0 ) > (M − M0 )τ . Moreover, from (P. 1)–(P. 2),

∂u z ≤ Mb + TMf and Lin ≤ Mb + TMf + Mτ for small τ . This gives Lin − ∂u z Lin ∂u z

≥

(M − M0 )τ (Mb + TMf )(Mb + TMf + Mτ )

.

(4.17)

Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

10

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

To simplify the analysis we note that

Lin −∂u z Lin ∂u z

−

2τ Υ 2 (1 Dm

)

–

+ Λ21 ) ≥ 0 is a sufficient condition for the positivity of the last term

on the left hand side of (4.15) . Inequality (4.17) implies that this is satisfied when

( 1−

2Υ 2 Dm

(1 + Λ21 )(Mb + TMf )τ

) M ≥ M0 +

2Υ 2 Dm

(1 + Λ21 )(Mb + TMf )2 .

Hence by defining

τ∗ =

Dm

4Υ

2 (1

+

Λ21 )(Mb

and G =

+ TMf )

4Υ 2 Dm

(1 + Λ21 )(Mb + TMf )2 ,

we get that for all τ ≤ τ ∗ and M ≥ M1 , Lin − ∂u z

∂

Lin u z

−

2τ Υ 2 Dm

(1 + Λ21 ) ≥ 0.

Now consider τ ≤ τ ∗ and M ≥ M1 . Inequality (4.15) is restated as:

∫ Ω

Lin

2 ein

 2 | | + τ Dm ∇ ein  +

∫ Ω

2 ein−1

∂u z(ζ0 )|

∫

| ≤ Ω

Lin−1

2 eni−1

|

∫

2

| + Ω

|Lin − Lin−1 ||ein−1 | .

(4.18)

Observe that Ω can be split into two disjoint sets defined as Ω1 , Ω2 such that Ω1 = {⃗ x ∈ Ω : |∂u z(ζ0 )| < Ω2 = {⃗x ∈ Ω : |∂u z(ζ0 )| ≥ 21 Mτ }. If ⃗x ∈ Ω1 , then from M ≥ M1 ≥ 4M0 we get

1 M 2

τ } and

∂u z(uin−2 ) ≤ ∂u z(ζ0 ) + max {∂uu z }|ζ0 − un | + max {∂uu z }|un − uni−2 | I (uin−1 ,un )

≤ ∂u z(ζ0 ) + M0 τ + M0 τ ≤

1 2

I (un ,uin−2 )

Mτ +

1 2

Mτ

which implies that Lin−1 = max{∂u z(uin−2 ) + Mτ , 2Mτ } = 2Mτ . By the same logic, as |∂u z(ζ0 )| < This means that

∫

1 M 2

τ , one has Lin = 2Mτ .

2

Ω1

|Lin − Lin−1 ||ein−1 | = 0.

If ⃗ x ∈ Ω2 then ∂u z(ζ0 ) ≥ 21 Mτ ≥ 2M0 τ . There can be four cases. If ∂u z(uin−1 ) ≤ Mτ and ∂u z(uin−2 ) > Mτ then i |Ln − Lin−1 | = ∂u z(uin−2 ) − Mτ ≤ ∂u z(uin−2 ) − ∂u z(uin−1 ). Similarly if ∂u z(uin−1 ) > Mτ and ∂u z(uni−2 ) ≤ Mτ then |Lin − Lin−1 | ≤ ∂u z(uin−1 ) − ∂u z(uin−2 ). Lin = Lin−1 if both ∂u z(uin−2 ), ∂u z(uni−1 ) ≤ Mτ and |Lin − Lin−1 | = |∂u z(uni−1 ) − ∂u z(uni−2 )| if both ∂u z(uni−2 ), ∂u z(uin−1 ) > Mτ . Combining everything we can state that

|Lin − Lin−1 | ≤ |∂u z(uin−1 ) − ∂u z(uin−2 )| ≤

{|∂uu z |}|uin−1 − uin−2 | ≤ 2M0 τ .

max

I (uin−1 ,uin−2 )

From the above, one obtains

∫

∫ ∫ 2 2 2 |Lin − Lin−1 ||ein−1 | = |Lin − Lin−1 ||ein−1 | ≤ 2M0 τ |ein−1 | Ω Ω2 Ω2 ∫ ∫ 2 2 ≤ ∂u z(ζ0 )|ein−1 | ≤ ∂u z(ζ0 )|ein−1 | . Ω2

Ω

Using this result in (4.18) one gets

∫

2

Ω

2

Lin |ein | + τ Dm ∇ ein  ≤



∫

2

Ω

Lin−1 |ein−1 | .

(4.19)

Taking the sum of (4.19) from i = 1 to i = p gives

∫

2

Ω

Lpn |epn | + τ Dm

∫ p ∑    i 2 2 ∇ e  ≤ L0n |e0n | ≤ L0n L1 (Ω ) Λ2 τ 2 , n i=1

(4.20)

Ω

∑∞ 

2

with L0n = max{∂u z(⃗ x, un−1 )  + Mτ, 2Mτ }. Similar estimates are given in [8,11,12]. In other words, the series i=1 ∇ ein  is convergent implying that ∇ ein  → 0 as i → ∞. This concludes the proof of the first part of Theorem 4.1. (b). To prove (b) we rearrange (4.14) as T1′ + τ T2 = T3′ + τ T5 + τ T6 , where Tj , j ∈ {2, 5, 6} have been defined before, and ′ T1 , T3′ are T1′ :=

∫

2

Ω

 2

Lin |ein | ≥ (m + Mτ )ein  ;

 2

2

T3 := ((Lin − ∂u z(ζ ))ein−1 , ein ) ≤ Mτ ein  + Mτ ein−1  ; ′



Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

11

the last inequality following from (3.6) and Young’s inequality. Moreover, a different estimate for T5 , T6 can be obtained as T5 ≤ T6 ≤

1 Dm

1 Dm

∥(Dni−1 − Dn )∇ un ∥2 + ∥F⃗ni−1 − F⃗n ∥2 +

Dm

4

Dm

4

∥∇ ein ∥2 ≤ dMF2

∥∇ ein ∥2 ≤

Dm

2 Λ21 DM

Dm

∥ein−1 ∥2 +

∥ein−1 ∥2 +

Dm

4

Dm

4

∥∇ ein ∥2 ;

∥∇ ein ∥2 .

Adding in the estimates for all remaining terms, one gets from (4.14),

[ ] 2  i 2   i 2 τ (dMF2 + DM Λ22 )    ei−1 2 ,   m en + Dm ∇ en ≤ τ M + n 2

Dm

which rewrites as

 i 2 τ D m  i 2 e  + ∇ e  ≤ τ n n 2m

[

m

M+

2 (dMF2 + DM Λ21 )

](

)  i−1 2 τ Dm  i−1 2 e  + ∇ e  . n n 2m

Dm

Taking

τ¯ =

√

mDm 2 Λ21 ) MDm + (dMF2 + DM

and α =

τ m

[ M+

2 Λ21 ) (dMF2 + DM

Dm

] ,

(4.21)

one observes that the iteration converges in the equivalent H 1 (Ω )-norm ∥u∥H 1 (Ω ) =

√

convergence rate α = O( τ ) < 1 if τ < τ¯ . □

√ ∥ u∥ 2 +

τ Dm 2m

∥∇ u∥2 and has the

5. Numerical results For the numerical discussions we consider the Richards equation, which is widely used in groundwater modelling. In terms of the capillary pressure p, the non-dimensional Richards equation reads:

ˆ ]+f, ∂t S(p) = ∇ · [k(S(p))(∇ p − g)

(5.1)

where gˆ is the unit vector along the direction of gravity and f is the source term. Richards equation involves nonlinearities in all the terms. The saturation function S is increasing and the permeability function k takes non-negative values. Without entering into the details, as the specific forms will be given later, we mention that in general one has k(0) = 0 and S ′ (p) → 0 as p → −∞. Further S ′ (p) = 0 for all p ≥ 0. However, if the flow does not become completely unsaturated, meaning that S(p) → 0 (this being the case when the initial and boundary conditions are taken accordingly), Assumption (P. 3) will be satisfied. Also (4.13) is satisfied as discussed in Remark 4.1. The theoretical results presented in the previous sections do not depend upon the spatial discretization. Hence, for the numerical results, different methods like finite difference, finite element or finite volume, can be used to discretize (5.1) in space. Here we have used a two point flux approximation finite volume scheme [27] defined on two dimensional triangular unstructured grids. We take

Ω = (0, 1) × (0, 1) and T = 1,

(5.2)

and use FVCA8 benchmark meshes of different sizes. For the triangulation T the mesh size is h = sup{diam(T ) : T ∈ T }. With S and k defined as 1

{ S(p) =

1

(1−p) 3

1

if p < 0 if p ≥ 0,

k(S) = S 3 ,

(5.3)

the first numerical example is constructed in such a way that p˜ (x, y, t) = 1 − (1 + t 2 )(1 + x2 + y2 ),

(5.4)

is the exact solution (see also [16]). The source term is f (x, y, t) =

2(1 − y2 ) (1 + y2 )2

− √ 3

3 (1 +

2t t 2 )4 (1

− +

y2 )

2x (1 + t 2 )(1 + x2 + y2 )2

.

(5.5)

The boundary and initial conditions are as in Table 1. gˆ points along the positive x-axis. Fig. 1 (left) shows the comparison of the analytical solution p˜ with the numerical solution. The numerical solution is p−˜p obtained from the MS with M = 1. The maximum relative error ∥ p˜ ∥L∞ (Ω ) is 1.38% and the L2 (Ω ) error is 0.0116 which is of the order of the discretization errors, implying that the computational results are accurate. To ensure correctness, such kind of profile match has been conducted for all the results shown afterwards. Fig. 1 (right) shows how the errors decrease for the same computation for different schemes. Both L2 and L∞ errors of the MS decrease monotonically. The decrease of the L∞ error shows that Assumption (A. 1) is valid and the linear profile of the Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

12

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

Table 1 Assumed initial and boundary conditions. IC

t=0

p(x, y, 0) = p˜ (x, y, 0)

on Ω

BC

x=0: y=0:

p(0, y, t) = p˜ (0, y, t), ∂y p = 0,

x=1: y=1:

p(1, y, t) = p˜ (1, y, t), k(S(p))∂y p = k(S(p˜ (x, 1, t))∂y p˜ (x, 1, t).

Fig. 1. (left) Analytical solution vs. numerical solution for S, k and f given in (5.3), (5.5) and for initial and boundary condition stated in Table 1. The computational parameters are M = 1, τ = 0.001 and h = 0.02. (right) The error decay of the same computation at t = 0.5 for different schemes. Both L2 error and L∞ error of the MS have been plotted and L = 0.2 has been used for the LS.

L2 and L∞ errors shows that the convergence is indeed linear as pointed out in Theorem 4.1 and Proposition 3.2. Observe that the error plot of the PS almost coincides with the error plot of the MS. This is because τ = 0.001 is quite small, and for reasons explained afterwards this is true in general for small values of τ . The NS converges faster than the PS and the MS. The error of the LS decreases linearly but the speed of convergence is considerably less than the other schemes shown. The values of M and L used are the optimal values that give fastest convergence in their respective cases. This is explained in detail later. To illustrate the strength of the MS compared to the other schemes mentioned above, we increase the complexity of the problem by using relations used for real-life simulations of such processes. For this purpose we replace the expressions in (5.3) with the standard van Genuchten relationship: for ℓ > 1 and q = 1 − 1ℓ ,

{ S(p) =

1 (1+(−p)ℓ )q

1,

,

if p < 0 if p ≥ 0,

√ k(S) =

S(1 − (1 − S 1/q )q )2 .

(5.6)

In the computations, ℓ = 3 has been used throughput. Other conditions and definitions remain as in Table 1. The problem is degenerate as inside the domain there are regions where p > 0 or S ′ (p) = 0. For the numerical results M = 10 has been used for the MS and L = 0.4 has been used for the LS, unless specified otherwise. These values give optimal convergence for both the schemes. The choice is motivated later. First a mesh study is conducted where the time step is fixed and the mesh size h is varied from h = 0.1 to h = 0.02. The results are shown in Fig. 2 for two time step sizes τ = 0.01 and τ = 0.001 for a fixed time t = 0.5. It follows from Fig. 2 that the MS, the LS and the PS show linear convergence. The MS is faster than the PS for all mesh sizes for τ = 0.01 (Fig. 2 (left)). The convergence rate of the LS is higher than that for other schemes and in consequence, LS converges slower. For h = 0.05 and h = 0.02 no further decrease in error is visible after a few iterations in case of the NS. The convergence rates for all the schemes vary with the mesh size. The dependence of the convergence rates of the MS and the LS on mesh size stems from the fact that M or L chosen in these cases is not the M0 or Lipschitz constant value (which in this case is L ≈ 0.65) but the values that give the optimal convergence properties. This is also seen in Fig. 6. In fact, taking L = 1 for the LS results in all the error plots for different mesh sizes being on top of each other. This error characteristic is presented by the curve labelled L = 1. Also, it appears that for a fixed time step size, decreasing h makes the convergence slower. This will be explained in detail later. Fig. 2 (right) shows the results for τ √ = 0.001. The convergence is much faster than the τ = 0.01 case for the MS as it was shown in Theorem 4.1 that α = O( τ ). The difference between the convergence rates of the PS and of the MS is very small. This is because τ ≪ 1 and M = 10 so that Mτ ≪ S ′ (p) ∼ O(1) in most of the domain Ω . So difference between the MS and the PS is small. Fig. 2 (right) also shows that convergence rate of the LS is stable with respect to the mesh size and is not impacted greatly by the change in time step size. This is probably due to Condition 16 of [7] being satisfied for this value of τ = 0.001. However, the convergence is slower when compared to the other schemes. As for the NS, for h = 0.1 it is marginally faster than the PS and the MS. But for finer meshes, NS becomes slower than the PS or the MS. Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

13

Fig. 2. Mesh study of the NS, the PS, the LS and the MS at t = 0.5. (left) Results for τ = 0.01; (right) results for τ = 0.001. The mesh sizes used are h = 0.1, 0.05, 0.02. The L = 1 curve on the left figure corresponds to the error of the LS for L = 1 and h = 0.05. In all other computations M = 10 and L = 0.4 and the PS corresponds to M = 0.

Fig. 3. Numerical study of different linearization techniques for τ = 0.1 at t = 0.5. The mesh sizes used are h = 0.1, 0.05, 0.02. In all the computations M = 10 and L = 1 and the PS corresponds to M = 0.

The quadratic convergence of the NS is not observed from Fig. 2. This might be due to the small time step sizes required for the NS to show quadratic convergence in degenerate cases (see [3] and Table 2). It could also be attributed to the errors of the linear solver as it is known from literature [4,7,16] that the stiffness matrices for the NS are relatively ill-conditioned. Indeed, the GMRES solver gives higher residual errors for the NS compared to other methods. The errors introduced in computing the Jacobian could possibly be another reason that causes this deviation. One of the advantages of the MS is its stability and for this reason we must look at larger time step sizes. Fig. 3 shows the results for τ = 0.1. The PS converges, albeit much slower than the MS, for h = 0.1. For h = 0.05 the convergence is very slow for the PS and the errors are not decreasing monotonically. In fact for h = 0.05 the iterations fail to converge starting from t = 0.7. For h = 0.02 the PS fails to converge even at t = 0.5. Similar behaviour is observed for the NS. For h = 0.05 the errors start to increase after a few iterations and the solution diverges at t = 0.8. The reason for this behaviour of the PS and the NS is the bound given in (1.2). As the mesh size decreases the time step size has to be reduced in order to guarantee the convergence for the PS and the NS [3]. However this constraint is not there for the MS. Although there is a threshold on the time step size for the MS’s convergence, this threshold does not depend on the mesh size and so in all the cases shown in Fig. 3, the errors for the MS are decreasing and the convergence is faster when compared to the NS or the PS. This shows that for numerical computations the MS is more stable than the NS or the PS in general. As the convergence rate of the LS is quite stable with respect to mesh sizes for L = 1, the error behaviour of the LS has been plotted only for h = 0.05 and L = 1 in Fig. 3. LS also has an upper bound on the time steps for linear convergence that does not depend upon the spatial discretization [7]. An interesting observation is that for larger time step sizes the convergence of the LS can actually be faster than the convergence of the MS. This is due to the fact that in the computations Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

14

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

Fig. 4. Convergence rate (α ) vs. time step size (τ ) for h = 0.05 at t = 0.5. The convergence rates are calculated as the geometric mean of ∥pin+1 − pin ∥/∥pin − pni−1 ∥ over first 10 iterations. The dotted black line shows the slope of 0.5.

Fig. 5. The dependence of the error behaviour on M for the MS. (left) Results for τ = 0.01; (right) results for τ = 0.1.

for τ = 0.1, Mτ ∼ L and so Mτ + S ′ (p) can be greater than L. So the apparent L value for the MS can be larger than the value of L required to be imposed on the LS. This might make the convergence of the MS slowerthan that of  the  LS.  Fig. 4 shows how the convergence rate α , calculated here as the geometric average of pin+1 − pin /pin − pin−1  over the first 10 iterations, changes with τ . From Theorem 4.1 one expects that for small enough τ the convergence rate should scale √ with τ which implies that the slope of the solid line in Fig. 4 should nearly be 12 for small τ values. This is indeed the case. For τ ≤ .01 the line is almost parallel to the reference dashed line representing a slope of 21 . It is to be noted though that Theorem 4.1 was proved for the non-degenerate case whereas this test case is degenerate. Next we study the effect of the choice of M on the rate of convergence. One expects from Theorem 4.1 and (4.21) that for small values of M the condition (Lin − ∂u z) > 0 will not be satisfied whereas large values of M would result in a slower convergence. This behaviour is observed precisely if one varies L in the case of the LS [7,16]. From Fig. 5 we see that it describes the MS as well. For both time step sizes τ = 0.01 and τ = 0.1 we get that the optimal M is close to 10. If one chooses M value away from the optimal M then α increases. As M → 0 the rates tend towards the rate of the PS. In the τ = 0.01 case the effect of an optimal M is less pronounced than in the τ = 0.1 case and from Fig. 2 (right) one can speculate that it is even less important in the τ = 0.001 case. This a priori knowledge can be useful when developing the linearisation scheme as one can decide what M to choose by running the computation for only the first time step with a coarser grid. In our case, we chose M = 10 for the results presented earlier for this reason. Finally, we present a comparison study with Example 1 of [7]. The functional relationships are taken exactly as stated in [7] but different meshes and linear solvers are used. The initial condition used in this example is p(x, y, 0) =

⎧ ⎨pvad ⎩−y +

1 4

in Ω ∩ {y <

1 4

},

in Ω ∩ {y >

1 4

}.

(5.7)

Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

15

Fig. 6. Iterations required by different schemes to have ∥pin − pin−1 ∥ < 10−5 , for a fixed τ = 1 and pv ad = −2. The parameters are taken from [7, Fig. 1]. Table 2 Iterations required by different schemes to have ∥pin − pin−1 ∥ < 10−5 , for a fixed h = pv ad = −2. The parameters are taken from [7, Table 1].

τ MS M = 0.01 PS M = 0 NS LS L = 0.25 LS L = 0.15

1 40

and

Number of iterations 1

0.1

0.01

0.001

18 19 – 54 35

22 22 – 50 33

12 12 – 39 26

7 7 7 154 99

Here pv ad < 0 is a constant. The gravity points towards negative y direction and the results are for the first time step, i.e. n = 1. As the initial condition is discontinuous, the NS fails to converge in all cases in our computations. However, the NS converged in most cases if the discontinuity was regularized, e.g. by considering a linear interpolation over a small interval of length 0.1. After this step the results obtained match closely the values given in [7]. In particular, we partly reproduce [7, Figure 1] with pv ad = −2 and [7, Table 1] with pv ad = −3 in Fig. 6 and Table 2 respectively. The result shows that the NS is the fastest when it converges but it requires smoother initial conditions and smaller time step sizes to converge. The PS and the MS have comparable stability properties. The MS, as before, is at least as fast as the PS in all cases. The LS is relatively slower and taking a smaller value of L (L = 0.15) compared to the Lipschitz constant (L ≈ 0.25) speeds up the convergence, but makes the convergence rate more susceptible to changes in mesh size. This was seen in Fig. 2 (left) too. The convergence rate decreases with the time step size for all the schemes except for the LS. A drastic increase in number of iterations required is seen for the LS at τ = 0.001, see Table 2. This is explained by the fact that the convergence order is L/(L + C τ ), which approaches 1 when τ goes to 0. M = 0.01 is used for the MS in Fig. 6 and Table 2 as it gives the optimal convergence rate. 6. Conclusion In this paper we propose a linearization scheme for a general class of nonlinear parabolic partial differential equations. We have given a rigorous convergence proof of the scheme and compared its behaviour with that of other linearization schemes: the Newton scheme (NS), the modified Picard scheme (PS) or the L-scheme (LS), are generally used to solve the sequence of elliptic equations obtained from time-discretization of such problems. The NS and the PS have the drawback that they converge only if the initial guess for the iterations is close enough to the solution of the elliptic problem. For the concerned sequence of elliptic equations, this leads to a severe restriction on time step size. On the other hand, the LS converges irrespective of the initial guess but is much slower than the schemes mentioned above. To resolve these issues, a combination of the LS and the PS, termed modified L-scheme (MS) in this context, is proposed. The MS uses local estimates and the solution of the previous time step as the initial guess to improve the convergence behaviour of the LS. This is shown first for quasilinear equations that have linear diffusivity and advection terms. It is proved that the scheme converges linearly irrespective of the spatial discretization and for any time step, even in degenerate cases. Moreover, for small time step sizes, the linear convergence rate is proportional to the time step size, implying that the scheme converges faster as the time step size is made smaller. Next, this result is generalized to the case when the diffusivity and the advection terms are nonlinear. It is proved that if the time step size is smaller than an upper-bound which is independent of the spatial discretization, the scheme Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

16

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

converges even for degenerate cases. Linear convergence is achieved in the non-degenerate case with a convergence rate that is proportional to the square root of the time step size for sufficiently small time steps. Finally, numerical results are presented for the Richards equation for all the schemes mentioned above. It is seen that the MS is faster than the PS and the LS. Moreover, the MS is more stable than the NS or the PS in the sense that the MS converges for larger time steps. Numerically it is observed that the MS indeed gives convergence rates proportional to the square root of the time step size. The final numerical results imply that the convergence rate of the MS can be controlled by tuning the parameter M. Acknowledgements K. Mitra is supported by Shell and the Netherlands Organisation for Scientific Research (NWO), Netherlands through the CSER programme (project 14CSER016) and by the Hasselt University, Belgium through the project BOF17BL04. The research of I.S. Pop is supported by the Research Foundation-Flanders (FWO), Belgium through the Odysseus programme (project G0G1316N). The authors thank the anonymous referees for their valuable comments and suggestions that helped improve the manuscript. Appendix. Proof of Propositions 3.1 and 3.2

Proof of Proposition 3.1. We consider the L-scheme for (3.4). Let uin be the solution of the ith iteration of the L-scheme (2.7) initiated by u0n = un−1 . We use the properties of the sequence {uin }i∈N to prove Proposition 3.1. From (2.7) and (3.1), uin satisfies (Luin , φ ) + τ (D∇ uin , ∇φ ) = (Luin−1 , φ ) − (b(uin−1 ) − b(un−1 ), φ ) + τ (F⃗ , ∇φ ) + τ (fni−1 , φ ).

(A.1)

Let ρ i denote the difference between consecutive iterates, i.e.

ρ i = uin+1 − uin . Observe that for i = 0, ρ 0 solves the problem L(ρ 0 , φ ) + τ (D∇ρ 0 , ∇φ ) = τ [(f (⃗ x, tn , un−1 ), φ ) + (F⃗ , ∇φ ) − (D∇ un−1 , ∇φ )]. As ∇ · (D∇ un−1 ) ∈ L∞ (Ω ) is assumed in Proposition 3.1, ρ 0 satisfies the equation Lρ 0 − τ ∇ · (D∇ρ 0 ) = τ [−∇ · F⃗ + f (⃗ x, tn , un−1 ) + ∇ · (D∇ un−1 )], in the classical sense (see [17][Section 6.3]). As the term f (⃗ x, tn , un−1 ) − ∇ · F⃗ + ∇ · (D∇ un−1 ) is bounded, the maximum principle applies [17][Section 6.4] and therefore a C0 > 0, independent of τ , exists such that

∥ρ 0 ∥L∞ (Ω ) = C0 τ . We claim that there exists a β ∈ (0, 1) such that

∥ρ k ∥L∞ ≤ β∥ρ k−1 ∥L∞ .

(A.2)

The value of β will be specified later. The proof of (A.2) is similar to the proof of Lemma 3.1. Subtracting (A.1) for i from the same equation written for (i + 1)th, we get an expression similar to (3.5), i.e. L(ρ i , φ ) + τ (D∇ρ i , ∇φ ) = ((L − ∂u z(ζ ))ρ i−1 , φ ), where ζ : Ω → R is a function satisfying ζ ∈ I (uin , uin−1 ). Let Ci−1 = ∥ρ i−1 ∥L∞ (Ω ) /τ . Taking φ = [ρ i − β Ci−1 τ ]+ as test function in the above and using the bounds for D we get

2

L(ρ i , [ρ i − β Ci−1 τ ]+ ) + τ Dm ∇[ρ i − β Ci−1 τ ]+  ≤



∫ Ω

(L − ∂u z(ζ ))ρ i−1 [ρ i − β Ci−1 τ ]+ ,

which implies that

2

L∥[ρ i − β Ci−1 τ ]+ ∥2 + τ Dm ∇[ρ i − β Ci−1 τ ]+  ≤



∫ Ω

((L − ∂u z(ζ ))ρ i−1 − Lβ Ci−1 τ )[ρ i − β Ci−1 τ ]+ .

With this, we analyse the sign of the expression (L − ∂u z(ζ ))ρ i−1 − Lβ Ci−1 τ . If ρ i−1 ≤ 0 then this is clearly negative. If ρ i−1 > 0, since ∥ρ i−1 ∥L∞ (Ω ) = Ci−1 τ and 0 ≤ L − ∂u z(ζ ) ≤ L − m, see (2.4), by taking

β = (L − m)/L < 1, one obtains (L − ∂u z(ζ ))ρ i−1 − Lβ Ci−1 τ ≤ [(L − m) − Lβ]Ci−1 τ = 0. This shows that ρ i ≤ β Ci−1 τ a.e. The inequality ρ i ≥ −β Ci−1 τ follows in a similar fashion. Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.

K. Mitra, I.S. Pop / Computers and Mathematics with Applications (

)

–

17

∑∞

∑∞

i i i+1 i ∞ − uin , it implies that uin → un in the Observe that i=0 β = C0 τ /(1 − β ) and as ρ = un i=0 ∥ρ ∥L (Ω ) ≤ C0 τ L (Ω )-norm as i → ∞. Note that, uin also converges to un in the H 1 (Ω ) norm as shown in [13] and due to the uniqueness of the weak solutions of (3.4) (see [20]), un is the unique solution of (3.4). Finally, defining Λ = C0 /(1 − β ) we get

∞

∥un − un−1 ∥L∞ (Ω ) ≤

∞ ∑

∥ρ i ∥L∞ (Ω ) ≤ Λτ .

i=0

Rewriting (3.4) as ( ∇ · ( D ∇ un ) , φ ) =

(

1

τ

(b(un ) − b(un−1 )) + ∇ · F − fn , φ

)

,

we see right away that ∇ · (D∇ un ) = τ1 (b(un ) − b(un−1 )) − fn + ∇ · F a.e. and as the terms on the right hand side are bounded in the L∞ (Ω )-norm, so is ∇ · (D∇ un ). □ Proof of Proposition 3.2. The proof is almost identical to the  proof of Proposition  3.1. We subtract (3.4) from (3.5) and follow the steps of the proof of Proposition 3.1 to obtain that uin − un L∞ (Ω ) ≤ β uin−1 − un L∞ (Ω ) if

β ≥ max

{

Lin − ∂u z(ζ )

}

Lin

,

ζ ∈ I (un , uin−1 ).

Observe that Lin ≥ m + Mτ . Moreover, from (3.6) we obtain, Lin − ∂u z(ζ ) ≤ 2Mτ for M ≥ M0 . Hence, Lin − ∂u z(ζ ) Lin

≤

2Mτ m + Mτ

≤

2Mτ m

.

(A.3)

By defining,

β=

2Mτ

,

m we observe that for τ <

(A.4) m , 2M

β < 1 and β = O(τ ). □

References [1] L. Bergamaschi, M. Putti, Mixed finite elements and Newton-type linearizations for the solution of Richards’ equation, Internat. J. Numer. Methods Engrg. 45 (8) (1999) 1025–1046. [2] F. Lehmann, P.H. Ackerer, Comparison of iterative methods for improved solutions of the fluid flow equation in partially saturated porous media, Transp. Porous Media 31 (3) (1998) 275–292. [3] F.A. Radu, I.S. Pop, P. Knabner, Newton-type methods for the mixed finite element discretization of some degenerate parabolic equations, in: Numerical Mathematics and Advanced Applications, Springer, 2006, pp. 1192–1200. [4] J.E. Jones, C.S. Woodward, Newton–Krylov-multigrid solvers for large-scale, highly heterogeneous, variably saturated flow problems, Adv. Water Resour. 24 (7) (2001) 763–774. [5] K. Brenner, C. Cances, Improving Newton’s method performance by parametrization: The case of the Richards equation, SIAM J. Numer. Anal. 55 (4) (2017) 1760–1785. [6] M.A. Celia, E.T. Bouloutas, R.L. Zarba, General mass-conservative numerical solution for the unsaturated flow equation, Water Resour. Res. 26 (7) (1990) 1483–1496. [7] F. List, F.A. Radu, A study on iterative methods for solving Richards’ equation, Comput. Geosci. 20 (2) (2016) 341–353. [8] W. Jäger, J. Kačur, Solution of porous medium type systems by linear approximation schemes, Numer. Math. 60 (1) (1991) 407–427. [9] W. Jäger, J. Kačur, Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes, ESAIM Math. Model. Numer. Anal. 29 (5) (1995) 605–627. [10] J. Kačur, Solution to strongly nonlinear parabolic problems by a linear approximation scheme, IMA J. Numer. Anal. 19 (1) (1999) 119–145. [11] E. Magenes, R.H. Nochetto, C. Verdi, Energy error estimates for a linear scheme to approximate nonlinear parabolic problems, ESAIM Math. Model. Numer. Anal. 21 (4) (1987) 655–678. [12] I.S. Pop, W.A. Yong, A numerical approach to degenerate parabolic equations, Numer. Math. 92 (2) (2002) 357–381. [13] I.S. Pop, F.A. Radu, P. Knabner, Mixed finite elements for the Richards equation: linearization procedure, J. Comput. Appl. Math. 168 (1) (2004) 365–373. [14] I.S. Pop, W.A. Yong, On the existence and uniqueness of a solution for an elliptic problem, Stud. Univ. Babeş-Bolyai Math. 45 (2000) 97–107. [15] F.A. Radu, J.M. Nordbotten, I.S. Pop, K. Kumar, A robust linearization scheme for finite volume based discretizations for simulation of two-phase flow in porous media, J. Comput. Appl. Math. 289 (2015) 134–141. [16] D. Seus, K. Mitra, I.S. Pop, F.A. Radu, C. Rohde, A linear domain decomposition method for partially saturated flow in porous media, Comput. Methods Appl. Mech. Engrg. 333 (2018) 331–355. [17] L.C. Evans, Partial Differential Equations, Wiley Online Library, 1988. [18] O.A. Ladyzhenskaya, V. Solonnikov, N. Uraltseva, Linear and quasilinear parabolic equations of second order, in: Translation of Mathematical Monographs, AMS, Rhode Island, 1968. [19] H.W. Alt, S. Luckhaus, Quasilinear elliptic-parabolic differential equations, Math. Z. 183 (3) (1983) 311–341. [20] F. Otto, L1 -contraction and uniqueness for quasilinear elliptic–parabolic equations, J. Differential Equations 131 (1) (1996) 20–38. [21] J. Kačur, Method of Rothe in evolution equations, Lecture Notes in Math. 1192 (1986) 23–34. [22] M. Slodicka, A robust and efficient linearization scheme for doubly nonlinear and degenerate parabolic problems arising in flow in porous media, SIAM J. Sci. Comput. 23 (5) (2002) 1593–1614. [23] F.A. Radu, K. Kumar, J.M. Nordbotten, I.S. Pop, A robust, mass conservative scheme for two-phase flow in porous media including Hölder continuous nonlinearities, IMA J. Numer. Anal. 38 (2) (2018) 884–920. [24] M. Borregales, F.A. Radu, K. Kumar, J.M. Nordbotten, Robust iterative schemes for non-linear poromechanics, 2017, arXiv preprint arXiv:1702.00328. [25] J. Wang, C.V. Pao, Finite difference reaction–diffusion equations with nonlinear diffusion coefficients, Numer. Math. 85 (3) (2000) 485–502. [26] O.A. Ladyzhenskaya, N.N. Uraltseva, L. Ehrenpreis, Linear and Quasilinear Elliptic Equations, Academic Press, 1968. [27] R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, Handb. Numer. Anal. 7 (2000) 713–1018.

Please cite this article in press as: K. Mitra, I.S. Pop, A modified L-scheme to solve nonlinear diffusion problems, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.09.042.