Convergence of the scheme with weights for the numerical solution of a heat conduction equation with delay for the case of variable coefficient of heat conductivity

Applied Mathematics and Computation 256 (2015) 83–93

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Convergence of the scheme with weights for the numerical solution of a heat conduction equation with delay for the case of variable coefficient of heat conductivity Andrey Lekomtsev, Vladimir Pimenov ⇑ Department of Computational Mathematics, Ural Federal University, 19 Mira street, Ekaterinburg 620002, Russia

a r t i c l e

i n f o

Keywords: Parabolic equations Delay Scheme with weights Variable coefficient of heat conductivity

a b s t r a c t One-dimensional parabolic equations with delay effects in the time component for the case of variable coefficient of heat conductivity are considered. The scheme with weights is constructed for the numerical solution of these equations. The order of approximation error for the constructed scheme, stability, and order of convergence are investigated. Ó 2015 Elsevier Inc. All rights reserved.

1. Introduction In the simulation of dynamics processes often simultaneously two effects are found: distribution of parameters in space and heredity in time [1,2]. Such equations include diffusion equation with delay of general type, constant or variable, concentrated or distributed. Numerical methods for solving such equations were considered in many papers, see for example [3–8]. Options of the method of lines, in which discretization is carried out only on state variables, reduce problem to the numerical solving systems of functional differential equations [4,5]. However at such discretization there are the stiff systems which rigidity increases at increase in number of points of discretization on a state. In a number of works grid methods were investigated, see for example [6] and the bibliography to this work. In these works the main idea consists in replacement of an initial problem its discrete analog and introduction of intermediate interpolation space. At such approach there are following problems: solvability of the received system of the difference equations, stability of algorithms and their convergence. Thus, since the effect of delay can enter nonlinearly, the nonlinear system of the difference equations of large dimension turns out. In the present paper we consider the quasilinear parabolic equation with the effect of delay entering only in inhomogeneous term. By linear part algorithms known for the parabolic equations are used (schemes with weights). For the accounting of delay along with idea of carrying out interpolation of discrete prehistory for creation of the differential model adequate to the initial equation the idea of extrapolation of prehistory is used. Extrapolation allows to avoid the nonlinear equations in spaces of large dimension. As a result on each temporary layer we obtain a linear system of the equations which can be effectively solved by the sweep method. This approach allows us to construct effective algorithms that can be the basis of the corresponding program packages. Justification of convergence of these algorithms is carried out by application both the general theory of differential schemes [9], and the general technique of research of difference schemes of the solution of functional differential equations [10,11]. The last technique uses the ideas of work [12] developed for the ordinary differential equations. In works [13,14] this technique was applied to research of numerical methods of the solution of the parabolic ⇑ Corresponding author. E-mail addresses: [email protected] (A. Lekomtsev), [email protected] (V. Pimenov). http://dx.doi.org/10.1016/j.amc.2014.12.149 0096-3003/Ó 2015 Elsevier Inc. All rights reserved.

84

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

equations in case of constant coefficients, with one and two spatial variables and with effect of delay. In the present paper existence of dependence from time of coefficient of heat conductivity demanded modification of this general technique of research of difference schemes of the solution of the functional differential equations. The structure of paper is the following. In the second section the formulation of the problem is provided and the main assumptions become. In the third section discretization of a task is carried out. Interpolation’s and extrapolation’s constructions are for this purpose given and the family of schemes with weights is described. The concept of the residual of a method without interpolation is entered and its order is studied. The fourth section contains the main theoretical base: the description of the general difference scheme with effect of heredity. Concepts of a temporary grid are entered; discrete model in each timepoint as element of finite-dimensional normed space, prehistory of discrete model, interpolation space, obvious step-by-step formula, starting values, function of exact values. In a step-by-step formula the operator of transition is allocated, its properties define stability of the scheme. The main statement the theorem of an order of convergence of the scheme which depends on an order residual with interpolation is proved. In the fifth section the embedding of the algorithm described in the third section is carried out to this scheme. Connection of an order of residual without interpolation, an order of interpolation and an order of residual with interpolation is established. The main attention is paid to a conclusion of the condition guaranteeing stability. As a result of the stated results, the theorem of orders of convergence of the methods described in this paper turns out. The order of convergence of the methods is quadratic in time and space steps. Results of numerical experiments are given in the last section: test example with the distributed delay and a model example from population dynamics with constant delay. The equations contain variable coefficients in both examples.

2. Problem statement and main assumptions Let us consider a one-dimensional heat conduction equation with delay in the form:

@u @2u ¼ a2 ðtÞ þ f ðx; t; uðx; tÞ; ut ðx; ÞÞ; @t @x

ð1Þ

here, x 2 ½0; X  R1 and t 2 ½t0 ; h  R1 are independent space and time variables, respectively; uðx; tÞ is the required function; ut ðx; Þ ¼ fuðx; t þ sÞ; s 6 s < 0g is the prehistory of the required function by the time t; and s is the value of delay. Let the initial conditions

uðx; tÞ ¼ uðx; tÞ;

x 2 ½0; X; t 2 ½t0  s; t 0 

ð2Þ

and the boundary conditions

