Two-grid economical algorithms for parabolic integro-differential equations with nonlinear memory

Accepted Manuscript Two-grid economical algorithms for parabolic integro-differential equations with nonlinear memory

Wansheng Wang, Qingguo Hong

PII: DOI: Reference:

S0168-9274(19)30038-8 https://doi.org/10.1016/j.apnum.2019.02.001 APNUM 3503

To appear in:

Applied Numerical Mathematics

Received date: Revised date: Accepted date:

18 June 2018 31 January 2019 8 February 2019

Please cite this article in press as: W. Wang, Q. Hong, Two-grid economical algorithms for parabolic integro-differential equations with nonlinear memory, Appl. Numer. Math. (2019), https://doi.org/10.1016/j.apnum.2019.02.001

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Two-grid economical algorithms for parabolic integro-differential equations with nonlinear memory Wansheng Wanga , Qingguo Hongb,∗ a

Department of Mathematics, Shanghai Normal University, Shanghai, 200234, China b Department of Mathematics, Pennsylvania State University, State College, PA 16802,U.S.A.

Abstract In this paper, several two-grid finite element algorithms for solving parabolic integro-differential equations (PIDEs) with nonlinear memory are presented. Analysis of these algorithms is given assuming a fully implicit time discretization. It is shown that these algorithms are as stable as the standard fully discrete finite element algorithm, and can achieve the same accuracy as the standard algorithm if the coarse grid size H and the fine grid size h satisfy r−1 H = O(h r ). Especially for PIDEs with nonlinear memory defined by a lower order nonlinear operator, our two-grid algorithm can save significant storage and computing time. Numerical experiments are given to confirm the theoretical results. Keywords: parabolic integro-differential equation, two-grid method, error estimate, finite element method, stability, backward Euler scheme 2008 MSC: 65M60, 65R20, 65L05 1. Introduction The main purpose of this paper is to present some discretization techniques based on two finite element subspaces for solving parabolic integro∗

Corresponding author Email addresses: [email protected] (Wansheng Wang), [email protected] (Qingguo Hong)

Preprint submitted to Applied Numerical Mathematics

February 19, 2019

differential equations (PIDEs) with nonlinear memory:  t K(t − s)Bu(s)ds = f (x, t), (x, t) ∈ Ω × (0, T ], (1.1) ut + Au + 0

u(x, t) = 0, (x, t) ∈ ∂Ω × (0, T ], u(x, 0) = u0 (x), x ∈ Ω,

(1.2) (1.3)

where Ω ⊂ Rd (d ≥ 1) is a bounded and polyhedral domain with a piecewise smooth boundary ∂Ω, K(t) is a smooth or nonsmooth memory kernel, and f is a known function, u is the solution we need to solve which is scalar function . A is a symmetric positive definite second-order elliptic operator with smooth coefficients in x and t, and B is a nonlinear operator of at most second order; that is, Bu = −∇ · (α(x, u)∇u + β(x, u)) + γ(x, u) · ∇u + g(x, u).

(1.4)

For brevity, we will drop the dependence of variable x in α(x, u), β(x, u), γ(x, u), and g(x, u) in the following exposition. We assume that the functions α(u) ∈ Rd×d is a tensor function, β(u)) ∈ Rd , γ(u) ∈ Rd are vector functions, and g(u) ∈ R1 is scalar function, respectively. And all the functions α(u), β(u), γ(u) and g(u) are smooth and bounded together with the Gateaux derivative. For the functions β(u) and g(u), we also assume that β(0) = 0 and g(0) = 0. Equations of the above type, or linear versions thereof, can arise from many physical processes in which it is necessary to take into account the effects of memory due to the deficiency of the usual diffusion equations [20, 33, 39]. For approximating the solution u of PIDEs, both finite difference and finite element methods have been investigated extensively in the past for both the linear and nonlinear problem (see, for example, [8, 9, 29, 31, 37, 12, 54]). Recently, several new numerical methods such as mixed finite element method, finite volume element method, and discontinuous Galerkin method for space discretization or time discretization have been proposed to solve PIDEs (see, for example, [19, 42, 36, 41, 6, 35]). The two-grid method based on two finite element spaces, one on a coarse grid and one on a fine grid, was first developed by Xu [47, 48, 49, 50] for nonsymmetric linear and nonlinear elliptic problems. Since then, the two-grid method for elliptic problems has been investigated further, e.g., Axelsson and Layton [3], Xu and Zhou [51], Li and Huang [28], and Bi and Ginting [4, 5]. 2

In these works, theoretical study and numerical experiments show that the combined use of the numerical method such as finite element method and finite difference method, and the two-grid technique is computationally more efficient than the original method. Due to this better practical performance, the two-grid method has been widely applied to the study of eigenvalue problems [52, 53, 24], steady Navier Stokes equations [27, 21, 23, 15], the timedependent Navier Stokes problem [22, 1, 2, 40, 44], the nonlinear parabolic problem [16, 32, 17, 10, 46, 14, 13, 38], and nonlinear hyperbolic equations [11]. Recently, Jin, Shu, and Xu [26] used this technique to solve decoupling systems of partial differential equations; Mu and Xu [34] and Cai, Mu, and Xu [7] employed it for the mixed Stokes-Darcy model. In [45], we proposed the two-grid algorithms based on the backward Euler scheme and finite element approximation for semi-linear PIDEs, and studied the long-time stability and error estimates of the two-grid algorithms. In this paper, we present some two-grid algorithms for PIDEs with nonlinear memory and perform theoretical analysis that demonstrates our methods’ ability to match the accuracy of the classic finite element method by (1) solving a nonlinear problem on a coarse space SH and (2) solving a symmetric positive definite linear problem on the fine space Sh . Thus, solving PIDEs with nonlinear memory is not much more difficult than solving one linear problem, as dim SH  dim Sh and the work involved in solving the nonlinear problem on the coarse grid is relatively limited. It is worth adding that when α ≡ 0, our algorithm significantly reduces computational memory and storage requirements. A practical difficulty of numerical methods for PIDEs is that all previous values must be stored, as they all enter subsequent equations. In order to reduce memory requirements, some economical schemes have been proposed (for example, see, [43, 25]). However, these schemes either require more regularities on the solution u [43], or they cannot be applied to nonlinear problem [25]. The remainder of this article is organized as follows: In Section 2, we present some conventions and notations that will be used throughout the article. In Section 3, the stability and error estimate of the classic fully discrete finite element method are discussed. The two-grid algorithms for PIDEs with nonlinear memory are presented and the stability and error estimates of these algorithms are discussed in Section 4. In Section 5, we offer some concluding remarks. Throughout this paper, we use the letters C and c (with and without subscripts) to denote a generic positive constant that stand for different values 3

depending on the context in different equations. When it is not important to keep track of these constants, we conceal the letter C or c in the notation  or , such that x  y means x ≤ Cy and x  y means x ≥ cy. 2. Preliminaries r,p (Ω) be the stanFor any non-negative integer r and number p ≥ 1, let W p dard Sobolev space with a norm · r,p given by v r,p = |κ|≤r Dκ v pLp (Ω) (with the usual modification if p = ∞). This Sobolev space is also equipped with the seminorm |v|pr,p = |κ|=r Dκ v pLp (Ω) . For p = 2, we denote Hr = W r,2 (Ω) and take H01 as the subspace of H1 consisting of functions with a vanishing trace on ∂Ω. For simplicity, we also use notations · r , · and · ∞ , and | · |r such that · r = · r,2 , · = · 0,2 and · ∞ = · 0,∞ , and | · |r = | · |r,2 . Let {Sh }0
inf { v − χ + h v − χ 1 }  hr v r ,

χ∈Sh

d r ≥ 1 + . (2.1) 2

v ∈ Hr ∩ H01 ,

We also assume that {Sh }0 0 independent of h such that ∇χ ∞ ≤ Ch−d/2 ∇χ ,

χ ∈ Sh .

(2.2)

The weak formulation of the problem (1.1), (1.3) is: Find u ∈ H01 (Ω) such that  t (ut , v) + A(u, v) + K(t − s)B(u(s), v)ds = (f, v), v ∈ H01 , (2.3) 0

u(0) = u0 ,

(2.4)

