Generalized collocation method for two-dimensional reaction-diffusion problems with homogeneous Neumann boundary conditions

Computers and Mathematics with Applications 56 (2008) 2360–2370

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Generalized collocation method for two-dimensional reaction-diffusion problems with homogeneous Neumann boundary conditions R. Revelli ∗ , L. Ridolfi Department of Hydraulics, Politecnico, Turin, Italy

article

info

Article history: Received 21 April 2008 Accepted 23 May 2008 Keywords: Collocation methods Nonlinear differential equations Reaction-diffusion models Two-dimensional problems

a b s t r a c t Nonlinear reaction-diffusion problems are ubiquitous in science, ranging from physics to engineering, from economy to biology. In particular, in this paper two-dimensional problems with generic Dirichlet boundary conditions or homogeneous Neumann boundary conditions are focused on. To simulate this class of mathematical problems, a quite simple generalized collocation method is proposed and developed. Moreover, a test case and an application to the Schnakenberg model are described, in order to show the efficiency of the simulation method and the goodness of the numerical results. © 2008 Elsevier Ltd. All rights reserved.

1. Introduction A well-known solution technique of nonlinear initial boundary value problem for nonlinear partial differential equations is the generalized collocation method, originally called the differential quadrature method. A concise account concerning the application of the method can be given with reference to the initial boundary value problem for models described by partial differential equations. The brief description given in what follows is analyzed in detail and generalized in the chapters which follow. Bearing this in mind, the method can be applied as follows: the method discretizes the original continuous model (and problem) into a discrete (in space) and continuous (in time) model, with a finite number of degrees of freedom, while the initial boundary value problem is transformed into an initial value problem for ordinary differential equations simply implementing the boundary conditions in the discretization points corresponding to the boundary and assigning initial conditions in all points of the collocation. This method is well documented in the literature on applied mathematics; it was first proposed by Bellman and Casti [1] and developed by several authors in deterministic [2] and stochastic frameworks [3]. Additional applications and developments are due, among others, to Chen [4,5] in the context of ordinary and partial differential equations, Shu and Richards [6] concerning the solution of Navier–Stokes Equations, Bert and Malick [7,8] concerning the analysis of the vibration of cylinder shells or rectangular plates. Let us also mention the papers by Karami and Malekzadeh [9,10], concerning the application of the method to various problems of vibration analysis, and by Artioli, Gould and Viola [11] and Artioli and Viola [12], concerning the solution of mathematical problems in structural mechanics. Finally, the reader can refer to the valuable monograph of Shu [6] that deals with various aspects of the generalized quadrature method in the context of fluid mechanics. Without discussing, at this stage, the validity of the method, with respect to alternative ones (which can only be done for well-defined problems), it is well documented that the method can provide a useful discretization of continuous models and efficiently deals with nonlinearities including the ones related to implicit boundary conditions. These features make the method interesting for engineering applications, as is documented in the review paper by Bert and Malik [7] and in the bibliography therein cited. Their paper provides an interesting and detailed report on the application of the original differential quadrature (collocation) method to several engineering problems.

∗

Corresponding author. E-mail addresses: [email protected] (R. Revelli), [email protected] (L. Ridolfi).

0898-1221/$ – see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2008.05.041

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

2361