uð0; tÞ ¼ uðX; tÞ ¼ 0;

t 2 ½t0 ; h

ð3Þ

be given. Problem (1)–(3) is the simplest boundary-value problem for the one-dimensional heat conduction equation with delay effect in the general form for the case of variable coefficient of heat conductivity. We assume that coefficient a2 ðtÞ satisfies the following condition

0 < c1 6 a2 ðtÞ 6 c2 ;

t 2 ½t 0 ; h:

ð4Þ

We assume also that the functions aðtÞ; uðx; tÞ and the functional f are such that problem (1)–(3) has a unique solution uðx; tÞ in the classical sense. Moreover, we assume that the function uðx; tÞ has the smoothness required in the reasonings below. We denote by Q ¼ Q ½s; 0Þ the set of functions uðsÞ that are piecewise continuous on the half-open interval ½s; 0Þ with a finite number of points of discontinuity of the first kind and right continuous at the points of discontinuity. In addition, the functions uðsÞ have a finite left-hand limit at zero. We define the norm of a function on Q by the relation kuðÞkQ ¼ sups6s<0 juðsÞj. We additionally assume that the functional f ðx; t; u; v Þ is given on ½0; X  ½t 0 ; h  R  Q and is Lipschitz with respect to the last two arguments; i.e., there exists a constant Lf such that, for all x 2 ½0; X; t 2 ½t0 ; h; u1 2 R; u2 2 R; v 1 ðÞ 2 Q , and v 2 ðÞ 2 Q , the following inequality holds:

jf ðx; t; u1 ; v 1 ðÞÞ  f ðx; t; u2 ; v 2 ðÞÞj 6 Lf ðju1  u2 j þ kv 1 ðÞ  v 2 ðÞkQ Þ:

3. Scheme with weights Let us divide the interval ½0; X of variation of the space variable into parts with step h ¼ X=N, introducing the points xi ¼ ih; i ¼ 0; . . . ; N, and the interval ½t 0 ; h of variation of the time variable into parts with step D > 0, introducing the points tk ¼ t 0 þ kD; k ¼ 0; . . . ; M. We assume that the value s=D ¼ m is an integer.

85

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

We denote by uik approximations of the functions uðxi ; t k Þ at the nodes. We use the notation ak ¼ aðtk Þ and akþ1=2 ¼ aðtk þ D=2Þ. For every fixed i ¼ 0; . . . ; N we introduce the discrete prehistory by the time tk ; k ¼ 0; . . . ; M: fuil gk ¼ fuil ; k  m 6 l 6 kg. An interpolation–extrapolation operator of the discrete prehistory is, by definition, a mapping I that takes each time tk ðk ¼ 0; . . . ; MÞ and each discrete prehistory fuil gk by the time tk to a function v ik ðÞ 2 Q ½s; D. In what follows, we omit the index k at the function v ik ðÞ. Definition 1. We will say that an interpolation–extrapolation operator has order of error Dp on the exact solution if there exist constants C 1 and C 2 such that, for all i; k, and t 2 ½t k  s; t kþ1  the following inequality holds:

kv i ðtÞ  uðxi ; tÞk 6 C 1 max juil  uðxi ; t l Þj þ C 2 Dp : km6l6k

For example, the piecewise linear interpolation

v i ðtk þ sÞ ¼ ððtl  tk  sÞuil1 þ ðtk þ s  tl1 Þuil Þ=D;

t l1 6 t k þ s 6 t l ;

with extrapolation by extension

v i ðtk þ sÞ ¼ ððsÞuik1 þ ðD þ sÞuik Þ=D;

tk 6 t k þ s 6 t kþ1 ;

has second order [11]. As an approximation of the functional f ðx; t; uðx; tÞ; ut ðx; ÞÞ in the grid nodes we will consider the functional Fðv ðÞÞ. We denote by F ik ðv i ðÞÞ value of the functional Fðv ðÞÞ at the node ðxi ; tk Þ. For F ik ðv i ðÞÞ in the simplest case, for example, we can take the functional f ðxi ; tkþ1=2 ; uikþ1=2 ; v itkþ1=2 ðÞÞ. In this case the functional F ik ðv i ðÞÞ is defined on Q ½s; D and is Lipschitz

v ðÞ with the constant Lf . Argument uikþ1=2 in the functional f ðxi ; tkþ1=2 ; uikþ1=2 ; v it ðÞÞ is calcuof extrapolation, but argument v it ðÞ is calculated explicitly by means of interpolation and

with respect to the variable lated explicitly by means

kþ1=2

kþ1=2

extrapolation. Functional f ðx; t; uðx; tÞ; ut ðx; ÞÞ in the grid nodes can be approximated in another way. Other ways of approximation the functional f ðx; t; uðx; tÞ; ut ðx; ÞÞ will be considered after the theorem about the order of the residual. For 0 6 s 6 1, consider the family of methods

uiþ1  2uikþ1 þ ui1 uikþ1  uik uiþ1  2uik þ uki1 kþ1 ¼ sa2kþ1 kþ1 þ ð1  sÞa2kþ1 k þ F ik ðv i ðÞÞ; i ¼ 1;. .. ;N  1; k ¼ 0; .. .; M  1; 2 2 2 2 D h h