where A(·, ·) is the bilinear form associated with the operator A on H01 × H01 and B(·, ·) is defined by B(u, v) = (α(u)∇u + β(u), ∇v) + (γ(u) · ∇u + g(u), v),

u, v ∈ W 1,∞ ∩ H01 .

(·, ·) denotes the inner product in L2 (Ω). We always assume that A is coercive and continuous with coercivity constant ν0 and continuity constant ν1 . That is, we have A(v, v) ≥ ν0 v 21 ∀v ∈ H01 , ∀u, v ∈ H01 . |A(u, v)| ≤ ν1 u 1 v 1 4

(2.5) (2.6)

In view of the assumptions on the functions α(u), β(u), γ(u), and g(u), it is easily verified that there exists a positive constant μ0 such that |B(u, v)| ≤ μ0 u 1 v 1 .

(2.7)

For the time discretization of (1.1)-(1.3) we will consider the backward Euler scheme. To analyze the discretization on a time interval (0, T ], let N be a positive integer, Δt = T /N , and let tn = nΔt. As the truncation error of the backward Euler scheme is O(Δt), we introduce a quadrature formula with a truncation error O(Δt),  tn n  ωni g(ti ) = K(tn − s)g(s)ds + O(Δt). (2.8) Δt i=1

0

Given our emphasis on two-grid discretization in space, we will not discuss how to obtain the numbers ωni , but only assume that there exists a positive constant K1 such that |ωni | ≤ K1 for any 1 ≤ n ≤ N, 1 ≤ i ≤ n and that ωnn = 0. Therefore, the problem considered in this paper must be discretized by a fully implicit scheme. Thus, the backward Euler fully discrete finite element approximation of problem (1.1), (1.3) is defined as a sequence {U n }N n=0 , such that 

n   n n ¯ ωni B(U i , v) = (f n , v), v ∈ Sh , n ≥ (2.9) 1, ∂U , v + A(U , v) + Δt i=1 0

U =

uh0 ,

(2.10)

¯ n = U n −U n−1 , uh is an appropriate approximation of u0 in Sh , where ∂U 0 Δt f n = f (tn ). We know that (2.9) will result in a truncation error O(Δt) in time. But for nonlinear problems considered in this paper (ωnn = 0), the solution of a nonlinear algebraic system is required at each time step. To decrease the amount of computational work, we propose using a two-grid technique to solve the PIDEs with nonlinear memory. With this technique, at each time step, solving a nonlinear problem on the fine space Sh is reduced by solving a nonlinear problem on the coarse space SH and solving a linear SPD problem on the fine space Sh . For functions that vanish on the boundary, we recall Poincare’s inequality: there exists a constant P such that ∀v ∈ H01 ,

v ≤ P|v|1 . 5

We make extensive use of the −type inequality 2ab ≤ a2 + b2 / , > 0, and of the inequality a2 + b2 ≤ (|a| + |b|)2 . The results of this paper are based on the identity 2(an+1 , an+1 − an ) = |an+1 |2 − |an |2 + |an+1 − an |2 ,

(2.11)

and the following Gronwall lemma proved in [18]. Lemma 2.1 (Discrete Gronwall lemma [18]) Let 0 ≤ λ < 1, and an , bn , cn , λn ≥ 0 with {cn } being monotonically increasing. Then an + bn ≤

n−1 

n = ,  + 1, · · ·

λj aj + λan + cn ,

(2.12)

j=

implies for n = ,  + 1, · · · n−1  cn λj cn  an + bn ≤ ≤ exp 1+ 1 − λ j= 1−λ 1−λ



n−1 1  λj 1 − λ j=

 .

3. Error estimate for the classic fully discrete finite element method In this section, we discuss the stability and error estimate of the standard fully discrete finite element method (2.9), (2.10). First, we prove the stability of the solution of (2.9) and (2.10). Theorem 3.1 Let U n be the solution obtained by (2.9) and (2.10). Then for all

1 7ν02 , (3.1) Δt ≤ min , 2 8μ20 K12 T we have

U n +

≤ En1/2

n 

1/2 U i − U i−1 2

i=1

U 0 2 + Δt

n 

√

ν0 + 2

1/2 f i 2

i=1 2

where En = 6 max{e2tn , e(2μ0 K1 tn /ν0 ) }. 6

,



n 

1/2 Δt U i 21

i=1

(3.2)

Proof. By taking v = 2ΔtU n in (2.9) and using (2.11), we obtain n 2

U − U

n−1 2

+ U − U n

n−1 2

+

2ν0 Δt U n 21

+ 2(Δt)

2

n  i=1

≤ 2Δt f U . n

ωni B(U i , U n )

n

(3.3)

Using (2.7), we have U n 2 − U n−1 2 + U n − U n−1 2 + 2ν0 Δt U n 21 n    1 2 ≤ μ0 (Δt) |ωni |( U i 21 + U n 21 ) + Δt U n 2 + f n 2 . (3.4) i=1 Choose = ν0 /(μ0 K1 tn ) to obtain U n 2 + U n − U n−1 2 + ν0 Δt U n 21 n  μ20 K12 tn n−1 2 2 ≤ U + (Δt) U i 21 + Δt f n 2 + Δt U n 2 . (3.5) ν0 i=1 By summation, we have n 2

U +

n 

U − U i

i−1 2

+ ν0 Δt

i=1 0 2

≤ U + Δt

n 

U i 21

i=1

n 

i 2

U + (Δt)

i=1

2

n  i=1

i n  μ20 K12 ti  j 2 U 1 + Δt f i(3.6) 2 , ν0 j=1 i=1

which implies that n  μ20 K12 tn U − U + ν0 − Δt Δt U i 21 (1 − Δt) U + ν 0 i=1 i=1



n−1 i n 2 2    4μ ν K t i 0 U i 2 + Δt ≤ U 0 2 + Δt max 1, 0 2 1 U j 21 + Δt f i 2 . ν 4 0 i=1 j=1 i=1 n 2

n 



i

i−1 2

(3.7) Since condition (3.1) implies that 1 − Δt ≥

7

1 2

and ν0 −

μ20 K12 tn Δt ν0

≥

ν0 , 8

with

the aid of discrete Gronwall lemma 2.1, we obtain n 2

U + 2

n 

U − U i

i=1 0 2

U + Δt

≤ En

n Δtν0  i 2 + U 1 4 i=1 

i−1 2

n 

f i 2

,

(3.8)

i=1