On the other hand, it is known that the method does not generally work in some circumstances. For instance, referring to the Dirichlet problem, the classical Lagrange interpolation is not useful to deal with problems in unbounded domains or with solutions that are oscillating, with high frequency, in the space variables. This problem can be overcome by a suitable use of Sinc functions in the analysis of nonlinear problems, [13–15]. The use of Sinc functions was introduced by Bellomo and Ridolfi [16] and subsequently applied to the solution of various problems in applied sciences, see [17–20]. It can be concluded that several studies have been developed to generalize and improve the mathematical method, which has subsequently been applied to interesting engineering problems. These improvements have generated a mathematical method, called the generalized collocation method, which is useful to solve a large class of nonlinear problems in applied sciences. The book by Bellomo et al. [21] reports about various generalizations of the method with applications to the solutions of several methods of interest in engineering sciences. This paper develops the application of the method to the computational solution of nonlinear diffusion and reactiondiffusion problems with generic Dirichlet boundary conditions or homogeneous Neumann boundary conditions. The aforementioned class of equations has attracted a lot of attention because of its possibility to generate spatial patterns, and this explains the widespread use of reaction-diffusion models in several scientific fields, e.g., biology, chemistry, fluid dynamics, etc. [22–24]. This explains why spatial reaction-diffusion models represent one of the most investigated forms of partial differential equations (e.g., [25,26]). In the following two sections, the mathematical model and the numerical solution method will be described. Afterwards, a test case (where the analytical solution is known) will be simulated in order to evaluate the error bounds. Finally, an application to a nonlinear reaction-diffusion model will be discussed and some conclusions will be drawn. 2. The mathematical model The computational method developed in this present paper refers to systems of equations. This section first describes the statement of the mathematical problems under consideration in the case of a scalar equation, while the technical generalization to systems of equations is subsequently illustrated. Specifically, let us consider the following class of equations:

∂u ∂ 2u ∂ 2u = f (u) + 2 + 2 , ∂t ∂x ∂y

(1)

where u = u(t , x) :

[0, T ] × D → R,

x = (x, y),

(2)

is the dependent variable, f = f (u) is a nonlinear function of u, and where the domain D ⊆ R of the space variables has a boundary ∂ D. The function f is assumed to be a regular function of its argument. After a suitable dimensional analysis, one can assume that u ∈ [0, 1] and that the domain D is bounded as follows: 2

D ⊆ [0, 1] × [0, 1].

(3)

The statement of the mathematical problems is proposed assuming boundary ∂ D of D regular so that a unit vector n, directed towards the interior of D, can be defined in any point of ∂ D. Bearing all above in mind consider the following problems: Problem 1 (Dirichlet Boundary Conditions). Consider the initial boundary value problem for Eq. (1) with initial condition u(0, x) = φ(x),

x = (x, y) ∈ D,

(4)

and time depending Dirichlet boundary condition

∀t ∈ [0, T ],

x ∈ ∂D :

u(t ; x ∈ ∂ D) = α(t , x),

(5)

given as smooth function of their argument consistent, for t = 0, with the initial condition (4). Problem 2 (Neumann Boundary Conditions). Consider the initial boundary value problem for Eq. (1) with initial condition (4) and homogeneous Neumann boundary condition

∀t ∈ [0, T ],

x ∈ ∂D :

∂u (t ; x ∈ ∂ D) = 0, ∂n

(6)

consistent, for t = 0, with the initial condition (4). The above statement can be technically particularized to specific problems when D is a rectangular domain in R2 . Up to suitable normalization of the independent variables, one can assume that D = [0, 1] × [0, 1], which means that D has been transformed in the unit square of R2 . In this case, the boundary ∂ D consists only in four distinct parts, namely,

∂ D = Γ1 ∪ Γ2 ∪ Γ3 ∪ Γ4 ,

(7)

2362

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

where

Γ1 = {(x, 0), x ∈ [0, 1]}, Γ3 = {(0, y), y ∈ [0, 1]},

Γ2 = {(x, 1), x ∈ [0, 1]}, Γ4 = {(1, y), y ∈ [0, 1]}.

(8)

In this case, one can restate Dirichlet boundary conditions as follows, for any t ∈ [0, T ]

∀x ∈ [0, 1] :

u(t ; x, 0) = α(t , x), u(t ; x, 1) = β(t , x),

(9)

∀y ∈ [0, 1] :

u(t ; 0, y) = γ (t , y), u(t ; 1, y) = δ(t , y),

(10)

and

where α , β , γ and δ are given smooth functions of their arguments consistent with the initial condition (4) and satisfying the consistency conditions

α(t , 0) = γ (t , 0),