ð5Þ

with the initial conditions

ui0 ¼ uðxi ; t0 Þ;

i ¼ 0; . . . ; N;

v i ðtÞ ¼ uðxi ; tÞ; t < t0 ;

i ¼ 0; . . . ; N

and the boundary conditions

u0k ¼ uNk ¼ 0;

k ¼ 0; . . . ; M:

For s ¼ 0, we obtain an explicit scheme; for other s (0 < s 6 1) and any fixed k, system (5) is a linear tridiagonal system with respect to uikþ1 with diagonal dominance and can be effectively solved by the sweep method. Note that the value of the functional F ik ðv i ðÞÞ is calculated explicitly by means of interpolation and extrapolation. Definition 2. The residual of the method is, by definition, the value

Wik ¼

uðxi ; t kþ1 Þ  uðxi ; tk Þ uðxiþ1 ; tkþ1 Þ  2uðxi ; t kþ1 Þ þ uðxi1 ; t kþ1 Þ  sa2kþ1  ð1  sÞa2kþ1 2 2 2 D h uðxiþ1 ; tk Þ  2uðxi ; tk Þ þ uðxi1 ; tk Þ i   F k ðutk ðxi ; ÞÞ: 2 h

ð6Þ

One can find the order of the residual of the method for specific s and F by using the Taylor expansion of the function uðx; tÞ (under the corresponding smoothness conditions). For example, the following statement is valid, which completely corresponds to the similar statement for the scheme with weights without delay [9]. Theorem 1. For F ik we take the functional f ðxi ; tkþ1=2 ; uikþ1=2 ; v itkþ1=2 ðÞÞ. Suppose that the exact solution uðx; tÞ of problem (1)–(3) is three times continuously differentiable in t and four times continuously differentiable in x; the second derivatives of the solution 2 with respect to x are twice continuously differentiable in t. Then, the residual of the scheme with weights has order Dbs þ h ; where bs ¼ 1 when s – 1=2, and bs ¼ 2 when s ¼ 1=2.

86

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

Proof. We expand the exact solution of problem (1)–(3) by Taylor’s formula in a neighborhood of the point ðxi ; t kþ1=2 Þ under the assumption that derivatives of the corresponding orders are bounded. Then,

Wik ¼ ¼

@u @2u @2u 2 ðxi ; tkþ1 Þ þ OðD2 Þ  sa2kþ1  2 ðxi ; tkþ1 Þ  ð1  sÞa2kþ1  2 ðxi ; t k Þ þ Oðh Þ  f ðxi ; tkþ1 ; uðxi ; t kþ1 Þ; utkþ1 ðxi ; ÞÞ 2 2 2 2 2 @t @x1 @x1 2 @u @2u D @3u 2 ðxi ; tkþ1 Þ  a2kþ1  2 ðxi ; t kþ1 Þ þ ð1  2sÞa2kþ1  ðxi ; t kþ1 Þ  f ðxi ; tkþ1 ; uðxi ; tkþ1 Þ; utkþ1 ðxi ; ÞÞ þ OðD2 þ h Þ: 2 2 2 2 2 2 2 @t 2 @x1 @t@x21 2

Since uðx; tÞ is the solution of the equation (1), we obtain

Wik ¼

D @3u 2 ð1  2sÞa2kþ1  ðxi ; tkþ1 Þ þ OðD2 þ h Þ: 2 2 2 @t@x21

This relation implies the conclusion of the theorem. h The selection of the functional F ik ðv i ðÞÞ is directly related to the order of the residual. Following are some ways to choose of the functional F ik ðv i ðÞÞ. If we consider F ik ðv i ðÞÞ ¼ f ðxi ; tk ; uik ; v itk ðÞÞ or F ik ðv i ðÞÞ ¼ f ðxi ; t kþ1 ; uikþ1 ; v itkþ1 ðÞÞ, then the residual of 2 the scheme with weights has order D þ h . 2 For F ik ðv i ðÞÞ ¼ 12 f ðxi ; t k ; uik ; v itk ðÞÞ þ 12 f ðxi ; t kþ1 ; uikþ1 ; v itkþ1 ðÞÞ we obtain that the residual has order D2 þ h when s ¼ 1=2. i 1 h2 i If we consider sk ¼ 2  12Da2 and take F k ðv ðÞÞ the following way k

F ik ð i ðÞÞ

v

5 1 i1 iþ1 iþ1 ¼ f ðxi ; tkþ1=2 ; uikþ1=2 ; v itkþ1=2 ðÞÞ þ ðf ðxi1 ; t kþ1=2 ; ui1 kþ1=2 ; v t kþ1=2 ðÞÞ þ f ðxiþ1 ; t kþ1=2 ; ukþ1=2 ; v t kþ1=2 ðÞÞÞ; 6 12 4

then the residual of the scheme with weights has order D2 þ h [13]. We denote by eik ¼ uðxi ; tk Þ  uik the value of error of the method at the nodes. q

