# Artificial boundary conditions for the Burgers equation on the plane

## Artificial boundary conditions for the Burgers equation on the plane

Applied Mathematics and Computation 286 (2016) 1–14 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepage: ...
Applied Mathematics and Computation 286 (2016) 1–14

Contents lists available at ScienceDirect

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

Artiﬁcial boundary conditions for the Burgers equation on the plane Nadaniela Egidi∗, Pierluigi Maponi Università di Camerino, Scuola di Scienze e Tecnologie, Camerino, MC 62032, Italy

a r t i c l e

i n f o

Keywords: Burgers equation Boundless domain Fourier transform Finite element method

a b s t r a c t The numerical solution of the initial value problem for the two-dimensional Burgers equation on the whole plane is considered. Usual techniques, like ﬁnite difference methods and ﬁnite element methods cannot be directly applied for the solution of this problem, because the corresponding domain is unbounded. We propose a new method to overcome this difﬁculty. The eﬃciency of the proposed method is tested by several numerical examples. © 2016 Elsevier Inc. All rights reserved.

1. Introduction The Burgers equation is a mathematical model employed in a large number of applications. It was introduced by Bateman [1] and it was later studied by Burgers [2,3], as a mathematical model for turbulence in ﬂuid dynamics. However, after these preliminary studies, it has been considered in various areas of applied mathematics [4], such as for examples: gas dynamics, waves in shallow water, hydromagnetic waves in cold plasma, ion-acoustic waves in cold plasma [5]; modeling of traﬃc ﬂow [6,7]; dynamics of soil water [8]; shock formation in inelastic gases [9]; cosmology [10]; turbulence in ﬂuid dynamics [11,12]. Besides its effective modeling features, it can be used to set-up computational methods in order to deal with more difﬁcult problems as the Navier–Stokes equation, whose analytical properties and approximation techniques are only partially known. The three-dimensional incompressible Navier–Stokes equation is one of the main open problems in Mathematics; the Clay Mathematics Institute [13] has identiﬁed the existence, smoothness and breakdown of the Navier–Stokes solution as one of the millennium problems. In Burgers equation, discontinuities may appear in ﬁnite time, even if the initial condition is smooth [14,15]. These discontinuities give rise to the phenomenon of shock waves with important applications in physics [16]. So that, Burgers equation is a proper model to test numerical algorithms for the simulation in ﬂows with severe gradients or shocks, and it provides a good training model to assess stability and eﬃciency of computational methods. Various numerical techniques have been proposed to solve Burgers equation, such as: ﬁnite difference methods [17–21]; ﬁnite element methods [18,22,23]; Adomian decomposition method [24–27]; integral equations methods [14,23,28,29]; cellular automata methods [30,31]. Finite difference methods and ﬁnite element methods are the most general and effective techniques for the numerical solution of differential equations. However, these methods are able to solve differential equations on bounded domains. In order to deal with unbounded domains, suitable numerical techniques must be considered. In [23,32] two different artiﬁcial boundary methods are proposed. The general idea of artiﬁcial boundary methods is to reduce

Corresponding author. Tel.: +39 347 6644705. E-mail addresses: [email protected], [email protected] (N. Egidi), [email protected] (P. Maponi).

http://dx.doi.org/10.1016/j.amc.2016.04.008 0 096-30 03/© 2016 Elsevier Inc. All rights reserved.

2

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

an initial value problem on an unbounded domain to an equivalent problem on a bounded artiﬁcial domain by introducing suitable conditions on the artiﬁcial boundary. This kind of methods, after the trailblazing paper [33] for the wave equation, has been proposed for the Helmholtz equation [34], for diffusion equations [35], for Maxwell equations [36], and for Navier–Stokes equation [37,38]. In the present paper, we develop the idea described in [23], where the integral formulation of the Burgers equation on the Fourier transform space is used to compute the artiﬁcial boundary condition, and the Galerkin method is used to solve the resulting initial-boundary value problem. In particular, an eﬃcient method for the computation of these boundary conditions is proposed; this method is based on the numerical solution of the above mentioned integral equation by using FFT algorithm. The main idea of the proposed method is to follow the solution of the problem in the original space as well as in the Fourier transform space. In particular, the numerical solutions of these two formulations are computed by using a truncation of the original space and of the conjugate space, respectively, so these solutions provide a complementary information of the original problem in the unbounded domain. This is a simple idea, but, from our experience, it produces an effective technique. Moreover it can be easily generalized to deal with Navier–Stokes equation, by using the corresponding integral formulation on the Fourier transform space [39]. A numerical experiment is used to test the eﬃciency and the accuracy of the proposed method on several numerical examples. In Section 2, the integral formulation of the Burgers equation is derived and the FFT algorithm for its solution is given. In Section 3, we illustrate the Galerkin algorithm used for the numerical solution of Burgers equation on a bounded domain. In Section 4 we give the artiﬁcial boundary condition method proposed to solve Burgers equation on the whole plane. In Section 5 some numerical experiments with the proposed method are presented. In Section 6 conclusions and future developments are given. 2. Integral equation Let N be the set of natural numbers, R the set of real numbers, R+ the set of positive real numbers, and C the set of complex numbers. We denote with ι the imaginary unit. Let d ∈ N, we denote with Rd , Cd the d-dimensional real Euclidean space and the d-dimensional complex Euclidean space, respectively. Let x, y ∈ Rd be column vectors, we denote with x T the transpose of x, with x T y the scalar product of x and y, and with | x | the Euclidean norm of x . Let L2 be the space of functions that are square integrable, and L1 be the space of functions whose absolute value is integrable. Let f : Rd → Ck , f ∈ L1 (Rd , Ck ) ∩ L2 (Rd , Ck ), the Fourier transform  f of f is deﬁned as

 f (ξ ) =

 1 d/2  T e−ιξ x f (x )dx, d 2π R

ξ ∈ Rd ,

(1)

and the inverse Fourier transform ˇf of f is

ˇf (x ) =

 1 d/2  T eιξ x f (ξ )dξ , d 2π R

x ∈ Rd .

(2)

We consider an initial boundary value problem for the Burgers equation on the plane:

  ∂ u(x, t ) − ν u(x, t ) + uT (x, t )∇ u(x, t ) = f (x, t ), x ∈ R2 , t > 0, ∂t

(3)

u(x, 0 ) = u0 (x ),

x ∈ R2 ,

(4)

t > 0,

(5)

lim u(x, t ) = 0,

|x|→∞

where ∇ is the gradient operator,  is the Laplace operator, ν > 0 is the viscosity, u0 : R2 → C2 is the known initial solution, f = ( f1 , f2 )T is the source term, and u = (u1 , u2 )T : R2 × R+ → C2 is the unknown function. We suppose that, for each t > 0, and for s = 1, 2,

u(·, t ),

∂u ∂u ∂u (·, t ), (·, t ), us (·, t ) ∈ L1 (R2 , C2 ) ∩ L2 (R2 , C2 ), ∂t ∂ xs ∂ xs

∂u ( x , t ) = 0. ∂ |x|→∞ xs lim

(6) (7)

Let  u(ξ , t ), ξ ∈ R2 , t > 0 be the Fourier transform of u(x, t ), x ∈ R2 , t > 0, with respect to variable x . Then, when s = 1, 2, ξ ∈ R2 , t > 0, we have

  T ∂u 1 ∂ (ξ , t ) = e−ιξ x u(x, t )dx = ∂t 2π R2 ∂t

∂ u ( ξ , t ), ∂t

