rjs-2008-lwsomathilake-final.dvi RUHUNA JOURNAL OF SCIENCE Vol. 4, September 2009, pp. 1–12 http://www.ruh.ac.lk/rjs/ I S S N 1800-279X © 2009 Faculty of Science University of Ruhuna. Numerical Solutions of Reaction-Diffusion systems with coupled diffusion terms L.W. Somathilake Department of Mathematics, University of Ruhuna, Matara, Sri Lanka. (Correspondence: sthilake@maths.ruh.ac.lk) Abstract. Numerical solutions of some initial-boundary value problem associated with a particular reaction-diffusion systems namely Gray-Scott model responsible for spatial pattern formation are considered. The aim of this paper is to numerically solve the above system with coupled diffu- sion terms. Firstly, using some linear transformations, a general form of diffusion coupled reaction- diffusion system is converted into reaction-diffusion system with uncoupled diffusion terms and then, some finite difference schemes (based on (Hoff 1978)) are constructed to obtain the solutions. Finally, the graphical representation of the numerical solutions are presented. Key words : Reaction-Diffusion equations, Pattern formations, Finite Difference Methods. 1. Introduction Properties of analytical solutions of diffusion-coupled reaction-diffusion systems have been reported for instance in (Kirane 1989, Badraoui 2002, P. Collet 1994, Kirane 1(1993). For example consider the diffusion-coupled reaction-diffusion system ∂u ∂t = a∆u − uh(v) x ∈ Ω t > 0, ∂v ∂t = b∆u + d∆v + uh(v) x ∈ Ω t > 0,      (1) with initial conditions u(x, 0) = u0(x), v(x, 0) = v0(x), x ∈ Ω (2) on a bounded domain Ω ⊂ Rn with Neumann boundary conditions, b > 0, a 6= d, v0 ≥ bu0 a − d ≥ 0, h(s) is a differentiable nonnegative function on R. The major result of this has been regarded in (Kirane 1989) while the existence of global solutions for the system (1) on unbounded domains has been reported in (Badraoui 2002). The existence of global solutions in Rn for (1) with h(s) = vm has been studied in (P. Collet 1994). The quasilinear system of reaction-diffusion equations ∂u ∂t = ∇.(a(u)∇u) − uh(u)v x ∈ Ω t > 0, ∂v ∂t = ∇.(b(v)∇v) + uh(u)v − λv x ∈ Ω t > 0,      (3) 1 LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... 2 Ruhuna Journal of Science 4, pp. 1–12, (2009) with initial conditions u(x, 0) = u0(x), v(x, 0) = v0(x), x ∈ Ω (4) and with Neumann or Dirichlet boundary conditions, is studied in (Kirane 1(1993) where in particular, the existence of a globally bounded solution has been shown. Moreover, the large time behavior of the solution has also been discussed. This type of mathematical models may arise when constructing mathematical models for the reaction-diffusion processes of substances (eg. chemicals, species, deceases etc.) such that diffusion of one substance depend on the diffusion of some other substance appear in precess. The Turing instability condition which give rise to pattern formation in reaction-diffusion system are well known (Murray 2003). Some of such patterns formation reaction diffusion systems are Gray-Scott model(Webpage ????a), Fitzhugh-Nagumo model(Webpage ????b). Many researchers have studied numerical solutions of these systems without coupled diffu- sion terms(Seaid 2001, Turk 1992, Grindrod 1991). Above investigations of analytical solutions of diffusion-coupled reaction-diffusion sys- tems and numerical solutions of pattern formation reaction diffusion systems motivated us to investigate numerical solutions of diffusion-coupled pattern formation reaction-diffusion systems. We aim in investigating numerical solutions of diffusion-coupled pattern formation reaction-diffusion systems using finite difference techniques. The paper is organized as fol- lows: In Section 2, we consider a general form of diffusion-coupled reaction-diffusion sys- tem and its transformation into a reaction-diffusion system with uncoupled diffusion terms. In Section 3, some finite difference schemes are constructed for a general form of reaction- diffusion system with uncoupled diffusion terms. In Section 4, numerical solutions, which are obtained using constructed semi-implicit finite difference method, of Gray-Scott type diffusion-coupled reaction diffusion system are presented. When semi-implicit finite differ- ence scheme is applied to a reaction diffusion system, a system of linear equations arise in each time step. In order to solve these linear systems conjugate-gradient method has been used in computer programs. 2. General form of a diffusion-coupled reaction-diffusion system Let v1(x, t), v2(x, t) ∈ Rm be variables explaining some characteristics of model system (eg. concentration, population size) in suitable units. Consider the following general form of diffusion-coupled reaction-diffusion system: ∂v1 ∂t = A∆v1 + F(v1, v2), x ∈ Ω, t > 0, ∂v2 ∂t = B∆v1 + D∆v2 + G(v1, v2), x ∈ Ω, t > 0,      (5) with initial conditions v1(x, 0) = v10(x) for x ∈ Ω v2(x, 0) = v20(x) for x ∈ Ω } (6) LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... Ruhuna Journal of Science 4, pp. 1–12, (2009) 3 and boundary conditions ∂v1 ∂n = 0, for x ∈ ∂Ω, t > 0. ∂v2 ∂n = 0, for x ∈ ∂Ω, t > 0.          (7) Here F = (F1, F2, ..., Fm) T , G = (G1, G2, ..., Gm) T , Ω ⊂ Rd is a connected, bounded open set with piecewise smooth boundary, and n is the outward unit normal vector to the boundary. A, B, and D are diagonal matrices of order m whose ith diagonal entries are respectively ai, bi, and di each of which is a real constants and ai 6= di for i = 1, 2, ..., m. 2.1. Transformation of a diffusion-coupled reaction-diffusion system to a reaction-diffusion system with uncoupled diffusion terms Let us define the Linear transformations L : (v1, v2)−→ ( v1, v2 − (A − D) −1Bv1 ) = (v, w). Under L, the system of equations (5) is transformed to ∂v ∂t = A∆v + F1 ( v, w + (A − D)−1Bv ) , x ∈ Ω, t > 0, ∂ ∂t (w + (A − D)−1Bv) = B∆v +D∆(w + (A − D)−1Bv) +G1 (v, w + (A − D)−1Bv) , x ∈ Ω, t > 0,                  (8) By simplifying have ∂v ∂t = A∆v + F1 ( v, w + (A − D)−1Bv ) , x ∈ Ω, t > 0, ∂w ∂t = −(A − D)−1(BA − DB)∆v + B∆v +D∆w + G1 ( v, w + (A − D)−1Bv ) −(A − D)−1BF1 ( v, w + (A − D)−1Bv ) , x ∈ Ω, t > 0,                  (9) This system can be written in the form: ∂v ∂t = A∆v + F2(v, w), x ∈ Ω, t > 0, ∂w ∂t = D∆w + G2(v, w), x ∈ Ω, t > 0,          (10) LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... 4 Ruhuna Journal of Science 4, pp. 1–12, (2009) where F2(v, w) = F1 (v, w + (A − D)−1Bv) , G2(v, w) = G1 (v, w + (A − D)−1Bv) −(A − D)−1BF1 (v, w + (A − D)−1Bv)        (11) The initial conditions are reduced to: v(x, 0) = v0(x) = v1(x, 0) = v10(x), for x ∈ Ω w(x, 0) = w0(x) = v2(x, 0) − (A − D) −1Bv1(x, 0) = v20(x) − (A − D) −1Bv10(x), for x ∈ Ω        (12) and boundary conditions are reduced to ∂v ∂n = 0, for x ∈ ∂Ω, t > 0. ∂w ∂n = 0, for x ∈ ∂Ω, t > 0.          (13) Now take u = (v, w)T , C = ( A 0 0 D ) , and F = (F2(v, w), G2(v, w)) T . Then the system is reduced to the form: ∂u ∂t = C∆u + F (u), x ∈ Ω, t > 0, (14) with initial data u(x, 0) = u0(x) for x ∈ Ω and boundary conditions ∂u ∂n = 0 for x ∈ ∂Ω t > 0. This system has 2m number of equations. The reaction-diffusion system (5) with coupled diffusion terms is now reduced to the reaction-diffusion system (14) which has no coupled diffusion terms. This system can be solved numerically for approximate solutions of u = (v, w)T . Finally, approximations for v1 and v2 can be obtained using L−1 transformation. 3. Finite Difference Schemes for reaction-diffusion systems In this section finite difference schemes for reaction diffusion system (14) are constructed based on (Hoff 1978). Let x = (x1, x2, ..., xd ) ∈ R d , ∆xi = hi (i = 1, 2, ..., d) be an increment in xi (i = 1, 2, ..., d) and τ be an increment in t. Also let xk = (k1h1, k2h2, ..., kd hd ) for k = (k1, k2, ..., kd ) ∈ Z d and tn = nτ for n ∈ Z. We shall approximate u(xk, tn) ≡ U nk for k ∈ I ; where I is an appropriate index set contained in Zd such that k ∈ S implies xk ∈ Ω. LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... Ruhuna Journal of Science 4, pp. 1–12, (2009) 5 Let M = 2m, S = d ∏ i=1 [ai, bi] and Y = {{Uk}k∈I : Uk ∈ R M} be a vector space of RM valued functions on S. Suppose that the second order accurate operator ∆2j on Y are constructed such that ∣ ∣ ∣ ∣ ( ∆2j u h2j ) k − ∂2u ∂x2j (xk) ∣ ∣ ∣ ∣ ∞ ≤ C1||u||4 h 2 j , j = 1, 2, ..., d ( ∆2j u h2j ) k = ∂2u ∂x2j (η) for some η ∈ Ω j = 1, 2, ..., d where C1 is a constant which is independent of u and h j. Assuming that I has been defined and operators ∆2j have been constructed, we replace the differential equation (14) with the finite difference equation U n+1k −U n k τ = d ∑ j=0 ∆2j h2j ( θU n+1k + (1 − θ)U n k ) + F (15) Where F is evaluated at (xk, tn,U n k ). Here 0 ≤ θ ≤ 1, and U n ∈ Y is defined by U nk = U (xk, tn) which is the corresponding finite difference approximation of u at the point (xk, tn). Let and β j = τ/h2j . Then (15) may be written as: ( I − θ d ∑ j=1 (β jC∆2j ) ) U n+1k = ( I + (1 − θ) d ∑ j=1 β jC∆2 j ) U nk + τF (16) Applying initial and boundary conditions the operator ∑dj=1 β j D(xk, t,Uk)∆ 2 j can be decom- posed. Let N be the cardinality of I and define F : [0, ∞) × SN−→Y by F (t,U )k = F(xk, t,Uk). Also let L : [0, ∞) × SN−→Y and Z : [0, ∞) × SN−→Y be two mappings such that L is linear and [L(t,U )u + Z(t,U )]k = { d ∑ j=1 [ β jC∆2j ] u } k , (17) where t ≥ 0, U ∈ SN , and u ∈ Y . Using these notations the difference scheme (16) can be written in the following simple form: (I − θL)U n+1 = (I + (1 − θ)L)U n + Z + τF . (18) Where L, Z and F are evaluated at (tn,U n). When θ = 0 in (18) we get U n+1 = (I + L)U n + Z + τF . Since L, Z, and F are evaluated at nth time level, U -values at (n+1)th time level can be evaluated explicitly by above equation. These types of finite difference schemes are called LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... 6 Ruhuna Journal of Science 4, pp. 1–12, (2009) fully explicit schemes. When θ 6= 0 in (18) U -values at (n + 1)th time level are expressed implicitly by U -value at nth time level. These types of finite difference schemes are called implicit schemes. When θ = 1/2 the scheme is called the Crank-Nicolsion finite difference scheme. The case θ = 1 is called semi-implicit scheme. This semi-implicit scheme is used in numerical simulations of this paper. In the following example the finite difference scheme is implemented on a domain in two dimensional space. EXAMPLE 1. Consider two dimensional case with Ω = (a1, b1) × (a2, b2). Let I = {(i, j) ; (i = 1, 2, ..., N1), ( j = 1, 2, ...N2)}, where N1 and N2 are such that (N1 + 1)h1 = b1 − a1, (N2 + 1)h2 = b2 − a2 and usual second order approximation for the Laplacian operator. Then ∆22 =      A A . . . A      N2×N2 +Z1 and ∆21 =        A1 B1 B1 A1 B1 . . . . . . . . . B1 A1 B1 B1 A1        N2×N2 +Z2; where A =        −2I I I −2I I . . . . . . . . . I −2I I I −2I        N1×N1 , A1 =        −2I −2I . . . −2I −2I        N1×N1 , B1 =        I I . . . I I        N1×N1 ; LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... Ruhuna Journal of Science 4, pp. 1–12, (2009) 7 here I is an M × M identity matrix. Z1 = (Z1,1,Z1,2, ...,Z1,N2 ) T and Z2 = (Z2,1,Z2,2, ...,Z2,N2 ) T . In the aboveZ1i = (1, 0, ..., 0, 1) T 1×N1 for i = 1, 2, ..., N2 Z2,i = { (1, 1, ..., 1, 1)T1×N1 if i = 1, N2 (0, 0, ..., 0, 0)T1×N1 otherwise. ¿From the equation (17) we get Lu + Z = (β1D∆21 + β2D∆ 2 2)u That is we get L(t,U ) =    L11 . . . L1N1 ... ... LN11 . . . LN1N1    N2×N2 where Li j =    β1DA1 + β1DB1 + β2DA if i = j β1DB1 if i = j ± 1 0 otherwise Also Z = β1DZ1 + β2DZ2. 4. Numerical Experiments In this section we consider pattern formation of diffusion coupled Gray-Scott model. The Gray-Scott model includes the following two irreversible reactions: U + 2V −→ 3V V −→ P where U and V are two reacting specimens and P an inert precipitate. The mathematical model for this reaction-diffusion process is of the form: ∂u ∂t = d1∆u + u2v − (η + ζ)u ∂v ∂t = d2∆v − u2v + η(1 − v).      (19) where u and v are the concentrations of U and V respectively and d1 and d2 are their respec- tive diffusion coefficients. η and ζ are dimensionless feed rates of first and second reaction respectively (Webpage ????a). LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... 8 Ruhuna Journal of Science 4, pp. 1–12, (2009) Figure 1 Time evolution of V : Surface plots of V at different time levels 4.1. Gray-Scott model with coupled diffusion terms We consider the following diffusion-coupled Gray-Scott model: ∂u ∂t = d1∆u + u2v − (η + ζ)u ∂v ∂t = d∆u + d2∆v − u2v + η(1 − v)      (20) LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... Ruhuna Journal of Science 4, pp. 1–12, (2009) 9 Figure 2 Time evolution of V : Surface plots of V at different time levels Here η, ζ, d, d1 and d2 are constants and d1 6= d2. We consider this reaction-diffusion system on the bounded domain [0, 1] × [0, 1] under no flux boundary conditions and under the parameters d1 = 1.0 × 10 −7, d2 = 8 × 10 −6, d = 2.5 × 10−6, ζ = 0.005, η = 0.0006. In this case initial conditions are: u(x, y, 0) = 0.101215; (x, y) ∈ (0, 1) × (0, 1) v(x, y, 0) = v0 + (rand(200)/100000.0)v0; (x, y) ∈ [0, 1] × [0, 1]; } (21) where v0 = 0.055328. LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... 10 Ruhuna Journal of Science 4, pp. 1–12, (2009) Figure 3 Time evolution of V : Density plots of V at different time levels It can be shown that above initial state and parameters satisfy the Turing instability condi- tions which are the conditions should satisfy by reaction diffusion systems in order to form spatial patterns(Murray 2003, Turing 1952). Numerical solutions of the reaction-diffusion system (20) subject to no-flux boundary conditions under initial data (21) are obtained using above introduced semi-implicit finite difference scheme. Surface plots of the v-component at different time levels of those numerical solutions are shown in Figures 1 and 2 (V denotes the numerical solutions of v). Density plots of the same are shown in Figures 3 and 4. LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... Ruhuna Journal of Science 4, pp. 1–12, (2009) 11 Figure 4 Time evolution of V : Density plots of V at different time levels According to these density plots it can be seen that these solutions form some spatial patterns. 5. Discussion The constructed implicit Finite Difference scheme can be directly applied to solve the trans- formed (system with coupled diffusion terms to system without coupled diffusion terms) reaction diffusion system. After that, the solutions of the coupled-reaction diffusion system are obtained by applying inverse of the transformation. LW Somathilake: Numerical Solutions of Reaction-Diffusion systems ... 12 Ruhuna Journal of Science 4, pp. 1–12, (2009) However when transforming a diffusion-coupled reaction diffusion system to a reaction- diffusion system with uncoupled diffusion terms the reaction terms become more compli- cated. This may affect to convergence of the finite difference scheme. In order to get rid of this problem step sizes of time have to be shorter. Again this may cause to increase computational errors and computational time. It is expected to estimate and compare com- putational errors and computational time of these two methods in my future work. References Badraoui, S. 2002. Existence of global solution for systems of reaction-diffusion equations on unbounded domain. Electronic Journal of Differential Equations 2002(2002) 1–10. Grindrod, P. 1991. Patterns and Waves- The theory and applications of Reaction-Diffusion equa- tions. Oxford Applied Mathematics and computing Science Series, Oxford. Hoff, David. 1978. Stability and convergence of finite difference methods for systems of nonlinear reaction-diffusioon equations. SIAM J. NUMER. ANAL. 15 1161–1177. Kirane, M. 1(1993). Global pointwise a priori bounds and large time behavior for a nonlinear system describing the spread of infectious disease. Applicationes Mathematicae 1–9. Kirane, M. 1989. Global bounds and asymptotics for a system of reaction-diffusion equations. Journal of Mathematical Analysis and Applications 138 328–342. Murray, J.D. 2003. Mathematical Biology: Spatial models and biomedical applications, vol. II. Springer-Verlag Berlin Heidelberg. P. Collet, J. Xin. 1994. Global existence and large time asymptotic bounds of l∞ solutions of thermal diffusive combustion systems on Rn. arXiv:chao-dyn/9412003 5. Seaid, M. 2001. Implicit-explicit approach for coupled systems of nonlineat reaction-diffusion equa- tions in pattern formation. Science Letters 3. Turing, A. M. 1952. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society of London 237 37–72. Turk, G. 1992. Texturing surfaces using reaction-diffusion. Phd theses, University of North Carolina, Chapel Hill. Webpage. ????a. http://www-swiss.ai.mit.edu/projects/amorphous/grayscott/. Webpage. ????b. http://www.scholarpedia.org/article/fitzhugh-nagumo−model.