Definition 3. We will say that the method converges with order h 1 þ Dq2 if there exists a constant C independent of eik ; h, q and D such that jeik j 6 Cðh 1 þ Dq2 Þ for all i ¼ 0; . . . ; N and k ¼ 0; . . . ; M. In view of the nonlinear dependence of the functional f (and, consequently, F) on the state and its prehistory, usual methods of investigating the order of convergence of difference schemes [9] are not applicable. However, to investigate the convergence of the schemes in this problem, as in other evolutionary problems with delay effect, we can apply the apparatus of abstract schemes with aftereffect developed earlier in [10] for the case of functional differential equations. Let us describe the main points of this apparatus as applied to our case. 4. General difference scheme with aftereffect and its order of convergence Let a segment ½t 0 ; h be given, and let a number s > 0 be the value of the delay. A step of a grid is, by definition, a number D > 0, such that s=D ¼ m is an integer; fDg is the set of steps. A (uniform) grid is, by definition, a finite set of numbers

RD ¼ ft i ¼ t0 þ iD 2 ½t0  s; h; i ¼ m; . . . ; Mg: þ We use the notation R D ¼ ft i 2 RD ; i 6 0g; RD ¼ ft i 2 RD ; i P 0g. A discrete model is defined as a grid function t i 2 RD ! yðt i Þ ¼ yi 2 Y; i ¼ m; . . . ; M, where Y is a q - dimensional normed space with norm k  kY . We will assume that the dimension q of the space Y depends on a number h > 0. For n P 0, the set fyi gn ¼ fyi 2 Y; i ¼ n  m; . . . ; ng is called the prehistory of the discrete model by the time tn . Let V be a linear normed space with norm k  kV (an interpolation space). A mapping I : Iðfyi gn Þ ¼ v 2 V is called an operator of the interpolation of the discrete prehistory of the model. We will say that the interpolation operator satisfies the Lipschitz condition if there exists a constant LI such that, for all prehistories fy1i gn and fy2i gn of the discrete model, the following inequality holds:

kv 1  v 2 kV 6 LI max ky1i  y2i kY : nm6i6n

Starting values of the model are defined by the function R D ! Y:

yðti Þ ¼ yi ;

i ¼ m; . . . ; 0:

ð7Þ

The formula of the advance of the model by a step is, by definition, the algorithm

ynþ1 ¼ Sn yn þ DUðt n ; Iðfyi gn Þ; DÞ; þ D

ð8Þ

where U : R  V  fDg ! Y is the function of advance by a step and the transition operator Sn : Y ! Y is a linear operator. Thus, a discrete model (a numerical method; in what follows, simply a method) is defined by starting values (7), formula of advance by a step (8) and an interpolation operator. We assume that the function Uðt n ; v ; DÞ in (8) is Lipschitz with respect

87

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

to the second argument; i.e., there exists a constant LU such that, for all t n 2 Rþ D ; D 2 fDg and inequality holds:

v1; v2 2 V n

the following

kUðt n ; v 1 ; DÞ  Uðtn ; v 2 ; DÞkY 6 LU kv 1  v 2 kV : The function of exact values is, by definition, the mapping

Zðti ; DÞ ¼ zi 2 Y;

i ¼ m; . . . ; M:

We assume that the specification of the function of exact values is a consequence of the specification of an exact solution to problem (1)–(3). p We will say that starting values of the model have order Dp1 þ h 2 if there exists a constant C such that p

kzi  yi kY 6 CðDp1 þ h 2 Þ;

i ¼ m; . . . ; 0: p2

We will say that the method converges with order Dp1 þ h

if there exists a constant C such that

p2

p1

kzn  yn kY 6 CðD þ h Þ for all n ¼ m; . . . ; M. Method (8) is called stable if kSn k 6 1 for all n ¼ 0; . . . ; M  1, where

kSn k ¼ sup y–0

kSn ykY ; kykY

n ¼ 0; . . . ; M  1:

An error of approximation with interpolation (a residual) is, by definition, the grid function

dn ¼ ðznþ1  Sn zn Þ=D  Uðtn ; Iðfzi gnnm Þ; DÞ;

n ¼ 0; . . . ; M  1:

ð9Þ p1

We will say that method (8) has order of error of approximation with interpolation D such that

p2

þh

if there exists a constant C

p

kdn kY 6 CðDp1 þ h 2 Þ for all n ¼ 1; . . . ; M. The following main theorem is valid. Theorem 2. Suppose that method (8) is stable, the function U satisfies the Lipschitz condition with respect to the second argument, p the interpolation operator I satisfies the Lipschitz condition, the starting values have order Dp1 þ h 2 , and the error of p4 p3 approximation with interpolation has order D þ h , where p1 > 0; p2 > 0; p3 > 0 and p4 > 0. Then, the method converges and minfp2 ;p4 g the order of the convergence is at least Dminfp1 ;p3 g þ h . Proof. We use the notation dn ¼ zn  yn ; n ¼ m; . . . ; M ; then, for n ¼ 0; . . . ; M  1, we have

d n þ Ddn ; dnþ1 ¼ Sn dn þ Db

ð10Þ

where

b d n ¼ Uðt n ; Iðfzi gn Þ; DÞ  Uðtn ; Iðfyi gn Þ; DÞ: The assumptions that the mappings U and I are Lipschitz imply that

d n kY 6 K max fkdi kY g; kb

ð11Þ

nm6i6n