α(t , 1) = δ(t , 0),

β(t , 0) = γ (t , 1),

β(t , 1) = δ(t , 1).

(11)

In the same way, Neumann boundary conditions are stated as follows, for any t ∈ [0, T ]:

∀x ∈ [0, 1] :

∂u ∂u (t ; x, 0) = (t ; x, 1) = 0, ∂y ∂y

(12)

∀y ∈ [0, 1] :

∂u ∂u (t ; 0, y) = (t ; 1, y) = 0, ∂x ∂x

(13)

and

consistent, for t = 0, with the initial condition (4). The above reasoning holds also when the space domain is convex with respect to both axes. Otherwise, additional analysis is needed. The above statement of problems refers to equations with second-order derivatives on both space variables. It requires boundary conditions on the closed contour ∂ D, while equations with the first-order derivatives need boundary conditions on open contours. Higher-order equations require additional conditions. We also point out that the statement may differ for each equation: for instance Dirichlet conditions for the first equation and Neumann condition for the second. Remark 2.1. The statement of problems can be straightforwardly generalized to systems of equations simply implementing, with the same rules, initial and boundary conditions to each equation consistently with its structure.  3. Solution method Let us consider the class of problems described in the preceding section. The solution method is developed through the following steps: The first step consists in the discretization of the spatial domain into the collocation Ixy = {x1 , . . . , xn ; y1 , . . . , ym },

(14)

which may be equally spaced xi = (i − 1)hx

hx =

yj = (j − 1)hy

hy =

1 n−1 1 m−1

i = 1, . . . , n, j = 1, . . . , m,

(15) (16)

or identified by Chebychev–type collocation with decreasing values of the measures |xi − xj | and |yi − yj | towards the borders: xi = yj =

1 2 1 2

− −

1 2 1 2

 cos

n−1

 cos

i−1

π

j−1 m−1



π



i = 1, . . . , n, j = 1, . . . , m.

(17)

(18)

The dependent variable u = u(t , x, y) is interpolated and approximated by the values uij = u(t , xi , yj ) as follows: u(t , x, y) ' unm (t , x, y) =

n X m X i=1 j=1

Li (x)Lj (y)uij (t ),

(19)

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

2363

where the expressions of the Lagrange polynomials Li (x) and Lj (y) are Li (x) =

(x − x1 ) . . . (x − xi−1 )(x − xi+1 ) . . . (x − xn ) , (xi − x1 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )

(20)

Lj (y) =

(y − y1 ) . . . (y − yj−1 )(y − yj+1 ) . . . (y − ym ) . (yj − y1 ) . . . (yj − yj−1 )(yj − yj+1 ) . . . (yj − ym )

(21)

and

This interpolation is then used to approximate the spatial partial derivatives in the nodal points of Ixy , namely n X ∂r u (r ) (t ; xi , yj ) ' ahi uhj (t ), r ∂x h=1

(22)

m X ∂r u (r ) ( t ; x , y ) ' bkj uik (t ), i j ∂ yr k=1

(23)

n X m X ∂ 2u (1) (1) ahi bkj uhk (t ), (t ; xi , yj ) ' ∂ x∂ y h=1 k=1

(24)

and

where (r )

ahi =

dr Lh dxr

(xi ),

(25)

(yj ).

(26)

and (r )

bkj =

dr Lk dyr

Technical calculations provide the following result (1)

ahi

Q (xi ) Q = , (xi − xh ) (xh )

(1)

aii =

1

X

x − xh h6=i i

,

h 6= i,

(27)

and (1)

aii =

1

X

x − xh h6=i i

,

h = i,

(28)

where

Y Y (xi − xp ), (xi ) =

Y

(xh ) =

p6=i

Y (xh − xp ),

(29)

p6=h

while for the y direction (1)

bkj =

(yj ) Q , (yj − yk ) (yk ) Q

(1)

bjj =

X

1

k6=j

yj − yk

,

k 6= j,

