Identification of a memory kernel in a nonlinear integrodifferential parabolic problem

Applied Numerical Mathematics 120 (2017) 305–323

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Identification of a memory kernel in a nonlinear integrodifferential parabolic problem K. Van Bockstal a,∗ , M. Slodiˇcka a , F. Gistelinck b a

Department of Mathematical Analysis, research group of Numerical Analysis and Mathematical Modeling (NaM2 ), Ghent University, Krijgslaan 281 – S8, Ghent 9000, Belgium b Department of Data Analysis, Ghent University, Henri Dunantlaan 1, Ghent 9000, Belgium

a r t i c l e

i n f o

Article history: Received 31 May 2016 Received in revised form 20 January 2017 Accepted 6 June 2017 Available online 9 June 2017 Keywords: Parabolic IBVP Integral overdetermination Convolution kernel Reconstruction Convergence Time discretization

a b s t r a c t In this contribution, the reconstruction of a solely time-dependent convolution kernel in an nonlinear parabolic equation is studied. The missing kernel is recovered from a global integral measurement. The existence, uniqueness and regularity of a weak solution is addressed. More specific, a numerical algorithm based on Rothe’s method is designed. Numerical experiments support the obtained results. © 2017 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction In this contribution, the following inverse problem for a nonlinear parabolic equation with memory is considered: determine the unknown couple  K , u  obeying