(8)

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

  T ∂u 1 ∂ (ξ , t ) = e−ιξ x u(x, t )dx = ιξs u ( ξ , t ), ∂ xs 2π R2 ∂ xs

3

(9)

 2u ∂ ∂u ( ξ , t ) = ιξ (ξ , t ) = −ξs2 u ( ξ , t ), s ∂ xs ∂ x2s    ∂u 1 ∂u ι us (ξ , t ) = us (ξ , t )  (ξ , t ) = ys us (ξ − y, t ) u(y, t )dy, ∂ xs 2π ∂ xs 2π R2

(10) (11)

where  denotes the convolution product. Let f (·, t ) ∈ L1 (R2 , C2 ) ∩ L2 (R2 , C2 ) for every t > 0, we can consider the Fourier transform of both sides of equation (3), from which we obtain that  u must satisfy

∂ ι u (ξ , t ) = − ∂t 2π

 R2

yT  u(ξ − y, t ) u(y, t )dy − ν|ξ |2 u (ξ , t ) +  f ( ξ , t ),

(12)

where ξ ∈ R2 , t > 0,  f (ξ , t ) is the Fourier transform of f( x, t) with respect to x . From condition (4), and by applying the Duhamel principle [40] to formula (12), we have that  u must satisfy the following integral equation:

 t ι 2 2  u(ξ , t ) = e−ν|ξ | t u0 ( ξ ) − e−ν|ξ | (t−τ ) 2π 0   · yT  u(ξ − y, τ ) u(y, τ )dy + 2π ι f (ξ , τ ) d τ , R2

ξ ∈ R2 , t > 0,

(13)

where  u0 is the Fourier transform of u 0 . Integral equation (13) is a useful tool in the numerical solution of problem (3)–(5). We propose a simple approach for the computation of the approximate solution of this integral equation. Let  u(·, t ) be the solution of (13) at a generic time t ≥ 0, let δ t > 0 be a time step, then

2  u(ξ , t + δt ) = e−ν|ξ | δt  u (ξ , t ) −

·

 R2

ι 2π



t+δt

t

2 e−ν|ξ | (t−τ )





yT  u(ξ − y, τ ) u(y, τ )dy + 2π ι f (ξ , τ ) d τ ,

ξ ∈ R2 .

(14)

The time integral in formula (14) can be approximated by a usual quadrature rule producing a formula similar to multistep equations for differential equations; for example, the use of the rectangular rule produces an explicit formula resembling the Euler method, in fact we have:

  2   u(ξ , t + δt ) = e−ν|ξ | δt  u(ξ , t ) + δt i(ξ , t ) +  f (ξ , t ) , where

i (ξ , t ) = −

ι 2π

 R2

yT  u(ξ − y, t ) u(y, t )dy,

ξ ∈ R2 .

ξ ∈ R2 ,

(15)

(16)

An approximate solution of integral equation (13) can be actually computed by formula (15) when this computation can be restricted to a ﬁnite domain R ⊂ R2 . This is usually achieved from the vanishing property of  u(·, t ) at inﬁnity. In this way, we can obtain a practical scheme for the numerical approximation U ( ·, t) of the solution u ( ·, t) of (13). We note that, formulas (15) and (16) require the approximation of i, that can be expressed in terms of four convolution integrals of complex valued functions. So the numerical approximation of these integrals can be computed eﬃciently by a quadrature rule with equidistant nodes and by the FFT algorithm, see [41] for a detailed discussion on computational and stability properties of FFT algorithms. Thus, formulas (15) and (16) can be used as a marching-time scheme in the numerical solution (·, t ) from U ( ·, t), and then for the computation of of problem (3)–(5), when an eﬃcient algorithm for the computation of U  U (·, t + δt ) from U (·, t + δt ) can be determined. Of course the same method can be used for both the computations, and it can be based on a discrete Fourier transform, that in turn can be performed at FFT speed. For simplicity, we brieﬂy describe this method for a function of one variable; its generalization to functions of two real variables is immediate and can be omitted. Let f ∈ L1 (R ) ∩ L2 (R ) be a generic function that is zero or negligible outside the interval (−2π M, 2π M ), where M > 0 is a generic real number. So, the Fourier transform of f can be approximated as follows

1  f (ξ ) ≈ √ 2π

 2π M −2π M

f (x )e−ιξ x dx,

ξ ∈ R,

(17)

4

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

and by the midpoint quadrature rule with nodes x j = −2π M + ( j + 12 )2π h, j = 0, . . . , N − 1, N ≥ 1, and h = ing approximation is obtained:

2M N ,

the follow-

1 1  f (ξ ) ≈ √ 2π h f (x j )e−ι(−M+( j+ 2 )h )2πξ 2π j=0 N−1

=

N−1

√ h 2π 2π heι(M− 2 )2πξ f ( x j ) e −ι N

j2Mξ

,

ξ ∈ R.

(18)

j=0

Finally, by considering the evaluation nodes ξi = − 4NM +

 f ( ξi ) ≈

1 2M ( i

+ 12 ), i = 0, 1, . . . , N − 1, we have

N−1

√ h N 1 1 2π N 1 2π 2π heι2π (M− 2 )(− 4M + 2M (i+ 2 )) f (x j )e−ι N j (− 2 + 2 ) e−ι N ji .

(19)

j=0

Note that the sum in formula (19) is the discrete Fourier transform of the vector having components 2π

f (x j )e−ι N j (− 2 + 2 ) , j = 0, 1, . . . , N − 1, and so it can be computed eﬃciently by the FFT algorithm. In summary, we have the following algorithm for the computation of a generic time step in the solution of problem (3)–(5), where the source term f and the solution u = (u1 , u2 )T are zero or negligible outside the rectangle R = [−a1 , a1 ] × [−a2 , a2 ] ⊂ R2 . Let N1 and N2 be the number of discretization points, along x1 , x2 axis, respectively, let N

1



xs,i = −as + i +



1 2as , i = 0, 1, . . . Ns − 1, 2 Ns

(20)

be the discretization points along xs , s = 1, 2. We denote with U s,0 ∈ RN1 N2 , s = 1, 2, the vector whose entries are us (x1, i , x2, j , 0), s = 1, 2, i = 0, 1, . . . N1 − 1, j = 0, 1, . . . N2 − 1, and with U s,k ∈ RN1 N2 , s = 1, 2, k = 1, 2, . . . , the vector whose entries are approximations of us (x1, i , x2, j , kδ t ), s = 1, 2, i = 0, 1, . . . N1 − 1, j = 0, 1, . . . N2 − 1, where δ t > 0 is the time step. s, k , s = 1, 2, k ≥ 0, f, δ t > 0, N1 and N2 as above. Compute U s,k+1 , s = 1, 2, as follows:  as the discrete Fourier transform of U for s = 1, 2, compute U s, k ; s,k   for s = 1, 2, compute U s,k+1 from U s,k by formulas (15) and (16);  for s = 1, 2, compute U s,k+1 as the inverse discrete Fourier transform of U s,k+1 .

Algorithm 1. Given U

Note that the ﬁrst and the third step of Algorithm 1 requires the computation of four discrete Fourier transforms of vectors having N1 N2 entries. This computation can be performed by the two dimensional version of (19), and needs O (4N1 N2 log2 (N1 N2 ) ) arithmetic operations, see [42] for a discussion of FFT complexity. The second step requires the computation of function (16) that contains four convolution integrals, each one of them requires three FFT of vectors having 4N1 N2 entries, so its computation cost is O (48N1 N2 log2 (4N1 N2 ) ). 3. Galerkin method with ﬁnite element basis Let ⊂ R2 be a bounded domain with boundary ∂ . Let u be the solution of problem (3)–(5), and let g(x, t ) = (g1 (x, t ), g2 (x, t ))T = u(x, t ), x ∈ ∂ , t > 0. Then u is also the solution of the following initial-boundary value problem:





T ∂ t u(x, t ) − ν u(x, t ) + u (x, t )∇ u(x, t ) = f (x, t ),

u(x, 0 ) = u0 (x ), u ( x, t ) = g ( x , t ) ,

x ∈ , t > 0 x ∈ , x ∈ ∂ , t > 0 .

(21)

This problem can be solved by standard approximation methods such as ﬁnite difference methods, ﬁnite elements methods. In the following, we describe the Galerkin method with a ﬁnite elements basis. In particular, let T be a triangulation of having NT triangles τ m , m = 1, 2, . . . NT and P nodes p n , n = 1, 2, . . . P; a node is either a vertex of a triangle or a middle point of an edge of a triangle. Let PI be the number of the nodes inside , that are supposed to be at the beginning of the node list, without loss of generality, so that pP +n , n = 1, 2, . . . , PB ≡ P − PI are nodes on the boundary ∂ . Let φn : → R, I

n = 1, 2, . . . , P, be the piecewise quadratic continuous function such that, for each triangle τ of T , φ n |τ ≡ 0 if pn ∈ τ , and φ n |τ is a quadratic polynomial if pn ∈ τ , with φn ( pn ) = 1 and φn ( p) = 0 when p is a node of τ different from pn . See [43] for the construction of φ n , n = 1, 2, . . . , P . Let w(x, t ) = (w1 (x, t ), w2 (x, t ))T with

w s ( x, t ) =

PI P

ws(n ) (t )φn (x ) + gs ( pn , t )φn (x ), n=1

x ∈ , t > 0 ,

(22)

n=PI +1

be the approximation of u, where, for n = 1, 2, . . . , PI , ws(n ) (t ), s = 1, 2, are unknown coeﬃcients that must be computed. To this aim, the Galerkin method requires that, for each t > 0, and for n = 1, 2, . . . PI , the following equation holds:



 T  ∂ w(x, t ) − ν w(x, t ) + w (x, t )∇ w(x, t ) − f (x, t ) φn (x )dx = 0, ∂t

(23)

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

5

so that, from Green’s formula and from the fact that φn (x ) = 0 when x ∈∂ , the following equations are obtained:



∂ T w (x, t ) + (w (x, t )∇ ws (x, t )) − fs (x, t ) φn (x )dx ∂t s



∇ T ws (x, t )∇φn (x )dx = 0,

s = 1, 2.

(24)

By approximating the derivative term in (24) with the forward difference quotient

ws (x, t + δt ) − ws (x, t ) ∂ w (x, t + δt ) ≈ , ∂t s δt

(25)

where δ t > 0 is the discretization step for the time variable, and by linearizing the non-linear term in (24) in this way

  (wT ∇ ws )(x, t + δt ) ≈ wT (x, t )∇ ws (x, t + δt ),

we have the following linear system:



1

δt



+ +ν

ws (x, t + δt )φn (x )dx −







1

δt



(26)

w s ( x, t ) φn ( x ) d x

wT (x, t )∇ ws (x, t + δt ) − fs (x, t + δt )

 φn ( x ) d x

∇ T ws (x, t + δt )∇φn (x )dx = 0, n = 1, 2, . . . , PI , s = 1, 2,

(27)

where the unknowns are the coeﬃcients ws(n ) (t + δt ) ∈ R2 , n = 1, 2, . . . , PI , of ws (·, t + δt ), s = 1, 2. For k = 0, 1, . . . and s = (P )

(1 ) (2 ) 1, 2, we denote with tk = kδt the discretization times, with ws,k = (ws,k , ws,k , . . . , ws,kI )T ∈ RPI the vectors of representation (n )

(P )

(n )