(30)

and (1)

bjj =

X

1

k6=j

yj − yk

k = j,

(31)

where

Y Y (yj ) = (yj − yp ),

Y Y (yk ) = (yk − yp ).

p6=j

p6=k

(32)

Higher-order coefficients may be computed exploiting the following recurrence formula (r )

ahi = r

(1) (r −1)

ahi aii

(r −1)

+

ahi

xh − xi

! ,

(r +1)

aii

=−

X h6=i

(r +1)

ahi

,

(33)

2364

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

and (r )

bkj = r

(r −1)

(1) (r −1)

+

bkj bjj

!

bkj

xk − xj

(r +1)

,

bjj

=−

X

(r +1)

bkj

.

(34)

k6=j

The initial boundary value problem is then transformed into an initial value problem for ordinary differential equations that describes the evolution of the values uij (t ) of u in the nodes. The boundary conditions have to be enforced in the nodal points on the boundaries x1j , xnj , xi1 , and xim (i = 1, . . . , n; j = 1, . . . , m) of the domain D. Finally, the solution of the initial boundary value problem is obtained through the numerical solution of the initial value problem for ordinary differential equations, and, then, the interpolation of the solution by using Eq. (19). If the previous procedure is applied to Problem 1, with conditions (9) and (10), then, one obtains the system of ordinary differential equations duij dt

= f (uij ) +

n X

(2)

ahi uhj (t ) +

h=1

m X

(2)

bkj uik (t ),

(35)

k=1

supplemented with the initial conditions: ui1 = α(t , xi ),

uim = β(t , xi ),

(36)

u1j = γ (t , yj ),

unj = δ(t , yj ),

(37)

and for any i = 1, . . . , n, j = 1, . . . , m and r = 1, 2. Remark 3.1. The equations concerning the boundary conditions can be written in differential form. For example, dui1 dt

=

dα dt

,

duim dt

=

dβ dt

. 

(38)

The solution method for the Problem 2 can be developed according to the same conceptual paths that have been described for the Problem 1. The only difference concerns the fact that Neumann boundary conditions have to imposed. This means that it is necessary to set the consistency with the boundary conditions that yield the following linear algebraic system

 n−1 X (1) (1) (1)   ak1 ukj a u + a u = − 1j nj  11 n1   k=2    n−1 X   (1) (1) (1)  akn ukj  a1n u1j + ann unj = − k=2

(39)

m −1 X   (1) (1) (1)  bk1 uik b u + b u = −  i1 im m1   11  k=2   m −1  X  (1) ( 1)  ( 1 ) b1m uim + bmm uim = − bkm uik , k=2

(with i = 2, . . . , n − 1; j = 1, . . . , m). Such a system can be algebraically solved with respect to u1j , unj , u1j , and u1j and leads to the following: u1j =

unj =

ui1 =

ann

n−1 X

an1 a1n − a11 ann k=2 a1n

n−1 X

an1 a1n − a11 ann k=2 bmm

ak1 ukj −

ak1 ukj +

an1 a1n − a11 ann k=2

bm1 b1m − b11 bmm k=2

n−1 X

a11

an1 a1n − a11 ann k=2

m−1

X

n−1 X

an1

bk1 uik −

bm1

akn ukj ,

(40)

akn ukj ,

(41)

m−1

X

bm1 b1m − b11 bmm k=2

bkm uik ,

(42)

and uim =

b1m

m−1

X

bm1 b1m − b11 bmm k=2

bk1 uik −

b11

m−1

X

bm1 b1m − b11 bmm k=2

bkm uik .

The above relations can be introduced into system (35), ensuring the self-consistency of the system itself.

(43)

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

2365

Remark 3.2. The technical generalization to systems of partial differential equations is immediate. In this case the dependent variable is a vector, say u. The solution technique leads to a system of equations for each component of u.  Remark 3.3. The system can be forwardly applied to more general class of equations. For instance:

∂u =f ∂t



 ∂ 2u ∂ 2u . t , x, u, 2 , ∂ x ∂ y2

(44)

In this case the corresponding discrete model is as follows: duij dt

t , xi , yj , uij ,

=f

n X

(2)

ahi uhj (t ),

h=1

m X

(2)

!

bkj uik (t ) . 

(45)

k=1

4. Test case: The heat equation over a square In this section the simulation of the heat equation is proposed as test case; in fact, the knowledge of the analytical solution allows us to evaluate the errors due to the numerical method. Let us consider the two-dimensional heat equation 1 ∂u = 2 ∂t π



 ∂ 2u ∂ 2u + 2 , ∂ x2 ∂y

(46)

over the square domain D = [0, 1] × [0, 1], where the dimensionless dependent variable is u = u(t , x) :

[0, T ] × D → [0, 1],

x = (x, y) ∈ D.

(47)

The partial differential equation (46) is linked to the following initial condition: u(t = 0, x, y) = u0 (x, y) = cos

π 2

(x + y) + sin π (x − y),

(48)

and boundary conditions

π t u(t , x, 0) = e− 2 cos x + e−2t sin π x, 2

(49)

π t u(t , x, 1) = e− 2 cos (x + 1) + e−2t sin π (x − 1), 2

(50)

π t u(t , 0, y) = e− 2 cos y + e−2t sin π y, 2

(51)

π t u(t , 1, y) = e− 2 cos (1 + y)e−2t sin π (1 − y). 2

(52)

and

The initial boundary problem (46)–(52) admits the following analytical solution:

π t u(t , x, y) = e− 2 cos (x + y) + e−2t sin π (x − y). 2

(53)

The knowledge of such an analytic solution allows the comparison with that obtained by the computational method and shall permit to evaluate the accuracy of the numerical simulation. Applying the procedure developed above, the corresponding system of ordinary differential equations is: duij dt

=

1

n X

π2

h=1

(2)

ahi uhj (t ) +

m X

(2)

!

bkj uik (t ) ,

(54)

k=1

with i = 1, . . . , n and j = 1, . . . , m; while the initial conditions, given by Eq. (48), are uij = cos

π 2

(xi + yj ) + sin π (xi − yj ),

(55)

for i = 1, . . . , n and j = 1, . . . , m. Finally, boundary conditions are implemented according to Eqs. (49)–(52). Simulations have been obtained with n = m = 9 nodes. Fig. 1 shows the numerical solution u(t , x, y) for t = 1 respectively. One notes the smoothness of the solution despite the small number of nodes. This accuracy of the solution method is emphasized in Fig. 2 that shows the error

E n = ku − un k∞ = sup |u − un |, x∈D

(56)

2366

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

Fig. 1. Heat equation over a square: u versus time and space: t = 1.

Fig. 2. Heat equation over a square: computational error E n versus time and space: t = 1.

at t = 1, with respect to the analytical solution (53). It can be seen that the error is lower than 4 × 10−7 . Finally Fig. 3 shows the errors E n and the average error identified by the L1 -distance

E1n =

1 meas(D)

Z

|u − un |(t , x)dx,

(57)

D

with D = [0, 1] and meas(D) = 1, versus the number of node n. The errors confirm the good result obtained by Lagrange interpolation and show the fast decrement of the errors when the nodes number is increasing: when n is up to 5–10, the errors practically tend to zero. 5. Application: A nonlinear reaction-diffusion equation As a quite severe application of the proposed computational method, we consider here a model introduced by Schnakenberg [27] to describe a series of trimolecular autocatalytic reactions. The model consists in the following system of reaction-diffusion equations in two space dimensions:

  2  ∂u ∂ u ∂ 2u   = κ(a − u + u2 v) + ε +  2 2 ∂t  2 ∂x 2 ∂y ∂v ∂ v ∂ v    = κ(b − u2 v) + + 2 , ∂t ∂ x2 ∂y

(58)