where K ¼ LU LI . We use the notation Sn;l ¼ Sn  Sn1  . . .  Sl for n P 0 and 0 6 l 6 n. If l > n, then Sn;l ¼ E is identity operator. It follows from (10) that

dnþ1 ¼ Sn;0 d0 þ D

n n X X d j þ D Sn;jþ1 dj : Sn;jþ1 b j¼0

ð12Þ

j¼0

From (11) and (12), and the definition of the stability of the method (8), we have

kdnþ1 kY 6 K D

n X j¼0

max fkdi kY g þ kd0 kY þ ðh  t 0 Þ max fkdi kY g:

jm6i6j

06i6M1

ð13Þ

We use the notation

R0 ¼ max fkdi kY g; m6i60

R ¼ max fkdi kY g; 06i6M1

D ¼ R0 þ ðh  t0 ÞR;

ð14Þ

88

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

Then, we can write estimate (13) as follows:

kdnþ1 kY 6 K D

n X j¼0

max fkdi kY g þ D:

ð15Þ

jm6i6j

Let us prove the estimate

kdn kY 6 Dð1 þ K DÞn

ð16Þ

by induction on n ¼ 1; . . . ; M. Induction base. If we set n ¼ 0 in (15), then

kd1 kY 6 K DR0 þ D 6 K DD þ D ¼ Dð1 þ K DÞ: Induction step. Let estimate (16) be valid for all indices from 1 to n. Let us show that the estimate is also valid for n þ 1. Let us fix j 6 n. Let i0 ¼ i0 ðjÞ be the index for which maxjm6i6j fkdi kY g is attained. The following two situations are possible: 1. i0 6 0; then, maxjm6i6j fkdi kY g ¼ kdi0 kY 6 R0 6 Dð1 þ K DÞj . 2. 1 6 i0 6 j; then, by the induction assumption,

max fkdi kY g ¼ kdi0 kY 6 Dð1 þ K DÞi0 6 Dð1 þ K DÞj :

jm6i6j

Thus, the following estimate is valid in any case:

max fkdi kY g 6 Dð1 þ K DÞj :

jm6i6j

The estimate obtained and (15) imply

kdnþ1 kY 6 K D

n n X X ð1 þ K DÞn  1 Dð1 þ K DÞj þ D ¼ D þ K DD þ K D Dð1 þ K DÞj ¼ Dð1 þ K DÞ þ K DDð1 þ K DÞ 1 þ KD  1 j¼0 j¼1

¼ Dð1 þ K DÞ þ Dð1 þ K DÞðð1 þ K DÞn  1Þ ¼ Dð1 þ K DÞnþ1 : Thus, estimate (16) is proved. From this estimate, we obtain the estimate

kdn kY 6 D expðKðh  t 0 ÞÞ:

ð17Þ

Since, by definition (14) of the value D, the inequality

D 6 CðDminfp1 ;p3 g þ h

minfp2 ;p4 g

Þ

holds, (17) implies the conclusion of the theorem. h

5. The embedding into the general difference scheme with aftereffect Let us embed scheme (5) into the general difference scheme with aftereffect. For every t k 2 RD , we define the values of the T discrete model by the vector yk ¼ ðu0k ; u1k ; . . . ; uNk Þ 2 Y, where Y is a vector space of dimension q ¼ N þ 1, and T is the transposition symbol. In the space Y, we introduce an operator A:

Ai ðtkþ1 Þuik ¼ a2 ðt kþ1 Þ 2

uiþ1  2uik þ ui1 k k

2

2

h

;

i ¼ 1; . . . ; N  1;

A0 ðtkþ1 Þu0k ¼ AN ðt kþ1 ÞuNk ¼ 0; 2

2

T

Aðt kþ1 Þyk ¼ ðA0 ðt kþ1 Þu0k ; A1 ðt kþ1 Þu1k ; . . . ; AN ðtkþ1 ÞuNk Þ : 2

2

2

2

Then, we can write system (5) in the form of the equation

ykþ1  yk þ sAðtkþ1 Þykþ1 þ ð1  sÞAðt kþ1 Þyk ¼ F k ðv ðÞÞ; 2 2 D T

ð18Þ

where F k ðv ðÞÞ ¼ ðF 0k ðv 0 ðÞÞ; F 1k ðv 1 ðÞÞ; . . . ; F Nk ðv N ðÞÞÞ , v ðÞ ¼ Iðful gk Þ 2 Q q ½s; D. Here, the interpolation space V ¼ Q q ½s; D is the space of q-dimensional vector functions every component of which belongs to the space Q ½s; D. Using the identity

ykþ1 ¼ yk þ D

ykþ1  yk D

89

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

and introducing the operator

Bðtkþ1 Þ ¼ E þ DsAðt kþ1 Þ 2

2

(E is the identity operator), we bring Eq. (18) to the canonical form [9]:

Bðtkþ1 Þ 2

ykþ1  yk þ Aðtkþ1 Þyk ¼ F k ðv ðÞÞ; 2 D

k ¼ 0; . . . ; M  1:

ð19Þ

