A finite difference scheme for partial integro-differential equations with a weakly singular kernel

Applied Numerical North-Holland

APNUM

Mathematics

11 (1993) 309-319

309

375

A finite difference scheme for partial integro-differential equations with a weakly singular kernel * Tao

Tang

Department

of Mathematics

and Statistics,

Simon Fraser Uniuersity, Burnaby, BC, Canada WA

IS6

Abstract Tang, T., A finite difference scheme for partial integro-differential Applied Numerical Mathematics 11 (1993) 309-319.

equations

with a weakly singular

kernel,

A finite difference method for the numerical solution of partial integro-differential equations is considered. the time direction, a Crank-Nicolson time-stepping is used to approximate the differential term and product trapezoidal method is employed to treat the integral term. An error bound is derived for numerical scheme. Due to lack of smoothness of the exact solution, the overall numerical procedure does achieve second-order convergence in time. But the convergence order in time is shown to be greater than which is confirmed by a numerical example. Keywords.

Partial

integro-differential

equations;

Crank-Nicolson

method;

product

trapezoidal

In the the not one,

method.

1. Introduction This paper equation

is concerned

u, =Bu + subject

with the Crank-Nicolson

j,i(l -s)-1’2u,,(x,

to the boundary U(0, t)=u(l,

method

for the partial

integro-differential

s) ds,

(1.1)

condition t)=O,

t>o,

(1.2)

and the initial condition U(X, 0) = U&),

0
(l-3)

Correspondence to: T. Tang, Department of Mathematics and Statistics, Simon Canada WA lS6. E-mail: [email protected]. * This work was partially supported under NSERC Canada Grant OGP0105545. 016%9274/93/$06.00

0 1993 - Elsevier

Science

Publishers

B.V. All rights reserved

Fraser

University,

Burnaby,

BC,

310

T. Tang / Partial integro-differential equations

where B in (1.1) is a linear or nonlinear differential operator. In the practical applications, we often have Bu = -uux,

(1.4)

or, with p 2 0 a constant, Bu = pu,,.

(1.5)