1 ) (2 ) coeﬃcients of the approximate solution ws (·, tk ) i.e. ws,k = ws (tk ), n = 1, 2, . . . , PI , and with gs,k = (g(s,k , gs,k , . . . , gs,kB )T ∈ n) RPB where g(s,k = gs ( pP +n , tk ), n = 1, 2, . . . , PB . Then ws,k ∈ RPI is the solution of the following linear system:

1

δt

I



MI + NI (tk−1 ) + ν SI ws,k =

1

δt



M wTs,k−1 , gTs,k−1

T

1

δt



MB + NB (tk−1 ) + ν SB gs,k + bs,k ,

(28)

where vectors bs,k ∈ RPI have components



f s ( x , t k ) φn ( x ) d x ,

n = 1, 2, . . . , PI ,

(29)

the entries of matrices M, S, N (tk ) ∈ RPI ×P are



Mn,m =

 Sn,m =

N (tk )n,m =

φm (x )φn (x )dx,

(30)

∇ T φm (x )∇φn (x )dx, PI

l=1 s=1,2

(l ) ws,k



φl ( x )

(31)

PB

∂φm (x ) l) φn ( x ) d x + g(s,k ∂ xs l=1 s=1,2



φPI +l (x )

∂φm (x ) φn ( x ) d x . ∂ xs

(32)

where n = 1, 2, . . . , PI and m = 1, 2, . . . , P, and these matrices are partitioned according to the indices of nodes inside and on ∂ , that is M = (MI , MB ), S = (SI , SB ), N (tk ) = (NI (tk ), NB (tk )), MI , SI , NI (tk ) ∈ RPI ×PI , MB , SB , NB (tk ) ∈ RPI ×PB For each time step k, we have to solve linear system (28), which has 2PI equations and 2PI unknowns, i.e. the components ws(n ) (tk ), n = 1, 2, . . . , PI of ws,k , s = 1, 2. From the solution of system (28), we have that formula (22) deﬁnes an approximation w(x, tk ) = (w1 (x, tk ), w2 (x, tk ))T of u ( x, tk ) solution of problem (21). 4. Numerical method Let u = (u1 , u2 )T be the solution of problem (3)–(5), let T > 0 be the time limit for the numerical solution of this problem, let R = [−a1 , a1 ] × [−a2 , a2 ] ⊂ R2 be the rectangle containing the relevant part of u and f, so suppose that u (·, t) and f (·, t), t ∈ [0, T], are zero or negligible outside R. We propose a numerical algorithm to solve problem (3)–(5) that, at each time step, uses the FFT algorithm to solve integral equation (13); this numerical solution is then used to compute function g on the artiﬁcial boundary ∂ of the bounded domain ⊂ R2 in problem (21), which is solved by the Galerkin method described in Section 3. Note that, even if u 0 and f (·, t), t ∈ [0, T], have the same compact support, the solution u (·, t) usually has support spreading out more and more as t increases. So, every method for the numerical solution of (3)–(5) can provide an accurate result only for a bounded time interval, that depends on the size of the ﬁnite spatial domain where the solution is

6

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

actually computed. In the proposed method the eﬃciency of the integral method is exploited to solve the Burgers equation on a large interval R ⊂ R2 , this rough approximation is used to produce the boundary condition for the problem (21), that is deﬁned on a smaller domain where a more accurate method is used. Thus, at each time step tk = kδt , k = 0, 1, . . . , K, δ t > 0, K δt = T , we consider the approximation of u ( ·, tk ) computed by the integral equation method on interval R, and the one computed by the Galerkin method on domain , where ⊂ R. For k = 0, 1, . . . , K, let U s, k , s = 1, 2, be the approximation computed by the integral equation method with notation deﬁned in Section 2, and let w(·, tk ) ∈ R be the approximation computed by the Galerkin method with notation deﬁned in Section 3. In the proposed algorithm we have to exchange information between the FFT method and the Galerkin method. Formula (22) is used to transfer on the FFT grid the information computed by Galerkin method. Instead a bilinear interpolation formula is used to transfer on the triangulation nodes the information computed by the FFT. More precisely, given a point x = (x1 , x2 )T ∈ R2 in the rectangle R ⊂ R having vertices q1 = (x1,i , x2, j )T , q2 = (x1,i+1 , x2, j )T , q3 = (x1,i+1 , x2, j+1 )T and q4 =

(x1,i , x2, j+1 )T , where 0 ≤ i ≤ N1 − 2, 0 ≤ j ≤ N2 − 2, and given v = (v1 , v2 , v3 , v4 )T the values of a function v : R2 → R at these four vertices, the following approximation v˜ of v is computed

v˜ (x ) = v1 R2, j R1,i + v2 R2, j L1,i + v3 L2, j L1,i + v4 L2, j R1,i , x ∈ R,

(33)

where for s = 1, 2,

Ls,k =

xs − xs,k xs − xs,k+1 ; Rs,k = . xs,k+1 − xs,k xs,k − xs,k+1

(34)

The proposed algorithm is given below. Algorithm 2. Given problem (3)–(5), and given: δ t > 0 the size of the time step; K ∈ N the number of time steps; as > 0 s = 1, 2, that deﬁne the interval R for the integral equation method; Ns > 0 s = 1, 2, the number of discretization points inside (P ) (1 ) (2 ) R along axis xs ; ⊂ R the domain for the Galerkin method; T a triangulation of . Compute ws,k = (ws,k , ws,k , . . . , ws,kI )T ∈ RPI , s = 1, 2, k = 1, 2, . . . K, as follows: for s = 1, 2, compute U s,0 ∈ RN1 N2 by the evaluation of us ( ·, 0) at the grid nodes given by (20); for s = 1, 2, compute ws,0 = (us ( p1 , 0 ), us ( p2 , 0 ), . . . , us ( pP , 0 ))T ∈ RPI ; I

for s = 1, 2, compute



gs,0 = us ( pP +1 , 0 ), us ( pP +2 , 0 ), . . . , us ( pP +P , 0 ) I

I

I

T

∈ RPB ;

B

for k = 1, 2, . . . , K, do – set tk = kδt ; – for s = 1, 2, compute U s,k ∈ RN1 N2 by Algorithm 1;

n) – for s = 1, 2, compute the approximation g(s,k of gs ( pP +n , tk ), n = 1, . . . PB , by evaluating at pP +n the bilinear apI

I

proximation formula (33) obtained on the rectangle R containing pP +n , with function values v given by the entries I

of U

s, k

that correspond to the vertices of R; (P )

1 ) (2 ) – set gs,k = (g(s,k , gs,k , . . . , gs,kB )T ∈ RPB , s = 1, 2;

RPI ,

– compute ws,k ∈ s = 1, 2, solution of linear system (28); – reﬁne, by formula (22), the entries of U s, k , s = 1, 2, that correspond to grid points xi j = (x1,i , x2, j )T given by (20), when x ij ∈ . Note that, for each time iteration k = 1, 2, . . . , K, Algorithm 2 ﬁrst computes the approximated solution of problem (3)– (5), restricted to interval R, by the integral equation method. This solution is used to deﬁne the boundary values for problem (21) on domain ⊂ R, and this problem is solved by the Galerkin method. Finally, this solution is used to reﬁne the approximated solution computed by the integral equation method; this reﬁnement involves only the values corresponding to the discretization points of R that fall inside . Other algorithms can be considered by using the two methods described in the previous sections. The following ones will be used in the numerical experiments described in the next section in order to test the new proposed procedure (Algorithm 2), where the artiﬁcial boundary condition on the domain is computed by using the above mentioned FFT solution method. Algorithm 3. Given problem (3)–(5), and given: δ t > 0 the size of the time step; K ∈ N the number of time steps; as > 0 s = 1, 2, that deﬁne the interval R for the integral equation method; Ns > 0 s = 1, 2, the number of discretization points in R along axis xs . Compute U s,k ∈ RN1 N2 approximation of us ( ·, tk ) at grid nodes given by (20), where s = 1, 2, k = 1, 2, . . . K, tk = kδt , and u is the solution of (3)–(5), as follows: for s = 1, 2, compute U s,0 ∈ RN1 N2 by the evaluation of us ( ·, 0) at the grid nodes given by (20); for k = 1, 2, . . . , K, do

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

7

– set tk = kδt ; – for s = 1, 2, compute U s,k ∈ RN1 N2 by Algorithm 1. Note that, for each time iteration k = 1, 2, . . . , K, Algorithm 3 computes the approximated solution of problem (3)–(5), restricted to interval R, by the integral equation method. Of course, the bilinear formula (33) can be used to evaluate this approximate solution at a point different from the grid nodes given by (20). Algorithm 4. Given problem (21), and given: δ t > 0 the size of the time step; K ∈ N the number of time steps; ⊂ R2 the (P ) (1 ) (2 ) domain for the Galerkin method; T a triangulation of . Compute ws,k = (ws,k , ws,k , . . . , ws,kI )T ∈ RPI , s = 1, 2, k = 1, 2, . . . K, as follows: for s = 1, 2, compute ws,0 = (us ( p1 , 0 ), us ( p2 , 0 ), . . . , us ( pP , 0 ))T ∈ RPI ; I

for s = 1, 2, compute

gs,0 = (us ( pP +1 , 0 ), us ( pP +2 , 0 ), . . . , us ( pP +P , 0 )) ∈ RPB ; I

I

I

B

for k = 1, 2, . . . , K, do – set tk = kδt ; n) – for s = 1, 2, compute g(s,k = gs ( pP +n , tk ), n = 1, . . . PB ; I

(P )

1 ) (2 ) – set gs,k = (g(s,k , gs,k , . . . , gs,kB ) ∈ RPB , s = 1, 2;

– compute ws,k ∈ RPI , s = 1, 2 solution of linear system (28). Note that, for each time iteration k, Algorithm 4 computes, by the Galerkin method, the approximated solution of problem (21), where the boundary values g are given by a priori information. In particular, in the following experiments we will choose either the exact boundary condition g = u|∂ or the null boundary condition g = 0; this last choice is reasonable when the solution u is negligible outside the domain . 5. The numerical experiments The algorithms described in the previous section are tested by using three examples of problem (3)–(5), where the solution is known. This numerical experiment is based on the comparison of the accuracy in the computed solutions and of the computational time necessary to obtain such solutions. The following examples are considered in the experiment. Example 1. Let α > 0 be a given constant, let

u0 ( x ) =

x |x|2

να + α 2 e 4να

, x ∈ R2 ,

(35)

f ( x, t ) = 0 , x ∈ R 2 , t > 0 .

(36)

The exact solution of problem (3)–(5) is

u ( x, t ) =

x

ν (t + α ) + (t + α )

|x|2

, x ∈ R2 , t > 0 .

(37)

2 e 4ν (α +t )

In the numerical results we used α = 0.4, ν = 1 and t ∈ [0, 0.1]. Example 2. Let β > 0 be a given constant, f (x, t ) = 0, x ∈ R2 , t > 0, and let u0 = (u0,1 , u0,2 )T , where



u0,s (x ) =

h ( xr , β ) xs h ( xs , β ) − 4 e

h(y, z ) = e

−β

sin(xs ) √



π

β ( h ( xs , β )h ( xr , β ) + 1 ) 2

− y4z

2 z + ιy Erf √ 2 z

s, r ∈ {1, 2}, s = r,

,

2 z − ιy + Erf √ 2 z

(38)

, y, z ∈ R, z > 0,

(39)

and Erf is the Error function [44], page 297. The exact solution of problem (3)–(5) with ν = 1 is

u s ( x, t ) =

xs h ( xs , t + β ) − 4

e−(t+β ) sin(xs ) √

π

h ( xr , t + β ) · , s, r ∈ {1, 2}, s = r. (h(xs , t + β )h(xr , t + β ) + 1 )(t + β ) In the numerical results we used β = 10 and t ∈ [0, 0.1].

(40)

8

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

a

b

c

Fig. 1. (a) The triangulation T1 of , which has NT = 276 triangles, PI = 517 internal nodes and PB = 72 boundary nodes. The edges length ranges between 0.12r and 0.20r . (b) The triangulation T2 of , which has NT = 582 triangles, PI = 1111 internal nodes and PB = 108 boundary nodes. The edges length ranges between 0.08r and 0.15r .(c) Computational domains used in the numerical experiments. Algorithm 2 is applied on the small disk and on the square R with as = a, s = 1, 2, that is used to compute the artiﬁcial boundary condition; Algorithm 3 is applied on the square R with as = a, s = 1, 2; Algorithm 4 is applied on the small disk , with r = r1 , by using exact boundary condition or on the large disk , with r = r2 , by using the zero boundary condition.

Example 3. Let γ , ζ , c > 0 be given constants, and let v : R → R2 be a given function. For x ∈ R2 and t > 0, let h(x, t ) = ||x||2 − ct − γ ,

f ( x, t ) =

⎧ ζ  h (x,t ) v (t ) + ⎪ ⎨e ·e ⎪ ⎩

ζ h (x,t )

ζ

(c + 4ν − 2xT v(t )·  x||2 − 4hν|| 2 (x,t ) (ζ + 2h (x, t ) ) )v (t ) , h2 (x,t )

0,

and

 u0 ( x ) =

ζ

||x||2 − ct − γ < 0, ||x||2 − ct − γ ≥ 0,

||x||2 − γ < 0, ||x||2 − γ ≥ 0.

e h(x,0) v(0 ), 0,

The exact solution of problem (3)–(5) is



u ( x, t ) =

ζ

e h(x,t ) v(t ), 0,

||x||2 − ct − γ < 0, ||x||2 − ct − γ ≥ 0.

(41)

In the numerical results we used ν = γ = ζ = 1, c = 0.1, v(t ) = 0.1(t + 1, 2t + 1 )T and t ∈ [0, 0.01]. The numerical results of these examples are computed by Algorithms 2–4, where the integral equation method and the Galerkin method are applied with different values for the corresponding approximation parameters. In the integral equation method the interval R = [−a1 , a1 ] × [−a2 , a2 ] is chosen in order to contain the set where the solution of the three examples is practically different from zero; parameters as , s = 1, 2, depend on the example considered and are reported in Table 2–4; the number of discretization points is Ns = 30, 60, 120, s = 1, 2. In the Galerkin method the domain is a circular disk with center the axis origin and radius r which depends on the example considered and on the boundary condition chosen in the Galerkin method; the value of r is reported in Table 2–4. For each circular domain, two triangulations T1 and T2 are considered. In particular, Algorithm 2 is applied always with triangulation T1 , instead Algorithm 4 is applied with both triangulations T1 and T2 . Fig. 1(a) shows the triangulation T1 , which has NT = 276 triangles, PI = 517 internal nodes, PB = 72 boundary nodes, and the edge length ranges between 0.12r and 0.20r . Fig. 1(b) shows the triangulation T2 , which has NT = 582 triangles, PI = 1111 internal nodes, PB = 108 boundary nodes, and the edge length ranges between 0.08r and 0.15r . The comparison of the algorithms is given in terms of the elapsed time in the computation, and the accuracy of the numerical solution, that is described by means of the following two error indices

NE

E (t ) =

EC (t ) =

n=1



 |U1 (ξ n , t ) − u1 (ξ n , t )| + |U2 (ξ n , t ) − u2 (ξ n , t )|  , NE  n=1 |u1 (ξ n , t )| + |u2 (ξ n , t )|

maxn=1,...,NE

 |U1 (ξ n , t ) − u1 (ξ n , t )|, |U2 (ξ n , t ) − u2 (ξ n , t )|   , maxn=1,...,NE |u1 (ξ n , t )|, |u2 (ξ n , t )|

(42)



(43)

where U and u are the approximated solution and the exact solution of problem (3)–(5), respectively, and ξ n , n = 1, 2, . . . , NE , are uniformly distributed points in a circular disk centered at the origin and having radius rE ; in particular, these points are chosen as vertices of a grid of squares.

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

9

Table 1 The L1 -norm of u, computed over 100 points uniformly distributed on ∂ at time t for the three considered examples. ∂ is a circumference centered at the axes origin with radius r . Example 1

Example 2

Example 3

r t

2.6

3.0

0.75

1.14

δt

2.9(−1) 3.4(−1) 3.9(−1)

2.0(−2 ) 2.1(−2 ) 2.1(−2 )

0 0 0

5δ t 10δ t

7.1 δt = 0.01 2.2(−12 ) 2.8(−11 ) 3.8(−10 )

8.1 δt = 0.01

9.2(−2 ) 9.2(−2 ) 9.1(−2 )

7.7(−2 ) 7.4(−2 ) 7.3(−2 )

δt = 0.001

Table 2 The numerical results obtained by Algorithms 2–4, for Example 1. It is shown the number of time iterations K, the time t, the errors E(t) and EC (t), and the computational time TC . The errors E(t) and EC (t) have been computed by (42) and (43), respectively, with N = 7709 uniformly distributed points inside a circular disk centered at the origin and having radius rE = 2.6. The results obtained with Algorithm 2 are relative to as = 5.9, s = 1, 2, r = 2.6, and the triangulation of is T1 . The results obtained with Algorithm 3 are relative to as = 5.9, s = 1, 2. The results obtained with Algorithm 4 are relative to different boundary conditions g and/or different triangulations and/or different computational domains. Algorithm 2 K 1 5 10

t 0.01 0.05 0.10

K 1 5 10

t 0.01 0.05 0.10

Algorithm 3 K 1 5 10

t 0.01 0.05 0.10

K 1 5 10

t 0.01 0.05 0.10

N1 = N2 = 30 E(t) 1.9(−3) 3.1(−3) 4.0(−3)

EC (t) 1.1(−2 ) 1.3(−2 ) 1.4(−2 )

N1 = N2 = 60 E(t) 1.3(−3) 2.1(−3) 3.1(−3)

TC 0.20 1.36 2.65 N1 = N2 = 30 E(t) 2.6(−2) 3.0(−2) 3.6(−2)

EC (t) 4.5(−2 ) 4.4(−2 ) 4.9(−2 )

EC (t) 3.0(−3 ) 4.2(−3 ) 5.5(−3 ) TC 0.39 2.40 4.93

N1 = N2 = 60 E(t) 7.9(−3) 1.7(−2) 2.7(−2)

EC (t) 1.2(−2 ) 2.4(−2 ) 4.0(−2 )

N1 = N2 = 120 E(t) 1.2(−3) 2.2(−3) 3.5(−3)

EC (t) 2.5(−3 ) 3.0(−3 ) 5.9(−3 )

TC 2.37 12.33 24.99 N1 = N2 = 120 E(t) 4.0(−3) 1.5(−2) 2.5(−2)

EC (t) 5.7(−3 ) 2.2(−2 ) 3.9(−2 )

TC 0.03 0.17 0.29

TC 0.21 1.06 2.12

TC 2.19 10.95 21.79

Algorithm 4

r = 2 . 6 , T 1 , g = u| ∂

r = 7.1, T 1 , g = 0

r = 7.1, T 2 , g = 0

K 1 5 10

t 0.01 0.05 0.10

E(t) 1.1(−3) 1.9(−3) 2.6(−3)

K 1 5 10

t 0.01 0.05 0.10

EC (t) 2.4(−3 ) 2.9(−3 ) 4.2(−3 ) TC 0.17 0.98 1.67

E(t) 2.1(−2) 1.9(−2) 1.8(−2)

EC (t) 3.6(−2 ) 3.6(−2 ) 3.6(−2 ) TC 0.15 0.80 1.60

E(t) 8.0(−3) 6.8(−3) 6.3(−3)

EC (t) 1.6(−2 ) 1.3(−2 ) 1.4(−2 ) TC 0.56 2.93 5.26

The numerical results are relative to K = 10 time iterations, with step size δt = 0.01 for Examples 1 and 2, and δt = 0.001 for Example 3. Fig. 1 (c) shows the computational domains R and chosen for the numerical results. In particular Algorithm 2 is applied on the square R with ai = a, i = 1, 2, that contains a disk with r = r1 , where the artiﬁcial boundary conditions are evaluated. Algorithm 3 is applied on the square R with ai = a, i = 1, 2. Algorithm 4 is applied on the disk with r = r1 by using the exact boundary conditions g = u|∂ , and on the disk with r = r2 by using the zero boundary conditions g = 0. Note that Table 1 shows the magnitude, at three different times, of the exact solution u on the boundary ∂ of the domains , used for each examples. These values are given by the L1 -norm of u evaluated by 100 points that are uniformly distributed on ∂ . The numerical results obtained by Algorithms 2–4 are reported in Table 2–4. These tables show the errors (42) and (43) in the numerical solution and the elapsed time for the Examples 1–3. From Table 2–4, we can observe that Algorithm 2 has a higher computational time than Algorithm 3 when the same Ni , i = 1, 2, are used. Moreover

10

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14 Table 3 The numerical results obtained by Algorithms 2–4, for Example 2. It is shown the number of time iterations K, the time t, the errors E(t) and EC (t), and the computational time TC . The errors E(t) and EC (t) have been computed by (42) and (43), respectively, with N = 7709 uniformly distributed points inside a circular disk centered at the origin and having radius rE = 3.0. The results obtained with Algorithm 2 are relative to as = 6.7, s = 1, 2, r = 3.0, and the triangulation of is T1 . The results obtained with Algorithm 3 are relative to as = 6.7, s = 1, 2. The results obtained with Algorithm 4 are relative to different boundary conditions g and/or different triangulations and/or different computational domains. Algorithm 2 K 1 5 10

t 0.01 0.05 0.10

K 1 5 10

t 0.01 0.05 0.10

Algorithm 3 K 1 5 10

t 0.01 0.05 0.10

K 1 5 10

t 0.01 0.05 0.10

N1 = N2 = 30 E(t) 2.8(−4) 6.9(−4) 1.1(−3)

EC (t) 3.0(−3 ) 3.9(−3 ) 4.5(−3 )

N1 = N2 = 60 E(t) 8.9(−5) 1.8(−4) 3.0(−4)

TC 0.22 1.43 3.07 N1 = N2 = 30 E(t) 2.4(−3) 2.3(−3) 2.4(−3)

EC (t) 3.4(−3 ) 3.3(−3 ) 3.2(−3 )

TC 0.39 2.53 5.06 N1 = N2 = 60 E(t) 5.8(−4) 6.9(−4) 1.0(−3)

TC 0.04 0.17 0.33

Algorithm 4

r = 3 . 0 , T 1 , g = u| ∂

K 1 5 10

t 0.01 0.05 0.10

E(t) 3.4(−5) 3.4(−5) 3.6(−5)

K 1 5 10

t 0.01 0.05 0.10

EC (t) 6.7(−5 ) 6.7(−5 ) 7.1(−5 ) TC 0.17 0.88 1.75

EC (t) 7.7(−4 ) 1.1(−3 ) 1.3(−3 )

EC (t) 8.4(−4 ) 8.5(−4 ) 1.2(−3 ) TC 0.22 1.09 2.19

r = 8.1, T 1 , g = 0 E(t) 7.3(−4) 6.7(−4) 6.9(−4)

EC (t) 1.2(−3 ) 1.2(−3 ) 1.2(−3 ) TC 0.18 0.90 1.72

N1 = N2 = 120 E(t) 4.8(−5) 9.4(−5) 1.8(−4)

EC (t) 2.0(−4 ) 4.5(−4 ) 8.0(−4 )

TC 2.43 12.73 25.38 N1 = N2 = 120 E(t) 1.6(−4) 4.5(−4) 8.8(−4)

EC (t) 2.0(−4 ) 6.2(−4 ) 1.2(−3 )

TC 2.19 10.84 21.79 r = 8.1, T 2 , g = 0 E(t) 2.5(−4) 2.4(−4) 2.4(−4)

EC (t) 5.6(−4 ) 5.5(−4 ) 5.4(−4 ) TC 0.60 2.95 5.80

Algorithm 2 has also a higher computational time than Algorithm 4 when the same triangulation T1 is used. Note that this is a natural consequence of the structure of Algorithm 2. From these tables we can also observe that Algorithm 4 provides the most accurate results, with respect to the other algorithms, when the exact boundary condition is used. On the contrary Algorithm 3 seems to be the most coarse-grained one among the three considered here. However, when the exact boundary condition is substituted with the zero boundary condition, on a larger domain, Algorithm 4 is no longer the most accurate one and/or its computational cost becomes higher than the ones of the other two algorithms. Figs. 2–4 show the error E obtained for Examples 1–3, respectively, by Algorithm 2 with N1 = N2 = 60 and triangulation T1 on , by Algorithm 3 with N1 = N2 = 60, by Algorithm 4 with triangulation T1 on and with exact boundary condition, and by Algorithm 4 with triangulation T1 and T2 , on a larger domain and with the zero boundary condition. Note that the graphics of the error EC , as it is shown in Table 2–4, have a similar beaviour of the ones shown in Figs. 2–4, so they are not reported for brevity. Fig. 5 reports the mean computational time of Algorithm 2 with N1 = N2 = 30, 60, 120, and the mean computational time of Algorithm 4 with triangulation T2 ; the average is calculated from the computational times TC in the proposed examples. From these ﬁgures we can observe that when the zero boundary condition is used, Algorithm 4, that is the Galerkin method on the big disc of Fig. 1(c), does not provide the more accurate results. More precisely, for Examples 1 and 2, Algorithm 2 with N1 = N2 = 60 is more accurate then Algorithm 4 with the zero boundary condition; instead for Example 3, Algorithm 2 with N1 = N2 = 60 has an accuracy similar to the one of Algorithm 4 with the zero boundary condition and triangulation T2 but Algorithm 2 has a lower computational time. Therefore, Algorithm 2, with N1 = N2 = 60, has accuracy comparable with the one of Algorithm 4 with triangulation T2 , but with a reduced computational time. Note that the artiﬁcial boundary conditions allow the use of a small domain for the application of the Galerkin method. Moreover, it can be observed that proper discretization grids for the integral equation method give optimal results, in fact the case N1 = N2 = 120 requires a higher computational time that the case N1 = N2 = 60, but it does not provide relevant improvements in the numerical solution.

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

11

Table 4 The numerical results obtained by Algorithms 2–4, for Example 3. It is shown the number of time iterations K, the time t, the errors E(t) and EC (t), and the computational time TC . The errors E(t) and EC (t) have been computed by (42) and (43), respectively, with N = 7709 uniformly distributed points inside a circular disk centered at the origin and having radius rE = .75. The results obtained with Algorithm 2 are relative to as = .94, s = 1, 2, r = .75, and the triangulation of is T1 . The results obtained with Algorithm 3 are relative to as = .94, s = 1, 2. The results obtained with Algorithm 4 are relative to different boundary conditions g and/or different triangulations and/or different computational domains. Algorithm 2 K 1 5 10

t 0.001 0.005 0.010

K 1 5 10

t 0.001 0.005 0.010

Algorithm 3 K 1 5 10

t 0.001 0.005 0.010

K 1 5 10

t 0.001 0.005 0.010

N1 = N2 = 30 E(t) 3.1(−4) 9.8(−4) 1.5(−3)

EC (t) 2.5(−3 ) 4.3(−3 ) 4.7(−3 )

N1 = N2 = 60 E(t) 1.5(−4) 5.4(−4) 7.4(−4)

TC 0.24 1.68 3.36 N1 = N2 = 30 E(t) 3.0(−3) 3.6(−3) 4.1(−3)

EC (t) 3.6(−3 ) 4.9(−3 ) 5.7(−3 )

EC (t) 1.1(−3 ) 2.2(−3 ) 2.2(−3 ) TC 0.41 2.67 5.50

N1 = N2 = 60 E(t) 8.8(−4) 1.5(−3) 2.1(−2)

EC (t) 1.2(−3 ) 2.7(−3 ) 3.4(−3 )

N1 = N2 = 120 E(t) 1.2(−4) 4.2(−4) 5.3(−4)

EC (t) 7.9(−4 ) 1.7(−3 ) 1.6(−3 )

TC 2.55 13.16 26.45 N1 = N2 = 120 E(t) 3.5(−4) 1.0(−3) 1.6(−3)

EC (t) 8.2(−4 ) 2.3(−3 ) 2.9(−3 )

TC 0.04 0.17 0.33

TC 0.21 1.10 2.20

TC 2.20 11.00 22.12

Algorithm 4

r = .75, T1 , g = u|∂

r = 1.14, T1 , g = 0

r = 1.14, T2 , g = 0

K 1 5 10

t 0.001 0.005 0.010

E(t) 1.1(−4) 1.2(−4) 1.3(−4)

E(t) 3.4(−4) 9.0(−4) 1.7(−3)

K 1 5 10

t 0.001 0.005 0.010

EC (t) 4.4(−4 ) 4.3(−4 ) 4.3(−4 ) TC 0.19 0.89 1.93

EC (t) 1.6(−3 ) 2.9(−3 ) 4.5(−3 ) TC 0.18 0.86 1.74

E(t) 1.7(−4) 1.8(−4) 3.2(−4)

EC (t) 8.1(−4 ) 9.8(−4 ) 2.1(−3 ) TC 0.60 3.02 6.01

Fig. 2. The error E for Example 1 at time iteration k, k = 1, 5, 10, computed by (42) with rE = 2.6. Alg.2: Algorithm 2 has been applied with Ni = 60, ai = 5.9, i = 1, 2, and triangulation T1 on with r = 2.6. Alg.3: Algorithm 3 has been applied with Ni = 60, ai = 5.9, i = 1, 2. Alg.4: Algorithm 4 has been applied with triangulation T1 on with r = 2.6 and with exact boundary value g = u|∂ . Alg.4τ 1 and Alg.4τ 2 : Algorithm 4 has been applied with triangulation T1 and T2 , respectively, on with r = 7.1 and with the zero boundary value g = 0.

12

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

Fig. 3. The error E for Example 2 at time iteration k, k = 1, 5, 10, computed by (42) with rE = 3.0. Alg.2: Algorithm 2 has been applied with Ni = 60, ai = 6.7, i = 1, 2, and triangulation T1 on with r = 3.0. Alg.3: Algorithm 3 has been applied with Ni = 60, ai = 6.7, i = 1, 2. Alg.4: Algorithm 4 has been applied with triangulation T1 on with r = 3.0 and with exact boundary value g = u|∂ . Alg.4τ 1 and Alg.4τ 2 : Algorithm 4 has been applied with triangulation T1 and T2 , respectively, on with r = 8.1 and with the zero boundary value g = 0.

Fig. 4. The error E for Example 3 at time iteration k, k = 1, 5, 10, computed by (42) with rE = 0.75. Alg.2: Algorithm 2 has been applied with Ni = 60, ai = 0.94, i = 1, 2, and triangulation T1 on with r = 0.75. Alg.3: Algorithm 3 has been applied with Ni = 60, ai = 0.94, i = 1, 2. Alg.4: Algorithm 4 has been applied with triangulation T1 on with r = 0.75 and with exact boundary value g = u|∂ . Alg.4τ 1 and Alg.4 τ 2 : Algorithm 4 has been applied with triangulation T1 and T2 , respectively, on with r = 1.14 and with the zero boundary value g = 0.

Fig. 5. The mean computational time of Algorithms 2 and 4 at time iteration k = 1, 5, 10 obtained from Examples 1–3. For Algorithm 2 triangulation T1 and N1 = N2 = N = 30, 60, 120 are used; the average is calculated from the computational times TC in the proposed examples. For Algorithm 4 triangulation T2 is used.

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

13

In summary, the integral equation method used in Algorithm 2 seems to be quite accurate for the computation of the condition on the artiﬁcial boundary and Algorithm 2 seems to be an effective method for the numerical solution of the Burgers equation on the plane. In fact, it takes advantage on the eﬃciency of the FFT algorithm for the evaluation of the artiﬁcial boundary condition and of the accuracy of the ﬁnite elements methods for the solution of the corresponding boundary value problem in the domain of interest. The numerical results have been obtained by Matlab code running on a Workstation 5 processors Intel(R) Xeon(R) CPU X5650 @ 2.67 GHz, and operative system Red Hat Enterprise Linux Workstation, Release 6.4. 6. Conclusions The initial-boundary value problem for the Burgers equation on the plane is considered, and a method for the numerical solution of this problem is proposed. Given a bounded domain , this method computes the artiﬁcial boundary condition on ∂ by using the integral formulation of the original problem in the Fourier space. This integral equation is solved eﬃciently by the FFT algorithm, and the resulting boundary value problem for the Burgers equation on the bounded domain is solved numerically by the Galerkin method. This is a quite general method to compute the artiﬁcial boundary condition on a domain ; in particular, it can be easily generalized to take into account other discretization procedures in place of the Galerkin method and of the time-marching scheme. However, the directness of the Galerkin method and the explicit scheme for time discretization allow us to focus the main contribution of the present paper, that is the integral equation method for the evaluation of artiﬁcial boundary conditions. The proposed method gives quite interesting results, in fact the accuracy of the computed solution is comparable with the one obtained by the Galerkin method on a larger domain, where the boundary value of the solution is practically equal to zero. However, due to the eﬃciency of integral equation method, its computational cost is lower than the Galerkin method on this larger domain. So, this method deserves further studies: for the improvement of the numerical solution obtained by the integral equation method; for the adaptive choice of the interval R, containing the ”relevant support” of the solution, where the integral equation method is applied; for the evaluation of the scheme stability that is related on the stability properties of Volterra integral equations. An interesting study is also the generalization of the proposed method to the three dimensional Burgers equation, and to the Navier–Stokes equation. References [1] H. Bateman, Some recent researches on the motion of ﬂuids, Mon. Weather Rev. 43 (1915) 163–170. [2] J.M. Burgers, Mathematical examples illustrating relations occurring in the theory of turbulent ﬂuid motion, Trans. R. Neth. Acad. Sci. Amst. 17 (1939) 1–53.. [3] J.M. Burgers, A Mathematical Model Illustrating the Theory of Turbulence, vol. 1, Academic Press, New York, 1948. [4] J.D. Logan, An introduction to Nonlinear Partial Differential Equations, Wiley-Interscience, New York, 1994. [5] C.H. Su, C.S. Gardner, Korteweg–deVries equation and generalizations. III. Derivation of the Korteweg–deVries equation and Burgers equation, J. Math. Phys. 10 (1969) 536–539. [6] T. Musha, H. Higuchi, Traﬃc current ﬂuctuation and the Burgers equation, Jpn. J. Appl. Phys. 17 (1978) 811–816. [7] T. Nagatani, Density waves in traﬃc ﬂow, Phys. Rev. E 61 (20 0 0) 3564–3570. [8] N. Su, J.P.C. Watt, K.W. Vincent, M.E. Close, R. Mao, Analysis of turbulent ﬂow patterns of soil water under ﬁeld conditions using Burgers’ equation and porous suction-cup samplers, Aust. J. Soil Res. 42 (2004) 9–16. [9] E. Ben-Naim, S.Y. Chen, G.D. Doolen, S. Redner, Shock-like dynamics of inelastic gases, Phys. Rev. Lett. 83 (1999) 4069–4072. [10] S.F. Shandarin, Y.B. Zeldovich, Turbulence, intermittency, structures in a self-gravitating mediums, Rev. Mod. Phys. 61 (1989) 185–220. [11] J. Bec, K. Khanin, Burgers turbulence, Phys. Rep. 447 (2007) 1–66. [12] J.P. Bouchaud, M. Mézard, Velocity ﬂuctuations in forced Burgers turbulence, Phys. Rev. E 54 (1996) 5116–5121. [13] http://www.claymath.org/millennium/(accessed 22.09.14). [14] C. Boldrighini, S. Frigio, P. Maponi, Exploding solutions of the complex two-dimensional Burgers equations: computer simulations, J. Math. Phys. 53 (2012), 083101 1–12, doi:10.1063/1.4746814. [15] D. Li, Y.G. Sinai, Singularities of complex-valued solutions of the two-dimensional Burgers system, J. Math. Phys. 51 (2010) 015205. [16] M. Kardar, G. Parisi, Y.C. Zhang, Dynamic scaling of growing interfaces, Phys. Rev. Lett. 56 (1986) 889–892. [17] A.R. Bahadir, A fully implicit ﬁnite-difference scheme for two-dimensional Burgers’ equations, Appl. Math. Comput. 137 (2003) 131–137. [18] C.A.J. Fletcher, A comparison of ﬁnite element and ﬁnite difference solutions of the one- and two-dimensional Burgers’ equations, J. Comput. Phys. 51 (1983) 159–188. [19] M. Gülsu, T. Özis, Numerical solution of Burgers’ equation with restrictive Taylor approximation, Appl. Math. Comput 171 (20 05) 1192–120 0. [20] P.C. Jain, D.N. Holla, Numerical solution of coupled Burgers’ equations, Int. J. Non-Linear Mech. 13 (1978) 213–222. [21] W. Liao, An implicit fourth-order compact ﬁnite difference scheme for one-dimensional Burgers’ equation, Appl. Math. Comput. 206 (2008) 755–764. [22] E.N. Aksan, A numerical solution of Burgers’ equation by ﬁnite element method constructed on the method of discretization in time, Appl. Math. Comput. 170 (2005) 895–904. [23] N. Egidi, P. Maponi, An absorbing boundary condition for the Burgers equation, in: Proceedings of MASCOT13-IMACS/ISGG Workshop, IMACS Series in Computational and Applied Mathematics (submitted for publication). [24] S. Abbasbandy, M.T. Darvishi, A numerical solution of Burgers’ equation by time discretization of Adomian’s decomposition method, Appl. Math. Comput. 170 (2005) 95–102. [25] M. Dehghan, A. Hamidi, M. Shakourifar, The solution of coupled Burgers’ equations using Adomian–Pade technique, Appl. Math. Comput. 189 (2007) 1034–1047. [26] S.M. El-Sayed, D. Kaya, On the numerical solution of the system of two-dimensional Burgers’s equations by the decomposition method, Appl. Math. Comput. 158 (2004) 101–109.

14

N. Egidi, P. Maponi / Applied Mathematics and Computation 286 (2016) 1–14

[27] A. Gorguis, A comparison between ColeHopf transformation and the decomposition method for solving Burgers’ equations, Appl. Math. Comput 173 (2006) 126–136. [28] C. Boldrighini, N. Egidi, S. Frigio, P. Maponi, An integral equation method for the solution of the Burgers equation, in: Proceedings of MASCOT12 & ISGG12 Meeting on Applied Scientiﬁc Computing and Tools, Grid Generation, Approximation and Visualization, IMACS Series in Computational and Applied Mathematics (in press). [29] N. Egidi, P. Maponi, M. Quadrini, Integral equation method for the numerical solution of the Burgers equation, Math. Comput. Simul. (submitted for publication). [30] J. Yepez, Quantum lattice-gas model for the Burgers equation, J. Stat. Phys. 107 (2002) 203–224. [31] J. Zhang, G. Yan, Lattice Boltzmann method for one and two-dimensional Burgers equation, Physica A 387 (2008) 4771–4786. [32] X. Wu, J. Zhang, Artiﬁcial boundary method for two-dimensional Burgers’ equation, Comput. Math. Appl. 56 (2008) 242–256. [33] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comput. 31 (1977) 629–651. [34] J.B. Keller, D. Givoli, Exact non-reﬂecting boundary conditions, J. Comput. Phys. 82 (1989) 172–192. [35] L. Halperland, J. Rauch, Absorbing boundary conditions for diffusion equations, Numer. Math. 71 (1995) 185–224. [36] J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114 (1994) 185–200. [37] C.H. Bruneau, Boundary conditions on artiﬁcial frontiers for incompressible and compressible Navier–Stokes equations, Math. Model. Numer. Anal. 34 (20 0 0) 303–314. [38] Y. Zhou, Z.J. Wang, Absorbing boundary conditions for the Euler and Navier–Stokes equations with the spectral difference method, J. Comput. Phys. 229 (2010) 8733–8749. [39] D. Li, Y. Sinai, Blow ups of complex solutions of the 3D Navier–Stokes system and renormalization group method, J. Eur. Math. Soc. 10 (2008) 267–313, doi:10.4171/JEMS/111. [40] P. DuChateau, D. Zachmann, Applied Partial Differential Equations, Harper & Row, New York, 1989. [41] E.O. Brigham, The Fast Fourier Transform and its Applications, Prentice Hall, Englewood Cliffs, New Jersey, 1988. [42] S.G. Johnson, M. Frigo, A modiﬁed split-radix FFT with fewer arithmetic operations, IEEE Trans. Signal Process. 55 (2007) 111–119. [43] C. Cuvelier, A. Segal, A.A. VanSteenhoven, Finite Element Methods and Navier–Stokes Equations, Springer, Netherlands, 1986. [44] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, Dover Publication, Inc., New York, 1968.