In view of condition (4) for all tkþ1=2 the operator Aðt kþ1=2 Þ is self-adjoint and positive [9] in the sense of the scalar product T T of vectors y ¼ ðy0 ; y1 ; . . . ; yN Þ 2 Y; u ¼ ðu0 ; u1 ; . . . ; uN Þ 2 Y:

ðy; uÞ ¼

N X yi ui h: i¼0

pffiffiffiffiffiffiffiffiffiffiffi The norm in the space Y is defined as kyk ¼ ðy; yÞ. E (identity operator) is self-adjoint and positive. Consequently, for all tkþ1 the operator Bðtkþ1 Þ is self-adjoint and positive. 2 2 Since Bðtkþ1=2 Þ is a positive operator in a finite-dimensional Hilbert space for all t kþ1=2 , there exists B1 ðtkþ1=2 Þ for all tkþ1=2 (see [9]). Consequently, Eq. (19) can be reduced to the explicit form

ykþ1 ¼ Sðtkþ1 Þyk þ DUðtk ; Iðfyl gk Þ; D; t kþ1 Þ; 2

2

k ¼ 0; . . . ; M  1;

ð20Þ

where the transition operator Sk ¼ Sðt kþ1=2 Þ is defined by the formula Sk ¼ E  DB1 ðtkþ1=2 ÞAðt kþ1=2 Þ and the function of advance by a step is defined by the formula Uðt k ; v ; D; t kþ1 Þ ¼ B1 ðtkþ1 ÞF k ðv ðÞÞ 2

2

Let us investigate the stability of the scheme obtained, i.e. the fulfillment of condition kSk k 6 1 for k ¼ 0; . . . ; M  1. To this end, along with Eqs. (19) and (20), consider the homogeneous difference scheme in the canonical and explicit forms

Bðtkþ1 Þ 2

ykþ1  yk þ Aðtkþ1 Þyk ¼ 0; 2 D

ykþ1 ¼ Sk yk ;

k ¼ 0; . . . ; M  1:

k ¼ 0; . . . ; M  1:

ð21Þ

The operators Aðtkþ1=2 Þ and Bðtkþ1=2 Þ are permutable for all t kþ1=2 [15]. To obtain the conditions of stability we use the fact 2 that for all t kþ1=2 inequality holds kAðtkþ1=2 Þk 6 4c . Let us consider h2 2

Bðtkþ1 Þ  2

D D 1 1 h Aðt kþ1 Þ ¼ E þ sDAðt kþ1 Þ  Aðt kþ1 Þ ¼ E þ ðs  ÞDAðtkþ1 Þ P 0 under s P  : 2 2 2 2 2 2 4c2 D 2 2

ð22Þ

In view of the self-adjoint, positive and permutable operators Aðt kþ1=2 Þ and Bðtkþ1=2 Þ, and condition (22) we can apply Corollary 1 of Theorem 15 (see [15]). We obtain

kykþ1 k 6 kyk k;

k ¼ 0; . . . ; M  1:

Then, for the equivalent homogeneous equation (21) we obtain that

kSk yk k 6 kyk k;

k ¼ 0; . . . ; M  1;

consequently

kSk yk k 6 1; yk –0 kyk k

kSk k ¼ sup

2

k ¼ 0; . . . ; M  1; under s P

1 h  : 2 4c2 D

2

Thus, under condition s P 12  4ch2 D the scheme (20) is stable. We define the function of exact values by the relations

zk ¼ ðuðx0 ; t k Þ; uðx1 ; t k Þ; . . . ; uðxN ; t k ÞÞT 2 Y;

k ¼ 0; . . . ; M  1:

Starting values of the model can be taken equal to the function of exact values:

yj ¼ zj ¼ ðuðx0 ; t j Þ; uðx1 ; t j Þ; . . . ; uðxN ; t j ÞÞT ;

j ¼ m; . . . ; 0:

The definition of the residual without interpolation (6) in the scheme with weights for the heat conduction equation and the definition of the residual with interpolation (9) in the general scheme are essentially different. However, the following theorem connects these two definitions. p

Theorem 3. Suppose that the residual in the sense of (6) has order Dp1 þ h 2 , the functions F ik are Lipschitz, and the interpolation– extrapolation operator I is Lipschitz and has order of error Dp0 on the exact solution. Then, the residual with interpolation has order Dminfp1 ;p0 g þ hp2 .

90

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

Proof. Let us write the definition of the residual in the sense of (9)

dk ¼ ðzkþ1  Sk zk Þ=D  Uðtk ; Iðfzl gk Þ; D; t kþ1 Þ;

k ¼ 0; . . . ; M  1:

2

Since Bðt kþ1=2 Þ is a linear operator and there exists B1 ðtkþ1=2 Þ, the operator B1 ðtkþ1=2 Þ is also linear. In a finite- dimensional space, a linear operator is bounded; consequently, Bðtkþ1=2 Þ and B1 ðtkþ1=2 Þ are bounded for all t kþ1=2 in view of condition (4). For every component of the vector z, we obtain i jdk j

 " #   1   1 uðxi ; t kþ1 Þ  Sik uðxi ; t k Þ   i  i i i i i i ¼ B ðt kþ1 ÞB ðtkþ1 Þdk  ¼ B ðt kþ1 Þ B ðt kþ1 Þ  B ðt kþ1 ÞU ðtk ; Iðfuðxi ; tl Þgk Þ; D; tkþ1 Þ  2 2 2 2 2 2   D   2 3 1   uðxi ; tkþ1 Þ  ðE  DBi ðt kþ1 ÞAi ðtkþ1 ÞÞuðxi ; t k Þ   1 1 2 2 ¼ Bi ðt kþ1 Þ4Bi ðt kþ1 Þ  Bi ðt kþ1 ÞBi ðt kþ1 ÞF ik ðIðfuðxi ; t l Þgk ÞÞ5; 2 2 2 2 D  