⎧ t ⎪ ⎪ ⎪ ⎪ u ( x , t ) − ∇ · u ( x , t ))) + ( K ∗ u ( x ))( t ) = f (u (x, s)) ds + F (x, t ), ∂ (∇β( ⎪ ⎨ t ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

0

−∇β(u (x, t )) · ν = g (x, t ), u (x, 0) = u 0 (x),

in  × I ,

on  × I ,

(1)

in ,

where  is a bounded Lipschitz domain in Rd , d ≥ 1, with ∂ =  and I = [0, T ], T > 0, is the time frame. The usual

t

convolution in time is denoted by K ∗ u, namely ( K ∗ u (x))(t ) =

K (t − s)u (x, s)ds. The goal of this paper is to develop 0

*

Corresponding author. E-mail addresses: [email protected] (K. Van Bockstal), [email protected] (M. Slodiˇcka), fi[email protected] (F. Gistelinck). URLs: http://cage.ugent.be/~kvb (K. Van Bockstal), http://cage.ugent.be/~ms (M. Slodiˇcka).

http://dx.doi.org/10.1016/j.apnum.2017.06.004 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.

306

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

a numerical scheme for recovering the missing time-convolution kernel K = K (t ). This scheme is based on the additional integral-type measurement



u (x, t ) dx = m(t ),

t ∈ [0, T ].

(2)



The identification of missing memory kernels in partial integrodifferential equations of parabolic type is studied in [7,8, 5,21,6]. The same subject for equations of hyperbolic type is studied in [3,19,18,4,26,4,27]. These papers give some local and global in time existence results for the recovery of solely time-dependent memory kernels in semilinear integrodifferential models. However, these contributions are not focusing on a numerical implementation. The development of a numerical algorithm for the kernel K based on the global measurement (2) in the following semilinear parabolic problem

∂t u (x, t ) − u + K (t )h + ( K ∗ u (x))(t ) = f (u , ∇ u ) with Neumann-type boundary condition is done in [10]. This numerical algorithm is based on a measurement of this equation by using (2). The error analysis of the scheme is performed in [11]. The interested reader is also referred to [45], where the authors studied the same partial differential equation with a convolution of the form K ∗ u instead of K ∗ u and source term f (u ) instead of f (u , ∇ u ). More recently, in [9], using semigroup theory, the authors investigated the reconstruction of an unknown memory kernel from a more general linear integral overdetermination (u (t )) = m(t ) in an abstract linear (of convolution type) evolution equation of parabolic type:

∂t u + Au + K (t )h + ( K ∗ Bu )(t ) = f . In this contribution, the governing equation is nonlinear and the term K (t )h is skipped out of the equation, which complicates the analysis. Crucial in the analysis is that the term depending on the solution u in the right-hand side of (1) is in integral form. Note that other parameter (time-dependent) identification problems can be found in [35,12,13,41,16,30, 37,25,24,31,22,44,23]. For linear β , the type of integro-differential problems under consideration arise in the theory of reactive contaminant transport [14] and in the modeling of phenomena in viscoelasticity [36]. An nonlinear application arises from the modeling of flow in variably saturated rigid soils. The water flow in porous media is usually mathematically described by the Richards equation, a nonlinear, possibly degenerate, parabolic partial differential equation. The authors in [39] derived a non-local Richards equation for the water content in the mobile domain that is characterized by a memory kernel that encodes the local mass transfer dynamics as well as the geometry of the slow zones (cf. equation [19] from [39]). As mentioned before, the main goal of this paper is to design a numerical scheme describing a way of retrieving the couple  K , u . This is achieved not by the minimization of a cost functional (which is typical for IPs) but by the semi-discretization in time by Rothe’s method. The paper is organized as follows. First, the most important notations are introduced in Section 2 and a suitable variational formulation is derived in Section 3. Section 4 deals with the numerical scheme. A priori estimates are derived in Section 5 and the convergence of the proposed numerical scheme is shown in Section 6. The uniqueness of a weak solution is addressed in Theorem 6.2. Afterwards, an error analysis is performed in Section 7. Finally, in Section 8, the theoretically obtained results are illustrated by numerical experiments. 2. Notations Denote by (·, ·) the standard inner product of L2 () and · its induced norm. A similar notation is used when working at the boundary  , namely (·, ·) , L2 () and · . Consider an abstract Banach space X with norm · X . Let k ∈ N ∪ {0}. The set of k-times continuously differentiable functions w : [0, T ] → X equipped with the usual norm k  j =0

 

 

max  w ( j ) 

t ∈[0, T ]

X

is denoted by Ck ([0, T ], X ). The space L p ((0, T ), X ) is furnished with the norm



T 0

p

· X

1

p

with p > 1, cf. [17]. The symbol

X ∗ stands for the dual space to X . Finally, as is usual in papers of this sort, C , ε and C ε denote generic positive constants depending only on a priori known quantities, where ε is small and C = C ε −1 is large. ε

3. Derivation of the variational problem First, the PDE in (1) is multiplied with a test function φ ∈ H1 () and integrated over  to obtain for a.e. t ∈ [0, T ] that

⎛ t ⎞  (∂t u (t ), φ) − (∇ · (∇β(u (t ))) , φ) + (( K ∗ u )(t ), φ) = ⎝ f (u (s)) ds, φ ⎠ + ( F (t ), φ) . 0

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

307

Secondly, using Green’s theorem implies that

⎛ (∂t u (t ), φ) + (∇β(u (t )), ∇φ) + (( K ∗ u )(t ), φ) = ⎝

t

⎞ f (u (s)) ds, φ ⎠ + ( F (t ), φ) − ( g (t ), φ) .

(P)

0



Now, setting φ = 1 in (P) and using the measurement  u (t ) = m(t ), it is clear that

⎞ ⎛ t  m (t ) + ( K ∗ m)(t ) = ⎝ f (u (s)) ds, 1⎠ + ( F (t ), 1) − ( g (t ), 1) . 0

Finally, taking the time derivative of this equality, it holds that





m (t ) + K (t )m(0) + ( K ∗ m )(t ) = ( f (u (t )), 1) + F  (t ), 1 − g  (t ), 1  , where F  (t )

(MP)

and g  (t )

:= ∂t F (t ) := ∂t g (t ). The relations (P) and (MP) represent the variational formulation of the problem (1)–(2). In fact, the inverse problem is reformulated into a direct one. From now on, the real-valued everywhere differentiable function β satisfies the properties

β(0) = 0,

(3)



0 < β0  β (s)  β1 ,

∀s ∈ R,

(4)

i.e. β is strict monotone and Lipschitz continuous. Then β −1 exists and has the following properties

 1  β −1 ( s )  , β1 β0 1

∀s ∈ R.

The potential of β is defined by

z β ( z ) =

z ∈ R.

β(s) ds, 0

Using the properties of β , it is clear that β is convex and (see [42, Lemma 3.1])

0

β0 z 2 2

 β ( z ) 

β1 z 2 2

∀ z ∈ R.

,

(5)

Moreover, via the integral form of β , one can easily verify by the monotone character of β that

β( z1 )( z2 − z1 )  β ( z2 ) − β ( z1 )  β(z2 )( z2 − z1 ),

∀ z1 , z2 ∈ R.

(6)

4. Numerical scheme The well-posedness of (P) and (MP) is studied by using Rothe’s method [28]. This means that a time-discrete scheme based on Backward Euler’s method is designed and the convergence of the approximations towards the unique weak solution is proved under appropriate conditions on the data. Rothe’s method starts with an equidistant time-partitioning of the time frame [0, T ] into n ∈ N intervals [t i −1 , t i ]. The time step is denoted by τ = T /n < 1 and the discrete time points by t i = i τ , i = 0, . . . , n. The notations

zi ≈ z(t i ),

0  i  n,

and

δ zi =

z i − z i −1

τ

,

1  i  n,

are introduced for any function z. Based on (P) and (MP), the following decoupled system for approximating the unknowns ( K , u ) at time t i , 1  i  n, is proposed

(δ u i , φ) + (∇β(u i ), ∇φ) +

 i 



K k u i −k τ , φ

=

k =1

 i −1 



f (uk )τ , φ

+ ( F i , φ) − ( g i , φ)

(DPi)

k =0

and

mi + K i m(0) +

i  k =1





K k mi −k τ = ( f (u i −1 ), 1) + F i , 1 − g i , 1  .

(DMPi)

308

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

Define v i = β(u i ) for i = 1, . . . , n, then equation (DPi) can also be written as  A ( v i ), φ =  F˜ i , φ, with the nonlinear operator A : H1 () → H1 ()∗ and the linear functional F˜ i : H1 () → R defined by

 A ( v i ), φ = and

1

τ

(β −1 ( v i ), φ) + (∇ v i , ∇φ)

(7)

 i −1   i    1 ˜  F i , φ = (u i −1 , φ) + f (uk )τ , φ + ( F i , φ) − ( g i , φ) − K k u i −k τ , φ ,

τ

k =0

(8)

k =1

for any φ ∈ H1 (). For given i ∈ {1, . . . , n}, first (DMPi) is solved and then (DPi). Further, the index i is increased to i + 1. Using monotone operator theory [43], the existence and uniqueness of a solution ( K i , u i = β −1 ( v i )) ∈ R × H1 () for i = 1, . . . , n obeying (DMPi) and (DPi) can be proved. on a single Let f : R → R be bounded. The function β satisfies (3)–(4). Moreover, Theorem 4.1 (Existence of a solution

time step). assume that u 0 ∈ L2 (), g ∈ C1 [0, T ], L2 () , F ∈ C1 [0, T ], L2 () , and m ∈ C2 ([0, T ]) with m(0) = 0. Then there exist C > 0 and τ0 > 0 such that for any τ < τ0 and each i ∈ {1, . . . , n} there exists an unique couple (K i , u i ) ∈ R × H1 () solving (DMPi) and (DPi). Moreover, there exists a positive constant C such that

max | K i |  C .

i =1,...,n





Proof. If m (0) = 0, then m(0) + m (0)τ  > 0. If m (0) = 0, then set



τ0 = min 1,

 |m(0)| . Then for any τ < τ0 , we get by 2 |m (0)|

the triangle inequality that

      m(0) + m (0)τ   |m(0)| − m (0) τ > |m(0)| − m (0) τ0  |m(0)| > 0. 2

For each i ∈ {1, . . . , n}, the following recursive deduction can be made. Let u 0 , . . . , u i −1 ∈ L2 () and K 1 , . . . , K i −1 ∈ R be given. Then, (DMPi) implies the existence of a unique K i ∈ R such that









K i m(0) + m (0)τ = ( f (u i −1 ), 1) + F i , 1 − g i , 1  −

i −1  k =1

K k mi −k τ − mi .

(9)

The operator A defined in (7) is strict monotone because for all u , v ∈ H1 () holds that

   1 −1 1 −1  A (u ) − A ( v ), u − v ) = , 1 u − v 2 1 β (u ) − β ( v ), u − v + (∇ u − ∇ v , ∇ u − ∇ v )  min

τ

β1 τ

H ()

.

From A (0) = 0, it follows that A is coercive. The continuity of β −1 implies the demicontinuity of A. Moreover, the functional F˜ i defined in (8) is bounded. Starting from u 0 ∈ L2 (), the operator equation A ( v i ) = F˜ i has a unique solution v i ∈ H1 () for every i = 1, . . . , n [43]. Therefore, for i = 1, . . . , n, u i = β −1 ( v i ) is uniquely determined in H1 (). Finally, the relation (9) yields



|Ki |  C 1 +

i −1 



|Kk | τ ,

k =1

which is valid for any i = 1, . . . , n. An application of the discrete Grönwall lemma [2] gives the uniform bound of | K i |.

2

5. A priori estimates The a priori estimates proved in this section serve as uniform bounds to prove convergence. Lemma 5.1. Let the conditions of Theorem 4.1 be satisfied. Then there exist positive constants C such that for any τ < τ0 holds that

 

2 max u j  +

1 j n

n 

∇β(u i )2 τ  C

i =1

and





max β(u j )  C

0 j n

and

n  i =1

u i 2H1 () τ  C .

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

309

Proof. We set φ = β(u i )τ in (DPi) and sum up for i = 1, . . . , j, we obtain that j 

(δ u i , β(u i )) τ +

i =1

 i −1  

j 

∇β(u i )2 τ

i =1



j

=

i =1

f (uk )τ , β(u i )

τ+

j 

( F i , β(u i )) τ −

i =1

k =0

j 

( g i , β(u i )) τ −

i =1

 i j   i =1

 K k u i −k τ , β(u i )

τ.

(10)

k =1

For the first term on the left-hand side (LHS), we get that j 

⎛ ⎞     j  ⎝ ⎠ β (u i ) − β (u i −1 ) = β (u j ) − β (u 0 ). (u i − u i −1 , β(u i ))  (6)

i =1

i =1









Moreover,



 β0  u j 2 β ( u j )  (5)

2

and



      (5)  β (u 0 )  β1 u 0 2  C .   2   

The first two terms on the right-hand side (RHS) of (10) can be estimated by the boundedness of f and the Lipschitz continuity of β as

    i −1  j   j      u i 2 τ f ( u ) τ , β( u ) τ C + C  i k    i =1 k =0  i =1

and

    j   j     u i 2 τ . F , β( u )) τ C + C  ( i i    i =1  i =1

The third therm can be estimated by using the Neˇcas inequality [38], i.e.

z2  ε ∇ z2 + C ε z2 ,

∀ z ∈ H1 (),

0 < ε < ε0 .

(11)

Due to the properties of β , it holds true that

   j  j j     2     ∇β(u i )2 τ . g , β( u )) τ C + C u τ + ε  ( ε ε i i   i   i =1  i =1 i =1 By Young’s inequality, Hölder’s inequality, Theorem 4.1 and u 0 ∈ L2 (), we obtain that

    i  i  i   j j j     2  j   1 2     β(u i )2 τ  C + C u i 2 τ . K k u i −k τ , β(u i ) τ   2 Kk τ u i −k τ τ + 12   i =1 k =1  i =1 k =1 i =1 i =1 k =1

Collecting the previous estimates and fixing

ε sufficiently small gives that

j j   2  u j  + ∇β(u i )2 τ  C + C u i 2 τ . i =1

i =1

An application of the discrete Grönwall lemma concludes the proof.

2

Lemma 5.2. Let the conditions of Theorem 4.1 be satisfied and u 0 ∈ H1 (). Then there exist positive constants C such that for any τ < τ0 holds that n 

j  2  δ u i 2 τ + max ∇β(u j ) + ∇β(u i ) − ∇β(u i −1 )2  C

i =1

1 j n

i =1

and

  max u j H1 ()  C

1 j n

and

j  i =1

∇ u i − ∇ u i −1 2  C .

310

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

Proof. First, we put φ = β(u i ) − β(u i −1 ) in (DPi) and sum up for i = 1, . . . , j, we get that j 

(δ u i , β(u i ) − β(u i −1 )) +

i =1

j 

(∇β(u i ), ∇β(u i ) − ∇β(u i −1 )) =

i =1

+

j 

( F i , β(u i ) − β(u i −1 )) −

i =1

 i −1 j   i =1

j 

( g i , β(u i ) − β(u i −1 )) −

i =1

 i j   i =1

 f (uk )τ , β(u i ) − β(u i −1 )

k =0



K k u i −k τ , β(u i ) − β(u i −1 ) .

(12)

k =1

The strict monotonicity of β implies that j 

(δ u i , β(u i ) − β(u i −1 ))  β0

i =1

j 

δ u i 2 τ .

i =1

For the second term on the LHS of (12), we apply Abel’s lemma [20], which states that

2

j 

ai (ai − ai −1 ) = a2j − a20 +

i =1

j  (ai − ai −1 )2 . i =1

Hence (consider ∇β(u i ) instead of ai and integrate over ), j 

j   2 ∇β(u i ) − ∇β(u i −1 )2 . (∇β(u i ), ∇β(u i ) − ∇β(u i −1 )) = 12 ∇β(u j ) − 12 ∇β(u 0 )2 + 12

i =1

i =1

On all the first term on the RHS of (12), we apply the following summation by parts formula: for any real sequences { zi }∞ i =0 and { w i }∞ i =0 holds that j 

z i ( w i − w i −1 ) = z j w j − z 0 w 0 −

i =1

j 

( z i − z i −1 ) w i −1 .

i =1

Therefore,

 i −1 j   i =1

⎛ ⎞ j −1 j   f (uk )τ , β(u i ) − β(u i −1 ) = ⎝ f (uk )τ , β(u j )⎠ − ( f (u i −1 ), β(u i −1 )) τ , 

k =0

k =0 j 



( F i , β(u i ) − β(u i −1 )) = F j , β(u j ) − ( F 0 , β(u 0 )) − (δ F i , β(u i −1 )) τ , j

i =1

 j

i =1

i =1

i =1



( g i , β(u i ) − β(u i −1 )) = g j , β(u j )  − ( g 0 , β(u 0 )) − (δ g i , β(u i −1 )) τ , j

i =1

⎛ ⎞   j j j i i −1      K k u i −k τ , β(u i ) − β(u i −1 ) = ⎝ K k u i −k τ , β(u j )⎠ − K i u0 + K k δ u i −k τ , β(u i −1 ) τ . 



i =1

k =1

k =1

i =1

k =1

Using Lemma 5.1, it is clear that

   j  i −1     f (uk )τ , β(u i ) − β(u i −1 )   C   i =1 k =0 

and

   j     ( F i , β(u i ) − β(u i −1 ))  C .   i =1 

Due to Neˇcas inequality (11), Lemma 5.1 and the assumption u 0 ∈ H1 (), we readily obtain that

    2   j   .    g , β( u ) − β( u )) ( i i i −1    C ε + ε ∇β(u j )   i =1 

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

311

For the last term on the RHS of (12), it is clear by Lemma 4.1, Lemma 5.1 and Hölder’s inequality that

  2   j  i j  i −1          K u τ , β( u ) − β( u ) C + C K δ u τ    τ i i − 1 k i − k k i − k      i =1 k =1  i =1 k =1     j i −1 i −1     2 2 δ u i −k  τ τ C +C K τ k

i =1

 j

C +C



i =1

Combining the previous estimates and fixing j 

k =1 i 

k =1



2

δ uk  τ τ .

k =1

ε sufficiently small gives that

 i  j j    2  2 2 δ u i  τ + ∇β(u j ) + ∇β(u i ) − ∇β(u i −1 )  C + C δ uk  τ τ . 2

i =1

i =1

i =1

An application of the discrete Grönwall lemma concludes the proof.

k =1

2

In the following lemma, the PDE in (1) needs to be satisfied at t = 0. Hence, the following compatibility conditions are defined:

δ u 0 := ∂t u (0) := F 0 + ∇ · (∇β(u 0 )) ∈ L2 ()

if

∇β(u 0 ) ∈ H(div; )

(13)

and

u −1 := u 0 − τ δ u 0 ∈ L2 (). Now, the discrete measured problem (DMPi) is valid at t = 0. It holds that





m0 + K 0 m0 = ( f (u −1 ), 1) + F 0 , 1 − g 0 , 1  ,

(DMP0)

which defines K 0 .



C [0, T ], L2 () , m ∈ C3 ([0, T ]) and f : R → R is global Lipschitz continuous. Then, there exist positive constants C and τ0 such that for all τ < τ0 holds that Lemma 5.3. Let the assumptions of Lemma 5.2 be fulfilled. Moreover, assume that ∇β(u 0 ) ∈ H(div; ), g ∈ C2 [0, T ], L2 () , F ∈

2

n 

|δ K i |2 τ  C .

i =1

Proof. First, lets subtract (DMP0) from (DMPi) for i = 1 and divide the result by





δ K 1m0 = ( f (u 0 ) − f (u −1 ), 1) + δ F 1 , 1 − δ g 1 , 1



τ . We obtain

− δm1 − K 1m0 .

By the mean value theorem and the Lipschitz continuity of f , we get for

τ < τ0 holds that

|δ K 1 |  C δ u 0  τ + C  C , as m(0) = 0. Next, we apply the δ -operator to (DMPi) for 2  i  n. This gives us i −1 



δ K i m(0) = (δ f (u i −1 ), 1) + δ F i , 1 − δ g i , 1  − δmi − K i m0 − K k δmi −k τ . k =1

Theorem 4.1 and the Lipschitz continuity of f lead to

|δ K i |  C (1 + δ u i −1 ) Therefore, for n  i =1

for

2  i  n.

τ < τ0 , Lemma 5.2 implies that 2

2

|δ K i | τ = |δ K 1 | τ +

n  i =2

2



|δ K i | τ  C 1 +

n  i =2

 2

δ u i −1  τ

 C.

2

312

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

6. Convergence and existence of a solution Let us introduce the following piecewise linear function in time

 2

un : [0, T ] → L () : t →

u0

t=0

u i −1 + (t − t i −1 )δ u i

t ∈ (t i −1 , t i ]

,

1  i  n,

and a step function

 2

un : [0, T ] → L () : t →

u0

t=0

ui

t ∈ (t i −1 , t i ]

1  i  n.

,

Similarly, define K n , F n , F  n , g n , g  n , m n and m n . Using these so-called Rothe’s functions, (DPi) and (DMPi) can be rewritten on the whole time frame as1

⎛ ⎞ t τ  K n (tk )un (t − tk )τ , φ ⎠ (∂t un (t ), φ) + ∇β(un (t )), ∇φ + ⎝

k =1

⎞ ⎛ t τ 



f (un (tk ))τ , φ ⎠ + F n (t ), φ − g n (t ), φ  =⎝

(DP)

k =0

and

m n (t ) + K n (t )m(0) +

t τ 









K n (tk )m n (t − tk )τ = f (un (t − τ )), 1 + F  n (t ), 1 − g  n (t ), 1

k =1

 

.

(DMP)

This puts us in a position to prove the existence of a weak solution to (P) and (MP). Theorem 6.1 (Existence). Suppose the conditions of Lemma 5.2 are fulfilled. Moreover, let f : R → R be global Lipschitz continuous. Then there exists a weak solution  K , u  to the problem (P)–(MP), where K ∈ L2 (0, T ) and u ∈ C [0, T ], L2 () ∩ L2 (0, T ), H1 ()



with ∂t u ∈ L2 (0, T ), L2 () . Proof. From Lemmas 5.1 and 5.2, we have for all n > 0 that





max un (t )H1 ()  C

t ∈[0, T ]

T ∂t un (t )2 dt  C .

for all t ∈ [0, T ], 0

The Rellich–Kondrachov theorem [15, §5.8.1] implies the compact embedding

H1 () → → L2 (). Then, the conditions of [28, Lemma 1.3.13] are satisfied. This implies the existence of an element









u ∈ C [0, T ], L2 () ∩ L∞ (0, T ), H1 () , and a subsequence (unl )k∈N of (un )n∈N such that

⎧ ⎪ unl → u , in ⎪ ⎪ ⎪ ⎨ u (t ) u (t ), in nl

⎪ unl (t ) u (t ), in ⎪ ⎪ ⎪ ⎩ ∂t unl ∂t u , in



C [0, T ], L2 () , H1 () for all t ∈ [0, T ], H1 () for all t ∈ [0, T ],



(14a)

L2 (0, T ), L2 () ,

which we denote again by un for ease of reading. Note that un (0) − un (0) = 0. For all t ∈ (t i −1 , t i ] with 1  i  n, we have that

  un (t ) − un (t ) = |u i −1 + (t − t i −1 )δ u i − u i | = |(t − t i −1 − τ )δ u i | = |(t − t i )δ u i |  τ |δ u i | .

1

t τ = i and t τ = i − 1 when t ∈ (t i −1 , t i ].

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

313

Employing Lemma 5.2 gives





2 lim un − un L2 (0, T ),L2 ()  lim

n→∞

n→∞

n 

τ2

δ u i 2 τ  lim

n→∞

i =1

C n2

= 0,



such that un and un have the same limit in L2 (0, T ), L2 () , i.e.

un → u





L2 (0, T ), L2 () .

in

(15)

By [29, p. 88], the exists a subsequence {u n } of {un } such that

un → u ,

(0, T ) × .

a.e. in

Therefore, by the Lipschitz continuity of β , we get that





in L2 (0, T ), L2 () .

β(un ) → β(u )

(16)

From Lemma 5.2 follows that

T

  β(un (t ))2 1

H ()

dt  C .

0

such that β(un ) z in L (0, T ), H1 () . Due to (16), we obtain that z = β(u ). Therefore,

By the reflexivity of the space L2 (0, T ), H1 () , we get the existence of a subsequence (again denoted with the same index)

2





∇β(un ) ∇β(u )

L2 (0, T ), L2 () .

in

T Using Theorem 4.1, we have that

(17)

   K n (t )2 d t  C , which means by the reflexivity of L2 (0, T ) that (in the sense of a

0

subsequence being indexed by n as well)

Kn K

in L2 (0, T ). m

(18)

  n (t ) = m (t ) and limn→∞ m n (t )

= m (t ) in L2 ([0, T ]), lim n→∞ g n (t ) = g (t ) and limn→∞ g  n (t ) = and limn→∞ F  n (t ) = F  (t ) in L2 (0, T ), L2 () because m, F and g are

It is clear that limn→∞ g  (t ) in L2 (0, T ), L2 () and limn→∞ F n (t ) = F (t ) prescribed. Now, we integrate (DP) in time over (0, η) ⊂ [0, T ] for the resulting subsequence to get

η

η (∂t un (t ), φ) dt +

0



∇β(un (t )), ∇φ dt +

0

⎛ ⎞ η  t τ ⎝ K n (tk )un (t − tk )τ , φ ⎠ dt 0

k =1

⎛ ⎞ η  η η t τ

⎝ ⎠ = f (un (tk ))τ , φ dt + F n (t ), φ dt − g n (t ), φ  dt . 0

k =0

0

This expression is valid for any have for n → ∞ that

η

(19)

0

η ∈ [0, T ]. We want to pass the limit n → ∞ in (19). Using the stability result (14a), we

η (∂t un , φ) →

0

(∂t u , φ) ,

∀φ ∈ H1 ().

0

Thanks to (17), we have that

η



η

∇β(un (t )), ∇φ dt →

0

(∇β(u (t )), ∇φ) dt ,

∀φ ∈ H1 ().

0

By Theorem 4.1 and Lemma 5.1, it holds for all t ∈ [0, T ] that t τ  k =1

τt τ

K n (tk )un (t − tk )τ = ( K n ∗ un )(t ) +

K n (s)un (t − s) ds = ( K n ∗ un )(t ) + O (τ ) t

(20)

314

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

and t τ 

t f (un (tk ))τ = f (u 0 )τ +

k =0

t f (un (s)) ds −

t f (un (s)) ds =

τ t τ

0

f (un (s)) ds + O (τ ) .

(21)

0

Hence, by (15), (18) and the Lipschitz continuity of f , we see that

⎛ ⎞ η  η  η t τ t τ

⎝ ⎠ lim K n (tk )un (t − tk )τ , φ dt = lim K n (tk ) un (t − tk ), φ τ = (( K ∗ u )(t ), φ) dt

n→∞

and

n→∞

k =1

0

0 k =1

0

⎛ ⎛ ⎞ ⎞ η  η  t t τ ⎝ lim f (un (tk ))τ , φ ⎠ dt = ⎝ f (u (s)) ds, φ ⎠ ds.

n→∞

k =0

0

0

0

Now, taking the limit n → ∞ in (19) results in

η

η (∂t u (t ), φ) dt +

0

η (∇β(u (t )), ∇φ) dt +

0

(( K ∗ u )(t ), φ) dt 0

⎛ ⎞ η  t η η = ⎝ f (u (s)) ds, φ ⎠ dt + ( F (t ), φ) dt − ( g (t ), φ) dt . 0

0

0

0

Taking the derivative with respect to η , we arrive at (P). In the same way as before, we integrate (DMP) in time and pass the limit for n → ∞. This follows the same line as passing the limit in (19), therefore we skip the details. Note that

T lim

n→∞

  un (t − τ ) − un (t )2 d t = 0.

0

Finally, we differentiate the result with respect to time and arrive at (MP).

2

In the previous theorem, to prove the existence of a weak solution to (P)–(MP), only the weak convergence of K n to K in L2 (0, T ) is used. However, assume that the conditions of Lemma 5.3 are satisfied. Then, from this lemma, it is clear that

T max | K n (t )|  C

t ∈[0, T ]

for all t in [0, T ]

|∂t K n (t )|2 dt  C .

and 0

From this follows that K n and K n have the same limit in L2 (0, T ). Moreover, by the Arzelà–Ascoli theorem [40, Theorem 11.28], there exists a subsequence { K n } (which we denote by the same symbol again) that converges uniformly on [0, T ] to K , i.e. K ∈ C([0, T ]). The reflexivity of the space L2 (0, T ) implies that ∂t K n ∂t K in L2 (0, T ). Corollary 6.1. Suppose the conditions of Lemma 5.3 solution  K , u  to the problem (P)–(MP),

are fulfilled. Then there exists a weak where K ∈ C([0, T ]) with K  ∈ L2 (0, T ) and u ∈ C [0, T ], L2 () ∩ L2 (0, T ), H1 () with ∂t u ∈ L2 (0, T ), L2 () . Now, it is possible to establish the uniqueness of a solution to (P)–(MP). Theorem 6.2 (Uniqueness). Let the assumptions of Lemma (P)–(MP) has at most one solution 

 6.1 be satisfied.

Then the problem  K , u  ∈ L2 (0, T ) × C [0, T ], L2 () ∩ L2 (0, T ), H1 () with ∂t u ∈ L2 (0, T ), L2 () . Proof. First, we integrate (P) over the time variable t ∈ (0, ξ ) ⊂ (0, T ). We obtain that

⎛

(u (ξ ) − u 0 , φ) + ⎝

ξ

⎞

⎛ ξ ⎞  ∇β(u (t ))dt , ∇φ ⎠ + ⎝ ( K ∗ u )(t )dt , φ ⎠

0

0

⎞ ⎛ ξ ⎞ ⎛ ξ ⎞ ⎛ ξ t     f (u (s)) dsdt , φ ⎠ + ⎝ F (t )dt , φ ⎠ − ⎝ g (t )dt , φ ⎠ . =⎝ 0

0

0

0



(22)

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

315

Suppose that there are two solutions  K 1 , u 1  and  K 2 , u 2  solving (P)–(MP). Denote the difference of the solutions by K = K 1 − K 2 and u = u 1 − u 2 . Subtracting the integrated variational formulation (22) for the different solution gives that

⎛ ξ ⎞ ⎛ ξ ⎞ ⎛ ξ ⎞    (u (ξ ), φ) + ⎝ [∇β(u 1 (t )) − ∇β(u 2 (t ))] dt , ∇φ ⎠ + ⎝ ( K 1 ∗ u )(t )dt , φ ⎠ + ⎝ ( K ∗ u 2 )(t )dt , φ ⎠ 0

0

⎞ ⎛ ξ t   =⎝ [ f (u 1 (s)) − f (u 2 (s))] dsdt , φ ⎠ . 0

0

0

Now, we put φ = β(u 1 (ξ )) − β(u 2 (ξ )) and integrate again over the time variable ξ ∈ (0, η) ⊂ (0, T ). This gives

η 0

⎛ ⎞ η ξ (u (ξ ), β(u 1 (ξ )) − β(u 2 (ξ ))) dξ + ⎝ [∇β(u 1 (t )) − ∇β(u 2 (t ))] dt , ∇β(u 1 (ξ )) − ∇β(u 2 (ξ ))⎠ dξ 0

0

⎛ ⎛ ⎞ ⎞ η ξ η ξ + ⎝ ( K 1 ∗ u )(t )dt , β(u 1 (ξ )) − β(u 2 (ξ ))⎠ dξ + ⎝ ( K ∗ u 2 )(t )dt , β(u 1 (ξ )) − β(u 2 (ξ ))⎠ dξ 0

0

0

0

⎛ ⎞ η ξ  t = ⎝ [ f (u 1 (s)) − f (u 2 (s))] dsdt , β(u 1 (ξ )) − β(u 2 (ξ ))⎠ dξ. 0

0

(23)

0

Subtracting (MP) for the different solutions gives

K (t )m(0) + ( K ∗ m )(t ) = ( f (u 1 (t )) − f (u 2 (t )), 1) . By the Lipschitz continuity of f and Grönwall’s lemma, it follows easily that

| K (t )|  C u (t ) ,

t ∈ [0, T ].

(24)

Now, we are able to estimate each term of (23). The strict monotonicity of β implies that

η

η u (ξ )2 dξ.

(u (ξ ), β(u 1 (ξ )) − β(u 2 (ξ ))) dξ  β0 0

0

For the second term on the LHS, we get by the strict monotonicity of β and integration by parts that

 η 2 ⎛ ⎞   η ξ  1 ⎝ [∇β(u 1 (t )) − ∇β(u 2 (t ))] dt , ∇β(u 1 (ξ )) − ∇β(u 2 (ξ ))⎠ dξ =  [∇β(u 1 (ξ )) − ∇β(u 2 (ξ ))] dξ   0.   2  0

0

0

By Young’s, Hölder’s and Jensen’s inequality, the Lipschitz continuity of β and K 1 ∈ L2 (0, T ), it follows for the third term on the LHS that

 η⎛ ξ ⎞     η ξ η   2  ⎝ ( K 1 ∗ u )(t )dt , β(u 1 (ξ )) − β(u 2 (ξ ))⎠ dξ   C ε u (t ) dt dξ + ε u (ξ )2 dξ.     0

0

0

0

0

For the last term on the LHS, we use the fact that u 2 ∈ C [0, T ], L () and (24) such that

2

 η⎛ ξ ⎞     η ξ η   2  ⎝ ( K ∗ u 2 )(t )dt , β(u 1 (ξ )) − β(u 2 (ξ ))⎠ dξ   C ε u (t ) dt dξ + ε u (ξ )2 dξ.     0

0

0

0

0

The same upper bound can be obtained for the term on the RHS of (23) by the Lipschitz continuity of f . Collecting all the estimates above and fixing ε sufficiently small imply that

η

η ξ u (ξ )2 dξ  C

0

u (t )2 dt dξ. 0

0

An application of Grönwall’s lemma gives that K in L2 (0, T ) follows from (24). 2

η 0

u (ξ )2 dξ = 0, i.e. u = 0 or u 1 = u 2 a.e. in  × [0, T ]. The uniqueness of

316

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

The convergence of Rothe’s functions towards the weak solution have been shown for a subsequence {un } and { K n } in Theorem 6.1. However, taking into account the uniqueness of a solution, it is clear that the whole Rothe’s sequences {un } and { K n } converge against the solution. 7. Error analysis For the error analysis, the assumptions of Lemma 5.3 needs to be satisfied. First, the following notations are introduced:

e u = un − u ,

e u = un − u ,

eK = K n − K .

Analogously, e g , e g  , e F , e F  , em and em are defined. Secondly, the following estimates are frequently used

    un (t − τ ) − un (t ) + un (t ) − un (t )  ∂t un (t ) τ ,

∀ t ∈ [0, T ].

Note that this inequality is also valid for t = 0 by compatibility condition (13). Based on the result of Lemma 5.2, this yields

T

     un (t − τ ) − un (t )2 + un (t ) − un (t )2 dt  C τ 2 .

(25)

0

Finally, the following theorem contains the error estimates. Optimal convergences rates for Rothe’s method (O (τ )) are obtained. Theorem 7.1. Let the conditions of Lemma 5.3 be fulfilled. Then there exist positive constants C and τ0 such that for all τ < τ0 holds that

T

   K n (t ) − K (t )2 dt +

0

T

  un (t ) − u (t )2  C τ 2 .

0

Proof. The proof exists out of two steps. First, a bound on e K in L2 (0, T ) is proved. Afterwards, the error estimates are established. For the first step, we subtract (MP) from (DMP) and we use an analogue of (20) to get τt τ

em (t ) + e K (t )m0 +







K n (s)m n (t − s) ds + e K ∗ m n (t ) + K ∗ em (t )

t





= f (un (t − τ )) − f (u (t )), 1 + e F  (t ), 1 − e g  (t ), 1 .

(26)



By Theorem 4.1 and the mean-value theorem, we obtain that

 e

  τ t     τ      

 

     +  K ∗ e  (t ) +  e  (t ), 1  +  e  (t ), 1   C τ .  ( t ) + K ( s ) m ( t − s ) ds  n n m m F g      t

For the fourth term on the LHS of (26), it is true that

t          e K ∗ m n (t )  C e K (s) ds. 0

The Lipschitz continuity of f implies for the first term on the RHS of (26) that

  

  f (un (t − τ )) − f (u (t )), 1   C un (t − τ ) − un (t ) + C e u (t ) .

Gathering all the estimates, we obtain from (26) (after division by m0 ) that

⎛ ⎞ t       e (t )  C ⎝τ + un (t − τ ) − un (t ) + e u (t ) + e (s) ds⎠ . K K 0

We take the second power and integrate the result with respect to the time variable. For each

η 0

⎛

  e (t )2 dt  C ⎝τ 2 + K

η 0

  un (t − τ ) − un (t )2 dt +

η

η  t e u (t )2 dt +

0

0

0

η ∈ [0, T ] holds that ⎞

  e (s)2 dsdt ⎠ . K

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

317

Applying estimate (25) and Grönwall’s argument imply that

η

⎛ ⎞ η 2  e (t ) dt  C ⎝τ 2 + e u (t )2 dt ⎠ . K

0

(27)

0

The second part of the proof starts with the subtraction of (P) from (DP). Next, the result is integrated in time over (0, η) ⊂

(0, T ). We obtain

⎛ ξ ⎞ ⎛ ξ⎡ ⎤ ⎞    t τ K n (tk )un (t − tk )τ − ( K ∗ u )(t )⎦ dt , φ ⎠ (un (ξ ) − u (ξ ), φ) + ⎝ [∇β(un (t )) − ∇β(u (t ))] dt , ∇φ ⎠ + ⎝ ⎣ 0

k =1

0

⎤ ⎞ ⎛ ξ ⎞ ⎛ ξ ⎞ ⎛ ξ⎡   t   t τ f (un (tk ))τ − f (u (s)) ds⎦ dt , φ ⎠ + ⎝ ( F n (t ) − F (t ))dt , φ ⎠ − ⎝ ( g n (t ) − g (t ))dt , φ ⎠ . =⎝ ⎣ k =0

0

0

0

0



Now, we put φ = β(un (ξ )) − β(u (ξ )) and integrate again over the time variable ξ ∈ (0, η) ⊂ (0, T ). This gives

η



⎛ ⎞ η ξ ⎝ [∇β(un (t )) − ∇β(u (t ))] dt , ∇β(un (ξ )) − ∇β(u (ξ ))⎠ dξ

un (ξ ) − u (ξ ), β(un (ξ )) − β(u (ξ )) dξ +

0

+

⎛ ⎡ η ξ  t τ ⎝

0

=

⎣

0

0

⎣

t f (un (tk ))τ −

k =0

0

0

⎤

⎞

K n (tk )un (t − tk )τ − ( K ∗ u )(t )⎦ dt , β(un (ξ )) − β(u (ξ ))⎠ dξ

k =1

⎛ ⎡ η ξ  t τ ⎝

0

⎤

⎞

f (u (s)) ds⎦ dt , β(un (ξ )) − β(u (ξ ))⎠ dξ

0

⎛ η ξ

⎛ ⎞ η ξ ⎝ ( F n (t ) − F (t ))dt , β(un (ξ )) − β(u (ξ ))⎠ dξ − ⎝ ( g n (t ) − g (t ))dt , β(un (ξ )) − β(u (ξ ))⎠ dξ.

+ 0

⎞

0

0

0



After some rearrangements in the terms and using (20) and (21), it holds that

η



un (ξ ) − u (ξ ), β(un (ξ )) − β(u (ξ )) dξ +

0

η



un (ξ ) − un (ξ ), β(un (ξ )) − β(u (ξ )) dξ

0

 η 2 ⎛ ⎛ ⎞ ⎞   η ξ τt τ    ⎝ ⎝ K n (s)un (t − s) ds⎠ dt , β(un (ξ )) − β(u (ξ ))⎠ dξ + 12   [∇β(un (t )) − ∇β(u (t ))] dt  +   0

0

t

0

⎛ ⎛ ⎞ ⎞ η ξ η ξ + ⎝ (e K ∗ un )(t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ + ⎝ ( K ∗ e u )(t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ 0

0

⎛ ⎡ η ξ ⎝

= 0

⎣ f (u 0 )τ +

τt τ

0

⎤

0

⎞

f (un (s)) ds⎦ dt , β(un (ξ )) − β(u (ξ ))⎠ dξ

t

0

⎛ ⎡ ⎤ ⎞ η ξ  t

+ ⎝ ⎣ f (un (s)) − f (u (s)) ds⎦ dt , β(un (ξ )) − β(u (ξ ))⎠ dξ 0

0

⎛ η ξ ⎝

+ 0

0

0

⎛ ⎞ η ξ e F (t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ − ⎝ e g (t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ. ⎞

0

0



(28)

318

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

By the strict monotonicity of β , we get that

η



un (ξ ) − u (ξ ), β(un (ξ )) − β(u (ξ )) dξ  β0

0

η e u (ξ )2 dξ. 0

For the second term on the LHS, by the Cauchy and Young inequalities, the Lipschitz continuity of β and (25), we obtain that  η    η

   





un (ξ ) − un (ξ ), β(un (ξ )) − β(u (ξ )) dξ   C ε τ 2 + ε



0

e u (ξ )2 dξ. 0

Using the Cauchy and Young inequalities, the uniform boundedness of K n and un and the Lipschitz continuity of β , it follows that

 η ⎛ ξ ⎛ τ t  ⎞ ⎞      τ η   2  ⎝ ⎝  ⎠ ⎠ K n (s)un (t − s) ds dt , β(un (ξ )) − β(u (ξ )) dξ   C ε τ + ε e u (t )2 dt .    0

t

0

0

For the fifth term on the LHS of (28), we use Young’s, Hölder’s and Jensen’s inequality together with the previously obtained estimate (27). We obtain that

 η⎛ ξ ⎞        ⎝ (e ∗ un )(t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ  K     0

0

 2  η ξ η   2   Cε   (e K ∗ un )(t )dt  dξ + ε e u (ξ ) dξ   0 0 0 ⎛ ⎞ η ξ  t η 2 ⎝ |e (s)| ds⎠ dt dξ + ε e u (ξ )2 dξ  Cε K

0

0

0

0

η ξ

(27)

 Cετ 2 + Cε

η e u (s)2 ds dξ + ε

0

0

e u (ξ )2 dξ. 0

Moreover, in an analogous way, it holds that

 η⎛ ξ ⎞     η ξ η   2  ⎝ ( K ∗ e u )(t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ   C ε e u (s) ds dξ + ε e u (ξ )2 dξ.     0

0

0

0

0

Now, we estimate the terms in the RHS of (28). The first term can be estimated as

 η⎛ ξ ⎡ ⎤ ⎞  τt τ    η   2  ⎝ ⎣ f (u 0 )τ +  ⎦ ⎠ f (un (s)) ds dt , β(un (ξ )) − β(u (ξ )) dξ   C ε τ + ε e u (ξ )2 dξ.    0

t

0

0

By the Lipschitz continuity of f , we calculate that

 η⎛ ξ ⎡ t ⎤ ⎞       

 ⎝ ⎣ ⎦ dt , β(un (ξ )) − β(u (ξ ))⎠ dξ  f ( u ( s )) − f ( u ( s )) ds n     0

0

0

η ξ

η 2

e u (ξ )2 dξ.

e u (s) ds dξ + ε

 Cε 0

0

0

For the third term on the RHS, we use the mean value theorem such that

 η⎛ ξ ⎞     η   2  ⎝ e (t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ   C ε τ + ε e u (ξ )2 dξ. F     0

0

0

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

319

We use the integration by parts formula for the last term on the RHS, i.e.

⎛ ⎞ ⎞ ⎛ η η ξ  η ⎝ e g (t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ = ⎝ e g (t )dt , β(un (t )) − β(u (t ))dt ⎠ 0

0



⎛ η

ξ

0

0

⎝e g (ξ ),

−

0

⎞

0



β(un (t )) − β(u (t ))dt ⎠ dξ. 

We employ Cauchy’s and Young’s inequalities, the mean value theorem, the trace inequality and the Lipschitz continuity of β such that

 η⎛ ξ  ⎞    η   2  ⎝ e g (t )dt , β(un (ξ )) − β(u (ξ ))⎠ dξ   C ε τ + ε e u (t )2 dt     0

0

0



 η 2  2    η ξ η ξ     2  +C  [∇β(un (t )) − ∇β(u (t ))] dt  dξ.  dt +ε ∇β( u ( t )) − ∇β( u ( t )) e ( s ) ds d ξ + C [ ] n u         0

0

0

0

0

Collecting all the previous estimates give

η e u (ξ )2 dξ +

(β0 − ε ) 0

1 2

2  η       dt −ε  ∇β( u ( t )) − ∇β( u ( t )) [ ] n     0

η ξ  Cετ 2 + C 0

Fixing



0

 2  η ξ   2  e u (s) ds dξ + C  [∇β(un (t )) − ∇β(u (t ))] dt   dξ.   0

0

ε sufficiently small, an application of Grönwall’s lemma concludes the proof. 2

Corollary 7.1. Let the conditions of Lemma 5.3 be fulfilled. Then there exist positive constants C and τ0 such that for all τ < τ0 holds that

T | K n (t ) − K (t )|2 dt  C τ 2 . 0

8. Numerical experiments In this section, it is assumed that the assumptions of Theorem 7.1 are satisfied. The presented numerical algorithm in the previous sections is summarized in the following pseudo code: Algorithm: Numerical scheme in pseudo code. input : T > 0, n ∈ N and functions F , F  , f , β , g, g  , m, m , m and u 0 output: kernel K and solution u at discrete time steps 1 τ ← T /n; 2 θ ← [0 : τ : T ]; 3 K ← zeros(n + 1); 4 u[0] ← u 0 ;





1 ( f (u 0 − τ δ u 0 ), 1) + F 0 , 1 − g0 , 1  − m0 ; m0 6 for i = 1 to n do 5 K[0] ←

1



 i −1 



K k mi −k τ − mi ; ( f (u i −1 ), 1) + F i , 1 − g i , 1  −

7

K[i ] ←

8

u[i ] ← solveEP( A (u i ), φ =  F˜ i , φ;

m0 + m0 τ

k=1

The aim of the simulations is to demonstrated the established error estimate on the convolution kernel K in Corollary 7.1. For these simulations, the finite element library DOLFIN [32,34] from the FEniCS project [33,1] is used.

320

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

In each experiment, the domain  = [0, 1]. The number of time discretization interval is chosen to be n = 2 j , j = 2, . . . , 8, such that the time step τ for the equidistant time partitioning equals respectively 2− j T , j = 2, . . . , 8. At each time-step, the resulting nonlinear elliptic problems (see Line 8 in the Algorithm) are solved numerically by the finite element method (FEM) using first order (P1-FEM) Lagrange polynomials for the space discretization. At each timestep, the nonlinearity ∇β(u i ) is approximated by β  (u i −1 )∇ u i . For the space discretization, a fixed uniform mesh of 100 intervals is used. The L2 -error between the numerical and exact kernel is approximated by Simpson’s rule for the several values of the timestep τ :

T | K n (t ) − K ex (t )|2 dt

E K ex (τ ) = 0

≈

n  τ i =1

6



 2

( K i −1 − K ex (t i −1 )) +

K i −1 + K i 2

 − K ex

t i −1 + t i 2

!

2

+ ( K i − K ex (t i ))

2

.

8.1. Experiment 1 The exact solution in the first experiment is prescribed as follows



u ex (x, t ) = 1 + t + t 2



1 + x2



and

K ex (t ) = exp(t ),

x ∈ [0, 1],

t ∈ [0, 1].

(29)

The functions f and β are given by f (s) = s + 1 and β(s) = s + 1. Note that β is linear. The exact kernel K ex is compared with the numerical solution for τ = 2−3 , 2−5 and τ = 2−7 in Fig. A.1(a). The errors E K ex (τ ) are depicted in Fig. A.1(b), where the error log2 E K ex is plotted as a function of log2 τ . The linear regression line through all data points is given by log2 E K ex = 2.0233 log2 τ − 0.9666. This is in accordance with the predicted convergence rate in Corollary 7.1. The linear regression line for the relative error is given by log2 E rK ex = 2.0233 log2 τ − 2.6422. For instance, the CPU time for τ = 2−8 is 8 s. 8.2. Experiment 2 In the second experiment the unknown kernel is sinusoidal, i.e.





u ex (x, t ) = 1 + t + t 2 + t 3 (1 + sin (π x))

and

K ex (t ) = sin(2π t ),

" x ∈ [0, 1],

t ∈ 0,

1 2

# .

(30)

The functions f and β are given by respectively f (s) = s + 5 and β(s) = s2 + s. The function β can be linearized outside the range of u ex . The CPU time for τ = 2−8 is 8 s. The results of the numerical experiment are depicted in Figs. A.2(a)–A.2(b). The linear regression line through the data points is given by log2 E K ex = 2.0282 log2 τ + 0.0166. The linear regression line for the relative error is given by log2 E rK ex = 2.0282 log2 τ + 2.0166. The CPU time for τ = 2−8 is 3 s. 8.3. Experiment 3 In the third experiment, the prescribed solution u is more general, i.e.





u ex (x, t ) = 1 + t + t 2 + t 3 (1 + sin (tx))

and

K ex (t ) = sin(2π t ),

x ∈ [1, 2],

t ∈ [0, 1] .

(31)

The functions f and β are the same as in Experiment 2. Again, accurate numerical approximations are obtained when the timestep is sufficiently small, see Figs. A.3(a)–A.3(b). The linear regression line through all data points is given by log2 E K ex = 2.0142 log2 τ + 1.8951. The linear regression line for the relative error is given by log2 E rK ex = 2.0142 log2 τ + 2.8951. The CPU time for

τ = 2−8 is 10 s.

9. Conclusion A nonlinear parabolic integro-differential problem of second order with an unknown solely time-dependent convolution kernel is considered. The missing information is compensated by an integral-type measurement over the domain. The existence and uniqueness of a weak solution for the problem is proved using Rothe’s method. This convergence of the numerical approximations to the exact solution is optimal in time. Numerical experiments support this theoretically obtained results. For instance, future research can concern the detemination of K when a local measurement  χ (x)u (x, t )dx = m(t ) is considered and the implementation of numerical experiments with noisy data. Acknowledgement The research was supported by the IAP P7/02-project of the Belgian Science Policy Office.

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

Appendix A. Figures

Fig. A.1. Kernel reconstruction in Experiment 1.

Fig. A.2. Kernel reconstruction in Experiment 2.

Fig. A.3. Kernel reconstruction in Experiment 3.

321

322

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

References [1] M.S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M.E. Rognes, G.N. Wells, The FEniCS project version 1.5, Arch. Numer. Softw. 3 (100) (2015) 9–23. [2] D. Bainov, P. Simeonov, Integral Inequalities and Applications, Mathematics and Its Applications. East European Series, vol. 57, Kluwer Academic Publishers, Dordrecht, 1992. [3] A.L. Bukhgeˇım, G.V. Dyatlov, Inverse problems for equations with memory, in: Inverse Problems and Related Topics, Proceedings of a Seminar, Kobe Institute, Kobe, Japan, February 6–10, 1998, Chapman & Hall/CRC, Boca Raton, FL, 2000, pp. 19–35. [4] C. Cavaterra, M. Grasselli, On an inverse problem for a model of linear viscoelastic Kirchhoff plate, J. Integral Equ. Appl. 9 (3) (1997) 179–218. [5] F. Colombo, D. Guidetti, A global in time existence and uniqueness result for a semilinear integrodifferential parabolic inverse problem in Sobolev spaces, Math. Models Methods Appl. Sci. 17 (4) (2007) 537–565. [6] F. Colombo, D. Guidetti, Some results on the identification of memory kernels, in: M. Ruzhansky, J. Wirth (Eds.), Modern Aspects of the Theory of Partial Differential Equations. Including Mainly Selected Papers Based on the Presentations at the 7th International ISAAC Congress, London, UK, July 13–18, 2009, in: Oper. Theory, Adv. Appl., vol. 216, Birkhäuser, Basel, 2011, pp. 121–138. [7] F. Colombo, D. Guidetti, A. Lorenzi, On applications of maximal regularity to inverse problems for integrodifferential equations of parabolic type, in: Gisèle Ruiz Goldstein, et al. (Eds.), Evolution Equations: Proceedings of the Conference in Honor of the 60th Birthdays of Philippe Bénilan, Jerome A. Goldstein and Rainer Nagel, Blaubeuren, Germany, June 11–17, 2001, in: Lect. Notes Pure Appl. Math., vol. 234, Marcel Dekker, New York, NY, 2003, pp. 77–89. [8] F. Colombo, D. Guidetti, V. Vespri, Some global in time results for integrodifferential parabolic inverse problems, in: Angelo Favini, et al. (Eds.), Differential Equations. Inverse and Direct Problems. Papers of the Meeting, Cortona, Italy, June 21–25, 2004, in: Lecture Notes in Pure and Applied Mathematics, vol. 251, CRC Press, Boca Raton, FL, 2006, pp. 35–58. [9] R.H. De Staelen Guidetti, On a finite difference scheme for an inverse integro-differential problem using semigroup theory: a functional analytic approach, Numer. Funct. Anal. Optim. 37 (7) (2016) 850–886. [10] R.H. De Staelen, M. Slodiˇcka, Reconstruction of a convolution kernel in a semilinear parabolic problem based on a global measurement, Nonlinear Anal. 112 (2015) 43–57. [11] R.H. De Staelen, K. Van Bockstal, M. Slodiˇcka, Error analysis in the reconstruction of a convolution kernel in a semilinear parabolic problem with integral overdetermination, J. Comput. Appl. Math. 275 (2015) 382–391. [12] M. Dehghan, Parameter determination in a partial differential equation from the overspecified data, Math. Comput. Model. 41 (2–3) (2005) 196–213. [13] M. Dehghan, Identification of a time-dependent coefficient in a partial differential equation subject to an extra measurement, Numer. Methods Partial Differ. Equ. 21 (3) (2005) 611–622. [14] J.W. Delleur, The Handbook of Groundwater Engineering, CRS Press, 1999. [15] L.C. Evans, Partial Differential Equations, Graduate Studies in Mathematics, vol. 19, American Mathematical Society, USA, Providence, RI, 1998. [16] U.M. Fedus’, Identification of the coefficient of the time derivative in a quasilinear parabolic equation, J. Math. Sci. 168 (4) (2010) 523–543. [17] H. Gajewski, K. Gröger, K. Zacharias, Nichtlineare Operatorgleichungen und Operatordifferentialgleichungen, Mathematische Lehrbücher und Monographien. II. Abteilung, vol. 38, Akademie-Verlag, Berlin, 1974. [18] M. Grasselli, An identification problem for an abstract linear hyperbolic integrodifferential equation with applications, J. Math. Anal. Appl. 171 (1) (1992) 27–60. [19] M. Grasselli, S. Kabanikhin, A. Lorenzi, An inverse hyperbolic integrodifferential problem arising in geophysics. I, Sib. Math. J. 33 (3) (1992) 415–426. [20] R.B. Guenther, J.W. Lee, Partial Differential Equations of Mathematical Physics and Integral Equations, Dover Books on Mathematics, Dover Publications, United States of America, 1988. [21] D. Guidetti, Convergence to a stationary state for solutions to parabolic inverse problems of reconstruction of convolution kernels, Differ. Integral Equ. 20 (9) (2007) 961–990. [22] A. Hazanee, D. Lesnic, Determination of a time-dependent coefficient in the bioheat equation, Int. J. Mech. Sci. 88 (0) (2014) 259–266. [23] M. Hussein, D. Lesnic, Identification of the time-dependent conductivity of an inhomogeneous diffusive material, Appl. Math. Comput. 269 (2015) 35–58. [24] M. Ismailov, F. Kanca, The inverse problem of finding the time-dependent diffusion coefficient of the heat equation from integral overdetermination data, Inverse Probl. Sci. Eng. 20 (2012) 463–476. [25] M. Ismailov, F. Kanca, D. Lesnic, Determination of a time-dependent heat source under nonlocal boundary and integral overdetermination conditions, Appl. Math. Comput. 218 (8) (2011) 4138–4146. [26] J. Janno, Discretization and regularization of an inverse problem related to a quasilinear hyperbolic integrodifferential equation, Inverse Probl. 13 (3) (1997) 711. [27] J. Janno, E. Kiss, L. von Wolfersdorf, On Tikhonov regularization for identifying memory kernels in heat conduction and viscoelasticity, Z. Angew. Math. Mech. 80 (4) (2000) 259–272. [28] J. Kaˇcur, Method of Rothe in Evolution Equations, Teubner Texte zur Mathematik, vol. 80, Teubner, Leipzig, 1985. [29] A. Kufner, O. John, S. Fuˇcík, Function Spaces, Monographs and Textbooks on Mechanics of Solids and Fluids, Noordhoff International Publishing, Leyden, 1977. [30] M. Lakestani, M. Dehghan, The use of Chebyshev cardinal functions for the solution of a partial differential equation with an unknown time-dependent coefficient subject to an extra measurement, J. Comput. Appl. Math. 235 (3) (2010) 669–678. [31] D. Lesnic, S.A. Yousefi, M. Ivanchov, Determination of a time-dependent diffusivity from nonlocal conditions, J. Appl. Math. Comput. 41 (2013) 301–320. [32] A. Logg, G.N. Wells, DOLFIN: automated finite element computing, ACM Trans. Math. Softw. 37 (2) (2010) 28. [33] A. Logg, K.-A. Mardal, G.N. Wells, et al., Automated Solution of Differential Equations by the Finite Element Method, Springer, Berlin, Heidelberg, 2012. [34] A. Logg, G. Wells, J. Hake, DOLFIN: A C++/Python Finite Element Library, Springer, Berlin, Heidelberg, 2012, Ch. 10. [35] A. Lorenzi, Determination of a time-dependent coefficient in a quasi-linear parabolic equation, Ric. Mat. 32 (2) (1983) 263–284. [36] R.C. MacCamy, A model for one-dimensional, nonlinear viscoelasticity, Q. Appl. Math. 35 (1977) 21–33. [37] A. Mohebbi, M. Dehghan, High-order scheme for determination of a control parameter in an inverse problem from the over-specified data, Comput. Phys. Commun. 181 (12) (2010) 1947–1954. [38] J. Neˇcas, Les méthodes directes en théorie des équations elliptiques (Direct Methods in the Theory of Elliptic Equations), Academia, Prague, 1967. [39] I. Neuweiler, D. Erdal, M. Dentz, A non-local Richards equation to model unsaturated flow in highly heterogeneous media under nonequilibrium pressure conditions, Vadose Zone J. 11 (3) (2012), http://dx.doi.org/10.2136/vzj2011.0132. [40] W. Rudin, Real and Complex Analysis, McGraw-Hill, 1987. [41] A. Saadatmandi, M. Dehghan, Computation of two time-dependent coefficients in a parabolic partial differential equation subject to additional specifications, Int. J. Comput. Math. 87 (5) (2010) 997–1008. [42] M. Slodiˇcka, S. Dehilis, A nonlinear parabolic equation with a nonlocal boundary term, J. Comput. Appl. Math. 233 (12) (2010) 3130–3138, finite Element Methods in Engineering and Science (FEMTEC 2009).

K. Van Bockstal et al. / Applied Numerical Mathematics 120 (2017) 305–323

323

[43] M.M. Vainberg, Variational Method and Method of Monotone Operators in the Theory of Nonlinear Equations, translated from Russian by A. Libin, translation edited by D. Louvish, A Halsted Press Book, John Wiley & Sons, New York–Toronto, 1973; Jerusalem–London: Israel Program for Scientific Translations. xi, 356 p. (1973). [44] K. Van Bockstal, M. Slodiˇcka, Determination of a time-dependent diffusivity in a nonlinear parabolic problem, Inverse Probl. Sci. Eng. 23 (2) (2015) 307–330. [45] K. Van Bockstal, R.H. De Staelen, M. Slodiˇcka, Identification of a memory kernel in a semilinear integrodifferential parabolic problem with applications in heat conduction with memory, J. Comput. Appl. Math. 289 (0) (2015) 196–207.