where v and u are nondimensionalized concentrations of, respectively, self-activating and self-inhibiting chemical substances, and κ , a, b, and ε are dimensionless parameters.

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

2367

Fig. 3. Heat equation over a square: error E n (continuous line) and E1n (dashed line) versus number of nodes: t = 1.

The interesting aspect of the above model is its capability to select peculiar spatial patterns where initial small perturbations are amplified and spread, leading to the formation of spots that slowly move and interact. If relatively large diffusion coefficients are used, the model is stiff and, for this reason, Schnakenberg model is a severe test for the numerical algorithms (e.g., [26]). The two dependent variables u and v are such that: u = u(t , x) :

[0, 1] × D → [0, 1],

x = (x, y) ∈ D,

(59)

v = v(t , x) :

[0, 1] × D → [0, 1],

x = (x, y) ∈ D,

(60)

and

where D = [0, 1] × [0, 1]. The initial boundary problem is stated by linking the above model to the following initial conditions: u(0, x, y) = a + b +

1 1000

1 2 1 2 e−100[(x− 3 ) +(y− 2 ) ] ,

(61)

and

v(0, x, y) =

b

(a + b)2

,

(62)

while the boundary conditions are given by

 ∂u   (t , 0, y) = ∂x ∂v   (t , 0, y) = ∂x

∂u (t , 1, y) = 0 ∂x ∂v (t , 1, y) = 0, ∂x

(63)

 ∂u   (t , x, 0) = ∂y ∂v   (t , x, 0) = ∂y

∂u (t , x, 1) = 0 ∂y , ∂v (t , x, 1) = 0 ∂y

(64)

and

which correspond to homogeneous Neumann conditions. The application of the method generates the following system of ordinary differential equations:

!  n m X X duij  ( 2 ) ( 2 )  = κ(a − uij + u2ij vij ) + ε aki ukj + bkj uik   dt k=1 k=1! n m X X  dvij  (2) (2)  = κ(b − u2 vij ) + a vkj + b vik ,  dt

ij

ki

k=1

(65)

kj

k=1

where i = 2, . . . , n − 1 and j = 2, . . . , m − 1, and with the initial conditions given by (61) and (62). The solution of the initial boundary value problem requires that we express the values of the two dependent variables u and v in the boundary

2368

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

Fig. 4. Reaction-diffusion equations: u versus time and space: t = 1.

nodes. Technical calculations yield

 1   u1,j = (1) (1)   (1) (1)  an1 a1n − a11 ann       1    un,j = a(1) a(1) − a(1) a(1) n1 1n 11 nn  1    v = (1) (1)  (1) (1)  1,j  a a − a11 ann  n1 1n     1   vn,j = (1) (1) (1) (1) an1 a1n − a11 ann

(1)

ann

n−1 X

(1)

(1)

ak1 ukj − an1

n−1 X

k=2

(1)

a1n

n−1 X

ann

n−1 X

(1)

(1)

ak1 ukj − a11

n−1 X

a1n

n−1 X

akn ukj (1)

!

akn ukj

k=2

(1)

(1)

ak1 vkj − an1

n−1 X

k=2

(1)

!

k=2

k=2

(1)

(1)

(1)

(66)

!

akn vkj

k=2

(1)

(1)

ak1 vkj − a11

n−1 X

k=2

(1)

!

akn vkj

,

k=2

for j = 1, . . . , m, and

 1   u1,j = (1) (1)   (1) (1)  bm1 b1m − b11 bmm       1    um,j = b(1) b(1) − b(1) b(1) m1 1m

11

mm

(1)

bmm

m−1

X

(1)

(1)

bk1 ukj − bm1

k=2

(1)

b1m

m−1

X k=2

m−1

X

(1)

!

bkm ukj

k=2

(1)

(1)

bk1 ukj − b11

m−1

X

(1)

!

bkm ukj

k=2