where k ¼ 0; . . . ; M  1. Let us add and subtract F ik ðutk ðxi ; ÞÞ. We obtain

    1  uðxi ; t kþ1 Þ  uðxi ; tk Þ i jdk j ¼ Bi ðt kþ1 Þ Bi ðt kþ1 Þ þ Ai ðt kþ1 Þuðxi ; t k Þ  F ik ðutk ðxi ; ÞÞ  F ik ðIðfuðxi ; tl Þgk ÞÞ þ F ik ðutk ðxi ; ÞÞ ; 2 2 2 D k ¼ 0; . . . ; M  1:

Using the definition of the order of the residual of the method in the sense of (6), the facts that F is Lipschitz with respect to uðÞ and that the operator B1 ðt kþ1=2 Þ is bounded, and the definition of the order of an interpolation–extrapolation operator, we obtain, for k ¼ 0; . . . ; M  1 i

p

p

jdk j 6 C B1 C 1 ðDp1 þ h 2 Þ þ C B1 Lf C 2 Dp0 6 C B1 ðC 1 þ Lf C 2 ÞðDminfp1 ;p0 g þ h 2 Þ: Consequently, the norm of the vector dk satisfies the inequality

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kdk k ¼ ðdk ; dk Þ ¼

!12 !12 N N X X pffiffiffiffiffiffi i 2 p2 p minfp1 ;p0 g ðdk Þ h 6 C B1 ðC 1 þ Lf C 2 ÞðD þh Þ h 6 2X C B1 ðC 1 þ Lf C 2 ÞðDminfp1 ;p0 g þ h 2 Þ: i¼0

i¼0

This relation implies the conclusion of the theorem. h The embedding of the scheme with weights for the heat conduction equation into the general scheme has been carried out. Theorem 2 implies the following statement. 2

p

Theorem 4. Suppose that condition s P 12  4ch2 D holds, residual in the sense of (6) has order Dp1 þ h 2 , the functions F ik are Lipschitz, and the interpolation–extrapolation operator I is Lipschitz and has order of error Dp0 on the exact solution. Then, the p method converges with order Dminfp1 ;p0 g þ h 2 . Using this theorem, we conclude that the scheme with weights with a piecewise linear interpolation and extrapolation by 2 2 2 expansion has order D2 þ h , if condition s ¼ 1=2 holds; and has order D þ h , if conditions s – 1=2 and s P 12  4ch2 D hold.

u(x,t)

15 10 5 0 −5 −10 −15 2.5 2 1.5 1 x

0.5

0

10

5

15

20

t

Fig. 1. The number of grid points in t: M ¼ 300, in x: N ¼ 40. s ¼ 1=2.

91

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

6. Numerical examples All examples we will solve by means the scheme with weights for the parameter s equal s ¼ 1=2 (like Crank–Nicolson scheme). However, the Example 1 will be considered for the parameter s ¼ 0 and s ¼ 1. Example 1. Let us consider the following equation for the case of variable coefficient of heat conductivity with distributed delay in the variable t:

@uðx; tÞ @ 2 uðx; tÞ ¼ a2 ðtÞ  sinðtÞex @t @x2 "  a2 ðtÞ uðx; tÞ þ uðx; t  t=2Þ  cosðt=2Þex þ

Z

#

0

ex uðx; t þ sÞds  e2x ðsinðtÞ  sinðt=2ÞÞ ;

ð23Þ

t=2

where x 2 ½0:5; 2:5; t 2 ½p; 6p; a2 ðtÞ ¼ 2 þ cosðtÞ. The initial conditions:

uðx; tÞ ¼ cosðtÞex ;

x 2 ½0:5; 2:5; t 2 ½p=2; p:

The boundary conditions:

uð0:5; tÞ ¼ e0:5 cosðtÞ;

uð2:5; tÞ ¼ e2:5 cosðtÞ; t 2 ½p; 6p:

The exact solution is uðx; tÞ ¼ cosðtÞex . Below, we give the results of the numerical experiment. The table contains the comparison of a maximum of the module of a difference between the matrices of the exact and approximate solutions for different steps h and D for the parameter s ¼ 1=2. The approximate solution of the Eq. (23) is shown in Fig. 1 for the parameter s ¼ 1=2 and in Fig. 3 for the parameter s ¼ 1. The difference between the approximate and exact solutions of the Eq. (23) is shown in Fig. 2 for the parameter s ¼ 1=2 and in Fig. 4 for the parameter s ¼ 1. When parameter s ¼ 0 then constructed scheme with weights does not give a numerical solution of the Eq. (23), since method becomes unstable.

