A variational approach to the sum splitting scheme
Monika Eisenmann Eskil Hansen


A variational approach to the sum splitting scheme]A variational approach 
to the sum splitting scheme

M. Eisenmann]Monika Eisenmann
Monika Eisenmann

  Centre for Mathematical Sciences

  Lund University

  P.O. Box 118

  221 00 Lund, Sweden
monika.eisenmann@math.lth.se

E. Hansen]Eskil Hansen
Eskil Hansen

  Centre for Mathematical Sciences

  Lund University

  P.O. Box 118

  221 00 Lund, Sweden
eskil.hansen@math.lth.se



	

	Nonlinear parabolic equations are frequently encountered in applications and 
	efficient approximating techniques for their solution are of great importance. In 
	order to provide an effective scheme for the temporal approximation of such 
	equations, we present a sum splitting scheme that comes with a straight forward 
	parallelization strategy. 
  The convergence analysis is carried out in a variational framework that allows for 
  a general setting and, in particular, nontrivial temporal coefficients. 
  The aim of this work is to  illustrate the significant advantages of a variational 
  framework for operator splittings and use this to extend semigroup based theory 
  for this type of  schemes. 




[
    [
    December 30, 2023
=====================







§ INTRODUCTION


Nonlinear parabolic equations, which we state as abstract evolution equation of 
the form 

    u'(t)+A(t)u(t)=f(t),   t∈(0,T)  and   u(0)=u_0,

are frequently encountered in applications appearing in physics, chemistry and 
biology; see <cit.> and <cit.>. A few 
standard examples of the diffusion operator A(t)v are 

    -∇·(α(t) |∇ v|^p-2∇ v), 
        
      -Δ(α(t) |v|^p-2v)
        and  
      -∑_i=1^d D_i(α(t) |D_i v|^p-2 D_i v).

Here, the first and second operator is referred to as the p-Laplacian and the 
porous medium operator, respectively. 

Due to the problems' significance, effective techniques for their approximations 
become crucial. As we consider parabolic equations, for stability reasons the 
temporal approximation schemes need to be implicit. 
For equations which in addition are given in several 
spatial dimensions the resulting spatial and temporal approximation schemes 
require large scale computations. This typically demands implementations in 
parallel on a distributed hardware. One possibility to design temporal 
approximation schemes that can directly be implemented in a parallel fashion is 
to utilize operator splitting; see, e.g., <cit.> for an 
introduction. 
Note that the solutions of nonlinear parabolic problems typically lack high-order 
spatial and temporal regularity. Thus, there is little use to consider high-order 
time integrators. 

In order to illustrate the splitting concept, consider the simplest implicit scheme, 
namely, the backward Euler method. For N temporal steps, a step size
k=T/N and the starting value ^0=u^0 the backward Euler 
approximation ^n of u(nk) is given by the recursion

    1/k(^n-^n-1)+k^n^n= 
      ^n,   n∈{1,…, N},

where (^n)_n and (^n)_n are suitable 
approximations for A and  f, respectively. Assuming that the nonlinear 
resolvent of ^n exists, we find the reformulation 

    ^n=(I+k^n)^-1(^n-1+k^n),  
       n∈{1,…, N}.

To implement one step of the backward Euler scheme in parallel, we split the Euler 
step into s independently solvable problems. To this end, we decompose 
^n and ^n as 

    ^n=∑_ℓ=1^s^n_ℓ  and  ^n=∑_ℓ=1^s^n_ℓ,
         n∈{1,…, N}.

With this abstract operator splitting, one can design various temporal 
approximation schemes. Two possibilities to split one single Euler step are given 
by formally multiplying or adding the operators (I+k^n_ℓ)^-1, ℓ∈{1,…,s}. 
A composition of such fractional step operators yields the Lie splitting scheme

    (I+k^n)^-1≈∏_ℓ=1^s(I+k^n_ℓ)^-1.

Thus, we obtain s possibly easier subproblems that are solved after each other. 
For a straightforward parallelization it is more convenient to choose a splitting, 
where the single steps can be computed at the same time. The sum splitting

    (I+k^n)^-1≈1/s∑_ℓ=1^s(I+k s 
      ^n_ℓ)^-1

offers this crucial advantage. The s fractional steps are solved at the same time 
and their average is used as an approximation. 

The decomposition (<ref>) can be utilized in many different ways. A 
first possible application is a source term splitting, where the high-order terms 
are split from the low-order terms. For example a source term splitting of the 
reaction-diffusion equation governed by A(t)v= 
-∇·(α(t)∇ v) +p(t,v) would have the form 

    _1^nv= -∇·(α(nk)∇ v) 
        and  _2^nv=p(nk,v).

Here, the actions of (I+k^n_1)^-1 can be evaluated by a standard 
fast linear elliptic solver and the actions of the nonlinear resolvent 
(I+k^n_2)^-1 can often be expressed in a closed form. Examples 
of studies dealing with various source term splittings can be found in 
<cit.>. 

Another possibility is a dimension splitting, where each spatial derivative is 
considered as a separate differential operator. For example, the dimension 
splitting of the nonlinear porous medium operator and the third operator in 
(<ref>) are given by

    _ℓ^nv=-D_ℓℓ(α(nk) |v|^p-2v)
        and  _ℓ^nv=-D_ℓ(α(nk) |D_ℓ v|^p-2 D_ℓ v),

respectively. This splitting yields that the action of each nonlinear resolvent 
(I+k^n_ℓ)^-1
can be separated into lower-dimensional subproblems that can be solved on their 
own. Note that the p-Laplacian lacks a natural dimension splitting. Examples of 
convergence results for the dimension splitting of the third equation in 
(<ref>) can be found in <cit.>, where the Lie 
scheme is used, and in <cit.>, where the sum, Lie and 
Peaceman–Rachford schemes are considered for the autonomous case.

A limitation of the dimension splitting approach is the rather large need of 
communication between the subproblems, which can impede an effective 
distributed implementation. Dimension splitting is also quite restrictive in terms 
of the spatial domains that can be considered.
A modern alternative to dimension splitting, which is applicable to a very general 
class of spatial domains, is the domain decomposition based splitting. 
Here, the subproblems are given on s spatial subdomains that share a small 
overlap. As an example consider the three nonlinear diffusion 
operators (<ref>) and introduce 
a partition of unity (χ_ℓ)_ℓ=1^s, where each weight function 
χ_ℓ vanishes outside its corresponding spatial subdomain. The domain 
decompositions _ℓ^nv are then 

    -∇·(χ_ℓα(nk) |∇ v|^p-2∇ v),
        
      -Δ(χ_ℓα(nk) |v|^p-2 v)
      and  
      -∑_i=1^dD_i(χ_ℓα(nk) |D_i v|^p-2 D_i v),

respectively. This approach is well suited for a parallel computation, as the 
actions of (I+k^n_ℓ)^-1 can be solved independently of each 
other and the communication required is small, due to the small overlaps between 
the subdomains. Studies regarding domain decomposition based splittings 
applied to linear and autonomous parabolic equations include 
<cit.>. Convergence for the Lie and sum splittings are given in 
<cit.> for the autonomous p-Laplace and porous 
medium equations. 

Operator splitting schemes are typically analyzed in a semigroup framework, 
which yields convergence for a wide range of temporal 
approximation schemes, including the Lie and sum schemes; see 
<cit.> for more details on the solution concept. However, there does 
not seem to be a straightforward way to extend the semigroup based convergence 
analysis to nonautonomous evolution equations. 
Furthermore, the semigroup framework requires some additional regularity 
conditions to relate the intersection of the domains D(A_ℓ), ℓ∈{1,…,s}, with the domain D(A) of the full operator. The latter, e.g., implies 
restrictions on the domain decomposition of the p-Laplace equation 
<cit.>.

In a variational setting this problem is avoided in a natural way while at 
the same time the analysis of nonautonomous problems is accessible. Also the 
structure of this approach is well suited to include a Galerkin scheme and 
therefore, in particular, the finite element method. 
However, the analysis typically needs to be tailored for each method.
The variational setting is a standard tool for existence theories 
<cit.> and has been used 
in several works in the context of “unsplit” time integrators 
<cit.>. However, in the context of temporal splitting 
schemes for nonlinear parabolic equations the only variational studies that we are 
aware of is <cit.>. 
Here, a variational analysis is employed when proving the convergence of the Lie 
scheme applied to nonautonomous evolution equations and, as already stated, is 
applied to the dimension splitting of the third equation in (<ref>). 

Hence, the aim of this paper is threefold. 
Firstly, we aim to generalize the previous semigroup based analysis for the sum 
scheme to nonautonomous evolution equations without any implicit regularity 
assumptions. The latter generalization will be applicable to splittings of 
reaction-diffusion, dimension and domain decomposition type. 
Secondly, we intend to extend the abstract variational convergence results for the 
Lie scheme to the sum splitting scheme. As this requires a tailored convergence 
proof, it is not a trivial implication.
Thirdly, we also strive to illustrate the advantages of a variational approach in the 
context of splitting analyses. 

This paper is organized as follows: In Section <ref>, we state the exact 
assumptions that are needed on the abstract variational framework considered in 
the paper. This section also contains an example that shows that the relevant 
application of domain decomposition integrators for the p-Laplacian operator 
fits into our abstract framework. This in mind, we prove the well-posedness of the 
sum scheme, as well as suitable a priori bounds in Section <ref>. 
The main convergence results are proven in Section <ref>; see 
Theorem <ref> and Theorem <ref>.



§ ABSTRACT SETTING
 

In this section, we introduce an abstract setting for the convergence 
analysis of the sum splitting scheme. We begin by presenting the exact 
assumptions made on the data and present the temporal discretization of 
the problem. This at hand, we can state the scheme that we will work with in this 
paper. The section ends with a more concrete setting that exemplifies the abstract 
framework.


	Let (H,[H]··,·_H) be a real, separable Hilbert space 
	and let (V, ·_V) be a real, separable, reflexive Banach space such 
	that V is continuously and densely embedded into H.  
	Further, there exist a seminorm |·|_V on V and c_V ∈ 
	(0,∞) such that ·_V ≤ c_V ( ·_H + 
	|·|_V) is fulfilled.
	
	Furthermore, for s ∈ let (V_ℓ, ·_V_ℓ), ℓ∈{1,
	…, s}, be real reflexive Banach spaces that are continuously and
	densely embedded into H, fulfill ⋂_ℓ = 1^s V_ℓ = V 
	and ∑_ℓ = 1^s·_V_ℓ is equivalent to ·_V.
	For every ℓ∈{1,…,s}, there exists a seminorm 
	|·|_V_ℓ and c_V_ℓ∈  (0,∞) such that 
	·_V_ℓ≤ c_V_ℓ( ·_H + 
	|·|_V_ℓ) and ∑_ℓ  = 1^s|·|_V_ℓ is 
	equivalent to |·|_V.


Identifying H with its dual space H^*, we obtain the Gelfand triples

    V d H ≅ H^* d  V^*
    	  and  
    	V_ℓd H ≅ H^* d V^*_ℓ,
    	  ℓ∈{1,…,s}.

The next assumption states the properties of the differential operator that 
are of importance.



	Let H and V be given as stated in Assumption <ref>.
	Furthermore, for T>0 and p > 1, let {A(t)}_t ∈ [0,T] be a 
	family of operators such that A(t)   V → V^* satisfy the 
	following conditions:
	
		
  (1) The mapping Av  [0,T] → V^*, v ∈ V, given by t 
		↦ A(t)v is continuous.
		
  (2) The operator A(t)  V → V^*, t ∈ [0,T], is
		radially continuous, i.e., the mapping τ↦[V]A(t)(u+τ v)w is continuous on [0,1] for u,v,w∈ V.
		
  (3) The operator A(t)  V → V^*, t ∈ [0,T], fulfills a 
		monotonicity condition such that  there exists η≥ 0 with
		
    [V]A(t)v - A(t)wv - w≥η |v-w|_V^p,    v,w ∈ V.

		
  (4) The operator A(t)  V → V^*, t ∈ [0,T], is uniformly
		bounded such that there exists β >0 with
		
    A(t) v _V^*≤β(1 + v_V^p-1),    v ∈ V.

		
  (5) The operator A(t)  V → V^*, t ∈ [0,T], fulfills a 
		coercivity condition such that there exist μ >0 and λ≥ 0 
		with
		
    [V]A(t) vv + λ≥μ |v|_V^p,    v ∈ V.

	


Now, we can combine Assumption <ref> and Assumption <ref> to state 
a decomposition of the operator family {A(t)}_t ∈ [0,T] that we employ in 
the analysis of the sum splitting scheme.


	For s ∈ let H, V and V_ℓ, ℓ∈{1,…,s},
	fulfill Assumption <ref>. For p > 1 and T >0  let the 
	operator family {A(t)}_t∈ [0,T] be given such that it fulfills
	Assumption <ref>. 
	Further, let {A_ℓ(t)}_t∈ [0,T], ℓ∈{1,…,s}, be given 
	such that A_ℓ(t)  V_ℓ→  V_ℓ^* fulfills  
	Assumption <ref>, with V replaced by V_ℓ for every ℓ∈{1,…,s}. Moreover, let the sum property
	
    ∑_ℓ = 1^s A_ℓ(t)v = A(t)v   in  V^*,    t ∈ 
    		[0,T], v ∈ V

	be fulfilled.



	Note that the optimal coefficients β, η, λ, μ for the operator 
	families {A(t)}_t ∈ [0,T], {A_ℓ(t)}_t∈ [0,T], ℓ∈{1,…, 
	s}, of operators do not necessarily have to be the same. For the sake of 
	simplicity, we assume that these coefficients coincide.


We also consider the differential operators of Assumption <ref> as 
Nemytskii operators acting on spaces of Bochner integrable functions. For an 
introduction to Bochner integrable functions we refer the reader to 
<cit.> or <cit.>. Some 
properties of such Nemytskii operators are collected in the next lemma. The 
proofs can be found in <cit.>.

 
	For p > 1, q = p/p-1 and T>0 let {A(t)}_t∈[0,T] 
	fulfill Assumption <ref>.
	Then the operator (A v)(t) = A(t)v(t) maps L^p(0,T;V) into L^q(0,T;V^*).
	The operator is radially continuous, i.e., the mapping τ↦⟨ 
	A(u+τ v), w ⟩_L^q(0,T;V^*) × L^p(0,T;V)
	is continuous on [0,1] for all u,v,w ∈ L^p(0,T;V).  Furthermore, it 
	fulfills a monotonicity, a boundedness and a coercivity condition such that
	
    ⟨ A v - A w, v -w ⟩_L^q(0,T;V^*) × L^p(0,T;V)≥η∫_0^T |v(t) - w(t) |_V^p t,
       A v _L^q(0,T;V^*)≤β(T^1/q + v_L^p(0,T;V)^p-1),
       ⟨ A v, v ⟩_L^q(0,T;V^*) × L^p(0,T;V) + λ T
    		≥μ∫_0^T |v(t)|_V^p t

	for all v,w ∈ L^p(0,T;V).


The Nemytskii operator of {A_ℓ(t)}_t∈[0,T], ℓ∈{1,…,s}, as introduced in Assumption <ref> also fulfills the same 
bounds with V replaced by V_ℓ. To make our setting complete, it remains to 
state the assumptions on f.


	Let V and V_ℓ, ℓ∈{1,…,s}, fulfill
	Assumption <ref>. Let p be the same value as in 
	Assumption <ref>  and q = p/p-1. Further, let f be in 
	L^q(0,T;V^*). Assume that there exist functions f_ℓ∈ 
	L^q(0,T;V^*_ℓ), ℓ∈{1,…,s}, such that
	
    ∑_ℓ = 1^sf_ℓ(t)= f(t)  in  V^*
    		  and  f_ℓ (t) _V_ℓ^*≤ f(t) _V^*,    a.e.  t ∈ 
    		(0,T).


Note that this assumption can be generalized to functions f ∈ L^q(0,T;V^*) + 
L^2(0,T; H), compare, for example, <cit.>. In 
order to keep the presentation more simple, we only consider the smaller space 
L^q(0,T;V^*).

We can now state the abstract evolution equation that we want to consider. 
In the following, let {A(t)}_t ∈ [0,T] be as stated in 
Assumption <ref>, let f fulfill Assumption <ref> and let u_0 
∈ H be given. It is our overall goal to find an approximation to the solution 
u of

    u' + A u = f      in  L^q(0,T;V^*),
    
    			u(0) = u_0    in H.

This evolution equation is uniquely solvable in a variational sense with a
solution u in C([0,T];H), where

    = { v ∈ L^p(0,T;V) : v' ∈ L^q(0,T;V^*) };

see <cit.> and 
<cit.> for further details. In the following 
analysis, we employ the sum splitting in 
order to obtain a temporal discretization of (<ref>). To this end, we
consider an equidistant grid on [0,T], where N ∈, k = T/N 
and t_n = n k for n ∈{ 0, …, N}. For ℓ∈{1, …,s} and n 
∈{1,…,N} we introduce

    _ℓ^n = A_ℓ(t_n)
    	  and  _ℓ^n = 1/k∫_t_n-1^t_n f_ℓ(t) t.

We use this to construct an approximation ^n ≈ u(t_n) of the
solution u of (<ref>) for n∈{0,…,N }.
This approximation  is given through a recursion

    _ℓ^n - ^n-1/k + s _ℓ^n _ℓ^n = s 
    		_ℓ^n      in  V_ℓ^*,
    		  ℓ∈{ 1,…,s},
    ^n= 1/s∑_ℓ = 1^s_ℓ^n     in  H

for n ∈{1,…,N} with ^0 = u_0.


 
	 
  A useful example that fits into our abstract setting is to approximate the 
  solution of the parabolic p-Laplace equation.
  Let Ω⊂^d, d≥ 1, be given, where Ω is a bounded
  domain and the boundary ∂Ω is Lipschitz. For p ≥ 2 we 
  consider the problem
  
    u_t(t,x) - ∇· (α(t) |∇ u(t,x)|^p-2∇ u(t,x))  = 
          g(t,x),
               (t,x) ∈ (0,T) ×Ω,
    α(t) |∇ u(t,x)|^p-2∇ u(t,x) · n = 0,
             (t,x) ∈ (0,T) ×∂Ω,
    
          u(0,x) = u_0(x),    x ∈Ω,

  where n denotes outer pointing normal vector. The function α 
  [0,T] → is an element of C([0,T]), u_0 ∈ L^2(Ω) and g  
  (0,T) ×Ω→ is a suitably chosen integrable function that we 
  explain in more detail at a later point. 
  Applications for this type of equation can be found in <cit.>.
  Our theory allows to solve (<ref>) with the help of a 
	domain decomposition scheme. A similar setting can be found in 
	<cit.> for p = 2. The case p ≥ 2 for an 
	autonomous problem with a more restrictive domain decomposition around the 
	boundary can be found in <cit.>.
	For s∈ let {Ω_ℓ}_ℓ =1^s be a family of overlapping 
	subsets of Ω such that ⋃_ℓ =1^s Ω_ℓ = Ω 
	is fulfilled. Furthermore, let each Ω_ℓ, ℓ∈{1,…,s}, be 
	either an open connected set with a Lipschitz boundary or a union of pairwise 
	disjoint open, connected sets Ω_ℓ,i such that ⋃_i=1^rΩ_ℓ,i = Ω_ℓ and Ω_ℓ,i has a Lipschitz 
	boundary for every ℓ∈{1,…,s} and i ∈{1,…,r}; see 
	Fig. <ref>. 
 
  On these subdomains let the partition of unity {χ_ℓ}_ℓ 
  =1^s⊂ W^1,∞(Ω) be given such that
	
    χ_ℓ (x)>0 for all x∈Ω_ℓ,
    		  χ_ℓ (x) = 0 for all x∈Ω∖Ω_ℓ,
    		  ∑_ℓ =1^sχ_ℓ= 1

	for ℓ∈{1,…,s}. For such a function χ_ℓ, ℓ∈{1,…,s}, the weighted Lebesgue space L^p(Ω_ℓ,χ_ℓ)^d 
	consists of all measurable functions v = (v_1,…,v_d) Ω_ℓ→^d such that
	
    (v_1,…,v_d)_L^p(Ω_ℓ,χ_ℓ)^d
    		= (∫_Ω_ℓχ_ℓ |(v_1,…,v_d)|^p 
    		x)^1/p

	is finite. The space L^p(Ω_ℓ,χ_ℓ)^d is a reflexive Banach
	space; see <cit.> and 
	<cit.>. Note that L^p(Ω_ℓ)^d 
	is a subspace of L^p(Ω_ℓ,χ_ℓ)^d and it holds true that 	
	v_L^p(Ω_ℓ,χ_ℓ)^d≤v _L^p(Ω_ℓ)^d for 
	every v ∈ L^p(Ω_ℓ)^d.
	
	
	For (H, [H]··, ·_H ) we use L^2(Ω) 
	the space of square integrable functions on Ω with the usual norm and 
	inner product. The energetic spaces V and V_ℓ are then given as
	
    V    = { u ∈ H: there exists a  v = (v_1,…,v_d) ∈ 
    		L^p(Ω)^d  such 
    		that 
               ∫_Ω u D_i φx = ∫_Ω 
    		v_i 
    		φx  for all φ∈ C^∞_0(Ω), i=1,…,d  }
        = W_0^1,p(Ω)

	and
	
    V_ℓ   = { u ∈ H: there exists a  v = (v_1,…,v_d) ∈ 
    		L^p(Ω_ℓ,χ_ℓ)^d  such that 
               ∫_Ω u D_i (χ_ℓφ) x = 
    		∫_Ω_ℓ v_i χ_ℓφx  for all φ∈ C^∞_0(Ω), i=1,…,d   },

	which are equipped with the norms 
	
    ·_V
    		= ·_H + ∇·_L^p(Ω)^d  and  ·_V_ℓ
    		= ·_H + ∇·_L^p(Ω_ℓ,χ_ℓ)^d.

	For t∈ [0,T] we introduce the operator A(t)  V → V^*
	
    A(t) uv = ∫_Ωα(t) |∇ u|^p-2∇ u ·∇ v x,
    		   u,v∈ V.

	Together with the partition of unity we define the decomposed energetic operators 
	A_ℓ(t)  V_ℓ→ V^*_ℓ, ℓ∈{ 1,…,s},
	
    A_ℓ(t) uv
    		= ∫_Ω_ℓχ_ℓα(t) |∇ u|^p-2∇ u ·∇ v x,
    		   u,v∈ V_ℓ,    t∈ [0,T].

	It is also possible to allow for more general coefficients α [0,T] 
	×Ω×^d →^d where α(t,·,·) 
	fulfills the condition stated in <cit.>.
	
	We assume that for g  (0,T) ×Ω→ the 
	abstract function [f(t)](x) = g(t,x), (t,x) ∈ (0,T) ×Ω, is an 
	element of L^q(0,T;V^*). We exploit that f(t) ∈ V^*, t ∈ (0,T), can be 
	represented by
	
    f(t)v = ∫_Ω f^0(t) v x
    		+ ∑_i = 1^d∫_Ω f^i(t) D_i v x,    v ∈ V

	where f^i(t) ∈ L^q(Ω) for i ∈{0,…,d}. These functions are not 
	necessarily unique unless we exchange V = W^1,p(Ω) by V = 
	W_0^1,p(Ω), compare <cit.>. 
	This in mind, we introduce f_ℓ(t) for a.e. t ∈ (0,T) as
	
    f_ℓ(t)v
    		= ∫_Ωχ_ℓ f^0(t) v x
    		+ ∑_i = 1^d∫_Ωχ_ℓ f^i(t) D_i v x,   v ∈ 
    		V_ℓ.

	
	Note that in this type of setting, we can also consider homogeniuous Direchlet 
	boundary conditions in (<ref>). Then an additional condition on 
	the partition of unity becomes necessary. In this case, we have to make the 
	further assumption that for every function χ_ℓ there exists 
	ε_0>0 such that for all ε∈ (0,ε_0)
	
    Ω_ℓ^ε
    		= { x∈Ω_ℓ : χ_ℓ(x) ≥ε}

	is a Lipschitz domain.


Further examples that fit our framework are a domain decomposition scheme for 
the porous medium equation as presented in 
<cit.> or a source term splitting 
as in <cit.>. An application to the third equation 
of (<ref>) is presented in <cit.>. Numerical experiments for  
this equation and the p-Laplace equation can be found in 
<cit.> and <cit.>, respectively.





§ SOLVABILITY AND A PRIORI BOUNDS FOR THE DISCRETE SCHEME



The abstract setting from the previous section in mind we are now well-prepared 
to state some properties of the solution of the numerical scheme 
(<ref>). Since the scheme is implicit, we start to verify that 
(<ref>) is uniquely 
solvable. Once this is at hand, we can provide a priori bounds of the solution. 
These bounds are a crucial part of the further analysis and allow for the 
convergence analysis in Section <ref>


	Let Assumptions <ref> and <ref> be fulfilled.
	Then the semidiscrete problem (<ref>) is uniquely
	solvable.



	In order to prove the existence of the elements (_ℓ^i)_i∈{0,…,N}, ℓ∈{1,…,s}, that fulfill (<ref>), 
	we argue inductively. Assuming that for i ∈{1,…,N} the 
	previous elements	(_ℓ^j)_j∈{0,…,i-1}, ℓ∈{1,…,s}, exist in the corresponding spaces, we prove the existence of 
	_ℓ^i ∈ V_ℓ for  every ℓ∈{1,…,s}.
	The operator I + s k _ℓ^i, ℓ∈{1,…,s}, is strictly 
	monotone due to (3) of Assumption <ref>, i.e., it holds true that
	
    (I + s k _ℓ^i)v - (I + s k _ℓ^i) w  v - w > 0, 
    		   v,w ∈ V_ℓ with  v ≠ w

	for every ℓ∈{1,…,s}.
	Furthermore, I + s k _ℓ^i is radially continuous as A_ℓ(t), ℓ∈{1,…,s}, is radially continuous for every t ∈ [0,T]. It remains to 
	verify that the operator is coercive. Using (5) of Assumption <ref> and 
	the norm bound of Assumption <ref>, it follows 
	
    (I + s k _ℓ^i) vv/v_V_ℓ   ≥v_H^2 + sk  μ|v|_V_ℓ^p
    			- sk λ/ c_V_ℓ(v_H + |v|_V_ℓ)  
       ≥min(1,sk μ)/c_V_ℓ·v_H^2 + 
    			|v|_V_ℓ^p
    		/v_H + |v|_V_ℓ
    		- sk λ/ c_V_ℓ(v_H + |v|_V_ℓ) →∞

  as v_V_ℓ→∞ for v ∈ V_ℓ and ℓ∈{1,…,s}.
	Thus, for ^i-1 = 1/s∑_ℓ = 1^s_ℓ^i-1∈ H,
	there exists a unique solution _ℓ^i ∈
	V_ℓ of
	
    ( I + s k _ℓ^i ) _ℓ^i = s k _ℓ^i + ^i-1

	for every ℓ∈{1,…,s} due Browder–Minty theorem; see
	<cit.> for further details.


We can now turn our attention to the a priori bounds. 


	Let Assumptions <ref> and <ref> be fulfilled.
	Then for the unique solution of (<ref>) there exist constants 
	M, M' < ∞ such that for every step size k = T/N the a priori 
	bounds
	
    max_n ∈{1,…,N}(1/s∑_ℓ = 1^s_ℓ^n^2_H )
    			+ 1/s∑_i = 1^N∑_ℓ = 1^s_ℓ^i -
    			^i-1^2_H
    			+ k ∑_i = 1^N∑_ℓ =1 ^s_ℓ^i_V_ℓ^p 
    			≤ M

	and
	
    1/k∑_i=1^N^i - ^i-1/k^q_V^* 
        = k^1-q∑_i=1^N^i - ^i-1^q_V^*≤ M'

	are fulfilled.



	In the following, let i ∈{1,…,N} and ℓ∈{1,…,s} be 
	arbitrary but fixed. Recall the identity
	
    [H]_ℓ^i - ^i-1_ℓ^i
    		= 1/2( _ℓ^i^2_H - ^i-1^2_H
    		+ _ℓ^i - ^i-1^2_H )

	and the inequality ·_V_ℓ≤ c_1 ( ·_H + 
	|· |_V_ℓ) with c_1 = max_ℓ∈{1,…,s} 
	c_V_ℓ stated in Assumption <ref>. Using the 
	weighted Young inequality, see <cit.>), we 
	obtain 
	
    1/2k( _ℓ^i^2_H - ^i-1^2_H
    		+ _ℓ^i - ^i-1^2_H )
    		+ s _ℓ^i _ℓ^i_ℓ^i
       = s _ℓ^i_ℓ^i
    		≤ s c_1  _ℓ^i_V_ℓ^*( _ℓ^i _H + 
    		|_ℓ^i |_V_ℓ)
       ≤ s c_1  _ℓ^i_V_ℓ^*_ℓ^i _H 
    		+ s c_2  _ℓ^i_V_ℓ^*^q 
    		+ s μ/2  |_ℓ^i |_V_ℓ^p

	with c_2 = c_1^q (pμ)^1-q/q 2^1-q. Thus, together 
	with the coercivity condition from Assumption <ref> (5) it follows 
	that
	
    _ℓ^i^2_H - ^i-1^2_H + _ℓ^i - ^i-1^2_H
          + k s μ  |_ℓ^i |_V_ℓ^p 
       ≤ 2k s c_1  _ℓ^i_V_ℓ^*_ℓ^i _H 
          + 2 k s c_2  _ℓ^i_V_ℓ^*^q + 2 k s λ.

	Employing the specific structure of ^i-1, we obtain
	
    ^i-1^2_H
             = 1/s∑_ℓ = 1^s_ℓ^i-1^2_H
          ≤1/s^2( ∑_ℓ = 1^s_ℓ^i-1_H
          )^2
       ≤1/s^2∑_ℓ = 1^s 1^2 ·∑_ℓ = 1^s_ℓ^i-1_H^2
          = 1/s∑_ℓ = 1^s_ℓ^i-1^2_H

	for i ∈{2,…,N} due to the Cauchy–Schwarz inequality for sums.
	Inserting this inequality in (<ref>), summing up from ℓ 
	= 1 to s as well as dividing by s, yields
	
    1/s∑_ℓ = 1^s( _ℓ^i^2_H -
    		_ℓ^i-1^2_H + _ℓ^i - ^i-1^2_H ) 
    		+ k  μ∑_ℓ = 1^s  |_ℓ^i |_V_ℓ^p 
       ≤ 2k c_1 ∑_ℓ = 1^s_ℓ^i_V_ℓ^*_ℓ^i _H 
    		+ 2 k c_2 ∑_ℓ = 1^s_ℓ^i_V_ℓ^*^q 
    		+ 2 k s λ

	for i ∈{2,…,N} and
	
    1/s∑_ℓ = 1^s( _ℓ^1^2_H 
    		+ _ℓ^i - u_0 ^2_H ) 
    		+ k  μ∑_ℓ = 1^s  |_ℓ^1 |_V_ℓ^p 
       ≤u_0_H^2 + 2k c_1  ∑_ℓ = 1^s_ℓ^1_V_ℓ^*_ℓ^1 _H 
    		+ 2 k c_2 ∑_ℓ = 1^s_ℓ^1_V_ℓ^*^q 
    		+ 2 k s λ.

	After a summation from i=1 to n ∈{1,…,N} and using the 
	telescopic structure, we obtain
	
    1/s∑_ℓ = 1^s_ℓ^n^2_H 
    		+ 1/s∑_i=1^n∑_ℓ = 1^s_ℓ^i - ^i-1^2_H
    		+ k  μ∑_i = 1^n∑_ℓ = 1^s  |_ℓ^i |_V_ℓ^p 
       ≤u_0_H^2 + 2k c_1 ∑_i = 1^n∑_ℓ = 1^s_ℓ^i_V_ℓ^*_ℓ^i _H 
    		+ 2 k c_2 ∑_i = 1^n∑_ℓ = 1^s_ℓ^i_V_ℓ^*^q 
    		+ 2 T s λ.

	For the right-hand side we can bound the summands using
	Assumption <ref> and Hölder's inequality
	
    k ∑_i = 1^n∑_ℓ = 1^s_ℓ^i_V_ℓ^*^q
          = k ∑_i = 1^n∑_ℓ = 1^s1/k∫_t_i-1^t_i f_ℓ(t) t_V_ℓ^*^q
       ≤∑_i = 1^n∑_ℓ = 1^s∫_t_i-1^t_i f_ℓ(t) _V_ℓ^*^qt≤ s  f _L^q(0,T;V^*)^q

	and
	
    k _ℓ^i_V_ℓ^*≤ k 1/k∫_t_i-1^t_i f_ℓ(t) t_V_ℓ^*≤∫_t_i-1^t_i f(t) _V^*t.

	Thus, we get that
	
    1/s∑_ℓ = 1^s_ℓ^n^2_H 
    			+ 1/s∑_i = 1^n∑_ℓ = 1^s_ℓ^i - ^i-1^2_H
    			+ k  μ∑_i = 1^n∑_ℓ = 1^s  |_ℓ^i |_V_ℓ^p 
       ≤u_0_H^2 + 2k c_1   ∑_i = 1^n∫_t_i-1^t_i f(t) _V^*t∑_ℓ = 1^s_ℓ^i 
    			_H 
    			+ 2 k s c_2  f _L^q(0,T;V^*)^q
    			+ 2 T s λ.

	As this is fulfilled for every n ∈{1,…,N}, it also follows that
	
    max_n ∈{1,…,N}(1/s∑_ℓ = 1^s_ℓ^n^2_H 
    		+ 1/s∑_i=1^n∑_ℓ = 1^s_ℓ^i - ^i-1^2_H
    		+ k  μ∑_i = 1^n∑_ℓ = 1^s  |_ℓ^i 
    		|_V_ℓ^p) 
       ≤u_0_H^2 + 2k s c_1  f _L^1(0,T;V^*)max_n ∈{1,…,N}(1/s∑_ℓ = 1^s_ℓ^n_H^2)^1/2
            + 2 k s c_2  f _L^q(0,T;V^*)^q + 2 T s λ.

	We abbreviate the terms 
	
    x^2    = max_n ∈{1,…,N}(1/s∑_ℓ = 1^s_ℓ^n^2_H 
    		+ 1/s∑_i=1^n∑_ℓ = 1^s_ℓ^i - ^i-1^2_H + k  μ∑_i = 1^n∑_ℓ = 1^s  |_ℓ^i 
    		|_V_ℓ^p)
    
    		a    = k s c_1  f _L^1(0,T;V^*)
    
    		b^2    = u_0_H^2 + 2 k s c_2  f _L^q(0,T;V^*)^q
    		+ 2 T s λ

	to obtain x^2 ≤ 2ax + b^2. This implies, in particular, that
  
    (x-a)^2 = x^2 -2ax +a^2 ≤ a^2 +b^2.

  Taking the square root on both sides, this yields
  
    |x-a| ≤√( a^2 +b^2)≤ a + b.

  As x-a ≤ |x-a| is fulfilled, we obtain x ≤ 2a + b after adding a to 
  both sides of the inequality. This shows that
	
    max_n ∈{1,…,N}(1/s∑_ℓ = 1^s_ℓ^n^2_H 
    		+ 1/s∑_i=1^n∑_ℓ = 1^s_ℓ^i - ^i-1^2_H + k  μ∑_i = 1^n∑_ℓ = 1^s  |_ℓ^i 
    		|_V_ℓ^p) ≤ M_1,
 
	where M_1 ≥ 0 is independent of k.  Using the norm inequality from 
	Assumption <ref>, this implies that there exists M_2 ≥ 0, 
	which does not depend on k, such that
	
    (k ∑_i = 1^N∑_ℓ = 1^s_ℓ^i 
    		_V_ℓ^p)^1/p≤ c_1 (k ∑_i = 1^N∑_ℓ = 1^s_ℓ^i 
    		_H^p)^1/p
    		+ c_1 (k ∑_i = 1^N∑_ℓ = 1^s  |_ℓ^i 
    		|_V_ℓ^p)^1/p≤ M_2.

	Altogether, we have proved the first a priori bound (<ref>).
	
	
	In order to prove (<ref>), we test (<ref>) 
	with v ∈ V and use Assumption <ref> (4) to see that
	
    [H]^i - ^i-1/kv   = 1/s∑_ℓ = 1^s[H]_ℓ^i - ^i-1/kv
       = ∑_ℓ = 1^s( _ℓ^iv - _ℓ^i
    			_ℓ^i v)
       ≤ c_3 v_V(∑_ℓ = 1^s_ℓ^i_V_ℓ^*
    		+ β∑_ℓ = 1^s(1 + _ℓ^i _V_ℓ^p-1))

	for i ∈{1,…,N}, where c_3 is the maximal embedding constant of 
	V into V_ℓ for ℓ∈{1,…,s}.
	Thus, we can estimate the V^*-norm by
	
    k^-1^i - ^i-1_V^*≤  c_3(∑_ℓ = 1^s_ℓ^i_V_ℓ^*
    		+ β∑_ℓ = 1^s(1 + _ℓ^i _V_ℓ^p-1)).

	This bound can be used to see that there exists M' ≥ 0 such that
	
    (k^1-q∑_i=1^N^i - ^i-1^q_V^*)^1/q
       ≤ c_3(k ∑_i=1^N(∑_ℓ = 1^s_ℓ^i_V_ℓ^* + β∑_ℓ = 1^s(1 + _ℓ^i 
    		_V_ℓ^p-1) )^q)^1/q
       ≤  c_3(k ∑_i=1^N∑_ℓ = 1^s_ℓ^i_V_ℓ^*^q )^1/q
    		+  c_3 β (Ts)^1/q
    		+  c_3 β(k ∑_i=1^N∑_ℓ = 1^s_ℓ^i 
    		_V_ℓ^p )^1/q≤ M'.

	Due to the first a priori bound (<ref>) and (<ref>) 
	the constant M' is independent of k.




§ CONVERGENCE ANALYSIS
 

In the following, we introduce prolongations of the solution of the discrete 
problem (<ref>) to the interval [0,T]. The main goal of this 
section is to prove that the sequence of such prolongations converges to the 
exact solution u of (<ref>). 
Corresponding to the grid 0 = t_0 < t_1 < … <t_N = T with k = 
T/N and t_n = n k, n ∈{0,…,N}, we construct piecewise 
constant and piecewise linear functions on the interval [0,T].
We consider the piecewise constant functions for t ∈ 
(t_n-1,t_n], n ∈{1,…,N}, and ℓ∈{1,…,s} given by

    U^k_ℓ(t) = _ℓ^n,   
    	U^k(t) = ^n,   
    	A^k_ℓ(t) = _ℓ^n,   and  
    	f^k_ℓ(t) = _ℓ^n

as well as the piecewise linear function

    Ũ^k(t) = ^n-1 + t - t_n-1/k
    	(^n - ^n-1)

with U^k_ℓ(0) = U^k(0) = Ũ^k(0) = u_0, A^k_ℓ(0) = 
_ℓ^1 and f^k_ℓ(0) = _ℓ^1.
As we consider step sizes k = T/N for N ∈, we denote
the sequences (U^T/N_ℓ)_N ∈ as
(U^k_ℓ )_k >0 for ℓ∈{1,…,s} in the following to keep the
notation more compact.
The same simplification in notation is used for the other functions
introduced above. Due to the a priori bound (<ref>) we see that

    U^k_ℓ∈ L^p(0,T;V_ℓ) ∩ L^∞(0,T;H), 	U^k, Ũ^k ∈ L^∞(0,T;H),  and	f^k_ℓ∈ L^q(0,T;V_ℓ^*).

Furthermore, due to Lemma <ref> the operator A_ℓ^k
maps the space L^p(0,T;V_ℓ) into L^q(0,T;V_ℓ^*).
Using the prolongations introduced above, we can state a discrete version of 
the differential equation. We first note that after summing up 
(<ref>) from 1 to s and dividing by s, we obtain

    1/ks∑_ℓ = 1^s( _ℓ^n - ^n-1)
    	+ ∑_ℓ = 1^s_ℓ^n _ℓ^n
    	= ∑_ℓ = 1^s_ℓ^n   in V^*.

Thus, we see that

    (Ũ^k)'(t)
    			+ ∑_ℓ = 1^s A^k_ℓ(t) U_ℓ^k(t)
    			= ∑_ℓ = 1^s f^k_ℓ(t)
    			     in V^*,
    			   t ∈ (0, T),
    
    			U^k(0) = Ũ^k(0) = u_0      in H,

where (Ũ^k)' is the weak derivative of Ũ^k. In the 
following, we will consider the limiting process of all the
appearing terms to connect to the original problem (<ref>) with 
(<ref>).


 
	Let Assumption <ref> be fulfilled and let W∈ L^p(0, T;V) be
	given.
	For ℓ∈{1,…,s} it follows that A_ℓ^k(t) W(t) → A_ℓ(t) 
	W(t) in V_ℓ^* as k → 0 for a.e. t∈ (0,T). Furthermore, it holds 
	true that A_ℓ^k W → A_ℓ W in L^q(0,T;V_ℓ^*) as k → 0.



	Let ℓ∈{1,…,s} and ε> 0 be arbitrary. Due to the 
	continuity condition on A_ℓ, for almost every t ∈ (0,T) we find 
	δ >0 such that for all k < δ it follows that
	
    A_ℓ^k(t) W(t) - A_ℓ(t) W(t) _V^*
    		=  A_ℓ(t_n) W(t) - A_ℓ(t) W(t) _V^*≤ε,

	where t is within an interval (t_n-1,t_n], n ∈{1,…,N}.
	The second assertion of the lemma is a consequence of Lebesgue's 
	theorem of dominated convergence and the boundedness condition (4) 
	from Assumption <ref>.


 
	Let Assumption <ref> be fulfilled. Then it follows that f_ℓ^k →
	f_ℓ in L^q(0,T;V_ℓ^*), ℓ∈{1,…,s}, as k → 0.


 
	The statement above can easily be verified for a function from the space 
	C([0,T] ;V_ℓ^*). As the space C([0,T] ;V_ℓ^*) is a dense subspace 
	of L^q(0,T;V_ℓ^*) for ℓ∈{1,…,s} a density argument can 
	be used to verify the claimed statement.



	Let Assumptions <ref> and <ref> be fulfilled.
	Then there exists a subsequence  (k_i)_i ∈ of step sizes
	k_i = T/N_i and U ∈ such that
	
    U_ℓ^k_i U   in  L^p(0, T; V_ℓ)
    		  and  
    		U_ℓ^k_i U,    U^k_i U   in  
    		L^∞(0,T; H),

	as well as
	
    Ũ^k_i U   in  L^∞(0, T; H)
    		  and  
    		(Ũ^k_i)'  U'   in  L^q(0, T;V^*)

	for every ℓ∈{ 1,…, s } as i →∞. Here, U' denotes 
	the weak derivative of U.




	In the following proof, we do not distinguish between subsequences by 
	notation. Using Lemma <ref>, we obtain that
	
    U_ℓ^k_L^∞(0,T;H)^2
    		≤ s M,
        U^k_L^∞(0,T;H)^2
    		= Ũ^k_L^∞(0,T;H)^2
    		≤ M,
        U_ℓ^k_L^p(0,T;V_ℓ)^p = k ∑_i = 1^N_ℓ^i_V_ℓ^p
    		≤  M,
        (Ũ^k)' _L^q(0,T;V^*)^q = k^1-q∑_i=1^N^n - ^n-1^q_V^*≤ M'.

	Therefore, the sequence (U_ℓ^k)_k >0 is bounded in
	L^p(0, T;V_ℓ) as well as L^∞(0, T;H),
	(Ũ^k)_k >0 is bounded in L^∞(0, T;H), and
	( (Ũ^k)' )_k >0 is bounded in L^q(0, T;V^*).
	Since L^p(0, T;V_ℓ) is a reflexive Banach space and
	L^∞(0, T;H) is the dual space of the separable Banach space
	L^1(0, T;H), there exists a subsequence of (U_ℓ^k)_k>0 and
	U_ℓ∈  L^p(0, T;V_ℓ) ∩ L^∞(0, T;H) such that
	
    U_ℓ^k U_ℓ  in  L^p(0,T;V_ℓ)
    		  and  
    		U_ℓ^k U_ℓ  in  L^∞(0, T;H)

	as k → 0. Analogously, there exist a suitable further subsequence, 
	Ũ∈ L^∞(0, T;H) and W ∈ L^q(0, T; V^*) such that
	
    Ũ^kŨ  in  L^∞(0, T;H)
    		  and  
    		(Ũ^k)'  W   in  L^q(0, T;V^*)

	as k → 0.
	In the following, we prove that U_1 = … = U_s =: U is
	fulfilled. As ⋂_ℓ = 1^s V_ℓ  = V and the norm ∑_ℓ = 
	1^s·_V_ℓ is equivalent to ·_V this implies that U 
	∈ L^p(0,T;V).
	We show that U_1 = U_2 in L^p(0,T;V_1∩ V_2) ∩ 
	L^∞(0,T;H), the 
	other equalities follow in an analogous manner. 
	We can write  
	
    U_1^k(t) - U_2^k(t)
    		   = _1^n - ^n-1 - (_2^n - ^n-1)
       = k s ( _1^n - _1^n _1^n )
    		- k s ( _2^n - _2^n _2^n )
       = s ∫_t_n-1^t_n(( f_1(τ) - A_1(τ)
    		U_1^k(τ) ) - (  f_2(τ) - A_2(τ)
    		U_2^k(τ)  )) τ

	for t∈ (t_n-1,t_n], n∈{1,…,N},  as U_1^k(t_n-1) = ^n-1 
	= U_2^k(t_n-1) holds true by the construction of our scheme.
	Therefore, we obtain
	
    U_1^k(t) - U_2^k(t) _V^*   ≤ s ∑_ℓ  = 1^2  ∫_t_n-1^t_n f_ℓ(τ) - 
    		A_ℓ(τ)
    		U_ℓ^k(τ)_V^*τ
       ≤ s k^1/p∑_ℓ  = 1^2 ( ∫_t_n-1^t_nf_ℓ(τ) - A_ℓ(τ) U_ℓ^k(τ)_V^*^qτ)^1/q,

	where we used Hölder's inequality in the last step. We can bound the
	integrals by
	
    (∫_t_n-1^t_nf_ℓ(τ) - A_ℓ(τ)
    		U_ℓ^k(τ)_V^*^qτ)^1/q
       ≤(∫_t_n-1^t_nf_ℓ(τ)_V^*^qτ)^1/q
    		+ (∫_t_n-1^t_nA_ℓ(t_n)
    		_ℓ^n_V^*^qτ)^1/q
       ≤f_ℓ_L^q(0,T;V_ℓ^*)
    		+ k^1/qβ(1 + _ℓ^n_V_ℓ^p-1),

	for ℓ∈{1,2} which is bounded independently of n and k due to 
	the a priori bound (<ref>).
	Thus, it follows that U_1^k(t) - U_2^k(t) _V^*→ 0 as k → 0 
	for every t ∈ [0,T] and it is bounded by a  constant independent of t. Using
	Lebesgue's theorem of dominated convergence and the fact that L^q(0,T;V^*)  
	L^1(0,T;V^*), it follows that
	
    U_1^k - U_2^k _L^1(0,T;V^*)≤ c_1  U_1^k - U_2^k _L^q(0,T;V^*)→ 0   as  k → 0

	for c_1 ∈ (0,∞).
	This shows that U_1 - U_2 = 0 in L^1(0,T;V^*). By assumption the 
	embedding L^p(0,T;V_1∩ V_2) ∩ L^∞(0,T;H)   
	L^1(0,T;V^*) is continuous.
	The injectivity of the embedding operator implies that both U_1 - U_2 = 
	0 in L^p(0,T;V_1 ∩ V_2) and U_1 - U_2 = 0 in L^∞(0,T;H).
	The limits U and Ũ coincide in L^1(0,T;V^*) since
	
    U_1^k - U_2^k _L^1(0,T;V^*)^q
    		   ≤ c_1^q ∫_0^T U^k(t) - Ũ^k(t) _V^*^q t
       = ∑_n=1^N∫_t_n-1^t_n^n - ^n-1 -
    		t - t_n-1/k (^n - ^n-1) _V^*^q t
       =  1/k^q∑_n=1^N^n - ^n-1_V^*^q
    		∫_t_n-1^t_n (t_n - t)^q  t
       =  k/q∑_n=1^N^n - ^n-1_V^*^q
    		≤k^q/q M' → 0
    		   as  k → 0,

	where we used the a priori bound (<ref>). Making use of 
	the continuous embedding L^∞(0,T;H)  L^1(0,T;V^*), it follows that 
	U and Ũ coincide in L^∞(0,T;H).
	
	
	Last, we prove that the limit W ∈ L^q(0,T;V^*) is the weak
	derivative of U. To this end, let v ∈ V and φ∈
	C_c^∞(0,T) be arbitrary. Using Ũ^k  U in 
	L^∞(0,T;H) as k → 0, it yields
	
    -∫_0^TW(t) vφ(t) t   = -lim_k → 0∫_0^T(Ũ^k)'(t) vφ(t) 
    		t
       = lim_k → 0∫_0^T[H]Ũ^k(t) vφ'(t) 
    		t
    		= ∫_0^T[H]U(t) vφ'(t) t.

	Applying <cit.>, we obtain W = U' in 
	L^q(0,T;V^*).


The next lemma is an auxiliary result to identify the limit from the previous
lemma with the solution of (<ref>).

 
	Let Assumption <ref> be fulfilled and let the operator A_ℓ^k 
	 L^p(0,T;V_ℓ) → L^q(0,T;V_ℓ^*) be given as in
	(<ref>). Then for a sequence ( W_ℓ^k
	)_k >0 in L^p(0,T ;V_ℓ), ℓ∈{1,…,s}, and an
	element W ∈ L^p(0,T ;V) such that
	
    W_ℓ^k  W
    		  in  L^p(0,T ;V_ℓ)
    		  and  
    		A_ℓ^k W_ℓ^k  B_ℓ  in  L^q(0,T ;V_ℓ^*)

	for every ℓ∈{1,…,s} as k → 0 and ∑_ℓ = 1^s
	B_ℓ = B ∈ L^q(0,T;V^*) with
	
    lim sup_k → 0∑_ℓ = 1^s∫_0^TA_ℓ^k(t) W_ℓ^k(t)W_ℓ^k(t)t≤∫_0^TB(t) W(t)t

	it follows that AW = ∑_ℓ = 1^s A_ℓ W= ∑_ℓ = 1^s 
	B_ℓ = B in L^q(0,T;V^*).



	Due to the monotonicity of A_ℓ(t), t ∈ [0,T] and ℓ∈{1,…,s},
	we can write for every X ∈ L^p(0,T;V)
	
    ∑_ℓ = 1^s∫_0^TA_ℓ^k(t) W_ℓ^k(t)
    			- A_ℓ^k(t) X(t) W_ℓ^k(t) - X(t) t≥ 0.

	Thus, using Lemma <ref> it follows that
	
    ∑_ℓ = 1^s∫_0^T
    			A_ℓ^k(t) W_ℓ^k(t)W_ℓ^k(t)t
       ≥∑_ℓ = 1^s∫_0^T(  A_ℓ^k(t) W_ℓ^k(t) X(t)  +
    		 A_ℓ^k(t) X(t) W_ℓ^k(t) - X(t)) t
       k → 0⟶∑_ℓ = 1^s∫_0^T( B_ℓ(t) X(t)  + A_ℓ(t) X(t)  W(t) - 
    		X(t) ) t
       = ∫_0^T( B(t) X(t) 
    		+A(t) X(t)  W(t) - X(t) ) t.

	This implies
	
    lim inf_k → 0∑_ℓ = 1^s∫_0^T A_ℓ^k(t)
    			W_ℓ^k(t)W_ℓ^k(t)t
       ≥∫_0^T(B(t) X(t) 
    		+  A(t) X(t)  W(t) - X(t)) t.

	Applying (<ref>), this yields
	
    ∫_0^TB(t) W(t)t
       ≥∫_0^T( B(t) X(t) 
    		+ A(t) X(t) W(t) - X(t) ) t.

	The assertion of the lemma follows by the Minty monotonicity trick, 
	where X =  W ±θX̃ for θ∈ [0,1] and 
	X̃∈ L^p(0,T; V) is inserted in the inequality. 
	Dividing by θ and considering θ→ 0 then yields A W= B 
	in L^q(0,T;V^*). See, e.g., <cit.> for further 
	details.


Combining the prior lemmas, we can now state one of the main results of
this section. We prove that the limit of the sequence of prolongations 
is the solution of (<ref>).


	Let Assumptions <ref> and <ref> be fulfilled. Let u be
	the solution of (<ref>). Then for step sizes k = T/N the
	sequences (U_ℓ^k)_k>0, (U^k)_k>0 and  (Ũ^k)_k>0 
	defined in (<ref>) and (<ref>), respectively, 
	fulfill
	
		
		U_ℓ^k    u         in 
		L^p(0,T;V_ℓ), 

		
		U^k(t)   u(t)         in  H, 

		
		U^k u ,  Ũ^k    u         in  L^∞(0,T;H), 

		
		(Ũ^k)'    u'         in  L^q(0,T;V^*),

		
		∑_ℓ= 1^s A_ℓ^k U_ℓ^k    A u         in  L^q(0,T;V^*)
	
	as k → 0 for t∈ [0,T] and ℓ∈{1,…,s}.



	In the following, we will not distinguish between different subsequences by 
	notation. 
	Due to Lemma <ref> as well as the a priori bound
	(<ref>) there exists a constant M̃ > 0 such that for every 
	k>0, t ∈ [0,T], and ℓ∈{1,…,s}
	
    A_ℓ^k U_ℓ^k _L^q(0, T;V_ℓ^*)≤M̃  and   U^k(t) _H
    		≤M̃

	is fulfilled. Therefore, we can extract a subsequence of step sizes such 
	that there exits B_ℓ∈ L^q(0, T;V_ℓ^*) and y_t ∈ H with
	
    A_ℓ^k U_ℓ^k  B_ℓ  in 
    		L^q(0, T;V_ℓ^*)
    		  and  
    		U^k(t)  y_t   in  H

	as k → 0 for ℓ∈{1,…,s}. In the following, we abbreviate B := 
	∑_ℓ  = 1^s B_ℓ.
	Next we identify the derivative of U with equation (<ref>). Due to 
	Lemma <ref> and Lemma <ref>, 
	the following equality holds true
	
    U'
    		= _k → 0 (Ũ^k)'
    		= _k → 0∑_ℓ = 1^s(f_ℓ^k - A_ℓ^k
    		U_ℓ^k)
    		= f - B   in  L^q(0,T;V^*),

	where  denotes the weak limit. The limit U obtained in 
	Lemma <ref> is an element of C([0, T];H). 
	Thus, we can work with the continuous representative of U in the following.
	This in mind, we prove y_t = U(t) and u_0 = U(0) for t ∈ [0,T]. To
	this end, let v∈ V and φ∈ C([0,T]) ∩ C^1(0,T) be arbitrary.
	Recalling the equation for the time discrete values (<ref>), we 
	can write
	
    [H]U(t)vφ(t) - [H]U(0)vφ(0)
    		- ∫_0^t[H]U(τ)vφ'(τ) τ
        = ∫_0^tU'(τ)vφ(τ) τ
        = ∑_ℓ = 1^s∫_0^tf_ℓ(τ) -
    			B_ℓ(τ)vφ(τ) τ
        = lim_k → 0(∫_0^t_n (Ũ^k)' (τ)
    			+ ∑_ℓ = 1^s( A_ℓ^k(τ) U_ℓ^k(τ)
    			- B_ℓ (τ)) vφ(τ) τ
          -∑_ℓ = 1^s∫_t^t_nf_ℓ^k(τ) -
    			B_ℓ(τ)vφ(τ) τ)

	for t ∈ (t_n-1,t_n], n ∈{1,…,N}.
	Applying integration by parts and the fact that the linear and the
	constant interpolations always coincide at the grid points then shows that
	
    ∫_0^t_n 
    			(Ũ^k)'(τ)vφ(τ)
    		τ
       = [H]U^k(t_n)vφ(t) - [H]U^k(0)vφ(0) -
    		∫_0^t_n[H]Ũ^k(τ)vφ'(τ)
    		τ.

	Recall that in Lemma <ref> we have proved that 
	Ũ^k  U in L^∞(0,T;H) as k → 0.   
	Therefore, (<ref>) and the fact that (f_ℓ^k - 
	B_ℓ)_k>0 is a bounded sequence in L^q(0,T;V_ℓ^*) for every ℓ∈{1,…,s} shows that
	
    [H]U(t)vφ(t) - [H]U(0)vφ(0)
    		- ∫_0^t[H]U(τ)vφ'(τ) τ
        = [H]y_tvφ(t) - [H]u_0vφ(0)
    		- ∫_0^t[H]U(τ)vφ'(τ) τ,

	is fulfilled. This implies U(t) = y_t and U(0) = u_0 for every t ∈ [0,T]. 
	
	It remains to prove that B  = A U is fulfilled. To this end, we use
	Lemma <ref>. Applying Lemma <ref>, it 
	follows that U_ℓ^k  U in L^p(0, T;V_ℓ) as k → 0 for 
	ℓ∈{1,…,s}. Further, in (<ref>), we have seen that 
	A_ℓ^k U_ℓ^k  B_ℓ in L^q(0, T;V_ℓ^*) as k → 0 
	for every ℓ∈{1,…,s}.  Therefore, we still have to verify
	
    lim sup_k → 0∑_ℓ = 1^s∫_0^TA_ℓ^k(t) U_ℓ^k(t)U_ℓ^k(t)t≤∫_0^TB(t) U(t)t

	in order to apply Lemma <ref>.
	To this end, we test the semidiscrete problem (<ref>) 
	with _ℓ^n for ℓ∈{1,…,s} and n∈{1,…,N} to 
	obtain that
	
    _ℓ^n - ^n-1 + s k _ℓ^n
    			_ℓ^n_ℓ^n = sk _ℓ^n_ℓ^n.

	Summing up the equation form ℓ = 1 to s, dividing by s, and
	applying the identity from (<ref>), it follows that
	
    1/2s∑_ℓ = 1^s( _ℓ^n
    		_H^2 - ^n-1_H^2)
    		+ k ∑_ℓ = 1^s_ℓ^n _ℓ^n_ℓ^n≤ k ∑_ℓ = 1^s_ℓ^n_ℓ^n.

	After another summation for n = 1 to N and an application of  and
	(<ref>), we can rewrite this inequality to
	
    1/2s∑_ℓ = 1^s_ℓ^N _H^2 -  u_0 _H^2
    		+ k ∑_n=1^N∑_ℓ = 1^s_ℓ^n
    			_ℓ^n_ℓ^n≤ k ∑_n=1^N∑_ℓ = 1^s_ℓ^n_ℓ^n,

	due to a telescopic sum structure. Inserting the definition of the 
	prolongations from (<ref>) and (<ref>), we 
	see that
	
    1/2s∑_ℓ = 1^s( U_ℓ^k(T) _H^2 -  u_0 _H^2 
    		)
    		+∑_ℓ = 1^s∫_0^T A_ℓ^k(t)
    			U_ℓ^k(t)U_ℓ^k(t)t
       ≤∑_ℓ = 1^s∫_0^T
    			f_ℓ(t) U_ℓ^k(t)t.

	Together with (<ref>) and the weak lower semicontinuity
	of the norm this yields
	
    lim sup_k → 0∑_ℓ = 1^s∫_0^TA_ℓ^k(t) U_ℓ^k(t)U_ℓ^k(t)t
       ≤lim sup_k → 0(∑_ℓ = 1^s∫_0^T
    			f_ℓ^k(t)U_ℓ^k(t)t
    		- 1/2s∑_ℓ = 1^s(U_ℓ^k(T) _H^2
    		-  u_0_H^2))
       ≤lim sup_k → 0(∑_ℓ = 1^s∫_0^T
    			f_ℓ^k(t)U_ℓ^k(t)t
    		- 1/2 U^k(T) _H^2
    		+ 1/2 u_0_H^2))
       ≤∑_ℓ = 1^s∫_0^Tf_ℓ(t) U(t)t
    		- 1/2 U(T) _H^2 +  1/2 u_0_H^2
       = ∫_0^Tf(t) U(t)t -
    		∫_0^TU'(t)U(t)t.

	Therefore, we have proved that
	
    lim sup_k → 0∑_ℓ = 1^s∫_0^TA_ℓ^k(t) U_ℓ^k(t)U_ℓ^k(t)t≤∫_0^TB(t) U(t)t.

	Applying Lemma <ref>, this verifies that B = AU is fulfilled in
	L^q(0,T;V^*). Thus, U is a variational solution to the original 
	problem (<ref>). 
	As this problem has a unique solution u ∈, it follows that U = u.
	
	Next, we argue that the original sequence (U_ℓ^k)_k>0
	converges weakly to the unique solution u of (<ref>) in
	L^p(0,T;V_ℓ) for every ℓ∈{1,…,s}.
	The arguments above show that every converging subsequence of the
	bounded sequence (U_ℓ^k)_k>0 has the limit u. Applying the
	subsequence principle, see, e.g, <cit.>,
	yields that the entire sequence converges to this limit which proves
	(<ref>).
	An analogous argumentation shows that
	(<ref>)–(<ref>) hold true for the original
	sequence. To prove (<ref>), we recall (<ref>) and 
	the statement of Lemma <ref>. Inserting these two 
	limiting process in (<ref>) yields (<ref>).




	Let Assumptions <ref> and <ref> be fulfilled.
	Then for step sizes k = T/N the sequence (U^k)_k>0 defined 
	in  (<ref>) fulfills
	
    U^k(t) → u(t)   in H
    		  as  k → 0,   for  t ∈ [0,T],

	where u is the solution (<ref>). If η in 
	Assumption <ref> (3) is strictly positive then the sequence 
	(U_ℓ^k)_k>0 converges strongly to u in L^p(0,T;V_ℓ) for 
	ℓ∈{1,…,s}.



	For the analysis we split up the terms as follows
	
    1/s∑_ℓ = 1^s u(t) - U_ℓ^k(t) _H^2 
       + 2 ∑_ℓ = 1^s∫_0^tA_ℓ^k(τ)u(τ)
    			- A_ℓ^k(τ) U_ℓ^k(τ)  u(τ) - U_ℓ^k(τ)τ
       = X_1^k(t) + X_2^k(t) + X_3^k(t)

	with
	
    X_1^k(t)
    		   =  u(t) _H^2 + 2 ∑_ℓ = 1^s∫_0^tA_ℓ^k(τ)
    			u(τ)u(τ)τ,
    
    		X_2^k(t)
    		   = - 2/s∑_ℓ = 1^s[H] u(t)  U_ℓ^k(t)  -
    		2∑_ℓ = 1^s∫_0^tA_ℓ^k(τ) u(τ)U_ℓ^k(τ)τ
         - 2∑_ℓ = 1^s∫_0^tA_ℓ^k(τ)
    			U_ℓ^k(τ)  u(τ)τ,
    
    		X_3^k(t)
    		   = 1/s∑_ℓ = 1^s U_ℓ^k(t) _H^2
    		+ 2 ∑_ℓ = 1^s∫_0^tA_ℓ^k(τ)
    			U_ℓ^k(τ) U_ℓ^k(τ)τ.

	We analyze X_1^k, X_2^k and X_3^k separately.
	For X_1^k we apply Lemma <ref> and obtain
	
    lim_k → 0 X_1^k(t)
    		=  u(t) _H^2 + 2 ∫_0^tA(τ) u(τ)u(τ)τ.

	We use Lemma <ref>, (<ref>), 
	(<ref>) and (<ref>) as well as the definition of 
	U^k from (<ref>) to see
	
    lim_k → 0 X_2^k(t)
    		   = lim_k → 0(- 2[H] u(t)  U^k(t)  -
    		2∑_ℓ = 1^s∫_0^tA_ℓ^k(τ) u(τ)U_ℓ^k(τ)τ
         - 2∑_ℓ = 1^s∫_0^tA_ℓ^k(τ)
    			U_ℓ^k(τ)  u(τ)τ)
       = - 2 u(t) _H^2 - 4 ∫_0^tA(τ) u(τ)u(τ)τ.

	The convergence of (X_3^k(t))_k >0 needs somewhat more
	attention.
	Here, we assume that t ∈ (t_n-1,t_n], n ∈{1,…,N}, and obtain
	
    X_3^k(t)
    		   = 1/s∑_ℓ =1 ^s U_ℓ^k(t) _H^2
    		+ 2 ∑_ℓ =1 ^s∫_0^tA_ℓ^k(τ)
    			U_ℓ^k(τ) U_ℓ^k(τ)τ
       ≤1/s∑_ℓ =1 ^s_ℓ^n _H^2
    		+ 2k ∑_i =1 ^n∑_ℓ =1 ^s_ℓ^i _ℓ^i - _ℓ^i_ℓ^i
          + 2∑_ℓ =1 ^s∫_0^t_nf_ℓ^k(τ)U_ℓ^k(τ)τ
          - 2 ∑_ℓ =1 ^s∫_t^t_nA_ℓ^k(τ)
    			U_ℓ^k(τ) U_ℓ^k(τ)τ.

	Inserting (<ref>) and the identity (<ref>) 
	as well as applying a telescopic sum argument, it follows that
	
    2k ∑_i=1^n∑_ℓ =1 ^s_ℓ^i _ℓ^i -
    			_ℓ^i_ℓ^i
       = - 2/s∑_i=1^n∑_ℓ =1 ^s[H]_ℓ^i -
    			^i-1_ℓ^i≤- 1/s∑_ℓ =1 ^s_ℓ^n_H^2 + u_0_H^2.

	Assumption <ref> (5) yields the bound
	
    - 2 ∑_ℓ =1 ^s∫_t^t_nA_ℓ^k(τ)
    			U_ℓ^k(τ) U_ℓ^k(τ)τ≤  2k s λ.

	Therefore, we obtain that
	
    X_3^k(t)
    		   ≤ u_0 _H^2
    		+ 2 ∑_ℓ =1 ^s∫_0^t_nf_ℓ^k(τ) U_ℓ^k(τ)τ + 2k s λ.

	Due to Lemma <ref> and (<ref>), it follows that
	
    lim sup_k → 0 X_3^k(t)
    		≤ u_0 _H^2
    		+ 2 ∫_0^tf(τ) u(τ)τ.

	Thus, we have proved that
	
    lim sup_k → 0( 1/s∑_ℓ = 1^s u(t) - U_ℓ^k(t) 		
    		_H^2 
          + 2∑_ℓ = 1^s∫_0^tA_ℓ^k(τ)u(τ) -
    			A_ℓ^k(τ) U_ℓ^k(τ)  u(τ) - U_ℓ^k(τ)τ)
       ≤  -  u(t) _H^2 +  u_0 _H^2
    		+ 2 ∫_0^tf(τ) - A(τ) u(τ)u(τ)τ
       =  -  u(t) _H^2 +  u_0 _H^2
    		+ 2 ∫_0^tu'(τ)u(τ)τ
       =  -  u(t) _H^2 +  u_0 _H^2
    		+  ∫_0^t/t u(τ) _H^2 τ
    		= 0.

	The monotonicity condition from Assumption <ref> (3) and 
	(<ref>) then imply that
	
    u(t) - U^k(t) _H^2
       ≤1/s∑_ℓ = 1^s u(t) - U_ℓ^k(t) _H^2
    		+ 2η∑_ℓ = 1^s∫_0^t |u(τ) - 
    		U_ℓ^k(τ)|_V_ℓ^p τ
       ≤( 1/s∑_ℓ = 1^s u(t) - U_ℓ^k(t) _H^2 
          + 2∑_ℓ = 1^s∫_0^tA_ℓ^k(τ)u(τ) -
          A_ℓ^k(τ) U_ℓ^k(τ)  u(τ) - U_ℓ^k(τ)τ)
    		→ 0

	as k → 0 for every t ∈ [0,T]. This proves that U^k(t) → u(t) in 
	H as k →0 for every t ∈ [0,T].
	Assuming that η in Assumption <ref> (3) is strictly positive and 
	applying the norm bound from Assumption <ref>,  
	it follows that
	
    u - U_ℓ^k _L^p(0,T;V_ℓ)   ≤ c_V_ℓ(∫_0^T(u(t) - U_ℓ^k(t)_H
    		+ |u(t) - U_ℓ^k(t)|_V_ℓ)^p t)^1/p→ 0

	as k → 0 for every ℓ∈{1,…,s}.



	If η from Assumption <ref> (3) is only strictly positive for 
	some of the operator families {A_ℓ(t)}_t∈[0,T], ℓ∈{1,…,s}, 
	then one can see that in these particular spaces we have U_ℓ^k → u in 
	L^p(0,T;V_ℓ).



10
  
  AdamsFournier.2003
  R. A. Adams and J. J. F. Fournier.
  Sobolev Spaces, volume 140 of Pure and Applied
    Mathematics (Amsterdam).
  Elsevier/Academic Press, Amsterdam, second edition, 2003.
  
  Aronsson.1996
  G. Aronsson, L. C. Evans, and Y. Wu.
  Fast/slow diffusion and growing sandpiles.
  J. Differential Equations, 131(2):304–335, 1996.
  
  ArrarasEtAl.2017
  A. Arrarás, K. J. in 't Hout, W. Hundsdorfer, and L. Portero.
  Modified Douglas splitting methods for reaction-diffusion
  equations.
  BIT, 57(2):261–285, 2017.
  
  Barbu.1976
  V. Barbu.
  Nonlinear semigroups and differential equations in Banach
    spaces.
  Editura Academiei Republicii Socialiste România, Bucharest;
  Noordhoff International Publishing, Leiden, 1976.
  Translated from the Romanian.
  
  DiestelUhl.1977
  J. Diestel and J. J. Uhl, Jr.
  Vector Measures.
  American Mathematical Society, Providence, R.I., 1977.
  
  DrabekEtAl.1997
  P. Drábek, A. Kufner, and F. Nicolosi.
  Quasilinear Elliptic Equations with Degenerations and
    Singularities, volume 5 of De Gruyter Series in Nonlinear Analysis
    and Applications.
  Walter de Gruyter & Co., Berlin, 1997.
  
  Eisenmann.2019
  M. Eisenmann.
  Methods for the Temporal Approximation of Nonlinear,
    Nonautonomous Evolution Equations.
  PhD thesis, TU Berlin, 2019.
  
  EisenmannHansen.2018
  M. Eisenmann and E. Hansen.
  Convergence analysis of domain decomposition based time 
  integrators
  for degenerate parabolic equations.
  Numer. Math., 140(4):913–938, 2018.
  
  Emmrich.2004
  E. Emmrich.
  Gewöhnliche und Operator-Differentialgleichungen: Eine
    integrierte Einfüh­rung in Randwertprobleme und
    Evolutionsgleichungen für Studierende.
  Vieweg-Studium : Mathematik. Vieweg+Teubner Verlag, 2004.
  
  Emmrich.2009b
  E. Emmrich.
  Convergence of the variable two-step BDF time discretisation of
  nonlinear evolution problems governed by a monotone potential operator.
  BIT, 49(2):297–323, 2009.
  
  Emmrich.2009
  E. Emmrich.
  Two-step BDF time discretisation of nonlinear evolution problems
  governed by monotone operators with strongly continuous perturbations.
  Comput. Methods Appl. Math., 9(1):37–62, 2009.
  
  Emmrich.2009c
  E. Emmrich.
  Variable time-step ϑ-scheme for nonlinear evolution
  equations governed by a monotone operator.
  Calcolo, 46(3):187–210, 2009.
  
  EmmrichThalhammer.2010
  E. Emmrich and M. Thalhammer.
  Stiffly accurate Runge-Kutta methods for nonlinear evolution
  problems governed by a monotone operator.
  Math. Comp., 79(270):785–806, 2010.
  
  Evans.1998
  L. C. Evans.
  Partial Differential Equations, volume 19 of Graduate
    Studies in Mathematics.
  American Mathematical Society, Providence, RI, 1998.
  
  GGZ.1974
  H. Gajewski, K. Gröger, and K. Zacharias.
  Nichtlineare Operatorgleichungen und
    Operatordifferentialgleichungen.
  Akademie-Verlag, Berlin, 1974.
  Mathematische Lehrbücher und Monographien, II. Abteilung,
  Mathematische Monographien, Band 38.
  
  HansenHenningsson.2016
  E. Hansen and E. Henningsson.
  Additive domain decomposition operator splittings—convergence
  analyses in a dissipative framework.
  IMA J. Numer. Anal., 37(3):1496–1519, 2017.
  
  HansenOstermann.2008
  E. Hansen and A. Ostermann.
  Dimension splitting for evolution equations.
  Numer. Math., 108(4):557–570, 2008.
  
  HansenStillfjord.2013
  E. Hansen and T. Stillfjord.
  Convergence of the implicit-explicit Euler scheme applied to
  perturbed dissipative evolution equations.
  Math. Comp., 82(284):1975–1985, 2013.
  
  HundsdorferVerwer.2003
  W. Hundsdorfer and J. Verwer.
  Numerical Solution of Time-Dependent
    Advection-Diffusion-Reaction Equations, volume 33 of Springer
    Series in Computational Mathematics.
  Springer-Verlag, Berlin, 2003.
  
  KochEtAl.2013
  O. Koch, Ch. Neuhauser, and M. Thalhammer.
  Embedded exponential operator splitting methods for the time
  integration of nonlinear evolution equations.
  Appl. Numer. Math., 63:14–24, 2013.
  
  Leoni.2009
  G. Leoni.
  A First Course in Sobolev Spaces, volume 105 of Graduate Studies in Mathematics.
  American Mathematical Society, Providence, RI, 2009.
  
  LionsStrauss.1965
  J.-L. Lions and W.A. Strauss.
  Some non-linear evolution equations.
  Bull. Soc. Math. France, 93:43–96, 1965.
  
  Mathew.1998
  T.P. Mathew, P.L. Polyakov, G. Russo, and J. Wang.
  Domain decomposition operator splittings for the solution of
  parabolic equations.
  SIAM J. Sci. Comput., 19(3):912–932, 1998.
  
  Winkert.2018
  N. S. Papageorgiou and P. Winkert.
  Applied Nonlinear Functional Analysis. An Introduction.
  De Gruyter, Berlin, 2018.
  
  Roubicek.2013
  T. Roubíček.
  Nonlinear Partial Differential Equations with
    Applications, volume 153 of International Series of Numerical
    Mathematics.
  Birkhäuser/Springer Basel AG, Basel, second edition, 2013.
  
  Temam.1968
  R. Temam.
  Sur la stabilité et la convergence de la méthode des pas
  fractionnaires.
  Ann. Mat. Pura Appl. (4), 79:191–379, 1968.
  
  Vabishchevich.2008
  P. N. Vabishchevich.
  Domain decomposition methods with overlapping subdomains for the
  time-dependent problems of mathematical physics.
  Comput. Methods Appl. Math., 8(4):393–405, 2008.
  
  Vazquez.2007
  J. L. Vázquez.
  The Porous Medium Equation.
  Oxford Mathematical Monographs. The Clarendon Press, Oxford
  University Press, Oxford, 2007.
  Mathematical theory.
  
  Zeidler.1986
  E. Zeidler.
  Nonlinear Functional Analysis and its Applications. I.
  Springer-Verlag, New York, 1986.
  
  Zeidler.1990B
  E. Zeidler.
  Nonlinear Functional Analysis and its Applications.
    II/B.
  Springer-Verlag, New York, 1990.
  


	
