Convergence analysis of domain decomposition based time integrators for degenerate parabolic equations

Nonlinear parabolic equations of the form

equipped with suitable boundary and initial conditions, are frequently encountered in applications. If the diffusion constant \(D(u,\nabla u)\) vanishes for some values of u and \(\nabla u\), i.e., the equation is degenerate, one obtains a quite different dynamics compared to the linear case. The two main nonlinear features are finite speed of propagation and the absence of parabolic smoothening of the solution. Concrete applications can, e.g., be found when modelling gas flow through porous media, phase transitions and population dynamics. A survey of such applications is given in [25, Section 1.3 and Chapter 2]. In order to keep the presentation as clear-cut as possible, we will mostly ignore the presence of lower-order advection and reactions terms.

Approximating the solution of a partial differential equation typically results in large-scale computations, which require the usage of parallel and distributed hardware. One possibility to design numerical schemes that make use of such hardware is to decompose the equation’s domain into a family of subdomains. The domain decomposition method then consists of an iterative procedure where, in every step, the equation is solved independently on each subdomain and the resulting solutions are thereafter communicated to the adjacent subdomains. This independence of the decomposed equations and the absence of global communication enables the parallel and distributed implementation of domain decomposition methods. For linear parabolic equations the common procedure is to first discretize the equation in time by a standard implicit integrator. Then an elliptic equation on \(\varOmega \) is obtained in every time step, which is iteratively solved by a domain decomposition based discretization. We refer to the monographs [19, 21, 24] for an in-depth treatment of this approach. Another possibility is to apply the domain decomposition method to the full space-time domain \(\varOmega \times (0,T)\), which leads to an iterative procedure over parabolic problems that can be parallelized both in space and time; see, e.g., [12, 13, 15].

When considering nonlinear parabolic problems one finds that there are hardly any results concerning the analysis of domain decomposition based schemes. Two exceptions are the papers [17, 18], where domain decomposition schemes are analyzed for non-degenerate quasilinear parabolic equations and the degenerate two-phase Stefan problem, respectively. The lack of results in the context of degenerate equations is rather surprising from a practical point of view, as the equations’ finite speed of propagation is ideal for applying domain decomposition strategies. For example, a solution that is initially zero in parts of the domain \(\varOmega \) will in each time step only propagate to a small number of neighboring subdomains, which limits the computational work considerably. However, from a theoretical perspective the lack of convergence results is less surprising. The issue is that the standard domain decomposition schemes all link together the equations on the subdomains via boundary conditions. As the solutions of degenerate parabolic equations typically lack higher-order regularity, making sense of such boundary linking is, at the very least, challenging.

Examples of overlapping domain decompositions \(\{ \varOmega _{\ell } \}_{\ell =1}^{s}\) of a domain \(\varOmega \subset \mathbb {R}^{2}\), with \(s=4\) subdomains (left) and \(s=2\) subdomains that are further decomposed into families of pairwise disjoint sets (right), respectively

Full size image

In order to remedy this, we propose to directly introduce the domain decomposition in the time integrator via an operator splitting procedure. More precisely, let \(\{ \varOmega _{\ell } \}_{\ell =1}^{s}\) be an overlapping decomposition of the spatial domain \(\varOmega \), as exemplified in Fig.  1. On these subdomains we introduce the partition of unity \(\{ \chi _{\ell } \}_{\ell =1}^{s}\) and the operator decomposition, or splitting,

Two possible (formally) first-order integrators are then the sum splitting

which represents a “quick and dirty” scheme that is straightforward to parallelize, and the Lie splitting

which is usually more accurate but requires a further partitioning of the subdomains \(\varOmega _\ell \) in order to enable parallelization, as illustrated in Fig. 1. In contrast to the earlier domain decomposition based schemes, where an iterative procedure is required with possibly many instances of boundary communications, one time step of either splitting scheme only needs the solution of s elliptic equations together with the communication of the data related to the overlaps. Similar splitting schemes have, e.g., been considered in the papers [2, 16, 20, 26] when applied to linear, and to some extent semilinear, parabolic problems. However, there does not seem to be any analysis applicable to degenerate, or even quasilinear, parabolic equations in the literature.

Hence, the goal of this paper is twofold. First, we aim to derive a new energetic, or variational, framework that allows a proper interpretation of the operator decomposition (2) for two commonly occurring families of degenerate parabolic equations. These are the p-Laplace type evolutions, where the prototypical example is given by \(D(u,\nabla u)=|\nabla u|^{p-2}\), and the porous medium type equations, where \(D(u,\nabla u)=(p-1)|u|^{p-2}\) in the simplest case. For the porous medium application we will use the strategic reformulation

of the decomposition (2), in order to enable an energetic interpretation.

Secondly, we will strive to obtain a general convergence analysis for the domain decomposition based time integrators, including the sum and Lie splitting schemes. The main idea of the convergence analysis is to introduce the nonlinear Friedrich extensions of the operators f and \(f_{\ell }\), via our new abstract energetic framework, and then to employ a Lax-type result from the nonlinear semigroup theory [5].

Throughout the analysis \(\varOmega \subset \mathbb {R}^d\), \(d\ge 1\), will be an open, connected and bounded set and the parameter \(p\in (1,\infty )\) is fixed. Next, let \(\{ \varOmega _{\ell } \}_{\ell =1}^{s}\) be a family of overlapping subsets of \(\varOmega \) such that \(\bigcup _{\ell =1}^s \varOmega _{\ell } = \varOmega \) holds. Here, each \(\varOmega _{\ell }\) is either an open connected set, or a union of pairwise disjoint open, connected sets \(\varOmega _{\ell ,k}\) such that \(\bigcup _{k=1}^{r}\varOmega _{\ell ,k} = \varOmega _{\ell }\). On \(\{ \varOmega _{\ell } \}_{\ell =1}^{s}\) we introduce the partition of unity \(\{ \chi _{\ell } \}_{\ell =1}^{s}\subset C^{\infty }(\varOmega )\) such that

For details on the construction of explicit domain decompositions \(\{ \varOmega _{\ell } \}_{\ell =1}^{s}\) and partitions of unity \(\{ \chi _{\ell } \}_{\ell =1}^{s}\) we refer to [2, Section 3.2] and [20, Section 4.1].

The related weighted Lebesgue space \(L^p(\varOmega _{\ell },\chi _{\ell })\) can now be defined as the set of all measurable functions u on \(\varOmega _{\ell }\) such that the norm