! m −1 m −1  X X 1  ( 1 ) ( 1 ) ( 1 ) ( 1 )   bmm bk1 vkj − bm1 bkm vkj v = (1) (1)  (1) (1)  1,j  bm1 b1m − b11 bmm k=2 k=2   !   m −1 m −1 X X  1  (1) (1) (1) (1)  b1m bk1 vkj − b11 bkm vkj , vm,j = (1) (1) (1) (1) bm1 b1m − b11 bmm k=2 k=2

(67)

for i = 2, . . . , n − 1. Figs. 4 and 5 show the evolution of the dependent variables u and v at t = 1, respectively. Lagrange polynomials are used with 51 nodes and the following parameters are adopted: κ = 100, a = 0.1648, b = 0.8352, and ε = 0.05 (see also [26]). The typical feature of patterns generated by the Schnakenberg model, that is the amplification and the migration of the bumps, is plainly recovered by the simulations. This characteristic dynamics is even more evident in the representation by contour lines shown in the corresponding Figs. 6 and 7, where such lines join points of domain having the same value of population, u and v . We wish to point out that the present numerical simulation agree very well with the similar results reported by Hundsdorfer and Verwer [26], where different and more involved numerical methods are used. 6. Critical analysis A generalized collocation method for the numerical simulation of two-dimensional reaction-diffusion problems with generic Dirichlet boundary conditions or homogeneous Neumann boundary conditions. Several important processes in many scientific fields are modelled according to this class of equations. An account has been given Section 1, while recent

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

Fig. 5. Reaction-diffusion equations: v versus time and space: t = 1.

Fig. 6. Reaction-diffusion equations: contour plot of u versus time and space: t = 1.

Fig. 7. Reaction-diffusion equations: contour plot of v versus time and space: t = 1.

2369

2370

R. Revelli, L. Ridolfi / Computers and Mathematics with Applications 56 (2008) 2360–2370