which implies (3.2). Thus the proof is completed. Remark. From (3.1), we find that for a given integral interval (0, T ] the stepsize Δt is determined by the ratio of ν0 to μ0 and increases as the value of coercivity constant ν0 increases. ¯ n in (2.9) to obtain Now let us take v = 2Δt∂U ¯ n 2 + A(U n , U n ) − A(U n−1 , U n−1 ) + A(U n − U n−1 , U n − U n−1 ) 2Δt ∂U n  2 ¯ n) +2(Δt) ωni B(U i , ∂U i=1

¯ n ). = 2Δt(f , ∂U n

(3.9)

Since ¯ n ) ≤ 1 Δt f n 2 + 2Δt ∂U ¯ n 2 2Δt(f n , ∂U 2

(3.10)

and 2(Δt)

2

n 

¯ n )| ≤ 2Δtμ0 |ωni B(U , ∂U i

i=1

n 

|ωni | U i 1 U n − U n−1 1

i=1

≤

tn μ20 K12 ν0

Δt

n 

U i 21 + ν0 U n − U n−1(3.11) 21 ,

i=1

(3.9) becomes ν0 U n 21

1 tn μ20 K12  i 2 Δt f n 2 + ≤ Δt U 1 + ν1 U n−1 21 . (3.12) 2 ν0 i=1 n

Then we have the following result.

8

Theorem 3.2 Let U n be the solution obtained by (2.9) and (2.10). Then for all Δt ≤

ν02 , 2μ20 K12 T

(3.13)

we have  U 1 ≤ C n

U 0 21

+ Δt

n 

i 2

f

1/2 .

(3.14)

i=1 μ2 K 2 t

Proof. It follows from (3.13) that ν0 − 0 ν01 n Δt ≥ ν20 . Then an application of discrete Gronwall lemma 2.1 to (3.12) leads to 

n  f i 2 , U n 21 ≤ C U 0 21 + Δt i=1

which implies (3.14). This completes the proof. To estimate the error of the fully discrete approximation (2.9), we define, for w, u, v ∈ W 1,∞ ∩ H01 (Ω), B1 (w; u, v) = (α(w)∇u, ∇v) + (γ(w) · ∇u, v). Due to the assumptions on α(u) and γ(u), there exist a constant σ such that |B1 (w; u, v)| ≤ σ u 1 v 1 .

(3.15)

As usual, we write the error en = u(tn ) − U n as en = u(tn ) − U n = (u(tn ) − Vh u(tn )) + (Vh u(tn ) − U n ) = ρn + θn , where Vh u is the Ritz-Volterra projection of the solution u and defined by [9]  t K(t − s)B1 (u(s); u(s) − Vh u(s), v)ds = 0, v ∈ Sh(3.16) . A(u − Vh u, v) + 0

For ρ(t) = u(t) − Vh u(t), following the line of Cannon and Lin [9], we show that there exists C0 > 0, independent of h and t, such that (see, also, [12, 30]) t ≥ 0, ρ(t) + h ρ(t) 1 ≤ C0 hr |||u(t)|||r , r ρt (t) ≤ C0 h (|||u(t)|||r + |||ut (t)|||r ) , ρ(t) ∞ ≤ C0 hr | ln h||||u(t)|||r,∞ , 9

(3.17) (3.18) (3.19)

where  |||u(t)|||r = u(t) r +

t 0

 u(τ ) r dτ,

|||u(t)|||r,∞ = u(t) r,∞ +

t 0

u(τ ) r,∞ dτ,

and there exists a positive constant C = C(u), independent of h, such that ∇Vh u ∞ + ∇(Vh u)t ∞ ≤ C.

(3.20)

Now we need to estimate the error θn = Vh u(tn ) − U n . Theorem 3.3 Let u and U n be the solutions of (2.3)-(2.4) and (2.9)-(2.10), respectively. If 16σ 2 K12 T < ν02 ,

(3.21)

then, for sufficiently small Δt, we have  θn + θn − θn−1 + ν0 Δt θn 1  hr + Δt.

(3.22)

Proof. Firstly, it follows from (2.3) and (2.9) that   ¯ n , v + A(u − U n , v) ut − ∂U  tn n  K(t − s)B(u(s), v)ds − Δt ωni B(U i , v) = 0, + 0

v ∈ Sh .

i=1

Then we find that θn satisfies  tn  n n n ¯ ∂θ , v + A(θ , v) + A(ρ , v) + K(t − s)B1 (u(s); u(s) − Vh u(s), v)ds 0  tn + K(t − s)B1 (u(s); Vh u(s), v)ds



0



tn

+ 0

 = −

K(t − s)[(β(u(s)), ∇v) + (g(u(s)), v)]ds − Δt

ρ −ρ Δt n

n−1





, v − ut −



u(tn ) − u(tn−1 ) ,v , Δt

10

n 

ωni B(U i , v)

i=1

v ∈ Sh .

(3.23)

Using (3.16), we have 





tn

¯ , v + A(θ , v) + ∂θ n

−Δt +Δt

i=1 n 

0

K(t − s)B1 (u(s); Vh u(s), v)ds

ωni B1 (u(ti ); Vh u(ti ), v) + Δt

n 

ωni B1 (u(ti ); θi , v)

i=1

ωni

      α(u(ti )) − α(U i ) ∇U i , ∇v + γ(u(ti )) − γ(U i ) · ∇U i , v

i=1 tn

 + 0

−Δt +Δt  = −

n 

n

K(t − s)[(β(u(s)), ∇v) + (g(u(s)), v)]ds

n  i=1 n 

ωni [(β(u(ti )), ∇v) + (g(u(ti )), v)] ωni



   β(u(ti )) − β(U i ), ∇v + g(u(ti )) − g(U i ), v

i=1  u(tn ) − u(tn−1 ) ρn − ρn−1 , v − ut − ,v , Δt Δt

v ∈ Sh .

(3.24)

Now, in view of (2.8), we have   n   tn    K(t − s)B1 (u(s); Vh u(s), v)ds − Δt ωni B1 (u(ti ); Vh u(ti ), v)  Δt v 1 .    0 i=1

and    

tn

K(t − s)[(β(u(s)), ∇v) + (g(u(s)), v)]ds  n    −Δt ωni [(β(u(ti )), ∇v) + (g(u(ti )), v)]  Δt v 1 .  0

(3.26)

i=1

Due to (3.15), the fifth term on the left-hand side in (3.24) can be bounded as   n n       i ωni B1 (u(ti ); θ , v) ≤ ΔtK1 σ θi 1 v 1 . (3.27) Δt   i=1

i=1

11

(3.25)

By virtue of the assumptions on α(u), β(u), γ(u) and g(u), we know α, β, γ and g satisfy Lipschitz conditions with Lipschitz constant CL , and thus the sixth and ninth terms on the left-hand side in (3.24) are estimated as follows:   n            i i i i ω )) − α(U ) ∇U , ∇v + γ(u(t )) − γ(U ) · ∇U , v α(u(t Δt   ni i i   i=1   n          ωni α(u(ti )) − α(U i ) ∇θi , ∇v + γ(u(ti )) − γ(U i ) · ∇θi , v  ≤ Δt    i=1 n       ωni α(u(ti )) − α(U i ) ∇Vh u(ti ), ∇v Δt  i=1   + γ(u(ti )) − γ(U i ) · ∇Vh u(ti ), v  n n   i ≤ ΔtK1 σ θ 1 v 1 + CL ΔtK1 u(ti ) − U i ∇Vh u(ti ) ∞ v 1 i=1

≤ ΔtK1

n 

i=1 n 

σ θi 1 v 1 + CCL ΔtK1

i=1

( ρi + θi ) v 1 ,

(3.28)

i=1

where the estimate (3.20) has been used, and   n          i i ωni β(u(ti )) − β(U ), ∇v + g(u(ti )) − g(U ), v  Δt   i=1

≤ CL ΔtK1

n 

u(ti ) − U i v 1 ≤ CL ΔtK1

i=1

n 

( ρi + θi ) v (3.29) 1.

i=1

The first term on the right-hand side in (3.24) can be bounded as  n   tn  ρ − ρn−1  1 1 n n−1   ρt (s) ds v ,v  ≤ ρ − ρ v ≤  Δt Δt Δt tn−1  C0 h r t n ( |u(s) |r + |ut (s) |r ) ds v ; (3.30) ≤ Δt tn−1 and the last term can be bounded as    tj    ut − u(tn ) − u(tn−1 ) , v  ≤ utt ds v .   Δt tj−1 12

(3.31)

Taking v = 2Δtθn and substituting all the above estimates (3.25)-(3.31) into (3.24), we obtain θn 2 − θn−1 2 + θn − θn−1 2 + 2ν0 Δt θn 21 n  2 n 2 r n 2 ≤ CΔt θ 1 + C(Δt + Δth ) θ + 4σΔt K1 θi 1 θn 1 +4CCL Δt2 K1

n 

i=1

( ρi + θi ) θn 1

i=1 2

≤ CΔt + CΔt

2

θn 21

+2 1 σΔtK1 tn θn 21 2 + CCL Δt2 K1 2

n 

n  1 i 2 + C(Δt + h ) + CΔt θ + 2σΔt K1 θ 1 1 i=1 r 2

+ 2CCL ΔtK1

2

n 

n 2

2

ρi 2 + 2CCL Δt2 K1 tn θn 21

i=1

θi 2 + 2 2 CCL ΔtK1 tn θn 21 ,

(3.32)

i=1

where we have used the inequality CΔt

2

n  i=1

ai bn = C

n 

Δt1/2 ai Δt3/2 bn

i=1

   C C C C Δt ≤ a2i + Δt3 b2n ≤ Δt a2i + Δt2 tn b2n . 2 2 2 2 i=1 i=1 i=1 n

n

Using the estimate (3.17) for ρi , and taking 1 = we have

n

ν0 4σK1 tn

and 2 =

ν0 , 4CCL K1 tn

θn 2 − θn−1 2 + θn − θn−1 2 + ν0 Δt θn 21  8 2 2 2 r 2 ≤ C(Δt + h ) + C + C CL K1 tn Δt2 θn 2 ν0  n−1 8 2 2 8 2 2 2  i 2 2 n 2 + C + σ K1 tn + 2CCL K1 tn Δt θ 1 + σ Δt K1 tn θ 1 ν0 ν0 i=1  8 + C 2 CL2 Δt2 K12 tn θi 2 , ν0 i=1 n−1

(3.33)

13

Noting the condition (3.21) and taking sufficiently small Δt such that   8 2 2 8 2 2 2 ν0 1 C + σ K1 tn + 2CCL K1 tn Δt ≤ and C + C CL K1 tn Δt ≤ , ν0 2 ν0 2 we obtain θn 2 + θn − θn−1 2 + ν0 Δt θn 21 n−1 1 2 i 2 1 n−1 2 ≤ θ + Δt θ + Δt θn 2 + C(Δt + hr )2 2 2 i=1  ν0 ν0 + Δt θn 21 + Δt2 θi 21 . 2 2 i=1 n−1

(3.34)

Applying discrete Gronwall lemma 2.1 to the above inequality yields θn 2 + θn − θn−1 2 + ν0 Δt θn 21  (hr + Δt)2 .

(3.35)

which implies (3.22). This proves the theorem. Note that the condition (3.21), which implies that the equation (1.1) is diffusion-dominant, is appropriate, since the system may be blowup if the integral term is dominant. Under the condition (3.21), we can not study the long time behaviour of the numerical solution. Of course, if we assume that there exist positive constants α0 , α1 > 0 such that α0 |ξ|2 ≤ ξ T α(u)ξ ≤ α1 |ξ|2 ,

∀u ∈ R,

ξ ∈ Rd ,

(3.36)

then following the approach of [45], we can study the long time behavior of the exact solution and the numerical solution. We now give the H1 estimate of the error θn . Theorem 3.4 Let u and U n be the solutions of (2.3)-(2.4) and (2.9)-(2.10), respectively. Then, for all Δt satisfying ν02 , 16σ 2 K12 T

(3.37)

θn 1  hr + Δt.

(3.38)

Δt < we have

14

¯ n in (3.24), and estimating every terms in a Proof. Taking v = 2Δt∂θ way similar to Theorem 3.3, we get ¯ n 2 + A(θn , θn ) − A(θn−1 , θn−1 ) + A(θn − θn−1 , θn − θn−1 ) 2Δt ∂θ n  2 ¯ n 2 r n 2 ¯ ¯ n 1 θi 1 ∂θ ≤ CΔt ∂θ 1 + C(Δt + Δth ) ∂θ + 4σΔt K1 +4CCL Δt2 K1 2

n 

i=1

¯ n 1 ( ρi + θi ) ∂θ

i=1

2

C Δt ν0 n C 2 Δt n−1 2 ¯ n 2 ≤ + θ − θ 1 + (Δt + hr )2 + 2Δt ∂θ ν0 4 8 n n   16 2 ν0 n 16 2 2 2 i 2 n−1 2 2 + σ ΔtK1 tn θ 1 + θ − θ 1 + C CL ΔtK1 tn ρi 2 ν0 4 ν 0 i=1 i=1

 ν0 16 ν0 + θn − θn−1 21 + C 2 CL2 ΔtK12 tn θi 2 + θn − θn−1 21 . (3.39) 4 ν0 4 i=1 n

Using (2.3), (2.4), (3.17) and (3.35) yields ν0 θn 21

≤

ν1 θn−1 21

 16 + C(Δt + h ) + σ 2 ΔtK12 tn θi 21 . (3.40) ν0 i=1 r 2

n

Then when Δt satisfies (3.37), an application of discrete Gronwall lemma 2.1 to the above inequality leads to (3.38). This completes the proof We observe that if (3.21) holds, then for any Δt < 1, the conclusion (3.38) is valid. In the next theorem, we will establish the error estimate for the solution computed by the standard fully discrete finite element method (2.9)-(2.10). Theorem 3.5 (Error estimate for classic FEM) Let u be the solution of (2.3)-(2.4) and U n be the solution of (2.9)-(2.10). Then, for sufficiently small Δt, we have, for all n ≥ 1, U n − u(tn )  hr + Δt,

U n − u(tn ) 1  hr−1 + Δt.

(3.41)

Proof.The first inequality is a direct result of Theorem 3.3 and (3.17). From Theorem 3.4 and (3.17), we can prove the second inequality in (3.41).

15

4. Two-grid algorithms for PIDEs with nonlinear memory In this section, we present three two-grid algorithms of the backward Euler finite element method for PIDEs with nonlinear memory. The basic mechanism in these algorithms is the construction of two regular triangulations of Ω: a coarse triangulation TH with mesh size H and a fine one Th with mesh size h (h  H). For practical purposes, Th is a refinement of TH . The corresponding finite element spaces are SH and Sh , which will be called coarse and fine space, respectively. To state the algorithms, we define, for w, u, v ∈ W 1,∞ ∩ H01 (Ω), ˜ B(w; u, v) = (α(w)∇u + β(w), ∇v) + (γ(w) · ∇u + g(w), v). Due to the assumptions on α(u), β(u), γ(u), and g(u), there exist two constants μ1 and μ2 such that ˜ |B(w; u, v)| ≤ μ1 u 1 v 1 + μ2 w v 1 .

(4.1)

Let us now present our first two-grid algorithm. Algorithm 4.1. Step one (nonlinear problem on coarse grid TH ): Given UHn−1 , find UHn ∈ SH such that   1  n ˜ Hi ; UHi , v) = (f n , v), UH − UHn−1 , v + A(UHn , v) + Δt ωni B(U Δt i=1 n

v ∈ SH , n ≥ 1, UH0 = u0H ,

(4.2) (4.3)

Step two (linear problem on fine grid Th ): Given UHn , find Uhn ∈ Sh such that   1  n ˜ i ; U i , v) = (f n , v), ωni B(U Uh − Uhn−1 , v + A(Uhn , v) + Δt H h Δt i=1 n

Uh0

v ∈ Sh , = u0h .

n ≥ 1,

(4.4) (4.5)

Firstly, we observe that for the solution of (4.4) and (4.5), our stability result is similar to the solution of (2.9) and (2.10). 16

Theorem 4.1 (Stability of two-grid FEM Algorithm 4.1) Let Uhn be the solution obtained by Algorithm 4.1. Then when Δt satisfies (3.1) and

1 3ν02 Δt ≤ min , , (4.6) 2 2μ21 K12 T we have

sup Uhi +

1≤i≤n

1/2 Uhi − Uhi−1 2

i=1



≤ C

n 

Uh0 2 + UH0 2 + Δt

n 

√

ν0 + 2

1/2 f i 2



n 

1/2 Δt Uhi 21

i=1

.

(4.7)

i=1

Proof. Similar to (3.4), using (4.1), we have Uhn 2 − Uhn−1 2 + Uhn − Uhn−1 2 + 2Δtν0 Uhn 21  n  μ1 μ2 2 i 2 n 2 i 2 n 2 ≤ (Δt) |ωni | U + μ1 1 Uh 1 + U + μ2 2 Uh 1 4 1 h 1 4 2 H i=1   (4.8) +Δt Uhn 2 + f n 2 . After choosing 1 =

ν0 2μ1 K1 tn

and 2 =

ν0 , 2μ2 K1 tn

(4.8) becomes

Uhn 2 + Uhn − Uhn−1 2 + Δtν0 Uhn 21 n  2 2  μ1 K1 tn i 2 μ22 K12 tn i 2 n−1 2 2 ≤ Uh + (Δt) Uh 1 + UH + Δt f n 2 + Δt Uhn 2 . 2ν 2ν 0 0 i=1 (4.9) With arguments similar to those in Theorem 3.1, we obtain n Δtν0  i 2 + − + U 4 i=1 h 1 i=1 

n  ≤ C Uh0 2 + sup UHi 2 + Δt f i 2 .

Uhn 2

n 

Uhi

Uhi−1 2

1≤i≤n

(4.10)

i=1

As UHi satisfies inequality (3.14), we can obtain (4.7). To establish the error estimate for the solution computed by Algorithm 4.1, we need the following lemmas. 17

Lemma 4.2 Let U n and Uhn be the solutions obtained by (2.9)-(2.10) and Algorithm 4.1, respectively. If Δt satisfies condition Δt <

ν02 , 8μ21 K12 T

(4.11)

then for any n ≥ 1, we have √

2 Whn − Whn−1 + Whn 1  H r + hr−1 + Δt, ν0 Δt

(4.12)

where Whn = Uhn − U n . Proof. It follows from (2.9) and (4.4) that  1 ˜ i ; U i , v) − B(U ˜ i ; U i , v)) = 0. ωni (B(U (Whn − Whn−1 , v) + A(Whn , v) + Δt H h Δt i=0 n

(4.13)

˜ i ; U i , v)|. Firstly, we split B(U ˜ i ; U i , v)− ˜ i ; U i , v)−B(U Now let us bound |B(U H H h h ˜ i ; U i , v) as follows: B(U ˜ i ; U i , v) ˜ Hi ; Uhi , v) − B(U B(U = (α(UHi )∇(Uhi − U i ), ∇v) + ((α(UHi ) − α(U i ))∇U i , ∇v) + (β(UHi ) − β(U i ), ∇v) +(γ(UHi ) · ∇(Uhi − U i ), v) + ((γ(UHi ) − γ(U i )) · ∇U i , v) + (g(UHi ) − g(U i ), v). (4.14) It follows that |(α(UHi )∇(Uhi − U i ), ∇v)| + |(β(UHi ) − β(U i ), ∇v)| +|(γ(UHi ) · ∇(Uhi − U i ), v)| + |(g(UHi ) − g(U i ), v)| ≤ μ1 Whi 1 v 1 + CL UHi − U i v 1 ≤ μ1 Whi 1 v 1 + CL ( u(ti ) − UHi + u(ti ) − U i ) v 1 ≤ μ1 Whi 1 v 1 + CL (H r + hr + Δt) v 1 .

(4.15)

On the other hand, due to the assumption on α and γ, which implies that α and γ are bounded and satisfy Lipschitz condition, we have   (α(UHi ) − α(U i ))∇U i , ∇v)     ≤ (α(UHi ) − α(U i ))∇(U i − u(ti )), ∇v) + (α(UHi ) − α(U i ))∇u(ti ), ∇v) ≤ C ∇(U i − u(ti )) ∇v + CL UHi − U i ∇u(ti ) ∞ ∇v (4.16) ≤ C(u)(H r + hr−1 + Δt) ∇v , 18