is finite. The space \(L^p(\varOmega _{\ell },\chi _{\ell })\) is a reflexive Banach space, which follows by observing that the map \(G : L^p(\varOmega _{\ell },\chi _{\ell }) \rightarrow L^p(\varOmega _{\ell }):u\mapsto \chi _{\ell }^{1/p} u\) is an isometric isomorphism [9, Chapter 1]. We will also make frequent use of the product space \(L^p(\varOmega _{\ell },\chi _{\ell })^k\), equipped with the norm

which is again a reflexive Banach space [1, Theorem 1.23].

The idea is now to introduce a suitable pivot space H and the energetic spaces \(V\subset H\) and \(V_{\ell }\subset H\), such that the action of the linear differential operator, denoted by \(\delta \), arising in the nonlinear vector field and its weighted counterpart, \(\chi _{\ell }\delta \), can be interpreted as elements in \(L^p(\varOmega )^{k}\) and \(L^p(\varOmega _{\ell },\chi _{\ell })^{k}\), respectively. As a concrete example, consider the p-Laplacian \(u\mapsto \nabla \cdot (|\nabla u|^{p-2}\nabla u)\). Here, \(\delta =\nabla \) is the distributional gradient, \(H=L^{2}(\varOmega )\) and \(V=W^{1,p}(\varOmega )\). The space \(V_{\ell }\) becomes the intersection between \(L^{2}(\varOmega )\) and a weighted version of \(W^{1,p}(\varOmega _{\ell })\). The p-Laplace and porous medium type equations require different pivot and energetic spaces, and we therefore start by developing an abstract setting fitting both equation families.