applications show that diffusion models with source terms can also be used to model complex systems in life and biological sciences, see [24,28]. The proposed simulation method is quite simple, robust, and computationally efficient. A test case has shown that low errors are obtainable with a small number of nodes, while the application to a strongly nonlinear reaction-diffusion problem has shown how the method gives very good results and is competitive with respect to other more usual numerical methods (such as finite volumes [29], fractional steps [30], finite elements [31], and so on) at least when the geometry of the system is sufficiently simple as in the case considered in this present paper. The advantage essentially consists in a rapid implementation of boundary condition and in the robust control of nonlinearities. When the space variables are defined on a complex geometry it may be useful decomposing the space domain into several sub-domains properly connected at the boundary. The development of the method is straightforward. References [1] R. Bellman, J. Casti, Differential quadrature and long term integration, Journal of Mathematical Analysis and Applications 34 (1971) 235–238. [2] A. Satofuka, A new explicit method for the solution of parabolic differential equation, in: Numerical Methodologies and Properties in Heat Transfer, Hemisphere, New York, 1983, pp. 97–108. [3] N. Bellomo, F. Flandoli, Stochastic partial-differential equations in continuum physics — on the foundations of the stochastic interpolation method for itos type equations, Mathematics and Computers in Simulation 31 (1–2) (1989) 3–17. [4] C.N. Chen, Generalization of differential quadrature discretization, Numerical Algorithms 22 (1999) 167–182. [5] C.N. Chen, Differential quadrature element analysis using extended differential quadrature, Computer Mathematics with Applications 39 (2000) 65–79. [6] C. Shu, B.E. Richards, Application of generalized differential quadrature to solve two-dimensional incompressible navier-stokes equations, International Journal for Numerical Methods in Fluids 15 (1992) 791–798. [7] C. Bert, M. Malik, Differential quadrature method in computational mechanics, ASME Review 49 (1996) 1–27. [8] C. Bert, M. Malik, Three-dimensional elasticity solutions for free vibrations of rectangular plates by the differential quadrature method, International Journal of Solids and Structures 35 (1998) 299–318. [9] G. Karami, P. Malekzadeh, A new differential quadrature methodology for beam analysis and the associated differential quadrature element method, Computational Methods Applied Mechanical Engineering 191 (2002) 3509–3526. [10] G. Karami, P. Malekzadeh, Vibration of non-uniform thick plates on elastic foundation by differential quadrature method, Engineering Structures 26 (2004) 1473–1482. [11] E. Artioli, P. Gould, E. Viola, A differential quadrature method solution for the shear-deformable shells of revolution, Engineering Structures 27 (2005) 1879–1892. [12] E. Artioli, E. Viola, Static analysis of shear-deformable shells of revolution via g.d.q. method, Structural Engineering and Mechanics 19 (2005) 459–475. [13] J. Lund, K. Bowers, Sinc Methods, SIAM, Philadelphia, 1992. [14] F. Stenger, Numerical methods based on whittaker cardinal or sinc functions, SIAM Review 23 (1983) 165–224. [15] F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, Springer, Berlin, 1993. [16] N. Bellomo, L. Ridolfi, Solution of nonlinear initial-boundary value problems by Sinc collocation-interpolation methods, Computers & Mathematics with Applications 29 (4) (1995) 15–28. [17] R. Revelli, L. Ridolfi, Sinc collocation-interpolation method for the simulation of nonlinear waves, Computers & Mathematics with Applications 46 (8/9) (2003) 1443–1453. [18] R. Revelli, L. Ridolfi, P. Massarotti, Nonlinear convection-dispersion models with a distributed pollutant source I - Direct initial boundary value problems, Mathematical and Computer Modelling 39 (9/10) (2004) 1023–1034. [19] N. Bellomo, Nonlinear models and problems in applied sciences: From differential quadrature to generalized collocation methods, Mathematical and Computer Modelling 26 (4) (1997) 13–34. [20] R. Revelli, L. Ridolfi, Nonlinear convection-dispersion models with a localized pollutant source - II On a class of inverse problems, Mathematical and Computer Modelling 42 (2005) 601–612. [21] N. Bellomo, B. Lods, R. Revelli, L. Ridolfi, Generalize Collocation Methods - Solutions to Nonlinear Problems, Birkhäuser, Boston, 2008. [22] M.C. Cross, P.C. Hohenberg, Pattern-formation outside of equilibrium, Reviews of Modern Physics 65 (1993) 851–1112. [23] J.D. Murray, Mathematical Biology, in: Biomathematics, Springer-Verlag, Berlin, 1993. [24] M.B. Short, M.R. D’orsogna, V.B. Pasteur, G.E. Tita, P.J. Bratingham, A.L. Bertozzi, L.B. Chayes, A statistical model of criminal behaviour, Mathematical Models & Methods in Applied Sciences 18 (2008). [25] P. Grindrod, The Theory and Application of Reaction-Diffusion Equations. Pattern and Waves, in: Oxford Applied Mathematics and Computing Sciences Series, Oxford University Press, 1996. [26] W. Hundsdorfer, J.G. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer, Berlin, 2003. [27] J. Schnakenberg, Simple chemical reaction system with limit cycle behaviour, Journal of Theoretical Biology 81 (1979) 389–400. [28] N. Bellomo, N.K. Li, P.K. Maini, On the foundation of cancer modelling: Selected topics, speculations and perspectives, Mathematical Models & Methods in Applied Sciences 18 (2008) 593–646. [29] E. Bertolazzi, G. Manzini, On vertex reconstructions for cell-centered finite volume approximations of 2D anisotropic diffusion problems, Mathematical Models & Methods in Applied Sciences 17 (2007) 1–32. [30] A. Quaini, A. Quarteroni, A semi-implicit approach for fluid structure interaction based on algebraic fractional step methods, Mathematical Models & Methods in Applied Sciences 17 (2007) 957–983. [31] J. Xu, Y. Zhu, Uniform convergent multigrid methods for elliptic problems with strongly discontinuous coefficients, Mathematical Models & Methods in Applied Sciences 18 (2008) 77–107.