and   (γ(UHi ) − γ(U i )) · ∇U i , v)     ≤ (γ(UHi ) − γ(U i )) · ∇(U i − u(ti )), v) + (γ(UHi ) − γ(U i )) · ∇u(ti ), v) ≤ C ∇(U i − u(ti )) v + CL UHi − U i ∇u(ti ) ∞ v (4.17) ≤ C(u)(H r + hr−1 + Δt) v . Take v = 2(Whn − Whn−1 ) in (4.13), and combine (4.13), (4.14), (4.15), (4.16) and (4.17) to get 2 Whn − Whn−1 2 + ν0 Whn 21 + ν0 Whn − Whn−1 21 Δt n    ≤ 2Δt |ωni | μ1 Whi 1 + C(H r + hr−1 + Δt) Whn − Whn−1 1 + ν1 Whn−1 21 . i=1

(4.18) The first term on the right-hand side of the above inequality can be bounded as 2Δt

n 

  |ωni | μ1 Whi 1 + C(H r + hr−1 + Δt) Whn − Whn−1 1

i=1

 4 4 ≤ ΔtK12 μ21 tn Whi 21 + C 2 K12 t2n (H r + hr−1 + Δt)2 ν0 ν0 i=1 ν0 + Whn − Whn−1 21 , 2 n