Let \(\left( H, \left( \cdot ,\cdot \right) _{H}\right) \) be a real Hilbert space space and denote the space of distributions on \(\varOmega \) by \(\mathscr {D}'(\varOmega )\). For a given \(k\ge 1\) we introduce the linear operator

which is assumed to be continuous in the following fashion.

Assumption 1

If \(\lim _{n\rightarrow \infty } u_{n}=u\) in H then, for \(j=1,\ldots ,k\),

As the regularity of the weights \(\chi _{\ell }\) implies that \(\chi _{\ell }\varphi \in C^{\infty }_0(\varOmega _\ell )\) for all \(\varphi \in C^{\infty }_0(\varOmega )\), we can define the product \(\chi _{\ell }\delta u\) by

The energetic spaces V and \(V_{\ell }\) are then given as

and

respectively. On the energetic spaces we consider the operators

where \(\delta _{p}\) maps \(u\in V\) to the corresponding \(L^p(\varOmega )\) functions that \(\delta u\) can be represented by, and \(\delta _{p,\ell }\) maps \(u\in V_{\ell }\) to the corresponding \(L^p(\varOmega _{\ell },\chi _{\ell })\) functions that \(\chi _{\ell }\delta u\) can be represented by, respectively.

Lemma 1

\(V = \bigcap _{\ell =1}^s V_{\ell }\).

Proof

For an arbitrary \(u \in V\) it follows, for \(\ell = 1,\dots ,s\), that

for every \(\varphi \in C_0^{\infty }(\varOmega )\) and \(j = 1,\dots ,k\). As \((\delta _p u)_j|_{\varOmega _{\ell }} \in L^p(\varOmega _{\ell }) \subseteq L^p(\varOmega _{\ell },\chi _{\ell } )\), we have a representation of \((\delta u)_j\) in \(L^p(\varOmega _{\ell },\chi _{\ell } )\), i.e., \(u \in V_{\ell }\) for every \(\ell = 1,\dots , s\). Hence, \(V \subseteq \bigcap _{\ell =1}^s V_{\ell }\).

Next, assume that \(u\in \bigcap _{\ell =1}^s V_{\ell }\). Then we can write

for every \(\varphi \in C_0^{\infty }(\varOmega )\) and \(j =1,\dots ,k\). Let \(w_{\ell ,j}\) be the zero extension of \((\delta _{p,\ell } u)_j\) to the whole of \(\varOmega \). We can then define the measurable function \(v_{j}\) on \(\varOmega \) as \(v_j = \sum _{\ell =1}^{s} \chi _{\ell } w_{\ell ,j}\), which satisfies

Furthermore, the \(L^{p}(\varOmega )\) norm of \(v_{j}\) can be bounded by

This yields that \((\delta _p u)_j = v_j \in L^p(\varOmega )\) for \(j =1,\dots ,k\), i.e., \(u\in V\) and we thereby have the identification \(V=\bigcap _{\ell =1}^s V_{\ell }\). \(\square \)

Lemma 2

If Assumption 1 holds, then the operators \(\delta _p\) and \(\delta _{p,\ell }\), \(\ell =1,\dots ,s\), are linear and closed.

Proof

The linearity of the operators is clear, since \(\delta \) is a linear operator. Let the sequence \(\{u_n\}_{n\in \mathbb {N}}\subset V_{\ell }\) satisfy

Assumption 1 then yields that

for every \(\varphi \in C_0^{\infty }(\varOmega )\) and \(j=1,\dots ,k\). Hence, \((\chi _{\ell }\delta u)_j\) can be represented by the \(L^p(\varOmega _{\ell },\chi _{\ell })\) function \(v_j\), i.e., \(\delta _{p,\ell } u = v\) holds and the operator \(\delta _{p, \ell }\) is therefore closed. The closedness of \(\delta _{p}\) follows by the same line of reasoning. \(\square \)

On the energetic spaces V and \(V_{\ell }\), \(\ell =1,\dots ,s\), we define the norms

respectively.

Lemma 3

If Assumption 1 holds, then the spaces \((V, \Vert \cdot \Vert _V)\) and \((V_{\ell }, \Vert \cdot \Vert _{V_{\ell }})\), \(\ell =1,\dots ,s\), are reflexive Banach spaces.

Proof

Consider the reflexive Banach space \(X=H\times L^p(\varOmega _{\ell },\chi _{\ell })^k\), equipped with the norm \(\Vert (u_{1},u_{2})\Vert _{X}=\Vert u_{1}\Vert _{H}+\Vert u_{2}\Vert _{L^p(\varOmega _{\ell },\chi _{\ell })^k}\), and introduce the linear and isometric operator

The graph of the closed operator \(\delta _{p,\ell }\) coincides with the image \(G(V_{\ell })\), which makes \(G(V_{\ell })\) a closed linear subset of X. Here, \((G(V_{\ell }),\Vert \cdot \Vert _{X})\) is a reflexive Banach space [1, Theorem 1.22] and, as G is isometric, it is isometrically isomorphic to \((V_{\ell },\Vert \cdot \Vert _{V_{\ell }})\). Hence, the latter is also a reflexive Banach space. The same line of argumentation yields that V is a reflexive Banach space. \(\square \)

Hereafter, we will assume the following.

Assumption 2

The set V is dense in H.

Under this assumption it also holds that \(V_{\ell }\) is a dense subsets of H. By the construction of the energetic norms, one then obtains that the reflexive Banach spaces \((V, \Vert \cdot \Vert _V)\) and \((V_{\ell }, \Vert \cdot \Vert _{V_{\ell }})\) are densely and continuously embedded in H and we have the following Gelfand triplets

Here, the density of \(H^*\) in \(V^*\) and \(V^*_{\ell }\), respectively, follows, e.g., by [14, Bemerkung I.5.14]. For future reference, we denote the dual pairing between a Banach space X and its dual \(X^*\) by \(\langle \cdot \, ,\cdot \rangle _{X^{*}\times X}\), and the Riesz isomorphism from H to \(H^{*}\) by

Here, the Riesz isomorphism satisfies the relations

for all \(u\in H\), \(v\in V\) and \(v_{\ell }\in V_{\ell }\).

Remark 1

Throughout the derivation of the energetic framework we have assumed that the partition of unity \(\{ \chi _{\ell } \}_{\ell =1}^{s}\) consists of elements in \(C^{\infty }(\varOmega )\). This is somewhat restrictive from a numerical point of view, but this regularity is required if nothing else is known about the operator \(\delta : H \rightarrow \mathscr {D}'(\varOmega )^k\). Fortunately, in concrete examples; see Sects. 6 and 7, one commonly has that \(\delta (H)\subseteq H^{-1}(\varOmega )^k\). If we then choose a partition of unity \(\{ \chi _{\ell } \}_{\ell =1}^{s}\) in \(W^{1,\infty }(\varOmega )\), we have the property that \(\chi _{\ell }\varphi \in H^{1}_{0}(\varOmega )\) for every \(\varphi \in H^{1}_{0}(\varOmega )\), and we can once more derive the above energetic setting by testing with functions \(\varphi \) in \(H^{1}_{0}(\varOmega )\), instead of in \(C^{\infty }_{0}(\varOmega )\).

With the function spaces in place, we are now able to define the general energetic extensions of our vector fields.

Assumption 3

For a fixed \(p\in (1,\infty )\), let \(\alpha : \varOmega \times \mathbb {R}^k \rightarrow \mathbb {R}^k\) fulfill the properties below.

\((\alpha _{1})\) :

The map \(\alpha : \varOmega \times \mathbb {R}^k \rightarrow \mathbb {R}^k\) fulfills the Carathéodory condition, i.e., \(z \mapsto \alpha (x,z)\) is continuous for a.e. \(x\in \varOmega \) and \(x \mapsto \alpha (x,z)\) is measurable for every \(z\in \mathbb {R}^k\).

\((\alpha _{2})\) :

The growth condition \(|\alpha (x,z)| \le c_1 |z|^{p-1} +c_2(x) \) holds for a.e. \(x\in \varOmega \) and every \(z\in \mathbb {R}^k\), where \(c_1>0\) and \(c_2\in L^{p/(p-1)}(\varOmega )\) is nonnegative.

\((\alpha _{3})\) :

The map \(\alpha \) is monotone, i.e., for every \(z,\tilde{z} \in \mathbb {R}^k\) and a.e. \(x\in \varOmega \) the inequality \((\alpha (x,z) - \alpha (x,\tilde{z}))\cdot (z - \tilde{z}) \ge 0 \) holds.

\((\alpha _{4})\) :

The map \(\alpha \) is coercive, i.e., there exists \(c_3>0\) and \(c_4\in L^1(\varOmega )\) such that for every \(z\in \mathbb {R}^k\) and a.e. \(x\in \varOmega \) the condition \(\alpha (x,z) \cdot z \ge c_3 |z|^p - c_4(x)\) holds.

Compare with [27, Section 26.3].

We introduce the full energetic operator \(F : V \rightarrow V^*\) as

The operator F is well defined, as \(\delta _p v \in L^p(\varOmega )^k\) for \(v\in V\) and by (\(\alpha _{2}\)) we obtain that \(\alpha (\delta _p v) \in L^{p/(p-1)}(\varOmega )^k \cong \left( L^{p}(\varOmega )^k\right) ^*\). Furthermore, we define the decomposed energetic operators \(F_{\ell } : V_{\ell } \rightarrow V^*_{\ell }\), \(\ell =1,\dots , s\), by

These operators are well defined, as

is finite for every \(u,v\in V_{\ell }\), due to (\(\alpha _{2}\)). This family of operators is a decomposition of F, as it fulfills

We can now derive the basic properties of the (perturbed) energetic operators, namely, strict monotonicity and coercivity of \(\gamma + hF\) and \(\gamma + hF_{\ell }\). These operators are also hemicontinuous, i.e., the maps \(t \mapsto \langle (\gamma + hF_{\ell })(u +tv), w \rangle _{V_{\ell }^*\times V_{\ell }}\) are continuous on [0, 1] for \(u,v,w\in V_{\ell }\).

Lemma 4

If the Assumptions 1–3 hold and \(h>0\) , then the operators \(\gamma + hF: V \rightarrow V^*\) and \(\gamma + hF_{\ell }: V_{\ell } \rightarrow V^*_{\ell }\), \(\ell =1,\ldots ,s\), are strictly monotone, coercive and hemicontinuous.

Proof

We will only derive the properties for \(\gamma + hF_{\ell }\), as the same argumentation holds for \(\gamma + hF\). The strict monotonicity of the operator follows using (\(\alpha _{3}\)), as

holds for all \(u,v\in V_{\ell }\) with \(u\ne v\).

Next, we prove the coercivity of \(\gamma + hF_{\ell }\). By assumption (\(\alpha _{4}\)), we have

for every \(u \in V_{\ell }\). Hence, we have the limit

as \(\Vert u\Vert _{V_{\ell }} \rightarrow \infty \), which implies the coercivity of \(\gamma + hF_{\ell }\).

Last, we prove that \(F_{\ell }\) is hemicontinuous. Consider a sequence \(\{t_n\}_{n\in \mathbb {N}}\) in [0, 1] with limit t and introduce

As \(\lim _{n\rightarrow \infty } g(t_{n},x)=g(t,x)\) holds for almost every \(x\in \varOmega _{\ell }\), due to (\(\alpha _{1}\)), and

where the right-hand side is an \(L^1(\varOmega _{\ell })\) element, we obtain that

by the dominated convergence theorem. This implies that \(F_{\ell }\) is hemicontinuous, and the same trivially holds for \(\gamma + hF_{\ell }\).\(\square \)

Corollary 1

If the Assumptions 1–3 hold and \(h>0\) , then the operators \(\gamma + hF: V \rightarrow V^*\) and \(\gamma + hF_{\ell }: V_{\ell } \rightarrow V^*_{\ell }\), \(\ell =1,\ldots ,s\), are all bijective.

Proof

As \(\gamma + hF: V \rightarrow V^*\) and \(\gamma + hF_{\ell }: V_{\ell } \rightarrow V^*_{\ell }\) are all, by Lemma 4, strictly monotone, coercive and hemicontinuous, their bijectivity follows by the Browder–Minty theorem; see, e.g., [27, Theorem 26.A].\(\square \)

The energetic setting is too general for the convergence analysis that we have in mind. We therefore introduce the nonlinear Friedrich extensions of our vector fields, i.e., we restrict the domains of the energetic operators such that they become (unbounded) operators on the pivot space H. More precisely, we define the Friedrich extension \(f: {D}(f) \subseteq H \rightarrow H\) of the full vector field by

Analogously, we introduce the Friedrich extensions \(f_{\ell }: {D}(f_{\ell }) \subseteq H \rightarrow H\), \(\ell =1,\dots ,s\), of the decomposed vector fields by

We can derive that these operators are in fact maximal dissipative, i.e.,

Lemma 5

If the Assumptions 1–3 hold, then the operators \(f:{D}(f) \subseteq H \rightarrow H\) and \(f_{\ell }:{D}(f_{\ell }) \subseteq H \rightarrow H\), \(\ell =1,\dots ,s\), are all maximal dissipative.

Proof

By (\(\alpha _{3}\)) of Assumption 3, we have that

for all \(u,v\in {D}(f_{\ell })\). Next, for given \(h>0\) and \(v\in H\) one has, in virtue of Corollary 1, that there exists a unique \(u\in V_{\ell }\) such that \((\gamma +hF_{\ell })u=\gamma v\), or equivalently

Hence, \(u\in {D}(f_{\ell })\) and \((I-hf)u=v\) in H, i.e., \({R}(I-hf_{\ell })=H\). The operators \(f_{\ell }\), \(\ell =1,\dots , s\), are therefore maximal dissipative. The same argumentation holds for f. \(\square \)

Before we continue with our analysis we recapitulate a few properties of a general maximal dissipative operator \(g:{D}(g)\subseteq H\rightarrow H\). The resolvent

is well defined, for every \(h>0\), and nonexpansive, i.e.,

The latter follows directly by the definition of dissipativity. Furthermore, the resolvent and the related Yosida approximation \(g(I - hg)^{-1}\) satisfies the following.

Lemma 6

If \(g:{D}(g)\subseteq H\rightarrow H\) is maximal dissipative, then

in H for every \(u\in \overline{{D}(g)}\) and \(v\in {D}(g)\), respectively.

The proof of Lemma 6 can, e.g., be found in [3, Proposition II. 3.6] or [7, Proposition 11.3]. Next, we will relate the full vector field f with its decomposition \(\sum _{\ell =1}^{s} f_{\ell }\).

Lemma 7

If the Assumptions 1–3 hold, then \(\bigcap _{\ell =1}^s {D}(f_{\ell }) \subseteq {D}(f)\) and \(fu = \sum _{\ell =1}^{s} f_{\ell }u\) for every \(u\in \bigcap _{\ell =1}^s {D}(f_{\ell })\).

Proof

Choose a \(u\in \bigcap _{\ell =1}^s {D}(f_{\ell })\), then \(u\in \bigcap _{\ell =1}^s V_{\ell } = V\) and the sum \(z=\sum _{\ell =1}^{s} f_{\ell }u\in H\) satisfies the relation

for all \(v\in V\). Hence, \(Fu\in H ^{*}\), which yields that \(u\in {D}(f)\) and \(fu=-\gamma ^{-1} Fu =z\). \(\square \)

Unfortunately, the set \({D}(f)\) is in general not equal to \(\bigcap _{\ell =1}^s {D}(f_{\ell })\), as \(u\in {D}(f)\) does not necessarily imply that \(F_{\ell } u \in H^*\) for every \(\ell =1,\dots ,s\). This issue is well known and we will encounter it when decomposing the p-Laplacian; compare with Sect. 6. We will therefore assume that the mild regularity property below holds.

Assumption 4

\(V \subseteq {R}\bigl ( I - h f |_{\bigcap _{\ell =1}^s {D}(f_{\ell })} \bigr )\quad \)for all \(h>0\).

Under this assumption one has the following identification, which is sufficient for our convergence analysis.

Lemma 8

If the Assumptions 1–4 hold, then the closure of \(f|_{\bigcap _{\ell =1}^s{D}(f_{\ell })}\) is f, i.e.,

Proof

By Lemma 7 and the fact that the maximal dissipative operator f is closed [3, Proposition II.3.4], we obtain that

Next, choose an arbitrary \((u,fu) \in {{\mathrm{graph}}}(f)\). Since

we can define \(v_h \in \bigcap _{\ell =1}^s{D}(f_{\ell })\) via

for every \(h>0\). By Lemma 6, we have the limits

Hence, the set \({{\mathrm{graph}}}\bigl ( f|_{\bigcap _{\ell =1}^s {D}(f_{\ell })} \bigr )\) is dense in \({{\mathrm{graph}}}(f)\), i.e., its closure in \(H\times H\) is equal to \({{\mathrm{graph}}}(f)\). \(\square \)

With the Friedrich formulation of our full vector field \(f:{D}(f) \subseteq H \rightarrow H\), the parabolic equations all take the form of an abstract evolution equations, i.e.,

on H. Furthermore, with the decomposition \(f=\sum _{\ell =1}^s f_{\ell }\), the splitting schemes (3) and (4) are given by the operators

respectively. Here, \(S^{n}_{h}\eta \) and \(P^{n}_{h}\eta \) are both approximations of the exact solution u at time \(t=nh\).

As the resolvent of a maximal dissipative operator is well defined and nonexpansive on H, it is a natural starting point for a solution concept. To this end, consider the operator family \(\{\mathrm {e}^{tf}\}_{t\ge 0}\) defined by

where the limit is well defined in H for every \(\eta \in \overline{{D}(f)}\) and \(t\ge 0\); see [6, Theorem I]. The operator family \(\{\mathrm {e}^{tf}\}_{t\ge 0}\) is in fact a (nonlinear) semigroup and each \(\mathrm {e}^{tf}:\overline{{D}(f)}\rightarrow \overline{{D}(f)}\) is a nonexpansive operator on H. The unique mild solution of the evolution Eq. (5) is then given by the function \(u:t\mapsto \mathrm {e}^{tf}\eta \), which is continuous on bounded time intervals. An extensive exposition of the nonlinear semigroup theory can, e.g., be found in [3].

There is a discrepancy between the domain of the solution operator, i.e., \({D}(\mathrm {e}^{tf})=\overline{{D}(f)}\), and the fact that the operators \(S_{h}\) and \(P_{h}\) are not necessarily invariant over it. In order to avoid several technicalities induced by this, we will assume the following.

Assumption 5

The domain \({D}(f)\) is dense in H.

As f is the closure of \(f|_{\bigcap _{\ell =1}^s {D}(f_{\ell })}\), one has the inclusions

which implies that \(\overline{{D}(f)}=\overline{{D}(f|_{\bigcap _{\ell =1}^s {D}(f_{\ell })})}\). Hence, \({D}(f|_{\bigcap _{\ell =1}^s {D}(f_{\ell })})\) is also dense in H when Assumption 5 holds.

We can now formulate the following simplified version of the Lax-type convergence result given in [5, Corollary 4.3].

Lemma 9

Consider an operator family \(\{G_{h}\}_{h>0}\), where each operator \(G_{h}:H\rightarrow H\) is nonexpansive on H and the operator family is consistent, i.e.,

If the Assumptions 1–5 hold, then

for every \(\eta \in H\) and \(t>0\).

Theorem 1

If the Assumptions 1–5 hold, then the sum splitting (3) is convergent in H to the mild solution of the abstract evolution Eq. (5), i.e.,

for every \(\eta \in H\) and \(t>0\).

Proof

As each resolvent \((I-hsf_{\ell })^{-1}\) is nonexpansive on H for all values of \(hs>0\), one has the bound

and \(S_{h}\) is therefore nonexpansive on H. To validate the consistency of \(\{S_{h}\}_{h>0}\), we first observe that

The consistency can then be formulated in terms of the Yosida approximation, i.e., for every \(u\in \cap _{\ell =1}^s {D}(f_{\ell })\) one has the limit

in H, as \(h\rightarrow 0\); compare with Lemma 6. The desired convergence is then proven as the hypotheses of Lemma 9 hold. \(\square \)

Theorem 2

If the Assumptions 1–5 hold, then the Lie splitting (4) is convergent in H to the mild solution of the abstract evolution Eq. (5), i.e.,

for every \(\eta \in H\) and \(t>0\).

Proof

We once more prove convergence by validating the hypotheses of Lemma 9. The nonexpansivity of the operator \(P_{h}\) on H follows trivially as every resolvent \((I-hf_{\ell })^{-1}\) has the same property. In order to validate the consistency of \(\{P_{h}\}_{h>0}\), let \(u\in \cap _{\ell =1}^s {D}(f_{\ell })\) and consider the telescopic expansion

where \(u_{1,h}=u\) and

As the arguments of the Yosida approximations in (6) are h dependent, we can not directly use Lemma 6. Instead, we assume for the time being that the limit

exists. By introducing the maximal dissipative operator

which satisfies \((I-hf_{\ell })^{-1}u_{\ell ,h}=(I-he_{\ell })^{-1}(u_{\ell ,h}+hz_{\ell })\), we have the reformulation

By Lemma 6 and the nonexpansivity of \((I-he_{\ell })^{-1}\), one then obtains the limit

Hence, if (7) exists then \( \lim _{h\rightarrow 0}f_{\ell }(I-hf_{\ell })^{-1}u_{\ell ,h}=f_{\ell }u\). Furthermore, if (7) exists for every \(\ell =1,\ldots ,s\), then \(\lim _{h\rightarrow 0}1/h\,(P_{h}-I)u=fu\) in H.

The limit (7) obviously exists for \(\ell =1\). If it exists for \(\ell =k\) then it also exists for \(\ell =k+1\), as

in H, as \(h\rightarrow 0\). By induction, the limit (7) exists for every \(\ell =1,\ldots ,s\), and \(\{P_{h}\}_{h>0}\) is therefore consistent.\(\square \)

Remark 2

The results can be extended to perturbed equations \(\dot{u}=(f+g)u\), e.g., arising if a lower-order advection or reaction term is added to the diffusion process. Here, g and \(f+g\) are both assumed to satisfy a shifted dissipativity condition of the form

with M being a nonnegative constant, and the range condition \({R}(I-hg)=H\) for \(h\in (0,1/M)\). This is, e.g., satisfied when \(g:H\rightarrow H\) is Lipschitz continuous. More elaborate perturbation examples are given in [3, Section II.3.2]. For these perturbed evolution equations, one has convergence for the modified splitting schemes, with a single step given by

respectively. If \(g:H\rightarrow H\) is in addition Lipschitz continuous, then convergence is also obtained for the semi-implicit schemes

The convergence of the modified schemes follow just as for the proof of Theorem 2 together with the fact that [5, Corollary 4.3] is valid for operators \(G_{h}\) that have Lipschitz constants of the form \(1+Ch\).

As a first problem class we consider the parabolic equations of p-Laplace type with homogeneous Neumann boundary conditions, i.e.,

The domain \(\varOmega \subset \mathbb {R}^d\) is assumed to have a locally Lipschitz boundary \( \partial \varOmega \), and the map \(\alpha :\varOmega \times \mathbb {R}^{d}\rightarrow \mathbb {R}^{d}\) satisfies Assumption 3 for a given \(p\ge 2\). The classical p-Laplacian is then given by

After multiplication with v and a subsequent integration by parts, the variational form of (8) and its decomposition is formally given by

Here, we have introduce a domain decomposition \(\{\varOmega _{\ell }\}_{\ell =1}^s\), where \(\bigcup _{\ell =1}^s \varOmega _{\ell }=\varOmega \), together with a partition of unity \(\{\chi _{\ell }\}_{\ell =1}^{s}\) chosen in \(W^{1,\infty }(\varOmega )\); compare with Remark 1.

In order to fit the variational form into the abstract setting of Sect. 3, we choose the pivot space \(H = L^2(\varOmega )\) and the operator \(\delta \) as the distributional gradient

This choice of \(\delta \) fulfills the continuity Assumption 1, since for a convergent sequence \(\{u_n\}_{n\in \mathbb {N}}\subset L^2(\varOmega )\) and an arbitrary \(\varphi \in C^{\infty }_0(\varOmega )\) one can write

for every\(j=1,\dots ,d\), where \(D_j\) is the j-th partial derivative in a distributional sense. The space V is then

A bootstrap argument using the Sobolev embedding theorem yields the identification \(V = W^{1,p}(\varOmega )\). Since \(W^{1,p}(\varOmega )\) is dense in \(L^2(\varOmega )\), Assumption 2 is also fulfilled.

With these choices, \(\delta _{p}u\) is simply the weak gradient of \(u\in W^{1,p}(\varOmega )\) and we obtain the standard energetic form \(F :V \rightarrow V^*\) of p-Laplace type vector fields, i.e.,

The domain of the corresponding Friedrich extension can be written as

and fu is given by the weak divergence of \(\alpha (\nabla u)\). The same characterization can be made for \(F_{\ell }\) and \(f_{\ell }\), respectively. Applying Lemma 5 the operators f and \(f_{\ell }\), \(\ell =1,\dots ,s\), are maximal dissipative and Lemma 7 yields that

Validation of Assumption 5 requires further structure of the map \(\alpha \). For the classical p-Laplacian the related \(\alpha \) is continuously differentiable and \(\alpha (0)=0\), which implies that \(C_0^{\infty }(\varOmega )\) is a subset of \({D}(f)\). Hence, \({D}(f)\) is dense in \(L^2(\varOmega )\) and Assumption 5 is valid in this context. Finally, if Assumption 4 holds then the convergence results from Sect. 5 can directly be applied.

An example of a domain decomposition \(\{\varOmega _{\ell }\}_{\ell =1}^3\) that fulfills (10)

Full size image

Apart from the special cases when \(d=1\) or \(p=2\), the domains \({D}(f)\) of p-Laplace type vector fields can not be expected to coincide with \(\bigcap _{\ell =1}^{s} {D}(f_{\ell })\). The issue is that for an element \(u \in {D}(f)\) one has

where the function \(\alpha (\nabla u)\) only lies in \(L^{p/(p-1)}(\varOmega )^{d}\), with \(p>2\). The term \(f_{\ell }u\) is therefore, in general, not an \(L^{2}(\varOmega )\) function. In order to give a possible setting for which Assumption 4 is valid, we assume that the domain decomposition \(\{\varOmega _{\ell }\}_{\ell =1}^s\) is chosen such that

That is, the subdomain \(\varOmega _{s}\) separates the boundary \(\partial \varOmega \) from the other subdomains; as illustrated in Fig. 2.

Lemma 10

Consider a domain decomposition \(\{\varOmega _{\ell }\}_{\ell =1}^s\) that satisfies (10) and with subdomains \(\varOmega _{\ell }\), \(\ell =1,\ldots ,s-1\), that all have the segment property. If \(p\ge 2\) in addition satisfies \(p>(d+1)/2\) and the map \(\alpha \) fulfills Assumption 3(\(\alpha _{2}\)) with \(c_{2}\in L^{2}(\varOmega )\), then the Friedrich extension f of a p-Laplace type vector field and its decomposition into the operators \(f_{\ell }\) fulfill Assumption 4.

Proof

For an arbitrary \(g \in V=W^{1,p}(\varOmega )\) there exists a unique \(u\in {D}(f)\) such that \(u-hfu = g\) and Assumption 4 is then valid if \(u\in \bigcup _{\ell =1}^s{D}(f_{\ell })\). To prove this, we first observe that \(fu=\nabla \cdot \alpha (\nabla u) = (u-g)/h \in W^{1,p}(\varOmega )\) and \(W^{1,p}(\varOmega )\hookrightarrow L^r(\varOmega )\) for some \(r > dp/(p-1)\), as \(p \ge 2\) and \(p > (d+1)/2\). Hence, [8, Theorem 2 and Remarks pp. 829–830] implies that \(\nabla u\) is locally Hölder continuous on \(\varOmega \) and we obtain that

for every open domain \(\varOmega _{int}\) such that \(\overline{\varOmega }_{int}\subset \varOmega \).

As \(u\in {D}(f)\), we have the integration by parts

for every \(w\in W^{1,p}(\varOmega )\). Due to the extra interior regularity of \(\alpha (\nabla u)\) we can, e.g., extend (11) to all \(w=w_{1}+w_{2}\), where \(w_{1}\in W^{1,p}(\varOmega )\) and \(w_{2}\in H^{1}(\varOmega )\) is a.e. zero on \(\varOmega \setminus \varOmega _{int}\) for some open subdomain \(\varOmega _{int}\) that has the segment property and fulfills \(\overline{\varOmega }_{int}\subset \varOmega \). The latter implies that \(w_{2}\) is the zero extension of \(w_{2}|_{\varOmega _{int}}\in H^{1}_0 (\varOmega _{int})\); see, e.g., [1, Theorem 5.29].

Next, let \(v\in V_{\ell }\subset L^{2} (\varOmega )\), for \(\ell =1,\ldots ,s\), and consider \(\chi _{\ell } v \in L^2(\varOmega )\). Here,

for every \(\varphi \in C^\infty _0 (\varOmega )\), i.e., \(\chi _{\ell } v \in H^1(\varOmega )\) and \(\chi _{\ell } v = 0\) a.e. on \(\varOmega \setminus \varOmega _{\ell }\). If \(\ell <s\) then \(\chi _{\ell } v|_{\varOmega _{\ell }}\in H^1_{0}(\varOmega _{\ell })\).

For \(\ell =1,\ldots ,s-1\), we can test with \(w=\chi _{\ell }v\) and integrate by parts (11). Writing out \(\nabla (\chi _{\ell }v)\) and rearranging the terms gives us

i.e., \(u\in \bigcap _{\ell =1}^{s-1} {D}(f_{\ell })\), as the integrand on the right-hand side is in \(L^{2}(\varOmega )\).

It remains to prove that u lies in \({D}(f_{s})\). As the closure of \(\bigcup _{\ell =1}^{s-1}\varOmega _{\ell }\) does not intersect the outer boundary \(\partial \varOmega \), we can choose an open subset \(\varOmega _{out}\subset \varOmega _{s}\) such that \(\chi _{s}=1\) on \(\varOmega _{out}\), its boundary \(\partial \varOmega _{out}\) is locally Lipschitz continuous and \(\partial \varOmega \subset \partial \varOmega _{s}\). Let \(v\in V_{s}\), then \(\chi _{s}\delta _{p,s}v=\nabla v\) a.e. on \(\varOmega _{out}\) and \(\chi _{s}v|_{\varOmega _{out}}=v|_{\varOmega _{out}}\in W^{1,p}(\varOmega _{out})\). The local Lipschitz continuity of \(\partial \varOmega _{out}\) implies, e.g., via [1, Theorem 5.24], that there exists an extension \(w_{1}\in W^{1,p}(\varOmega )\) such that \(w_{1}=\chi _{s}v\) a.e. on \(\varOmega _{out}\). Furthermore, \(w_{2}=\chi _{s}v-w_{1}\in H^{1}(\varOmega )\) is zero a.e. on \(\varOmega _{out}\), i.e., it is a zero extension of an \(H^{1}_0 (\varOmega _{int})\) function on some subdomain \(\varOmega _{int}\), with \(\partial \varOmega _{int}\subset \varOmega _{out}\). For every \(v\in V_{s}\) we therefore have a partitioning of the form \(w=\chi _{s}v=w_{1}+w_{2}\) and the integration by parts (11) is well defined for \(\ell =s\). By the same argumentation as for \(\ell <s\), one obtains that u lies in \({D}(f_{s})\). \(\square \)

Remark 3

From a numerical perspective the construction (10) with a separating subdomain \(\varOmega _{s}\) is suboptimal for general time dependent PDEs, as it may increase the amount of communication in the implementation of scheme. However, as discussed in Sect. 1, we are foremost interested in the approximation of solutions with compact support in \(\varOmega \). Hence, for sufficiently short time intervals (0, T) there is obviously no communication related to \(\varOmega _{s}\); as exemplified in Fig. 2.

A second problem class that fits into our abstract setting is the parabolic equations of porous medium type with homogeneous Dirichlet boundary conditions, i.e.,

Here, the domain \(\varOmega \subset \mathbb {R}^d\) is assumed to have a locally Lipschitz boundary \( \partial \varOmega \), and the map \(\alpha :\varOmega \times \mathbb {R}\rightarrow \mathbb {R}\) fulfills Assumption 3 for a given p that satisfies

This restriction on p is made in order to assure the embedding

which is central in our forthcoming analysis. The standard porous medium equation is then given by

and the fast diffusion equation is obtained for the same \(\alpha \), but with \(1<p<2\); see [25]. The two-phase Stefan problem [11, Section 5.10] follows by choosing

where \(a,b>0\), and Assumption 3 is then valid for \(p=2\).

After multiplying (12) by w, where \(-\varDelta w = v\) in \(\varOmega \) and \(w=0\) on \(\partial \varOmega \), and integrating by parts twice, the variational form of (12) and its decomposition is formally

Above, we have once more introduced a domain decomposition \(\{\varOmega _{\ell }\}_{\ell =1}^s\) of \(\varOmega \) together with a partition of unity \(\{\chi _{\ell }\}_{\ell =1}^{s}\).

With the proper interpretation, the left-hand side of (14) is given by the inner product on \(H^{-1}(\varOmega )\); compare with [14, Bemerkung III.1.13]. The formal variational formulation (14) therefore leads us to choosing the pivot space \(H=H^{-1}(\varOmega )\) and the operator

The operator \(\delta \) obviously fulfills the continuity Assumption 1. The space V is now

and as before \(\delta _p u = v\), where v is the unique function stated in the definition of V. By the embedding (13) and [14, Bemerkung I.5.14], we obtain that

i.e., Assumption 2 is fulfilled. With these choices, we have the energetic form \(F :V \rightarrow V^*\) given by

In order to characterize the Friedrich operator f, we introduce the Dirichlet Laplacian \(-\varDelta :H^{1}_{0}(\varOmega )\rightarrow H^{-1}(\varOmega )\), where

As \(-\varDelta \) is the Riesz isomorphism from \(H^{1}_{0}(\varOmega )\) to \(H^{-1}(\varOmega )\), the inner product on \(H^{-1}(\varOmega )\) satisfies

for all \(u,v\in H^{-1}(\varOmega )\); compare with [10]. Next, for \(u\in {D}(f)\) there exists a \(z\in H^{-1}(\varOmega )\) such that

for all \(v \in \bigl (L^{p/(p-1)}(\varOmega )\bigr )^*\), or equivalently

Hence, \(-\,\alpha (\delta _{p} u)=(-\varDelta )^{-1}z\in H^{1}_{0}(\varOmega )\); see, e.g., [1, Lemma 3.31], and we obtain the characterization

and \(fu=\varDelta \alpha (\delta _{p} u)\) for \(u\in {D}(f)\).

Analogously to Sect. 6, we have \({R}(\delta )= H^{-1}(\varOmega ) \subset \mathscr {D}'(\varOmega )\) and we can therefore allow a partition of unity \(\{\chi _{\ell }\}\) in \(W^{1,\infty }(\varOmega )\). The spaces \(V_{ \ell }\), \(\ell =1,\dots ,s\), are then

Again, we write \(\delta _{p, \ell } u\) for the unique \(L^p(\varOmega _{\ell },\chi _{\ell })\) function v from this definition.

After introducing \(F_{\ell }\) and \(f_{\ell }\), as described in Sect. 3, we have by Lemmas 5 and 7 that the operators f and \(f_{\ell }\), \(\ell =1,\dots ,s\), are maximal dissipative and

Instead of Assumption 4 we can prove the stronger condition

To prove the equality take an arbitrary \(u\in {D}(f)\). Since \(\alpha (\delta _p u )\in H_0^1(\varOmega )\), we also have that \(\chi _{\ell }\alpha (\delta _p u ) \in H_0^1(\varOmega )\) for every weight function \(\chi _{\ell } \in W^{1,\infty }(\varOmega )\) and

That is, u also lies in \({D}(f_{\ell })\) for \(\ell =1,\dots ,s\).

Assumption 5 requires some further regularity of the map \(\alpha \) and the validation that \(\alpha (\delta _{p} u)\) vanishes on the boundary \(\partial \varOmega \). For the porous medium equation and the two-phase Stefan problem one has that \(\alpha (\varphi )\in H^{1}_{0}(\varOmega )\) for every \(\varphi \in C_0^{\infty }(\varOmega )\). The set of functionals of the form \(v\mapsto \int _{\varOmega }u v\,\mathrm {d}x\), where \(u\in C_0^{\infty }(\varOmega )\) and \(v\in H^{1}_{0}(\varOmega )\), is therefore a subset of \({D}(f)\). It is also a dense subset of \(H^{-1}(\varOmega )\), as \(C_0^{\infty }(\varOmega )\) is dense in \(L^{2}(\varOmega )\) and \(L^{2}(\varOmega )^*\) is dense in \(H^{-1}(\varOmega )\). Hence, Assumption 5 is valid for these two prototypical examples, and the convergence results of Sect. 5 hold.

Remark 4

The variational setting of porous medium type equations, with \(H^{-1}(\varOmega )\) as pivot space, is by no means standard. However, it enables a clear-cut way of introducing the related Friedrich operator. The variational setting has, e.g., been proposed in [14, Bemerkung I.5.14]. It has also been employed in [10] when proving convergence of finite element/implicit Euler approximations for the porous medium equation, on its very weak form. Note that the standard approach to prove that \(\varDelta \alpha \) is a maximal dissipative operator on \(H^{-1}(\varOmega )\) is to directly observe that it is the gradient of a convex function; see [4, Example 3].

We conclude by illustrating the convergence of the sum (3) and Lie (4) splitting schemes. In order to obtain a numerical example where the spatial error is negligible, we will consider a nonlinear parabolic equation with a known solution and a fine spatial grid. To this end, we choose the classical one-dimensional parabolic p-Laplace equation, i.e., Eq. (8) with

If the initial value is chosen as the Dirac delta then the solution of the parabolic p-Laplace equation in \(\mathbb {R}\times (0,T)\) is given by the Barenblatt solution

where \(\lambda =1/(2p-2)\), \(\kappa =\lambda ^{1/(p-1)}(p-2)/p\) and \([x]_{+}=\max \{0,x\}\). As we require the initial value \(\eta \) in (8) to be an \(L^{2}(\varOmega )\)-element, we choose \(\eta (x)=B_{p}(x,t_{0})\) for a \(t_{0}>0\). The solution is then given by \(u(x,t)=B_{p}(x,t+t_{0})\) for \(t\in (0,T)\), where the end time T is assumed to satisfy the condition \({\mathrm {supp}}(u(T))\subset \varOmega \).

Next, observe that one time step of either splitting scheme consists of solving equations of the form \( v_{\ell }=z+h f_{\ell } v_{\ell }\), or equivalently

A possible finite element discretization is obtained by equidistantly partitioning \(\bar{\varOmega }\) into M subintervals \([-L,-L+\varDelta x], (-L+\varDelta x,-L+2\varDelta x],\ldots ,(L-\varDelta x,L]\), with \(\varDelta x=2L/M\), and replacing \(V_{\ell }\) with the finite dimensional space \(S_{M}\), consisting of continuous functions on \(\bar{\varOmega }\) that are linear in each subinterval. Let \(\{\varphi _{j}\}\subset S_{M}\) denote the standard basis of hat functions and assume that \(\bar{\varOmega }_{\ell }\) is the union of the m subintervals \([x_{0},x_{1}],(x_{1},x_{2}],\ldots ,(x_{m-1},x_{m}]\). For \(z_{M}\in S_{M}\) being a known approximation of the element z from (15), with \(z_{M}=\sum _{j=0}^{m} d_{j}\varphi _{j}\) on \(\bar{\varOmega }_{\ell }\), the finite element approximation of \(v_{\ell }\) then has the form \(v_{\ell ,M}=\sum _{j=0}^{m} c_{j}\varphi _{j}\) on \(\bar{\varOmega }_{\ell }\) and \(v_{\ell ,M}=z_{M}\) otherwise. By employing mass lumping [23, Chapter 15], the unknown coefficients \(\{c_{j}\}_{j=0}^{m}\) satisfy the algebraic equation system

where \(\chi _{j+1/2}=1/\varDelta x\int ^{x_{j+1}}_{x_{j}}\chi _{\ell }\,\mathrm {d}x\). By utilizing Newton’s method for the algebraic equation systems, we obtain computable space time discretizations based on domaindecomposition splittings.

For our numerical experiments we set \(L=3\), \(t_{0}=0.1\), \(T=5\) and decompose \(\varOmega \) into eight consecutive subdomains \(\varOmega _{A},\ldots ,\varOmega _{H}\), all of the same width and with equally sized overlaps. We also introduce a fine equidistant partitioning of \(\bar{\varOmega }\) with \(M=8700\) subintervals, such that each of the seven overlaps contain 100 subintervals. In the implementation of the Lie splitting we set \(s=2\),

When the same domain decomposition is employed for the sum splitting, we refer to it as the Sum2 scheme. In order to study the influence of the number of operators \(f_{\ell }\) in the splitting schemes, we also consider the sum splitting with a \(s=4\) decomposition, where

This is referred to as the Sum4 scheme. In all cases, the weight functions \(\{\chi _{\ell }\}_{\ell =1}^{s}\subset W^{1,\infty }(\varOmega )\) are chosen to be continuous and piecewise linear. We finally estimate the full (space time) error at time \(t=5\) in the norm of \(L^{2}(\varOmega )\), by taking the difference between the numerical approximation and the linear interpolation of the exact solution.

The \(L^{2}\)-errors at time \(t=5\) with respect to the time step h for various time integrators and \(p=6\) (left) and for the Sum2 splitting scheme and various p values (right), respectively

Full size image

The \(L^{2}\)-errors for various time steps h and \(p=6\) are presented in Fig. 3 for the Lie, Sum2 and Sum4 schemes. As a reference, we also give the errors for the implicit Euler (IE) scheme. From the results it is clear that the approximation schemes converge and the temporal errors are dominating the full error. If the latter was not the case, the error curves would level out for small time steps h. In this somewhat academic example, we even obtain the classical first-order convergence for all the schemes. This is most likely due to the regularity of the solution, which is smooth away from the interface between \(u=0\) and \(u>0\). Note that such regularity is not to be expected in general and neither is first-order convergence.

Comparing the results of Sum2 and Sum4 reveals that the error constant dubbles when passing from \(s=2\) to \(s=4\) operators, which comes as no surprise as the time step h is scaled by the number of operators in the sum splitting (3). Hence, the number of operators should be kept at a minimum to assure a moderate error constant. Note that for complex geometries in \(d>1\) spatial dimensions, one is typically forced to use \(s>2\) operators in order to obtain a reasonable spatial accuracy.

We conclude by investigating the influence of the parameter p on the convergence of the sum splitting schemes (similar results hold for the Lie splitting). In Fig. 3 the convergence for the Sum2 method is illustrated for various time steps h and \(p=6,9,30\). As seen from the results, the convergence rate does not seem to be affected by the parameter p. The only observed impact on the splitting integrators is a mild time step restriction for extreme p values. This is in stark contrast to classical domain decomposition schemes applied to the elliptic p-Laplace equation. Here, the the convergence rate, with respect to the number of iterations, is predicted to deteriorate for large values of p [22, Section 5].