N = 10 N = 20 N = 40 N = 80

M = 150

M = 300

M = 600

0.0624 0.0501 0.0457 0.0446

0.0034 0.0020 0.0020 0.0020

0.0050 0.0020 0.0020 0.0020

The maximums of the modules of differences between the exact and approximate solutions of Eq. (23) obtained for different steps h and D for s ¼ 1=2. Numerical experiments confirm the theoretical results about what for s ¼ 1=2 the absolute error of the scheme with weights is less than for s ¼ 1.

error x 10

−3

2 1 0 −1 −2 −3 2.5 2 1.5 1 x

0.5

0

10

5

15

20

t

Fig. 2. The number of grid points in t: M ¼ 300, in x: N ¼ 40. s ¼ 1=2.

92

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93 u(x,t)

15 10 5 0 −5 −10 −15 2.5 2 1.5 1 0.5

x

10

5

0

15

20

t

Fig. 3. The number of grid points in t: M ¼ 300, in x: N ¼ 40. s ¼ 1.

error

0.06 0.04 0.02 0 −0.02 −0.04 2.5 2 1.5 1 0.5

x

0

15

10

5

20

t

Fig. 4. The number of grid points in t: M ¼ 300, in x: N ¼ 40. s ¼ 1.

u(x,t) 4500 4000 3500 3000 2500 2000 1500 1000 500 10 5 x

0

0

5

10

15

t

Fig. 5. The number of grid points in t: M ¼ 100, in x: N ¼ 50.

20

A. Lekomtsev, V. Pimenov / Applied Mathematics and Computation 256 (2015) 83–93

93

Example 2. Let us consider the following equation

@uðx; tÞ @ 2 uðx; tÞ ¼ a2 ðtÞ  qðtÞuðx; t  sÞ; @t @x2

ð24Þ

where x 2 ½0; 10; t 2 ½0; 20; s ¼ 0:1; aðtÞ ¼ sinðtÞ þ 2; qðtÞ ¼ cosðtÞ. The initial conditions: 2

uðx; tÞ ¼ 50eðx5Þ xð10  xÞ þ 1000;

x 2 ½0; 10; t 2 ½0:1; 0:

The boundary conditions:

uð0; tÞ ¼ 1000;

uð10; tÞ ¼ 1000; t 2 ½0; 20:

Such Eqs. (24) can be used to model problems of population dynamics with spatial migration [2]. The approximate solution of the Eq. (24) is shown below in Fig. 5. Acknowledgments This work was supported by the Russian Foundation for Basic Research (project No. 13-01-00089) and Act 211 Government of the Russian Federation No. 02.A03.21.0006 of 27 August 2013. We thank the anonymous reviewers for comments which have improved the clarity of the paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

J. Wu, Theory and Applications of Partial Functional Differential Equations, Springer-Verlag, New York, 1996. B. Zhang, Y. Zhou, Qualitative Analysis of Delay Partial Difference Equations, vol. 4, Hindawi Publishing Corporation, New York, 2007. L. Tavernini, Finite difference approximations for a class of semilinear Volterra evolution problems, SIAM J. Numer. Anal. 14 (5) (1977) 931–949. P.J. van der Houwen, B.P. Sommeijer, C.T.H. Baker, On the stability of predictor–corrector methods for parabolic equations with delay, IMA J. Numer. Anal. 6 (1986) 1–23. B. Zubik-Kowal, The method of lines for parabolic differential-functional equations, IMA J. Numer. Anal. 17 (1997) 103–123. K. Kropielnicka, Convergence of implicit difference methods for parabolic functional differential equations, Int. J. Mater. Anal. 1 (6) (2007) 257–277. P. Garcia, M.A. Castro, J.A. Martin, A. Sirvent, Numerical solutions of diffusion mathematical models with delay, Math. Comput. Modell. 50 (2013) 860–868. M.A. Castro, F. Rodriguez, J. Cabrera, J.A. Martin, Difference schemes for time dependent heat conduction models with delay, Int. J. Comput. Math. (2013) 257–277 www.tandfonline.com/doi/abs/10.1080/00207160.2013.77937. A.A. Samarskii, The Theory of Difference Schemes, Nauka, Moscow, 1989 (in Russian). V.G. Pimenov, General linear methods for numerical solving functional-differential equations, Differ. Equ. 37 (2001) 116–127. A.V. Kim, V.G. Pimenov, i-Smooth analysis and numerical methods for solving functional differential equations, Regular and Chaotic Dynamics, Moscow, Izhevsk, 2004 (in Russian). R.D. Skeel, Analysis of fixed-stepsize methods, SIAM J. Numer. Anal. 13 (1976) 664–685. V.G. Pimenov, A.B. Lozhnikov, Difference schemes for the numerical solution of the heat conduction equation with aftereffect, Proc. Steklov Inst. Math. 275 (2011) 137–148. A.V. Lekomtsev, V.G. Pimenov, Convergence of the alternating direction methods for the numerical solution of a heat conduction equation with delay, Proc. Steklov Inst. Math. 272 (1) (2011) 101–118. A.A. Samarskii, A.V. Gulin, Stability of Difference Schemes, URSS, Moscow, 2009 (in Russian).