(4.19)

where we have used 2Δt

n  i=1

 4tn  2 ν0 4tn  2 ν0 2 Δt ai + Δt b2 = Δt ai + b . ν0 4t ν 4 n 0 i=1 i=1 i=1 n

ai b ≤

n

n

Substituting (4.19) into (4.18), we get 2 ν0 Whn − Whn−1 2 + Whn 21 Δt 2 n 2 2 4μ K tn  ≤ ν1 |Whn−1 21 + 1 1 Δt Whi 1 + C(H r + hr−1 + Δt)2(4.20) . ν0 i=1 19

In view of (4.11), application of discrete Gronwall lemma 2.1 to the above inequality yields 4 Whn − Whn−1 2 + Whn 21 ≤ C(H r + hr−1 + Δt)2 . ν0 Δt

(4.21)

Then we arrive at (4.12). Combining Theorem 3.5 and Lemma 4.2 immediately yields the following theorem. Theorem 4.3 (Error estimate for two-grid FEM Algorithm 4.1) Let u be the solution of (2.3)-(2.4) and Uhn be the solution of Algorithm 4.1. Then, for sufficiently small Δt, we have, for all n ≥ 1, Uhn − u(tn ) 1  H r + hr−1 + Δt.

(4.22)

Proof. Using the triangular inequality Uhn − u(tn ) 1 ≤ U n − u(tn ) 1 + the second inequality in (3.41), and (4.12), we can obtain (4.22). From (4.22), it is easy to find that when the mesh sizes satisfy H = r−1 O(h r ) the two-grid Algorithm 4.1 achieves the same approximation for PIDEs with nonlinear memory as the classic finite element method does. Next we will present an algorithm that reduces a nonlinear problem to a symmetric positive definite (SPD) linear problem and a nonlinear system of smaller size. Algorithm 4.2. Step one (nonlinear problem on coarse grid TH ): Given UHn−1 , find UHn ∈ SH such that Uhn −U n 1 ,

  1  n ˜ i ; U i , v) = (f n , v), UH − UHn−1 , v + A(UHn , v) + Δt ωni B(U H H Δt i=1 n

v ∈ SH , n = 1, 2, · · · ,

UH0

=

(4.23) (4.24)

u0H .

Step two (SPD linear problem on fine grid Th ): Given UHn , find Uhn ∈ Sh such that   1  n ˜ hi ; Uhi , v) + Δtωnn B(U ˜ Hn ; UHn , v) ωni B(U Uh − Uhn−1 , v + A(Uhn , v) + Δt Δt i=1 n−1

= (f n , v), v ∈ Sh , n = 1, 2, · · · . Uh0 = u0h ,

(4.25) (4.26) 20

Obviously, this algorithm can also be applied to the nonsymmetric linear problem. Theorem 4.4 (Stability of two-grid FEM Algorithm 4.2) Let Uhn be the solution obtained by Algorithm 4.2. If Δt satisfies (3.1), then we have

Uhn

+

n 

1/2 Uhi

−

Uhi−1 2

i=1



Uh0 2 + Δt UH0 21 + Δt

≤ C

n 

n 1/2 ν0  i 2 + Δt Uh 1 2 i=1 1/2 √

f i 2

.

(4.27)

i=1

for any n ≥ 1. Proof. Similar to (3.4), using (4.1), we have Uhn 2 − Uhn−1 2 + Uhn − Uhn−1 2 + 2Δtν0 Uhn 21 n−1 μ μ    0 0 2 i 2 n 2 2 n 2 n 2 ≤ (Δt) |ωni | Uh 1 + μ0 Uh 1 + (Δt) |ωnn | UH 1 + μ0 Uh 1 i=1   (4.28) +Δt Uhn 2 + f n 2 . After choosing =

ν0 , μ0 K1 tn

the above inequality becomes

Uhn 2 + Uhn − Uhn−1 2 + Δtν0 Uhn 21

n−1  2 2  μ2 K 2 t n K t μ n 0 1 ≤ Uhn−1 2 + (Δt)2 Uhi 21 + 0 1 UHn 21 + Δt f n 2 + Δt Uhn 2 . ν ν0 0 i=1 (4.29) With arguments similar to those in Theorem 3.1, we obtain n Δtν0  i 2 Uh 1 4 i=1 i=1 

n n   ≤ C Uh0 2 + (Δt)2 UHi 21 + Δt f i 2 ,

Uhn 2 +

n 

Uhi − Uhi−1 2 +

i=1

21

i=1

in view of (3.14), therefore, n Δtν0  i 2 + − + U 4 i=1 h 1 i=1 

n  ≤ C Uh0 2 + Δt UH0 21 + Δt f i 2 ,

Uhn 2

n 

Uhi

Uhi−1 2

(4.30)

i=1

which implies (4.27). This completes the proof. Theorem 4.5 (Error estimate for two-grid FEM Algorithm 4.2) Let Uhn be the solution obtained by Algorithm 4.2. Then for sufficient small Δt, we have √ Uhn − u(tn )  ΔtH r−1 + hr + Δt, (4.31) for any n ≥ 1. Proof. As in Theorem 4.4, Whn = Uhn − U n satisfies the following error equation:  1 ˜ i ; U i , v) − B(U ˜ i ; U i , v)) (Whn − Whn−1 , v) + A(Whn , v) + Δt ωni (B(U h h Δt i=0 n−1

˜ n ; U n , v) − B(U ˜ n ; U n , v)) = 0. +Δtωnn (B(U H H

(4.32)

In view of the assumption on the coefficients of B, there exists a constants μB such that ˜ u, v) − B(w; ˜ |B(u; w, v)| ≤ μB u − w 1 v 1 . Then we have ˜ n ; U n , v) − B(U ˜ n ; U n , v))| |Δtωnn (B(U H H ≤ μB Δt|ωnn | UHn − U n 1 v 1 ≤ μB K1 Δt( UHn − u(tn ) 1 + u(tn ) − U n 1 ) v 1 ≤ μB K1 Δt(H r−1 + hr−1 + Δt) v 1 .

22

(4.33)

Take v = 2ΔtWhn in (4.32) to obtain Whn 2 − Whn−1 2 + Whn − Whn−1 2 + 2Δtν0 Whn 21 n−1  1 2 ≤ μB (Δt) |ωni |( Whi 21 + Whn 21 ) i=1  1 n 2 r−1 2 + Δt) . +μB K1 Δt Δt Wh 1 + Δt(H By choosing =

ν0 , tn μB K1

(4.34)

we get

Whn 2 − Whn−1 2 + Whn − Whn−1 2 + Δtν0 Whn 21 n−1  tn μ2B K12 tn μ2B K12 ≤ (Δt)2 Whi 21 + (Δt)2 (H r−1 + Δt)2 . (4.35) ν0 ν 0 i=1 Sum from 1 up to n to obtain Whn 2 ≤ (Δt)2

−

n  i=1

≤ (Δt)

2

n−1  i=0

≤ (Δt)

2

Wh0 2

+

n 

Whi

i=1 i−1 ti μ2B K12 

ν0

j=1

ν0 B

ν0

n  i=1

i ti+1 μ2B K12 

Whj 21

+

j=1 1

+ Δtν0

n 

Whi 21

i=1

Whj 21 +

i n−1  ti+1 μ2 K 2  i=0

−

Whi−1 2

ti μ2B K12 (Δt)2 (H r−1 + Δt)2 ν0

n  t i μ2 K 2 B

ν0

i=1

Whj 21 +

j=1

1

(Δt)2 (H r−1 + Δt)2

t2n μ2B K12 Δt(H r−1 + Δt)2 . ν0

(4.36)

An application of discrete Gronwall Lemma 2.1 yields Whn 2 +

n 

Whi − Whi−1 2 + Δtν0

i=1

t 2 μ2 K 2 ≤ n B 1 Δt(H r−1 + Δt)2 exp ν0



n 

Whn 21

i=1

t2n μ2B K12 ν02

.

(4.37)

Finally, (4.31) follows readily from this result when a triangular inequality is also applied. 23

Next we will present an algorithm that significantly reduces computational memory and storage requirements when B gathers lower-order spatial derivatives and nonlinear terms. To state the algorithm, we define ˜s (w; u, v) = (α(w)∇u, ∇v), B and N (w; u, v) = (β(w), ∇v) + (γ(w) · ∇u + g(w), v). In view of the assumptions on α(u), β(u), γ(u), and g(u), we find that there exist two constants μ3 and μ4 such that ˜s (w; u, v)| ≤ μ3 u 1 v 1 |B |N (w; u, v)| ≤ μ4 u v 1 .

(4.38) (4.39)

Then the algorithm can be stated as follows. Algorithm 4.3. Step one (nonlinear problem on coarse grid TH ): Given UHn−1 , find UHn ∈ SH such that   1  n ˜ Hi ; UHi , v) = (f n , v), ωni B(U UH − UHn−1 , v + A(UHn , v) + Δt Δt i=1 n

v ∈ SH , n ≥ 1,

UH0

=

(4.40) (4.41)

u0H .

Step two (linear problem on fine grid Th ): Given UHn , find Uhn ∈ Sh such that   1  n ˜s (U i ; U i , v) + N (U i ; U i , v)) ωni (B Uh − Uhn−1 , v + A(Uhn , v) + Δt H h H H Δt i=1 n

= (f n , v), Uh0 = u0h ,

v ∈ Sh ,

n ≥ 1,

(4.42) (4.43)

The stability of Algorithm 4.3 can be obtained by the same argument for Theorem 4.1. Theorem 4.6 (Stability of two-grid FEM Algorithm 4.3) Let Uhn be the solution obtained by Algorithm 4.3. Then when

3ν02 1 , , (4.44) Δt ≤ min 2 2μ23 K12 T 24

we have

sup Uhi +

1≤i≤n

1/2 Uhi − Uhi−1 2

i=1



≤ C

n 

Uh0 2

+

UH0 2

+ Δt

n 

√

ν0 + 2

1/2 f i 2



n 

1/2 Δt Uhi 21

i=1

.

(4.45)

i=1

Proof. Similar to (3.4), using (4.38) and (4.39), we have Uhn 2 − Uhn−1 2 + Uhn − Uhn−1 2 + 2Δtν0 Uhn 21  n  μ3 μ4 2 i 2 n 2 i 2 n 2 ≤ (Δt) |ωni | U + μ3 1 Uh 1 + U + μ4 2 Uh 1 4 1 h 1 4 2 H i=1   (4.46) +Δt Uhn 2 + f n 2 . After choosing 1 =

ν0 2μ3 K1 tn

and 2 =

ν0 , 2μ4 K1 tn

(4.46) becomes

Uhn 2 + Uhn − Uhn−1 2 + Δtν Uhn 21 n  2 2  μ3 K1 tn i 2 μ24 K12 tn i 2 n−1 2 2 ≤ Uh + (Δt) Uh 1 + UH + Δt f n 2 + Δt Uhn 2 . 2ν 2ν 0 0 i=1 (4.47) With arguments similar to those in Theorem 3.1, we obtain n Δtν0  i 2 + − + Uh 1 4 i=1 i=1 

n  ≤ C Uh0 2 + sup UHi 2 + Δt f i 2 .

Uhn 2

n 

Uhi

Uhi−1 2

1≤i≤n

As

(4.48)

i=1

UHi

satisfies inequality (3.2), we can obtain (4.45). To get an idea of the accuracy of Algorithm 4.3, we have the following theorem. Theorem 4.7 (Error estimate for two-grid FEM Algorithm 4.3) Let Uhn be the solutions obtained by Algorithm 4.3. Then for sufficient small Δt, we have √ Uhn − u(tn )  hr + Δt + ΔtH r , Uhn − u(tn ) 1  H r + hr−1 + Δt, (4.49) for any n ≥ 1. 25

Proof. Set Whn = Uhn − U n to get  1 ˜s (U i ; U i , v) − B ˜s (U i ; U i , v)) (Whn − Whn−1 , v) + A(Whn , v) + Δt ωni (B H h Δt i=1 n

+Δt

n 

ωni (N (UHi ; UHi , v) − N (U i ; U i , v)) = 0.

(4.50)

i=1

Similar to the proof of Lemma 4.2, we have N (UHi ; UHi , v) − N (U i ; U i , v) = (β(UHi ) − β(U i ), ∇v) + (γ(UHi ) · ∇(UHi − U i ), v) +((γ(UHi ) − γ(U i )) · ∇U i , v) + (g(UHi ) − g(U i ), ∇v)

(4.51)

and   (γ(UHi ) · ∇(UHi − U i ), v) ≤ C UHi − U i v 1 ≤ CL ( u(ti ) − UHi + u(ti ) − U i ) v 1 ≤ CL (H r + hr + Δt) v 1 . (4.52) The desired estimate can then be obtained in a way similar to proofs of Theorem 4.3 and Lemma 4.2. Remark. Observe that when α ≡ 0, the approximation of the integral term on the fine grid is identical to the approximation of the integral term on the coarse grid. This means that when we solve Uhn , all Uhi (i < n − 1) do not need to be stored on a fine grid. It also means that once the approximation of the integral term has been computed on the coarse grid it does not need to be computed on the fine grid. This significantly reduces computational memory and storage requirements. This result is novel and interesting even for linear problem. 5. Numerical experiments In this section, we show some experiments to confirm the effectiveness and theoretical analysis for Algorithm 4.3. We set the domain as [0, 1] × [0, 1] and T = 1.00. Noting that when α ≡ 0 in Algorithm 4.3, the algorithm does not need to store Uhi (i < n − 1), hence in order to confirm the efficiency and advantage of Algorithm 4.3, we set K(t) = e−t , α(u) = 0, β(u) = 26

(sin u, 1 − cos u)T , γ(u) = (1 − cos u, sin u)T , g(u) = sin u in (1.1) and we solve the following problem  t   ut − Δu + e−(t−s) − ∇ · β(u(s)) + γ(u(s)) · ∇u + g(u(s)) ds 0

= f (x1 , x2 ; t),

(5.53)

u(x1 , x2 ; t) = 0, (x1 , x2 ; t) ∈ ∂Ω × (0, T ], u(x1 , x2 ; 0) = u0 (x1 , x2 ), (x1 , x2 ) ∈ Ω. We further set u0 (x1 , x2 ) = x1 (1 − x1 )x2 (1 − x2 ) and f (x1 , x2 ; t)   = 2x1 (1 − x1 ) − x1 (1 − x1 )x2 (1 − x2 ) + 2x2 (1 − x2 ) + (1 − 2x1 )x2 (1 − x2 )t e−t  t   −t − 2(1 − 2x1 )x2 (1 − x2 )e cos x1 (1 − x1 )x2 (1 − x2 )e−s ds 0  t   + e−t es sin x1 (1 − x1 )x2 (1 − x2 )e−s ds. 0

Then we can verify that u(x1 , x2 ; t) = x1 (1 − x1 )x2 (1 − x2 )e−t is the true solution. We use linear finite element for the space discretization. The convergence rate and effectiveness of Algorithm 4.3 in H 1 norm given √ by Theorem 4.7 are confirmed in Table 1 with h = 21l , l = 2, · · · , 9; H = 12 h and Δt = 21l , l = 1, · · · , 8. H 1/4 1/6 1/8 1/12 1/16 1/23 1/32 1/46

h 1/4 1/8 1/16 1/32 1/64 1/128 1/256 1/512

Δt 1/2 1/4 1/8 1/16 1/32 1/64 1/128 1/256

UhT − u(T ) 1 2.17236 × 10−2 1.11164 × 10−2 5.59226 × 10−3 2.80089 × 10−3 1.40136 × 10−3 7.00760 × 10−4 3.50427 × 10−4 1.75207 × 10−4

h order −− 0.95 0.99 0.99 1.00 0.99 1.00 1.00

H order −− 1.95 1.99 1.99 2.00 1.99 2.00 2.00

Δt order −− 0.95 0.99 0.99 1.00 0.99 1.00 1.00

Table 1: Convergence rate and accuracy of Algorithm 4.3.

Following the Algorithm 4.3, in the numerical experiments, we do not store Uhi (i < n − 1) and save a lot of storege. Further, the method is much 27

more efficient than the standard fully discrete finite element algorithm √ since 1 we only need to solve a nonlinear problem with mesh-size H = 2 h and then solve the linear problem with mesh-size h. Using standard fully discrete finite element algorithm to solve the problem (5.53) by solving the nonlinear problem directly with mesh-size h and the convergence rate and error in H 1 norm are shown in Table 2. Comparing Table 1 and Table 2, we can clearly see that the effectiveness and accuracy of Algorithm 4.3 are the same as standard fully discrete finite element algorithm. The error estimate in L2 norm given by Theorem 4.7 can also be confirmed similarly, for simplicity, we omitted listing the tables here. h 1/4 1/8 1/16 1/32 1/64 1/128 1/256 1/512

Δt 1/2 1/4 1/8 1/16 1/32 1/64 1/128 1/256

UhT − u(T ) 1 2.17183 × 10−2 1.11115 × 10−2 5.58847 × 10−3 2.79844 × 10−3 1.39977 × 10−3 6.99958 × 10−4 3.49990 × 10−4 1.74996 × 10−4

h order −− 0.95 0.99 1.00 1.00 1.00 1.00 1.00

Δt order −− 0.95 0.99 1.00 1.00 1.00 1.00 1.00

Table 2: Error and convergence rate for standard fully discrete finite element algorithm.

6. Concluding remarks We have presented and derived error estimates for several two-grid finite element algorithms for PIDEs with nonlinear memory. With the backward Euler scheme, the two-grid strategy consists of two steps: (1) discretizing the fully nonlinear problem in space on a coarse grid with mesh-size H and time step-size Δt and (2) discretizing the linearized problem in space on a fine grid with mesh-size h and the same time step-size as in step (1). It is shown that these algorithms are as stable as the standard fully discrete finite element algorithm. We also present the error estimate at each time step. Compared with standard finite element methods, our algorithm not only keep good accuracy but also saves a lot of computational cost. As a byproduct of these results, we found that one of these algorithms, Algorithm 4.3, significantly reduces computational memory and storage requirements if the 28

nonlinear memory is defined by a first-order or zero-order nonlinear differential operator. Thus, the two-grid methods studied in this paper provide a new approach that takes advantage of some of the nice properties hidden in a complex problem. Numerical experiments for Algorithm 4.3 are provided to confirm the theoretical results and show that the two-grid method has the same effectiveness and accuracy as the standard fully discrete finite element algorithm. The analysis herein was carried out for an implicit Euler discretization in time. However, the results could be extended to the second-order accuracy backward differentiation formula (BDF) scheme. Moreover, the analysis is valid for a state-dependent forcing term f that satisfies certain conditions, e.g., ∂ ∂2 | f (x, t, u)| + | 2 f (x, t, u)| ≤ M, u ∈ R, ∂u ∂u where M is a positive constant. Acknowledgments The first author thanks Professor Jinchao Xu for suggesting this problem and for many stimulating and inspiring discussions. This paper was written at the School of Mathematical Sciences, Peking University, where the first author spent time as a visiting scholar. This work was partially supported by the National Natural Science Foundation of China [grant numers 11771060,11371074]. [1] H. Abboud and T. Sayah, A full discretization of the time-dependent Navier-Stokes equations by a two-grid scheme, M2AN Math. Model. Numer. Anal., 42 (2008), 141-174. [2] H. Abboud, V. Girault and T. Sayah, A second order accuracy for a fully discretized time-dependent Navier-Stokes equations by a two-grid scheme, Numer. Math., 114 (2009), 189-231. [3] O. Axelsson and W. Layton, A two-level method for the discretization of nonlinear boundary value problems, SIAM J. Numer. Anal. 33 (1996), 2359-2374. [4] C. Bi and V. Ginting, Two-grid finite volume element method for linear and nonlinear elliptic problems, Numer. Math., 108 (2007) 177198. 29

[5] C. Bi and V. Ginting, Two-grid discontinuous Galerkin method for quasi-linear elliptic problems, J. Sci. Comput., DOI: 10.1007/s10915011-9463-9. [6] I. H. Biswas, E. R. Jakobsen, and K. H. Karlsen, Differencequadrature schemes for nonlinear degenerate parabolic integro-PDE, SIAM J. Numer. Anal. 48 (2010), 1110-1135. [7] M. Cai, M. Mu and J. Xu, Numerical solution to a mixed NavierStokes/Darcy model by the two-grid approach, SIAM J. Numer. Anal., 47 (2009), 3325-3338. [8] J. R. Cannon and Y. Lin, Non-classical H 1 projection and Galerkin methods for nonlinear parabolic integro-differential equation, Calcolo 25 (1988) 187-201. [9] J. R. Cannon and Y. Lin, A priori L2 error estimates for finite element methods for nonlinear diffusion equations with memory, SIAM J. Numer. Anal. 27 (1990) 595-607. [10] C. Chen and W. Liu, Two-grid finite volume element methods for semilinear parabolic problems, Appl. Numer. Math., 60 (2010), 10-18. [11] C. Chen and W. Liu, A two-grid method for finite volume element approximations of second-order nonlinear hyperbolic equations, J. Comput. Appl. Math., 233 (2010), 2975-2984. [12] C. Chen and T. Shih, Finite element methods for integro-differential equations, Singapore: World Scientifi Pub. Co., 1998. [13] Y. Chen and L. Li, Lp error estimates of two-grid schemes of expanded mixed finite element methods, Appl. Math. Comput., 209 (2009), 197205. [14] Y. Chen, Y. Huang and D. Yu, A two-grid method for expanded mixed finite-element solution of semilinear reaction-diffusion equations, Int. J. Numer. Methods Eng. 57 (2003), 193-209. [15] X. Dai and X. Cheng, A two-grid method based on Newton iteration for the Navier-Stokes equtions, J. Comput. Appl. Math., 220 (2008), 566-573. 30

[16] C. N. Dawson and M. F. Wheeler, Two-grid method for mixed finite difference approximations fo non-linear parabolic equations, Contemp. Math., 180 (1994), 191-203. [17] C. N. Dawson, M. F. Wheeler and C. S. Woodward, A twogrid finite difference scheme for nonlinear parabolic equations, SIAM J. Numer. Anal., 35 (1998), 435-452 [18] E. Emmrich, Stability and error of the variable two-step BDF for semilinear parabolic problems, J. Appl. Math. Computing, 19 (2005), 33-55. [19] R. E. Ewing, Y. Lin, T. Sun, J. Wang, and S. Zhang, Sharp L2 error estimates and superconvergence of mixed finite element methods for non-Fickian flows in porous media, SIAM J. Numer. Anal., 40 (2002), 1538-1560. [20] M. Gurtin and A. Pipkin, A general theory of heat conduction with finite wave speeds, Arch. Rational Mech. Anal., 31 (1968), 113-126. [21] V. Girault and J. L. Lions, Two-grid finite-element schemes for the steady Navier-Stokes problem in polyhedra. Portugal Math., 58 (2001), 25-57. [22] V. Girault and J. L. Lions, Two-grid finite-element schemes for the transient Navier-Stokes equations, M2AN 35 (2001), 945-980. [23] Y. He and K. Li, Two-level stabilized finite element methods for the Steady Navier-Stokes problem, Computing, 74 (2005), 337-351. [24] X. Hu and X. Cheng, Accleration of a two-grid method for eigenvalue problems, Math. Comput. 80 (2011), 1287-1301. [25] Y. Q. Huang, Time discretization scheme for an integro-differential equation of parabolic type, J. Comput. Math., 3 (1994), 259-264. [26] J. Jin, S. Shu and J. Xu, A two-grid discretization method for decoupling systems of partial differential equations, Math. Comp. 75 (2006), 1617-1626. [27] W. Layton and L. Tobiska, A two-level method with backtracking for the Navier-Stokes equations, SIAM J. Numer. Anal., 35 (1998), 20352054. 31

[28] S. Li and Z. Huang, Two-grid algorithms for some linear and nonlinear elliptic systems, Computing, 89 (2010), 69-86. [29] Y. Lin, Galerkin methods for nonlinear parabolic integrodifferential equations with nonlinear boundary conditions, SIAM J. Numer. Anal., 27 (1990), 608-621. [30] Y. Lin, V. Thom´ ee and L. Wahlbin, Ritz-Volterra projections onto finite element spaces and applications to integro-differential and related equations, SIAM J. Numer. Anal. 28 (1991), 1047-1070. [31] J. C. Lopez-Marcos, A difference scheme for a nonlinear partial integrodifferential equation, SIAM J. Numer. Anal. 27 (1990), 20-31. [32] M. Marion and J. Xu, Error estimates on a new nonlinear Galerkin method based on two-grid finite elements, SIAM J. Numer. Anal. 32 (1995), 1170-1184. [33] R. K. Miller, An integro-differential equation for rigid heat conductions with memory, J. Math. Anal. Appl., 66 (1978), 313-332. [34] M. Mu and J. Xu, A two-grid method of a mixed Stokes-Darcy model for coupling fluid flow with porous media flow, SIAM J. Numer. Anal., 45 (2007), 1801-1813. [35] K. Mustapha, H. Brunner, H. Mustapha, and D. Schotzau, An hp-version discontinuous Galerkin method for integro-differential equations of parabolic type, SIAM J. Numer. Anal., 49 (2011), 1369-1396. [36] A. K. Pani, G. Fairweather, and R. I. Fernandes, Alternating direction implicit orthogonal spline collocation methods for an evolution equation with a positive-type memory term, SIAM J. Numer. Anal., 46 (2008), 344-364. [37] A. K. Pani and T. E. Peterson, Finite element methods with numerical quadrature for parabolic integrodifferential equations, SIAM J. Numer. Anal., 33 (1996), 1084-1105. [38] X. Qin and Y. Ma, Two-grid scheme for characteristics finite-element solution of nonlinear convection diffusion problems, Appl. Math. Comput., 165 (2005) 419-431. 32

[39] M. Raynal, On some nonlinear problems of diffusion, in Volterra Equations, S. London and 0. Staffans, eds., Lecture Notes in Math., 737, Springer-Verlag, Berlin, New York, 1979, pp. 251-266. [40] Y. Shang and K. Wang, Local and parallel finite element algorithms based on two-grid discretization for the transient Stokes equations, Numer. Algor. 54 (2010), 195-218. [41] R. K. Sinha, R. E. Ewing, and R. D. Lazarov, Mixed finite element approximations of parabolic integro-differential equations with nonsmooth initial data, SIAM J. Numer. Anal., 47 (2009), 3269-3292. [42] R. K. Sinha, R. E. Ewing, and R. D. Lazarov, Some new error estimates of a semidiscrete finite volume element method for a parabolic integro-differential equation with nonsmooth initial data, SIAM J. Numer. Anal., 43 (2006), 2320-2343. [43] I. H. Sloan and V. Thomee, Time discretization of an integrodifferential equation of parabolic type, SIAM J. Numer. Math., 23 (1986), 1052-1061. [44] T. Tachim Medjo and R. Temam, A two-grid finite difference method for the primitive equations of the ocean, Nonlinear Anal., 69 (2008), 1034-1056. [45] W. S. Wang, Long-time behaviour of two-grid finite element method for the fully discrete semilinear evolution equations with positive memory, J. Comput. Appl. Math., 250 (2013), 161-174. [46] L. Wu and M. B. Allen, A two-grid method for mixed finite element solution of reaction-diffusion equations, Numer. Methods Partial Differ. Equ. 15 (199), 317-332. [47] J. Xu, A new class of iterative methods for nonselfadjoint or indefinite elliptic problems, SIAM J. Numer. Anal. 29 (1992), 303-319. [48] J. Xu, Some Two-Grid Finite Element Methods, Tech. Report, P.S.U, 1992. [49] J. Xu, A novel two-grid method for semi-linear equations, SIAM J. Sci. Comput. 15 (1994), 231-237. 33

[50] J. Xu, Two-grid finite element discretization techniques for linear and nonlinear PDE, SIAM J. Numer. Anal. 33 (1996), 1759-1777. [51] J. Xu and A. Zhou, Local and parallel finite element algorithms based on two-grid discretization for nonlinear problems, Adv. Comput. Math., 14 (2001), 293-327. [52] J. Xu and A. Zhou, A two-grid discretization scheme for eigenvalue problems, Math. Comput., 70 (2001), 17-25. [53] J. Xu and A. Zhou, Local and parallel finite element algorithms for eigenvalue problems, Acta Math. Appl. Sin. Engl. Ser. 18 (2002), 185200. [54] T. Zhang, Finite element methods for partial differenio-integral equations, Beijing: Science Press, 2009 (in Chinese).

34