The problem (l.l)-(1.4) (hereafter referred to as Problem I> can be found in the modeling of physical phenomena involving viscoelastic forces, e.g., [5,10]. It has been pointed out by Sanz-Serna [14] that (1.1) and (1.4) give a simple model equation that combines the Eulerian derivative u, + uu, with a viscoelastic effect, just as Burgers’ equation provides a simple model for the study of more realistic situations involving Eulerian derivatives and viscous forces. Problem I has been recently considered numerically by several researchers [1,3,6,14]. The method implemented by Christie [3] treats the weakly singular integral kernel in (1.1) by means of the product trapezoidal technique (see, e.g., [ll, Chapter 81). The Crank-Nicolson method has been employed by this author for time-stepping but it has been observed that the overall procedure does not achieve second-order convergence in time. The problem (l.l)-(1.3) and (1.5) (hereafter referred to as Problem II) can be found in the modeling of heat flow in materials with memory, e.g. [4,7], and in linear viscoelastic mechanics, e.g., [2,12]. In linear viscoelastic problems, the integral term in (1.1) represents the viscosity part of the equation, and p in (1.5) is a Newtonian contribution to the viscosity. It can be seen that in (1.1) the kernel function has a weak singularity at the origin. This is particularly interesting in viscoelasticity, because it might smooth the solution when the boundary data is discontinuous [12]. Numerical investigations for Problem II have been given by several authors (see, e.g., [8,9,15]), but most of them considered smooth integral kernels only. In this paper, we shall construct a finite-difference scheme for (1.11, based on the CrankNicolson method and the product trapezoidal technique, and present convergence analysis for the scheme. Throughout this paper, we assume that u0 in (1.3) is such that the problem (l.G(l.3) has a unique solution in [O, 11 x [O, T]. Furthermore, we suppose that uXxXXand u, are continuous in [0, l] x [0, T]. We also assume that ult, uttt, and uXXtt are continuous for 0 G x G 1 and 0 < t G T, and that there exists a positive constant C, such that for 0 < x < 1 and O
I u,,(x, t) I G Cfp2,

Iu xxtt(X,t) I =GGp2

I U&,

t> I G C”t-3’2,

(1.6)

(see [6] for these assumptions). The remainder of the paper is organized as follows. In Section 2 we introduce a numerical scheme for solving Problems I and II, which is based on the Crank-Nicolson method in time and central differences in space. The product trapezoidal technique is used to approximate the integral term. The convergence properties of the numerical scheme are investigated in Section 3. Numerical results are presented in the final section. 2. Numerical scheme We introduce a grid xi =jh, j = 0, l,..., J, with h = l/J and J a positive integer. The steplength in time is denoted by k, k = T/N with N a positive integer, and a subscript IZrefers

T. Tang / Partial integro-differential equations

311

to the time level t, = nk, n = 0, 1,. . . , N. Moreover, we let t,,,,, = (n + i)k, 0 G IZG N - I. We first consider the approximation of an integral term with a similar type as that in (1.1). For any f~ C’([O, Tl) n C3((0, Tl) satisfying f”(t) = Ott-‘I*) and f”‘(t) = 0(tc3/*) as t + o+, we will approximate I(f, t) := /i(t - ~)-‘/~f(s> d s numerically. It is easy to see that I(.L t,+1,z) = 3[W>

tn> +1(f7 &+*)I + O(k2t,3’2),

12> 1

(2.1)

and I(.L t,,,) Since ti3/*

<

= 31(f, t,> + O(k”*).

2 3’2t,:(2

(2.2)

for II > 1, (2.1) and (2.2) yield that

I(f? L,,,)=

3[W7

t,> +q"f, L+,)] + O(k2t;Z2),

n 20.

P-3)

now use the product trapezoidal technique to approximate the integrals I(f, t,), 1 < ~1< N. For f(t, - 0) with 0 E [tj, tj+lI, 0
We

e,= ‘j”,-

f(tn -

ef(tn-j)

+

(2.4)

+R,l,

yf(tn_j_l)

and with 8 E [tn_l, t,l, we obtain

I3-in-‘p(t,)

t, - 8

fk - 0) = --f(tJ

+

o(k3’*)7

+

(2.5)

where in (2.4) the remainder term R,, can be bounded by O(k*t;?/! 1) = O( k3/*), 0 1. Further, using (2.4) and (2.3, we have z(f, t,) = /‘“eP/*f(

t, - e) de =

0

j=O

n-l

j=O

t, - e> de

f,

“+~-Bf(t,-j) +

cl+1 -I,*

C1

=

ncl jr’+‘s-1/2f(

0

since

~i(tn_j_~)]

de + Rn2,

P-6)

tj

with n-l I

I < c /tJ+1e-1/20(k3/2) de = O(W).

R,,

j=O It

follows by using integration n-1

C j=O

by parts that

c/+1 8 _I,2 fi+~Aef(tn_j)

/

2

de

+ E+(tn_j_J]

[

l, n-l

=

P-7)

t,

c

n-1

[t,!!ff(fn_j_l)

- t,!/2f(tn_j)] +

j=O

=&f(to)

2

c j=O

+

2 cpf(tn-_p)7

p=o

/tf+1fj1/2/(f”i) -[(“-j-‘)

de

t,

(2.8)

312

T. Tang / Partial integro-differential equations

where A,

=

2

I

t;/2

$/tn+1fj1,2

_

d()

tn

2 C” = - “81/2 de, k /o 2 tp+‘@/2 &I c, = k

I

/

I ,

(2.9) (2.10)

tP tp-I1, p21. - 1’”

e1j2

de

(2.11)

Combining (2.61, (2.71, and (2.81, we obtain that I(f,

t,) =A,f(t,)

+ 2 cpf(tn_p) p=o

+ 0(k3’2),

1
(2.12)

P222,

(2.13)

where A,, and c, are given by (2.9)-(2.11). Let 4& p0 = C” + ~ 3 P;

4& p1 = Cl - - 3 P;

Pp=cp,

where p is a nonnegative constant and is independent of k and h, i.e., p > 0 and p = O(1). The sequence {p,),“,O reduces to (c~}~=~ if the parameter /3 = 0. Since 4&E TP[f(Q

-f(t,-111= w3’2),

we obtain from (2.12) that t,) =A&,)

W

Now an application .y+l

+ 2 P,f(t,-,) p=o

of the standard Crank-Nicolson

k-u7 4

“7”:“‘) + 0(k2t,:(2

(2.14)

+ 0(k3’2).

method for (1.1) gives

+~+1,i(tn+l,2_S)1/2S2U~’ ‘) ds

+

h’),

(2.15)

for 0 G 12G N - 1 and 1
q+;q)

= Bu(xj,

tn+1,2)

+ 0(k2t;2f2

+

II”).

Also in (2.15) we have used the notation

operator

which is a discretized

(2.16)

a2y. = q+l - 25 + r/;.-I. The remainder term in (2.15) has been obtained by using Taylor’s expansion and in particular we have noticed that

= 0( k2t,-,3i2),

T. Tang / Partial integro-differential equations

313

where we have made use of (1.6). Combining (2.3), (2.141, and (2.1% we obtain that .;+I ku:

=B( “‘“:“‘)

+i,62U;+

~$+-;;Y-lj

+O(k%,;{2+k3’2+h2),

(2.17)

n>,O,

where k

= P, +A, +A, + 4&P/3 0

2h2

P n+l

>

2,

=

+A+1

+A,

2h2

(n > 1).

(2.18)

Our numerical scheme is based on the formula (2.17) and takes the form (2.19) for 0 < n < N - 1 and 1
O
(2.20)

q” = UO(Xj),

l
(2.21)

For Problem I, we choose the difference operator B’ such that (2.22) to approximate the differential AlI$ = y+1 - Ql. Since

a

tn+d +a 2

operator

4%)

- +,

B given in (1.41, where we have used the notation

t,,,,,)

= O(kQ;:<2)

< O(k%,-,3{2),

it can be verified that B’ given by (2.22) and B given by (1.4) satisfy the requirement For Problem II, we choose g such that

(2.16).

(2.23) and it can be shown that this operator and the one given by (1.5) satisfy (2.16). For later use we collect the following notations and results, as given in [6]. Throughout this paper, we shall use the notation U”, 0
314

T. Tang / Partial integro-differential equations

W=(W,, Iv,,...,

IV_,) are real vectors, then we define J-l

A+?

= y+1 - I$

IIv IIm=

max

I”; I,

l
IIVII1=

ChlF

I,

j=l

J-l (V,

W)

=

c

II~l12=W,

hvyq,

v>.

j=l

Furthermore,

the following identities hold (see, e.g., [13]):

(I$tl + y”+ (VAV+ Q2V,

3. Convergence

y.“l)Av;.n = ~“A~“+

A(V)‘,

(2.24)

A(yn)2,

(2.25)

V) = 0, A+W)

W) = -(A+V,

= -

J-1 c hA+F(A+lJ$. j=l

(2.26)

analysis

From now on, C denotes a positive constant which is independent of k, h, j, and II (with and O
O
Lemma 1. Let {a,}tzO be a sequence of real numbers with the properties a,

a n+1-

> 0,

a,

G

0,

a n+l-2a,+a,_,>0.

Then for any positive integer M, and real vector (VI, V,, . . . , V,) M-l

n=O +

(3.1)

with M real entries,

n

p=o

a,K+,-,

I

(3.2)

K+ld-

Now we will check if (c~}~=~defined by (2.10) and (2.11) satisfies (3.1). Lemma 2. The sequence (c~}~=~ defined by (2.10) and (2.11) satisfies cP > 0 and c~+~ < c,. Proof. It is obvious that c0 > 0. It can be shown from (2.11) that c,=/_X*(t,-R)-1’2(l-~)

de,

pal,

(3.3)

T. Tang / Partial integro-differential equations

which implies that c, 2 0 for fact that &k

p > 1.

- 8)-‘/2
Further, a direct calculation shows that c1 < cO. Using the

for r-21 and 06(-k,

we can obtain that cptl Q cP for

p 2

315

1.

k)

0

Lemma 3. The sequence {c~]~=~ satisfies c2 - 2c, + co < and c~+~ - 2c, + cP_* > 0 furp

> 2.

Proof. The first part of this lemma follows from c2 - 25 + co = ?[3&8fi+6]

(34


The second part follows from (3.3) and the fact that d2 #‘-

e>- 1’2>0

for r-21 and OE(---k,

k).

•!

It can be seen from Lemma 3 that {cJ=~ defined by (2.10) and (2.11) does not satisfy (3.1). In order to give a convergence result for the numerical scheme (2.19)-(2.21) we need to choose the parameter p 2 0 such that the sequence {P,}F=a in (2.13) satisfies the property (3.1). By recalling the proofs of Lemmas 2 and 3, we require that

which are equivalent to c2 - 2c, + co ZGC, -c2,

z p.

z> -

3

’

z
(3.6)

Since

c2 - 2c, + cg < 0,

c1 > ci - c2 2 c3 - 2c, + cl 2 0,

(3.7)

the inequalities (3.6) lead to -

c2 - 2c, + cg 3

< z < cg - 2c, + Cl,

P-8)


(3.9)

which is equivalent to -3&+8fi-6 3

The above discussion yields the following result. Lemma 4. If p satisfies (3.91, then the sequence {p,}~=o defined by (2.13) satisfies (3.1).

Now we begin to derive the error bounds for the numerical scheme (2.19)-(2.21). Let ey = q.” - u;, with ur = u(xj, t,) and C$F is the solution of (2.19M2.21). Subtraction of (2.17)

316

T. Tang / Partial integro-differential equations

from (2.19) gives e;+’ k

-ey

=l?(y)

-I+;)

+

k p, p=O

where Vj” := (CJn+l + Lsn)/2, uj” := (~7” + $)/2 bounded by 1rj” 1< 0( k*t;:{*

+ k3’*

lP(ent1-p

+ q-q

2h2

(3.10)

try,

and r,!’ is the residual term which can be

+ 15’).

(3.11)

Multiplying both sides of (3.10) by hk(er+’

+ ey> and summing in j, we obtain

IIe n+l II* - IIen II* =k(i(V”)

- -$

-E(L’“),

en+’ +e”)

p,( A+(e”+lmP + enep), A+(e”+’ + en))+ k( rn, en+’ + en),

$

(3.12)

P-O

where we have used (2.25) to get the second term of the right-hand side of (3.12). For Problem I, following Lopez-Marcos [6, p. 1281, we can show that the modulus of the first term of the right-hand side of (3.12) can be bounded by Ck IIV” - un I(* < Ck( IIe”+l 1)* + IIen )I‘).

(3.13)

For Problem II, it follows from (2.26) that k(B(V”)

-B(u”),

en+’ + e”) = 2k( l?(Vn) - B(u”),

I/” - u”)

(3.14) For the last term of (3.12) we have I k(r”,

en+’ +e”>l

~kllr”llm(Ile”+l

II 1+ II enII 1) G Ck II r” II ,A,

(3.15)

where A = maxo,.G, )Ien ]I, and we have made use of the fact that )IW 1)1 G d% (I W II = II W II. If we use the estymates (3.13)-(3.15) to (3.12) and sum over it (and note that 1)e” (1= O>,then we obtain, for both Problems I and II, that

II en+11(2
2

((Iem+’

II*+ Ile”‘ll*)+CA 2 kllrmllm

m=O

+ 2

m=O

y<

5

I-lm=O

( E

P,d+(e~+l-P

+ ey-p))A+(ey+l

+ ey),

(3.16)

p=o

for 0 G n G N - 1. It follows from Lemmas 1 and 4 that the last term of (3.16) is nonpositive. Moreover, it follows from (3.11) that 2 m=O

k )Irm IIm< C 2

( k3t,3/: + k5’* + kh*)

m=O

< Ck3/*

k m=O

(m +

1)-3’2 + Ck3’* + Ch* + O(k3’* + h*).

(3.17)

T. Tang / Partial integro-differential equations

Therefore,

317

(3.16) leads to n+l II e n+1112
II em )I 2 +

c

C(k3/* j- h*)A.

(3.18)

m=O An application of the discrete Gronwall lemma for (3.18) yields that

II en+q*
O,
(3.19)

The above inequality implies that A* 6 C(k 3/2 + h2)A, which is equivalent to A = O(k3/* + h*). Hence, we have obtained the following convergence result. Theorem 1. Assume that the solution of Problem I (Problem III) satisfies the smoothness requirements stated in the Introduction, and that (U’, . . . , UN> are solutions of (2.19)-(2.21) with the discretized operator given by (2.22) ((2.23) for Problem II). Zf /? in (2.13) satisfies (3.91, then as h and k tend to zero independently,

max (IUn--~“II

l
=O(k3/*+h2).

(3.20)

4. Numerical experiment For numerical verification of the above theorem we consider the following example, u, = {(t-s) u(0, t)=u(l,

-1’2ux,(x,

s)

t)=O,

u(x, 0) = sin(nx),

ds,

O,
(4.1) (4.2) (4.3)

The solution of this problem is [14] u(x, t> =M(TT 5/2t3/2>sin(~~), where M denotes the entire function M(Z) = E (- l)nT($z n=O

+ 1))‘z”.

(4.4)

It is well known that the central difference approximations yield a convergence order of (at most) 2, so the estimate of Theorem 1 for the space discretization is the best possible. The main purpose of this section is to verify that the error bound for the time discretization is the optimal one. In other words, we shall show that the predicted convergence order i in time is the best possible. For the numerical scheme (2.19)-(2.21) we obtain a linear algebraic system which takes the form MU n+i =F” ,

O
(4.5)

T. Tang / Partial integro-differential

318 Table 1 Errors and convergence

rates

N

Error

5 10 20 40

7.01D-2 2.49D-2 8.66D-3 3.05D-3

Rate (7.89D-2) (2.71D-2) (9.80D-3) (3.52D-3)

where F” E U&-l is independent

M=

‘a b 0 0

b a b 0

. .

. . .

\*

equations

0 b a b

1.49 (1.54) 1.52 (1.47) 1.51 (1.48)

of Un+’ and M is an (J - 1) x 0 0 b a

0 0 0 b

(J -

1) matrix given by

.** 0.. a** -0.

.

where a = 1 + k&/h2,

b = - ;k&/h2.

(4.6)

It can be shown that the solution Un+’ of (4.5) satisfies u n+l = PHPF”,

n > 0,

where P = ( pii) is a (J - 1) x (J - 1) symmetrical matrix with diagonal matrix which takes the form H =

diag

(4.7) pij =

\/2/Jsin(zjv/J)

and H is a

a + 2b cos y

Equation (4.7) provides an explicit form for the numerical solutions of problem (4.G(4.3). In Table 1 we list the errors (max, 4 n QN 11e” 11)and computed rates of convergence when uniform stepsizes h = k = T/N are used. In the calculations we have set /3 = 0.1 (see (3.9)) and T = 0.5. The numerical results reflect a convergence rate = : in time, which is in good agreement with the theoretical prediction of Theorem 1. Also in Table 1 we list the errors and rates for numerical solutions with p = 0; the values appear between parentheses. Although the convergence rate in time is seen to be about 5, it is unclear whether or not Theorem 1 can be extended to the case when /3 = 0.

Acknowledgement

I am grateful to Professor Zhen-huan Teng and the referees for many helpful comments and suggestions which have improved the presentation and content of this paper.

T. Tang / Partial integro-differential equations

319

References [l] P. Camino, A numerical method for a partial differential equation with memory, in: Proceedings 9th CEDYA, Valladolid (1987) 107-112. [2] R.M. Christensen, Theory of Viscoelasticity (Academic Press, New York, 1971). [3] I. Christie, Private communications to J.M. Sanz-Serna (1986). [4] M.E. Gurtin and A.C. Pipkin, A general theory of heat conduction with finite wave speed, Arch. Rational Mech. Anal. 31 (1968) 113-126. [5] A.S. Lodge, M. Renardy and J.A. Nohel, Viscoelasticity and Rheology (Academy Press, New York, 1985). [6] J.C. Lopez-Marcos, A difference scheme for a nonlinear partial integro-differential equation, SOlMJ. Numer. Anal. 27 (19901 20-31. [7] R.K. Miller, An integro-differential equation for rigid heat conductors with memory, J. Math. Anal. Appl. 66 (1978) 313-332. [8] B. Neta, Numerical solution of a nonlinear integro-differential equation, J. Math. Anal. Appl. 89 (1982) 598-611. [9] B. Neta and J.O. Igwe, Finite differences versus finite elements for solving nonlinear integro-differential equations, J. Math. Anal. Appl. 112 (1985) 607-618. [lo] J.M. Ortega, S.H. Davis, S. Rosemblat and W.L. Kath, Bifurcation with memory, SIAMJ. Appl. Math. 46 (1986) 171-188. [ll] P. Linz, Analytical and Numerical Methods for Volterra Equations, SIAM Studies in Applied Mathematics 7 (SIAM, Philadelphia, PA, 1985). [12] M. Renardy, Mathematical analysis of viscoelastic flows, Ann. Rec. Fluid Mech. 21 (1989) 21-36. [13] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Vulue Problems (Interscience Publishers, New York, 1967). [14] J.M. Sanz-Serna, A numerical method for a partial integro-differential equation, SLAM J. Numer. Anal. 25 (1988) 319-327. [15] Y. Yan and G. Fairweather, Orthogonal spline collocation methods for some partial integro-differential equations, SUM J. Numer. Anal. 29 (1992) 755-768.