HI": W‘Wld‘flfifilfl‘l Wmums!» ‘Hwé u CID—\(fl ZZZ". C This is to certify that the dissertation entitled Multigrid Methods for Solving Reaction-Diffusion Systems presented by Hsiu-Chuan Wei has been accepted towards fulfillment of the requirements for Ph .0. degree in W \ Major professor Date August 15 , 2000 MS U is an Affirmative Action/Squad! Opportunity Institution 0- 12771 PMCE IN RETURN Box to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE odfitézéoé 6/01 cJCIRC/DateDuopGS-pJS Multigrid Methods for Solving Reaction-Diffusion Systems By Hsiu—Chuan Wei A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics August 15, 2000 ABSTRACT Multigrid methods for solving reaction-diffusion systems By Hsiu—Chuan Wei A Multigrid method as an iterative method is pr0posed to solve a reaction-diffusion system. The model we are concerned in this thesis is a system of parabolic partial differential questions for two chemical species. A fully implicit finite difference scheme is used to discretize the differential equations. A V—cycle scheme with one smoothing step per grid is then applied to the Helmholtz equations [54] arising at each time step. The convergence of the V-cycle scheme for solving the linear system at each time step is obtained. The stability and convergence of the fully implicit scheme along with V-cycle as the iterative solver are proved. Numerical results for a reaction-diffusion system are presented. ACKNOWLEDGMENTS I would like to express my sincere appreciation to Dr. ChiChia Chiu, my thesis advisor, for her kind guidance during each phase of this project. Our many and varies discussions are vital to my development. With each obstacle encountered, she offered me countless suggestions that help me overcome each hurdle and move forward to succeed. I also wish to thank Dr. Chang Yi Wang, Dr. Charles MacCluer, Dr. Tien-Yien Li, Dr. Baisheng Yan and Dr. Zhengfang Zhou on my committee for their kind guidance and invaluable suggestions. Without their expertise, I could not have finished this project. iii TABLE OF CONTENTS LIST OF FIGURES LIST OF TABLES INTRODUCTION 1 Fully implicit discretization 2 Two-grid algorithm and its convergence 3 V-cycle and its convergence 4 Stability and Convergence 5 Experimental Results 5.1 The Thomas System ............................. 5.2 The Schnakenberg System .......................... 5.3 Comparison .................................. BIBLIOGRAPHY iv vi 13 25 34 43 44 49 52 56 2.1 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 LIST OF FIGURES Grids for 9k and Qk_1 ............................. Possible patterns estimated by linear stability analysis. .......... Turing’s patterns obtained by numerical simulation ............. Numerical simulation of Turing pattern of ( 5.1) and residual plot with '7 = 8. ................................... Numerical simulation of Turing pattern of ( 5.1) and residual plot with '7 = 48 .................................... Numerical simulation of Turing pattern of ( 5.1) and residual plot on a cone surface ................................. Numerical simulation of Turing pattern of ( 5.3) and residual plot of with 7 = 15 .................................... Numerical simulation of Turing pattern of ( 5.3) and residual plot with 'y = 30 .................................... Numerical simulation of Turing pattern of ( 5.3) and residual plot with 7 = 50 .................................... Numerical simulation of Turing pattern on a square ............. 14 47 48 48 50 51 51 53 5.1 5.2 5.3 5.4 LIST OF TABLES Comparison between single grid and multigrid with h = 0.05 and At = 0.01. 53 Comparison between single grid and multigrid with h = 0.025 and At = 0.0025 .................................... 53 Comparison between single grid and multigrid with h = 0.025 and At = 0.005. ................................... 53 Comparison between single grid and multigrid with h = 0.0125 and At = 0.0025 .................................... 54 vi Introduction Reaction-diffusion systems are systems of parabolic partial differential equations. They have been used successfully to model chemical and biological processes which involve pattern formation in morphogenesis. Since all models for spatial pattern generation are necessarily nonlinear, analytical solutions are generally not available. Numerical solutions of such equations are thus important for computer simulation of patterns and data analysis. In this thesis, general two-dimensional reaction diffusion systems are considered. A multigrid method is pr0posed for solving the system. Analysis and applications will be discussed. Fascinating spatial patterns, such as animal coat patterns and butterfly wing pat- terns, are exhibited on living creatures. The development of pattern and form in embryology is known as morphogenesis. It begins from a more or less homogeneous egg. The final complexity of pattern and form are generated during its development. How the patterns are established is still unknown. Although the process of the devel- Opment is genetically determined, the genes themselves cannot generate the pattern since the genetic material in most cells of an organism can be assumed the same [39]. Much research has tried to determine the mechanisms that can explain the gener- ation of patterns developed from less structured tissues. Therefore any mechanism assumed to explain this development must be capable of generating spatial patterns from almost homogeneous initial conditions. Several mechanisms have been proposed as possible pattern formation mechanisms in morphogenetic situations. Among these, the most powerful mechanism for forming coat patterns is a reaction diffusion mecha- nism proposed by Turing [3]. A. M. Turing first suggested in 1952 that some patterns that occur in chemistry result from interaction between chemical reaction and diffu- sion. Since then a substantial amount of research has been done on this subject. See [41] for a survey of related developments. According to Turing’s chemical prepattern approach, the first step of development in the morphogenetic process is the creation of a morphogen concentration spatial pattern. The underlying prepattern is believed to be laid down in the very early stages of the embryogenesis. Consider the case of the zebra. The gestation is about 360 days with the prepattern laid down about 3-5 weeks. After the prepattern is once established, the cells differentiate accordingly to produce melanin to reflect the spatial pattern of morphogen concentration. For example, the pattern of the giraffe could be generated when the morphogen concentration level is greater than some threshold level. When applying Thring’s theory of morphogenesis, an important step is to identify the chemical elements for the morphogens. For example, calcium has been identified as a morphogen for hair initiation in Acetabularia. In this thesis, we will assume that the pattern is generated when one of the morphogen concentrations is above some threshold level. Now let us consider the following two-dimensional reaction-diffusion system: at = DuAu+f(u,v) (0.1) vt = DvA’U +g(u,v), 92 _ a -0 an an _ an an — ’ where u(:1:,0) and v(a:,0) are given. Here we take the rectangular domain : Q = [0, L1] x [0, L2] C R2, 09 is the boundary of Q, and u and v are considered as two morphogen concentrations; f (u,v) and 2 g(u,v) are reaction kinetics which describe the interrelation between morphogens. The diffusion terms, DuAu and DvAv, reflect that each molecule can move around randomly with diffusivities Du > 0 and D” > 0. Reaction-diffusion mechanisms suggest that spatial patterns are evolved from diffusion-driven instability. Turing also demonstrated that, under certain conditions, two interactive morphogens could form a stable inhomogeneous concentration pat- tern. He suggested that without diffusion, u and v tend to a uniform steady state — no pattern will form. In the presence of diffusion, the system is unstable to small disturbances but the instability will be bounded by nonlinear reaction terms. Thus the inhomogeneous steady state is obtained and spatial patterns are generated [52]. All that is required for the creation of pattern is some sort of nonlinear activator- inhibitor mechanism [26]. The analysis of whether or not the systems are able to generate spatial patterns and how the pattern and mode is selected can be found in [41]. To formulate the mathematical model, we need initial and boundary conditions. Since the formation of the underlying prepattern could arise from instability to small perturbations in the featureless tissue, a natural way to simulate this situation in numerical computations is by introducing a small random perturbation about the uniform steady state. In the system given above, we assume initial conditions are given. As for the choice of boundary conditions, we take zero flux boundary conditions because it implies no external forces. Under certain conditions, the reaction-diffusion system described above has a unique solution. The proof of the existence and uniqueness of the solution can be found in [11]. The system is useful for modeling patterns in chemistry and biology. A few examples can be found in [13]. In these examples, the concentration of nutrient and the concentration of buffer satisfy a reaction-diffusion system as given in (0.1). Since the reaction term, f and g, in (0.1) is nonlinear, one can not hope in general 3 to obtain analytical solutions. Numerical solutions with preliminary linear analyses are required in studying the patterns generated by a given reaction-difl'usion system. Attempts have been made to demonstrate the patterns formed by such mechanisms and compare them with the mammalian coat patterns [3] [41]. Since the spatial patterns are generated by diffusion-driven instability, the system has to be unstable. Because of this nature, highly stable numerical methods are necessary for computer simulations of these patterns to ensure that the patterns obtained by computer simula- tions are formed by the instability of the original system, not by numerical instability. One approach is an ADI (alternating direction implicit) type of scheme [13] [14]. Other numerical methods that have been used to solve reaction-diffusion systems are finite element techniques, monotone iterative methods and explicit finite difference schemes [21] [27] [40] [46] [56]. Finite difference schemes have been used due to their simplicity. Explicit finite difl'erence schemes are especially easy to implement but the conditional stability forces small step size in time in numerical experiments. Consequently, the explicit schemes require long runtime. Another intuitive approach is using a fully implicit scheme which is unconditionally stable. The discrete problem arising from a fully implicit finite difference scheme is a large sparse system of linear equations. Existing methods for solving linear systems can be grouped into two branches: direct methods and iterative methods. Although direct methods produce exact solutions in a finite number of steps (regardless of roundoff errors), they can not be applied here because the complexity needed to invert the matrix is of cubic order, and because the inverse of a sparse matrix may not be sparse. The computation time and computer memory for a solution obtained by direct methods are much too great. We therefore consider using iterative methods. Our concern about the numerical scheme now turns to its convergence rate in practice. In order to obtain the inhomogeneous solutions, the computations have to be carried on until equilibrium is reached. Even though the unconditional stability of fully implicit finite difference schemes enable us to take a much greater step size in t, it still requires thousands of steps in t to reach the steady state. Thus in these numerical computations, we need to solve a pair of large linear systems at each time step for at least thousands of time steps. Therefore, the running time mainly depends on the efficiency of the iterative method used for solving the linear system at each time step. Unfortunately, the classic iterative methods, namely Jacobi, weighted Jacobi and Gauss-Seidel methods, have unsatisfactory convergence rates. They have been known to be able to remove the high frequency modes in an error efficiently, but are unable to damp the low frequency modes. Finding a fast solver for a large linear system becomes an important issue. In this thesis, a multigrid method is proposed for solving this problem. Consider that the smooth components which are efficiently approximated on coarse grids and the oscilatory components are fast to converge on fine grids. Multigrid techniques could eliminate all frequency components. Multigrid methods as iterative methods have been known to be a fast solver for linear systems arising from the discretization of partial differential boundary-value problems [4] [35]. Numerous works about numerical experiments and theoretical understanding of the convergence prOperties of these methods have appeared in the past three decades [2] [4]-[8] [20] [29]-[38]. Reported numerical experiments suggest that these methods are very efficient for a wide range of practical problems [10] [54]. There are many convergence proofs to multigrid algorithms. One approach to these methods is local Fourier analysis [22] [28] [54]. It is not generally rigorous. Local Fourier analysis gives realistic quantitative results on the convergence behavior yet assumes an unbounded domain in space. Thus it can be regarded as an analysis only for problems with periodic boundary conditions. With Neumann boundary condition in our model problem (0.1), we can view the multigrid as a single operator and study the norm of the operator. Some of the proofs have been given in [16]-[18] yet these proofs require a sufficiently large number of relaxation sweeps. In practice, multigrid is used only with one or a few relaxation sweeps. In our numerical computations (Chapter 5), only one presmoothing and one postsmoothing are needed in almost every example. Other proofs are based on the approach for variational formulation. Intensive re- search into the convergence of multigrid methods for variational problems can be seen in [4] [7]-[9] [29]-[38]. These proofs usually require a ‘regularity and approximation’ assumption. Convergence rates have been guaranteed with any amount of smoothing for solving linear systems An: = b, with A symmetric positive definite. The linear system arising from discretization of (0.1) is unfortunately nonsymmetric. Then re- defining an inner product and its corresponding norm by the midpoint rule [13] is proposed to formulate a variational-like problem and construct the proof of the V- cycle algorithm. Variational framework is a natural formulation when finite element discretization is used. In this thesis, since finite difference discretization is used and a new inner product is defined, more attention needs to be paid to the constructions of the transformations between grids. For the choice of the smoother, a red-black Gauss-Seidel method is applied because it has been used extensively as a smoother in multigrid methods [54] [60]. The convergence of this smoother under the norm used in this thesis is discussed in Chapter 3. The rigorous proof of the convergence of a V-cycle algorithm is also given with any amount of smoothing. This thesis focuses on the construction of the multigrid algorithm that can be ap- plied to our model problem and the proof of the convergence of the algorithm followed by the demonstration of the efficiency of the algorithm. In the numerical computa- tions, two reaction systems have been used to demonstrate the V-cycle algorithm. They are able to generate spots and stripes on rectangular domains. We should note that even the analysis has been discussed on a general rectangular domain due to its simplicity. We believe the V-cycle algorithm constructed in this thesis is also valid for domains in other shapes. Thus the pattern on the surface of a cone has been simulated in numerical experiments. Finally, we compare the running time of the V-cycle algorithm with that of the computations carried out only on the fine grid to demonstrate the high efficiency of the former. The outline of this thesis is as follows. In Chapter 1, the finite difference scheme is constructed. In Chapter 2, the two-grid algorithm is developed and its convergence is proved. In Chapter 3, we extend the two-grid algorithm to a V-cycle algorithm and its convergence is obtained. In Chapter 4, we discuss the stability and convergence of the fully implicit finite difference scheme along with the V-cycle algorithm described in Chapter 1 and 3. Finally in Chapter 5, the experimental results for the numerical solutions to two reaction-diffusion systems are presented. CHAPTER 1 Fully implicit discretization In this section we begin our study of the fully implicit finite difference dis- cretization of the reaction-diffusion system (0.1). For simplicity, consider a uni- form rectangular mesh on a rectangular domain 9 = [0, L1] x [0, L2] with mesh size h = Ll/Nl = L2/N2. Let (33,-,yj) E Q be a grid point then any-+1 = as, + h and 3/,“ = y]- + h. The increment in t, t"+1 — t", will be denoted by At, and we adopt the standard notation ufj z u(:r,-, yj, t"). The basic idea of the finite difference discretization is to approximate the deriva- tives in a differential equation by the difference quotients. For example, n+1 n —($iayj9t ) — 3t At can be seen from the derivative formula Bu u(:r, y, t + e) — u(a:, y, t) — :c, ,t = lim . 8t( 3/ ) HO 5 Similarly, the second order derivative %(z,, y], t") can be approximated as follows: 8% it ~ 71 _ n n 791:2 (313,5,yj,t )— ui+1j 21% + ui—l j. Here we propose to use a fully implicit scheme. We approximate (0.1) by the backward-time central-space scheme [53]. We have u".+l—u?. Arh 11" A u? —J——LAt : Du( h 211— + 14,12 __1__+1)+ f(un U2], Uznj) (11) v?.+1—v?. Arhv‘zH Ayh’vu :Jn+l ’Un . At : DU( )3 + 1,2” ) + 9(un uijav zj’) where A1,, and Ayh are the centered second order difference operators such that (Ax), + Ayh)21?j = (21?+1]— 221"]-z + 21, 1j)+(21[3+1—221,'-;+ 21,]; 1). The zero flux boundary condition is approximated by ”fix/2+1 = 213,24, 21:24) = 21,21, for all 2,71, 213,11,”- = 217,114,, 24:1”. = 21?], for all j,n The boundary treatment for 21 can be defined in the same way. Let 21" =[21? j] and 21" =g[v ] be the vectors obtained by the usual ordering. Then (1.1) can be written in the following matrix forms: (1 + 711w“ = 21" + At f" (1.2) (I + T'A)2)"+l = 21" + Atg", (1.3) where T = DuAt/h2,T' = DvAt/hz, A is the standard matrix resulting from the discretization operator —-Axh — Ayh, f" = [f(21§3,v{‘j)], and g" = [g(21,-j,v,-’; )]. Let A = I + TA and A’ = I + T’A. Then (1.2) and (1.3) become A21"+1 = 21"+Atf" (1.4) A'v"+1 = 21"+Atg". (1.5) A can be represented by the matrix F W —2TI —7'I W —TI —TI W —TI —TI W —TI —2TI W L. with I the (N1 + 1) x (N1 + 1) identity matrix and )- -1 1+ 47' —2T —7' 1+ 47' —T —7' 1+4T —1' —'r 1+47' —T —2T 1 + 47' J a (N1 +1) x (N1 +1) matrix. Note that the operators A and A’ are five-point stencils. The following lemma establishes the consistency of this fully implicit finite differ- ence discretization. Lemma 1.1 Let 21(23,, 3),, tn“) and 27(1ri, yj, tn“) be the exact solutions of (0.1). As- sume that 21 and 27 are of class C4. If the partial derivatives fu, fv, gu, and 9,, are continuous and uniformly bounded, then ~ttv+l_~n. ~n+l u‘ At"1 : Du(AIhU?1+1+‘Aflyhfi—W)+f(uij,@$)+rn'+l fifth-n. A "1+ A h-"tl 111751 = 1911(4",,321--+”—h‘L—)+g('u,-,-,2?'?,)+8’-‘,+1 10 where 0(At + h?) if ij is an interior index, r11.“ 31!.“ -_— 23 ’ 1.7 0(At + h) if ij is a boundary index. Proof: Use the Taylor series expansion {1;} = 22".“ —At(u)"-+1+O(At2) to obtain 713'“ _ {fl}. +1 1] _ ~ n . . T — (ut),j + 0(At) for all z] . Consider the Taylor series expansions h2 h3 n+1 + (axxx)%+l+0(h4) 71?:111' = 11:3“ + h(aat)?j+l + 3(firxlz‘j "'6— and ~ 1 ~ 1 ~ 1 hz ~ 1 ha 1 4 uff1j= 113+ -h(ux)?j+ + 3w”)? — Emma)? + 0(h ). Add the equations above and use the approximation for the zero flux boundary con- dition to obtain (ii-an)”+1 + 0(h2) ifi is an interior index, Axhflg+1 : (u$1)”*ln + 0(h) ifi is a boundary index. Similar expansions also hold for Ayhun+1 as well as the other concentratlon func- tion 17. Also, f(fl?j+1’fiilj+1) : f(flij’ fiij) + fu( qu whose stencil is l1 1—6 2 and 17; 0 0 4 O 1 2 4 0 q d for an interior point, for a boundary point with index j = 0, for the corner point with index i, j = 0. This restriction takes fine grid vectors and produces the coarse grid vectors according to the rule Pf‘lu = v for u E 0k and u E Qk_1, where 1 vij = filu2i+u2j+l + U21+1,2j-1 + u2i-l,2j+1 + U21—1,2j—1 + 2(u2i,2j+1 + ”21,2341 + U21+1,2j + u?i—1,2j)+ 4U2i,2j]a 1 S i S (mlk—l — 1) and 1 3.7 S (m2k—1_1)a 1 003' = I6[2u1,2j+1 + 2Ul,2j—1 + 2u0,2j+1 + 2”0,2341 + 4U1,2j + 4710,21], 14 1S j S (m2k—1 _1)1 and 1 ’Uoo = 1—6[4u00 + 4u01 + 4u10 + 4u11]. The values um, um,,_,,j, u,,m2k_,, ”mu—1,0, vomqu, and ”m1k_1,m21-1 are defined analo- gously. This operator is represented by the (mug/2 + 1)(m2k/2 + 1) x (m1,c + 1)(m2,c + 1) matrix _ l 21",, 21",c PIC-1 : i ik 2i]: A; k 16 , 2i",c 21",c J where lk is the (mm/2 +1) x (m1,c +1) matrix 2 2 ~ 1 2 1 1,, — 2 2 For the prolongation, we use the bilinear prolongation QLI : Qk_1 —> 9k. Its 15 stencil notation is given by and 1 2 1 i 2 4 2 for an interior point, 1 2 1 - r 1 2 1 1 . . . . _ :4- 2 4 2 for a boundary pomt w1th Index ] — 0, j 0 O 0 0 2 1 l i O 4 2 for the corner point with index i, j = 0. 0 0 0 If u E (2;, and u E Qk_1, then Qle = u with the components of u given by U2z',2j U21+1,2j u2i,2j+1 and u2i+l,2j+1 vij, 0 S Z S mlk—l and 0 S j S m2lc—la 1 . . 5(7111' + Ui+1,j)a 0 S 3 S 7nlk—1_1 and 0 S J S m2k—1, 1 . . 5(013' + vi,j+1)a 0 S 2 S mlk—l and 0 S J S m2k—1_11 1 . 10% + ’Uz‘+1,j + vi,j+1+ vi+l,j+l)1 0 S 2 S mlk—l — 1 and Ogjgm2k_1-1. 16 The operator is represented by the (muc + 1) (mg;c +1) x (Talk/2 +1)(m2k/2 +1) matrix 1- -| 2.];c J], J], 2i,c Q2—1=i jk jk 1 Jk Jk 2.7,c J where L, is the (m1;c + 1) x (mm/2 + 1) matrix - : 2 1 1 2 jk = 1 1 1 l . 2 .1 From the matrix forms of I), and Pf”, we see that the rank of Pf'l is (m1),c /2 + 1)(m2k/2 +1) and so is that of QLI. Next define the coarse grid operator Ak_1: Qk_1 ——> Qk_1, by Ak_1 = Pf—lAing—i- 17 If 5k and then Ak__1 has the same form as that of A, with dimension the number of points on - Gk 5k 2T;c Sk 218k ale .6}: 2’71: fit 71: Sic ilk 0k ’71: 3k Tk T, 5,, 2T k file file Gk Zflk 7k ‘71: file 271: 18 Sk 3k 011: ’71: a. . the grid Qk_1 and with components 01H = (36m, + 96/3,c + 6470/64, ,8k_1 = (60k + 32fik + 327k)/64, 7k—1 2 (Ok '1' 8flk + 167k)/64. Recall that the fine grid Operator A J = A is five-point. The transformation leads to a nine-point stencil for A J_.1. Fortunately, the coarser grids continue to have the nine-point stencils as shown above. This increases the computation cost somewhat. However the choices of the operators defined previously lead to nice properties by which the convergence of the multigrid process is guaranteed. Now A in (1.4) is symmetric and nonnegative definite with respect to the discrete L2 inner product (-, .)J on L2(QJ) determined by the midpoint rule [12] k,l=N (u9v)J : [2371”? Z uklvkl’ k,l=0 where 7,- : 1 / 2 if i = O or i = N, and 7,- : 1 otherwise, and h J = h. The correspond- ing norm is H u “,2 ‘/(u,u)J. We also define the inner product on each (2,, in the same way. The matrix A = I + TA is SPD (symmetric positive definite) with respect to (-, )J. Furthermore we have (Pf—lu, v)k_1 = (u,Q',§_.1v),c for every u E Q], and v E Qk_1. It will be shown in the following lemma that each A, is also SPD with respect to (., ~)k. Lemma 2.1 If Ak is SPD with respect to (~, -)k, then Ak_1 is SPD with respect to (-,-)k_1. Therefore (A1,_1)‘l exists. Proof: Let x E Qk_1 and x aé 0. Then (‘4k—133a-le—l '—' (Pf—lAkQLflJfi—l = (Asz—lxa Q:_1$)k ( positive ) : (Qt—lxiAkQIt—lxlk 19 = ($2 Pt—lAsz—lxht—l = (x, Ak_1x)k_1. (symmetric) Now Ak_1 is SPD with respect to (-, -)k_1. Consider the linear system Ak_1x = 0 with x 6 91,4. Note (Ak_1x,x)k_1 = (0,x)k_1 = 0. But (Ak_1x,x)k_1 > 0 unless x = 0. But x = O is the only solution to Ak_1x = 0. Thus (Ak_1)‘1 exists. Cl Since AJ is SPD, from Lemma 2.1, Ah is SPD for 1 g k 3 J and (A,,)‘1 exists. Now we can define an inner product by A(u,u) = (Aku, u)k, with corresponding norm |||u||| :2 (/A(u, u) for u, u E 5%. We recall the following well known properties of real symmetric matrices. The proof can be obtained by following the same procedure as for the standard 2-norm [1]- Lemma 2.2 (a) Let B be a (m1k+1)(m2k+1) x (m1k+1)(m2k+1) matrix, symmetric with respect to (~, -)k. Then II B “k: p(B), where p(B) is the spectral radius of B. (b) Let C be a (m1;c +1)(m2;c +1) x (m1,C +1)(m2;c +1) matrix, symmetric with respect to A(-, ). Then |||C|H = p(C). For the purpose of analysis, we define an auxiliary operator Sf‘l : (2,, —> Qk_1 by Sf'l = (Ak_1)”lP(°-1Ak. So, we have 21,431;—1 = Pf‘lAk. Next, we shall construct a two-grid scheme for the linear system in (1.4). Let D J be the diagonal matrix which consists of the diagonal elements of A = A J, f} = u" + Atf", and x"+1 be the exact solution of (1.4). So (1.4) becomes AJx"+1 = f} (2-1) The Gauss-Jacobi method can be expressed as follows. Let A J = B J + D J, DJUMI) = _BJu(i) + f?, 20 Dyna“) = (DJ—AJ)U(i)+fy, u = (DJ)_1(DJ—Aj)u(i)+(DJl—lf91 21““) = u“) — (DJ)"AJu“’ + (DJ)"f5‘, um'l) = u(i)+(DJ)—l(fy—AJu(i)). The two-grid algorithm for solving (2.1) is given below. Algorithm 2.1 Step]. W = u"+(DJ)-1(fy—Aju"). Step2. um = u(0)+Qj_1q, AJ_1q = Pf‘1(f}‘ — AJu‘”). Step3. u"+1 um + (DJ)_1(f}' — AJU(1)). Where Step1 is a pre—relaxation sweep using the Gauss—Jacobi iteration on the fine grid. Step2 is coarse grid correction which solves the coarse grid problem exactly. Step3 is a post-relaxation sweep on the fine grid. To establish the convergence of the two-grid algorithm, we need the following lemma. Lemma 2.3 (a)Sk+1Q',:+l——— (b) QLISS’I is symmetric with respect to A(-,-), i.e. for u,v E 0k, A(Qk— 15 41” v): A(u,Q’,:_ISf’1v). (c) Q’;_15,f_1 is idempotent, i.e. (Q’,§_IS,':"1)2 =Qk_ 18k 1, so its eigenvalues are 0and1. Proof: (a) The result follows directly from the definitions of Ak and S}: +1. Sk+1Qi§+1=(Akl’1Pf+1Ak+1Qi+l = (AU—1M0 = 1k- 21 (b) Let u, v E Qk, we have A(QLle'lvw) = (AkQ:_IS,’:_1v,u)k = (Sf‘lv,P,:‘—1Aku)k_1 = (Sf’lv, Ak_lS,’§_1u)k_1 = (Ak_1Sf_1v,S,’:—1u)k_1 = (Pf—lAkv,S,f_1u)k_1 = (Akvaz—ISf-lu)k = A(v,Ql§_13f‘1U)- (c) It follows from part (a) that ((2524354)? = Q£_le—1Q’,§_ISL°_1 = (Q,’§_11k_lS,’§‘1) = ( £4554). An idempotent matrix is similar to a diagonal matrix of the form diag(1, ..... , 1,0, ..... 0) [59]. CI Let x"+1 be the exact solution of the linear system (2.1) and an“ be the numerical solution obtained by the two-grid algorithm described above. Let the initial error be e0 2 x"+1 — u" and the error after the two-grid scheme be el = x"+1 — un“. Theorem 2.1 The error el 2 K80, and lllellll _<_ |||K||| . |||eo||| where K is the error operator with |||K||| < 1. The two-grid scheme is convergent. Proof: Since 61 = x — u = 2:“ — — (Do-1v: — Aw‘”) = x"+1— um — (D1)“1AJ(x"+1 —— um) = (1 -(DJ)-1AJ)(z"+l — u‘”) = (I — (DJ)“AJ)(x”+‘ — — 625-1(AJ-1)-1Pj-‘AJ(x"+l — u‘°’)> = (I — (Db-Wm — Q5_1(AJ—1)"Pf“AJ)(x"“ — u“”) 22 = (I — (DJ)‘1AJ)(I — 625-.(AJ-1>-1AJ_ISJ-1)u§.“"+g., (Rs-Inf." = (R.>-‘u§.“”+g.—A.u§:‘”, i i— i—l U56) = u), 1)+Rk(gk—Aku§c )). The V-cycle algorithm is given below. Algorithm 3.1 If k = 1, solve Alul = g1 exactly. We may assume that Bl = (All—l- For k > 1, solve Aku;c = gk by the algorithm described below with g = fJand Bk is defined in terms of Bk_1. Step1. ulcl) u?) + Rk(gk — AWE)”, with 115,“) = 0 if]: < J and 111°): u" z‘fk = J. Step2. ”53) = ”is” 'l' Q§_1Bk_1gk-1, where g.-. = Pt—‘(gk — Anti"). 25 Step3. Bkg,c = u)? + Rk(g,c — Aku(2)). For the smoother Rk, it is sufficient to choose a convergent iterative method. Here we will use the red-black Gauss-Seidel method because of its good smoothing rate for elliptic problems [54] [60]. Let ((Rk)‘1)‘ be the adjoint of (1‘2,,)‘1 with respect to (-, -)k. It can be verified that ((Rk)‘1)t = Ak + Dk — (RU-1, where D, is the diagonal matrix which consists of the diagonal elements of Ak. In the following lemma, we will show that (R,,)’1 + ((Rk)‘1)‘ — Ak = D, is positive definite. Thus red-black Gauss-Seidel method is convergent with respect to I” ' |||. Lemma 3.1 (a) The diagonal elements of Ak are positive for 1 g k S J, so D1, is positive definite. (5) ”III: - RkAklll <1 - Proof: (a) Since A1, is SPD, (Akeh e), is positive and so is the ith diagonal element of Ak for 0 g i g mk, where e,- is the vector with ‘1’ on the ith position and ‘0’ at all other entries. Thus D, is also positive definite. (b) For u, v E Ok, we need first to show that A(RkAku,v) = A(u, (Rk)‘Akv) and ((Rkltl—l =((Rk)'1)‘- NOW A(RkAku,v) = (RkAku,Akv);c = (Aku, (Rk)tAkv);c 1' A(u, (Rk)tAkU). But (Uwh = (Rk(Rk)‘1u,v)k = ((Rkl—lua(Rk)tv)k = (u,((Rk)’l)‘(Rk)‘v)k- Since ((Rk)“)‘(Rk)‘ = 1k, we have ((Rk)‘1)‘= ((Rk)‘)‘l- 26 Next we will prove the inequality in part (b). 0 S A((I;c — RkAk)u, (1k — RkAk)u) = A((Ik — (Rot/1km — RkAk)u, u) = A([Ik —— (Rut/4,. — RkAk + (Rk)‘AkRkAk]u, u) = A(v, u) — A([(Rk)‘Ak + RkAk — (Rk)‘AkRkAk]u, u) = A(v,u) — A([mky(rm-Inna,c + (Bowman-lam. — (Rk)‘AkRkAk]u,u) = A(U, U) — f1((1‘?Ic)tl(Rk)_1+((Rk)t)—1 — AkleAIcu, U) = A(Uau) - ((Rk)‘[(Rk)‘1+((Rk)‘)‘l — AkleAku,AkU)k = A(Uflt) — ([(RkY1 + (midi)—1 - AkleAku, Ric/11:10:: = A(u, u) — (DkRkAku, MAW». < A(u,u). Since DC is positive definite, (DkRkAku,RkAku)k > 0 if u is not zero. So |||I,c — RkAklll < 1 and the red-black Gauss-Seidel smoother is convergent for each k. [II For analysis purposes we will define the following Operators: 55 = Sf+15f121”°5§‘1, k_ k Ic+1 J-l PJ—Pk+1Pk+2”°PJ ’ and J J J—l k 1 Qk = QJ—IQJ—2 ' " k+ - Then we have the following lemma. Lemma 3.2 (a) (Pffvyw),c = (v,Q,{w)J ifv 6 OJ and w E 9k. (b) Ak -—- PfAJQi. 27 (C) 14ka = PfAJ. (d)35Q1{= Ik- Proof: To prove part (a), (vawh (P:+1P:—:21 "°Pf"1v,w)k = (131321 ' ° 'Pj—1ani+lw)k+1 = (v,Q.‘i_1QiI$'“Qi§+IW)J (’U, in)J For part (b), I: J _ k k+1 J—l J J—l k+1 PJAJQI: — Pk+1Pk+2 "'PJ AJQJ—lQJ—z ' " k k J—2 J—l k 1 : Pk+1"'PJ—1AJ—1QJ—2 ’ " k+ = PIf+1Ak+1QZ+l = Ak. For part (c), AkS§ = AkS£+1S£izl - - '55” = AkAzlPt+1AK+1$£ié - - ' 85‘1 = Pt+1Ak+IS£i%-~Sf“ = Pf+1PtfiAk+1-~Sf“ = lawn-vim = 105.4,. For part (d), sin = sis...---si:3s5-1Qi_.c25:;-~ = Sf+1°--Si:iIJ—1Q5:é~- Z“ 28 _ k J—2 J—l k+1 — Sk+1"’SJ—1QJ—2”' k _ k k+1 _ Sk+1°'°IJ-2"° k _ k Ic+1 — Ska [:1 Next we want to find the relation between the errors before and after the V-cycle. Let 2:], be the solution of Akxk = 9;, and uh be the numerical solution obtained at Step3 in the V-cycle algorithm. Then (I): — Bic/1k)“ = $1: — Bkgk = 3,, — uk = 33,, — v22) — RIC/luau, — u?) (2) = (Ik—RkAkXxk—u ) (1) = (1k - RkAk)(Ik — QLin—1Pf—1Akx33k — “k ) (1) = (1k — RkAkXI. — Qt-lBk_1Ak_ISt‘—‘)(xk — u. ) (0) = (1,. — Raoul. — Q£_.Bk-1Ak-ls,’:-‘)(Ik — RkAkxxk — uk ). (3.1) ForkS§ = I. — «2.185 +Qt<1k — R..A.)(I.. — Q2-.Bk_1A._ls,’:-1)(I. — RkAk)S’j. (3.4) Here we used (3.2) to rewrite (I,c — BkAk) . Next, using 5"} i = Ik, we have Similarly QiUk — RkAk) = Q: " QiRkAk = Q: — QiRkAkSSQi = (IJ — QiRkAkS§)Qi' (1k — RkAk)S§ = S"; — RkAkS'j = 5.]; — SSQiRkAkSS = 550'. — QngAksj). Rewrite the last term of (3.4) to obtain QiUk — RkAk)(Ik - Qi-lBk-lAk—ISIf—lxlk — 31.14055 (1.1 — QiRkAkS§)QI{(Ik — Qz—lBk—lAk—IS£_1)S§(IJ — QiRkAkSS) (IJ — QiRKAkSiWQiSi — Q,{_,Bk_1Ak_IS’j‘1)(IJ — QZRkAka). 30 Thus we have [J — QinAkS'JC = (I. — 155) +(IJ — Qi’RkAkSiXQiSIJC — Qi_1Bk—1Ak—15.’i_l)(1J — QiRkAkSil (3-5) Again, since S’jQfl = 1k, (IJ — QiS§)QiRkAkS§ = 0 and QlcleAkS§(IJ - Q1155) = 0. with (IJ — Qisile — QiRkAksi) = (IJ — Q1155) and (IJ - QiRkAkSSXIJ — Q1135) = (IJ - Qisil Thus (IJ — QiRkAksile — QiS§)(1J — QiRkAkS‘JC) = (IJ — Qisil (3-6) Then from (3.4), (3.5), and (3.6), [J — QinAksj = (IJ — QszAkSEXIJ — QLlBk—iAk—155—1XIJ — QiRkAkS’JC) Let TJ = IJ—RJAJ, T}, = (11— ,{RkAkS’fl forzgng—l, T1 = IJ—Qisi- Then I. — QgBkAks’j .—. m1. — Q,{_,B,,_,A,,_ls§-1)Tk, 2 g 1.: 3 J. (3.7) 31 Now we can write sun“ — an“ as a product of operators. From (3.3) and (3.7), 1"“ — an“ = (I. — RJAJ)(1. — Q5_,B._1AJ_lsj-1)(IJ — R,A.)(.;"+1— u") = TJ(IJ — Q5_1RJ_1AJ_ISj—1)Tj(xn+l — u") = TJTJ_,(IJ — Qj_,RJ_2AJ_Qsj-2)TJ_1TJ(x"+1 - u") = TJTJ—l "'T2T1Tz'"73—17301?"+1 — U") = K(:z:"+1 — u”), where K = TJTJ_1 ° ' ' T2T1T2 ° ° ' TJ_1TJ. Since TJ = IJ — RJAJ and since we have shown that IHTJHI < 1 in Lemma 3.1, the convergence will be obtained if we can show |||Tk||| = Hle — QiRkAkSflll S 1 . Lemma 3.3 (a) For 1 S k S J, QiSf is symmetric with respect to A(-, ). (b) The eigenvalues of QflSf are 0’s and 1’s, for 1 S k S J. (c) For u, v E 9k, A(u,v) = A(Qiu,Q,{v). ((1) ”Milk - RkA/clll <1. thenlllb - QngAksjlll S 1- Proof: The proofs of (a) and (b) and are similar to those of Lemma 2.3, parts (b) and (c). To prove part (c), let u, v E 9),. Then A(Qiani’U) = (QiuyAJinlJ = (U, PfAJQi’Ulk = (U. 14141)]: = A(u,v). For part ((1), let u 6 (L. Then 0 < A((IJ — QiRkAkSS)U, (IJ — QszAk55)u) = A(v, u) — A(v, QngAksju) — A(ankAksfiu, u) 32 +A(Q{RkAksju, QngAksy.) = A(u, u) — (AJU,QiRkAkS§U)J - (QiRkAiju, AJU)J +A(R,.A,,sf;u, RkAksju) = A(u, u) — (PfAJu, Rk.4kS§v)k — (RkAkS’ju, PjAJU)k +A(R,,A,,s§u, RkAksju) = A(u, u) — (AkS'ju, rah/1,352.), — (RkAkS'ju, Aij‘u)k +A(R,,A,,s§u, RkAkS’ju) = A(u, u) — A(S’jv, RkAkS’ju) - A(RkAksju, 5’52.) +A(R,,A,,s§u, RkAksfu) = A(u, v) — A(Sfu, S’ju) + A(Sfu, Sfu) —A(S’ju, RkAkau) — A(RkAkS’ju, 35..) + A(RkAksju, RkAksju) = A(u, u) — A(S'ju, Sfu) + A((I,c — RkAk)S’jv, (I,c — RkAk)S'ju) g A(u, u) — A(ssu, 352.) + A(S’ju, $5..) = A(u, it). Since 0 S A((IJ — QiRkAkSfim, (IJ —- QiRkAkS’jM) S A(u, v), IllUJ-QiRk/lksfllll S 1- D Theorem 3.1 Let co = 23"“ — u" be the initial error and el = 13"“ — an“ be the error after the V-cycle algorithm. Then el = Keo with |||K||| < 1. Therefore the V-cycle scheme is convergent. Proof: Since |||I — RkAklll < 1 for 1 S k S J, lllTkIHZIHIJ — QiRkAkS’jlll 31 by Lemma 3.3 (d). We also have HITJIH 2 III] -- RJAJHI < 1 and |||T1|H = p(Tl) = p(I — QIJSD = 1. Thus we have |||K||| < 1. The V-cycle scheme is convergent. C] 33 CHAPTER 4 Stability and Convergence Consider the linear systems (1.4) and (1.5). For each time step n, we solve each of (1.4) and (1.5) with m V-cycles. Let |||u|||A = A(u,u) and |||v|||AI = A’(u,u). For a V-cycle applied to (1.4) and (1.5), we have the convergent rates |||K|||A < 1 and |||K’|||Ar < 1 respectively. Consider w = (a) on the vector space 0 x Q with v u E Q and v E {2. We can define a norm for w E Q x Q by |||w||| = lllvlllA + HIleAI. It can be verified that I” - ||| satisfies the following prOperties: lll’w|||>O z'f waéO. |||w|||=0 Z'f w=0, lllcw|||=|c||||w||| for any complex number c, |||w1+w2||lSlllw1|||+|||w2|||- To show the stability and convergence of the scheme we need the following lemma. B1 Lemma 4.1 Let B = , where 81 and 32 have the same size as that of 0 B2 A-Thenll|B|||=max{ mam... mam. }. Proof: Note that Bw Bu + By lllB||l=maxm I”: MI” 1 III. III 2 III... we lllwlll w... Illummllvm. 34 But _ max IIIUIIIA+IIIUIHAI HIUIIIA’ lllvlllA' max{ Imam... mam. }. lllBIUHIA+|||B2v|||A' < IIIBIUIIIA |||B2’UIIIA' } |/\ Therefore IIIBIIISmax{ mam... |||B2IIIA' }. We want to study the stability of the fully implicit discretization with V-cycle scheme as a solver. The following theorem shows that the scheme is unconditionally stable. Then any instability exhibited in computation should come from the original reaction-diffusion system. Theorem 4.1 Let u" and v" be the numerical solution of (1.4) and {1.5) with the 0 ii of (1.4) and (1.5) with the initial values a9, and 27?]. Let w = (u) and w = (If). If v v initial values u = u0(2:,-j) and v9,- = v0(:1:,-j). Let u" and .7" be the numerical solution the partial derivatives fa, fv, gu, and gv are continuous and uniformly bounded, then lllw" - IVIII S (1+ CAI5)"|||w° - w°|||, where C > 0 is a constant independent of h and At . Remark: The general definition of stability for finite difference schemes can be found in [51]. Proof: Let K and K’ be the error operator for a V-cycle (see Chapter 3) performed on (1.4) and (1.5) respectively. For each n, let x", 5:", y", and g" satisfy Ax" 2 un‘1 +Atf""1, AIL” ___ 2171-1 + Atfn—l, 35 Alyn : Un_1+Atgn-l, A’g" -.—. in—1+Atg"-1, and u", a", v", and v" be the numerical solutions after m V-cycles are applied to the above linear systems respectively. So gn __ {In = (Kl)m(gn _ fin—l). In Chapter 3, we have shown that |||K|||A <1. So |||K’|||Ar < 1 as well. Since and we have (96" - 11") u" — it" u" — a" u" -— 21" Since “ally ij III" ___ A—lun—l + AtA—lfn—l in : A—lan—l + AtA—lfn—l, (1 — Km)(.:" — 5:") + Km(u"-1— a“), (1 — Km)A-1(u"-1— an~1)+ At(I — Km)A-1(f"-1 — fn-l) +K'"(u" — an“), (A-1+ Kmu — A-1))(u"-1- ran-1) +At(1 — Km)A-1(f"-1— fn-l). (4.1) — M222}. 27:3) = minim; - 22:3) + f.(C.-’;. 773)(v?j — 6:3,) 36 n n ijvv' where ( 3,773) is between (u 21') and (113,173.), (4.1) can be expressed as u" — a" = (A-1 + Km(1 — A-1))(u"-1— a"—1)+ At(1 — Km)A-IB;',(u"-1— a“) + At(1 — Km)A—IB;:,(U"-1 — a“), (4.2) where 8?, and Big are diagonal matrices depending on the values of fu and f,,. Similarly, 22“ — = (-1+ 0 fixed and any n such that nAt g T, we have IHEZHIHA S IIIESIIIA +AtZ(|||fj - fjHlA +63“) i=0 and n |||E3+1|||AI S IIIESIIIAI + AtZflllb’ - g’lllw +53“), i=0 where ej = 0(At + h3/2). 38 Proof: Again, we only establish the estimates for u, since the estimates for v follow similarly. Lemma 1.1 shows that Ax?“ = a" + Atf" + AtRn“, n+1 where R"+1 is a vector whose components are r,,- as in Lemma 1.1. Thus 21"“ = AM" + AtA-lf" + AtA‘an“. Set Assn“ = u" + Atf". Then 27"“ = A'lu" + AtA"1A‘lf". Let K be the error operator defined in Chapter 3 and assume m V—cycles are performed at each time step. Then 1”“ — u“1 = K"‘(:z:"+1 — u"). Moreover un+1 : (1_Km)xn+l +Kmun = (1 — Km)(A-1u" + AtA"1f") + Kmu” = Km(I — A‘1)u" + A'lu" + (I — Km)AA‘1f” and an+1— u"+1 = A-1(a" — u") — K’"(I — A-1)u" + AtA“(f" — f") + KmAtA-lf" +AtA‘1R"+1 = A-Iw” — u") + Km[AtA“lf" — (1 — A—1)u"] + AtA“1(f" — f”) +AtA‘1Rn“. Recall that MIA-1|“); S 1, |||I— A‘llllA S 1, and |||K|||A <1. Thus we may assume (HIKIHA)m S (At)2- But then Illa"+‘—u"+‘IHA s Illa"—u"III.+AtHIf"—f"IIIA+(HIKIHMHIMIH. +(IIIK1II.)'"AtIIIf"IIIA+AtIHR"+‘IHA < Illa"—u"IIIA+At|IIf"—I"IHA+(Atfllluvll.+(At)3lllf"lll.. +AtlllR"“H|A = Illil" - unlllA + A73(lllf‘.“ - f"|||A + AItIIIU’illlA + 5'“). (4.4) 39 where e"+1=(At)2lllf"|||A +|||Rn+1|llA + Atlllu"|||A. We have n 1 A(Rn+lu, Rn+1u) IIIR + IIIA = I535: A(u,u) (ARfl't-lu, Rn+1u)J u¢0 (AU,U)J x (AR"+1u,R"+lu)J(R”+1u,R"+lu)J u¢0 (R"+1u,R"+1u)J (Au, U)J Since A 2 AJ are SPD with respect to (-, ')J, (Au,u)J Z Am,n(A)(u, '11)] and (ARn+lu,Rn+IU)J < A . (R"+1u,R"+1u)J —” ”J Therefore |||R"+1IIIA S (lll A “J H R"+1 ”J and 2 HR" ll3= Z h2(.;;)2+ Z tweak (0(At+h3/2)). 1:069 mean since there are only 0(1 / h) boundary points where r?)- = 0(At + h). Thus 6"“ = AtlllntHA +HIRn+1H|A + AtllluanA = 0(At + h3/2). This together with (4.4) gives Illa"+1 - un'HlllA S III210 - UOIHA + NZUIIP - fjlllA + 6j“) 1:0 and lllEZHHIA S IIIEBIHA + AtZUllfj - fjHlA + 6’“)- i=0 Tl Corollary 4.1 Let E" = (g?!) with E3 and E3 defined in Theorem 4.2. If the U partial derivatives fu, fv, gu, and g,, are continuous and uniformly bounded, then for T > 0 such that nAt = T, we have lllE"||| S BCTUIIEOIII +8"). where e" = T0(At + h3/2). 40 Proof: From the properties of fu, f,,, g,,, and 9,, and the proof in Theorem 4.1, we obtain fj_. j j 1 f. f = B11 Bl? (1.13.1 _ wj), g] — 9] Bil Biz with Bj Bj 11 12 SEC). Bi». Biz Here Bu, Ba, 82), and 822 are defined in Theorem 4.1 and the constant C is inde- pendent of h and At. Therefore fi‘fi le—wjmmejm gJ__gJ and [HF - fjHlA + I”? ‘9jHlA’ S ClllEjlll- Hence "—1 n—l lllE"||| S HIEOIH +01“): HIEJIH +At 26’“- i=0 J—O Let e" 2 At 2?:1 61 = TO(At + h3/2). Then ej S e" for 0 g j g n and n—1 IIIE"||| S |||E°||| +6" +046: IIIE’III i=0 n-2 = |||E°||| +6" +CAtlHE"“III +0415: IIIE’III j=0 n—2 n—2 S IIIEOIH+6"+CAt(|l|E°H|+6”"1+CAt |||EJ|||)+CAtZ|||E’|H j=0 j=0 n—2 . n-2 . S |||E°|||+6"+CAt(|||E°||l+6"+CAtZ|||EJ|||)+CAtZ|HE’||| j=0 j=0 n—2 = (1+CAt)(lllE°lll+e")+(1+CAt)(CAtZIIIE’III) j=0 n—2 _ = (1+CAt)(|||E°|H+8"+XIIIEJIII) j=0 S (1+ CAt)”‘1(HIE°|H + 8" + CAtIIIEOIII) 41 S (1+ CAt)"(HIEOIII + e") |/\ eXpC"A‘(|||E°||| + e") = expCT(|||E°|l| +8")- [I Remark: The number m of V-cycles applied to (2.1) is assumed to be large in the proof above in order to obtain the stability and convergence theory. However, in the practical computations, one or two V-cycles performed at each time step are sufficient for the desired accuracy. This efficient scheme will be demonstrated in the next chapter. 42 CHAPTER 5 Experimental Results In this chapter, two reaction-diffusion systems are used to demonstrate the ef- ficiency of the fully implicit finite difference discretization with V-cycle as a solver. The first example is a substrate-inhibition reaction-diffusion mechanism experimen— tally studied by Thomas(1976). It has been shown to be able to form patterns on animal coats [42]. First, we do a linear stability analysis on a rectangular domain to predict the possible pattern on this domain. Then a numerical simulation is carried out to verify this estimated pattern. Finally, we use this system to generate spots and stripes by numerical simulations. To demonstrate that this numerical scheme is valid for other shapes of domains in practice, the pattern on a cone surface is also carried out. There are many other systems which has been used to study spatial patterns [3] [26] [42]. All these systems are capable of generating animal coat markings. The second example we use is the Schnakenberg( 1979) reaction [41]. It is used to generate spots, horizontal and vertical stripes in our numerical computations. For each numerical simulation, the computer program terminates when the nu- merical solutions u and v reach their equilibrium. In the simulation of the pattern on a tail, it requires two V-cycles at each time step for the numerical solution to con- verge to the equilibrium. Each of other simulations requires only one V-cycle at each 43 time step for the numerical solution to converge. Thus at the end of this chapter, we compare the running time needed in using the multigrid technique with that needed in using a classical iterative method on the fine grid alone. The results show that multigrid technique is much more efficient. 5.1 The Thomas System The example we consider here is taken from [41], the Thomas system. It is a substrate inhibition system: U, = DUAU+F(U, V) Vt = DVAV+G(U, V), where ksUV FUV :1. —k — (’ ) 1 2U k6+k7U+k8U2 and k5UV V :1. —k — . Gw’ ) 3 4V k6+k7U+k8U2 Here U and V are the concentrations of two chemical species with U a substrate and V an inhibitor. The constants DU and DV are positive diffusion coefficients and F (U, V) and G (U, V) are the kinetics with the positive rate constants ks. The nondimensional system is given by u, 2 Au + 7f(u,v) 2). = dAv+79(U.v). where f(u.v) = 0““ fl and (5.1) 9(u.v) = 0(1) — ’U) — Fifi—17 . . . Bu 8v . . . . . w1th zero flux boundary condltlon -—— = — = 0 and 11111313] condltlon u(2:, y, 0) an an an an and v(a:, y,0) given on Q. 44 Let us choose 0 as [0, 2.1] x [0,1.2] and parameter values d = 10, ’y = 20, a = 1.5, K = 0.1, p = 18.5, a = 92, and b = 64. The parameter values determine a steady state u, = 9.934 and v3 = 9.289. To form inhomogeneous spatial patterns, a reaction diffusion system must undergo diffusion driven instability or Turing instability. The uniform steady state is stable to small disturbance when the diffusion terms are absent but unstable to small spatial disturbance in the presence of diffusion. With the above parameter values, linear theory gives the range of modes k2 that are driven unstable: 2 2 1 . < k2 = 2 m— n— < . 05_ 7r(2.12+1.22)_1025, which means m2 n2 1.0637 < — — < 1.9507, — 2.12 + 1.22 _ where m and n are integers. Details in deriving the range of unstable modes can be found in [41]. The above range of unstable wavenumbers admits only the wavenumber m = 2, n 2 1. The solution which involves exponentially growing modes about the uniform steady state us and v, is given as 2 27m 7r 02,1810: )‘ cos —— cos i .2 2.1 1.2’ (5 ) where Cm is determined by initial conditions. Figure 5.1 is the pattern obtained from ( 5.2) to predict the pattern that could be formed. The values on the dark regions are greater than zero. The top graph is obtained with 02,1 < 0 and the bottom with 02,1 > 0. Which of these two solutions is obtained thus depends on the bias in the initial conditions. Turing(1952) suggested that, under certain conditions, the homogeneous steady state was unstable to small spatial perturbations and the stable nonuniform spatial patterns could evolve by diffusion driven instability. In the numerical experiments, the initial conditions are taken as random perturbations confined to a small region about the steady state. Figure 5.2 shows the concentrations in the morphogen u. 45 0.5 1 V is 2 Figure 5.1: Possible patterns estimated by linear stability analysis. The dark regions represent u < us. There are two patterns formed by using different seeds in the function calls of random number generator. They agree with the patterns suggested by linear stability analysis. The next three figures are obtained by numerical simulations of different patterns generated with different parameter values of 7 and different geometry. Each example is computed until the equilibrium is reached. The computing details of the simulation are as follows. In Figure 5.3, the pattern (top) of morphogen concentration u < uo = 10 in system ( 5.1) is computed on a grid of 65 x 257 points with 'y = 8. The graph of the residual (bottom) shows that morphogen u reaches the equilibrium at t z 400. In Figure 5.4, the pattern of morphogen concentration u < uo = 10 in system ( 5.1) is computed on a grid of 65 x 225 points with 7 = 48. The graph of the residual shows that morphogen it reaches the equilibrium at t z 60. In Figure 5.5, the pattern of morphogen concentration u < no 2 10 in system 46 Figure 5.2: Turing’s patterns obtained by numerical simulation. -2 h=o.01875 k=0.005 .4 1 V-que at each time step Boooorimeslepe 100100) 100 200 NO‘ 4W 500 Figure 5.3: Numerical simulation of Turing pattern of ( 5.1) and residual plot with 'y = 8. 47 _LQ h=0.01875 k=0.005 -5 . 1 V-cycle at each time step 12000 time steps I0910(7) -10 . -15 ‘ . . - . 0 Figure 5.4: Numerical simulation of Turing pattern of ( 5.1) and residual plot with 7 = 48. #12005 k=0.01 2 V-cycies at each time step 12000 time steps -6» 100100) 20 40 60 80 100 120 140 Figure 5.5: Numerical simulation of Turing pattern of ( 5.1) and residual plot on a cone surface. 48 ( 5.1) is computed on a grid of 65 x 225 points with 7 = 35. The graph of the residual shows that morphogen u reaches the equilibrium at t z 120. Rings of pattern typical of many spotted animal tails are generated when the reaction diffusion domains are tapering cylinders. Rings are at the tip and spots are obtained as the circumference increases. Although the domain is not rectangular, V-cycle algorithm could still be applied analogously. 5.2 The Schnakenberg System The main forms of patterns generated by a reaction-diffusion system are spots (Figure 5.4 and Figure 5.7), rings (Figure 5.5), horizontal stripes (Figure 5.8), and vertical stripes(Figure 5.3 and Figure 5.6). In this section, we choose another reac- tion kinetics. It is the simplest class of two—species reaction mechanism studied by Schnakenberg(1979) in limit cycle solutions of two-species reaction systems. It has also been shown to be able to generate animal coat markings [41]. Its nondimension- alization gives U: = A7~t+’7f(u,’U) v; = dAv+7g(u,v), where f(u,v) = a—u+u2v and (5-3) g(u,v) = b—u2v, with parameter values (1 = 50, a = 0.2, and b = 2. Thus the uniform steady state is uo = 2.2 and v0 = 0.413. Next we demonstrate the V-cycle algorithm with this reaction-diffusion system. We change the scale factor 7 and the domain to obtain different forms of pattern. The computing details follow. 49 h=-0.01875 k=0.005 1 V—cycle at each time step -5 . 12000 time steps 100100) -10’ -15 0 Figure 5.6: Numerical simulation of Turing pattern of ( 5.3) and residual plot of with 7 = 15. In Figure 5.6, the pattern of morphogen concentration u < uo = 2.2 in system ( 5.3) is computed on a grid of 65 x 225 points with 7 = 15. The graph of the residual shows that morphogen u reaches the equilibrium at t z 54. In Figure 5.7, the pattern of morphogen concentration u < uo = 2.2 in system ( 5.3) is computed on a grid of 65 x 217 points with 7 = 30. The graph of the residual shows that morphogen u reaches the equilibrium at t z 100. In Figure 5.8, the pattern (top) of morphogen concentration u < uo = 2.2 in system ( 5.3) is computed on a grid of 65 x 225 points with 7 = 50. The graph of the residual (bottom) also shows that morphogen u reaches the equilibrium at t z 26. 50 h=0.03750 10-0005 1 V—cycle at each time step 20000 time steps Figure 5.7: Numerical simulation of Turing pattern of ( 5.3) and residual plot with 7 = 30. h=0 03750 k-O. 005 1 V-cycle at each time step 1 5400 time steps -‘5 A L 0 Figure 5.8: Numerical simulation of Turing pattern of ( 5.3) and residual plot with 7 = 50. 51 5.3 Comparison We end this thesis by demonstrating the robustness of the V-cycle algorithm. We compute the solution of system ( 5.1) with 7 : 25 on the domain [0, 1.6] x [0, 1.6]. We discretize ( 5.1) with different spatial and time step sizes. Both the V-cycle algorithm and direct iterative method using red-black Gauss-Seidel iteration are applied to each case. In each case, the solution is obtained when the morphogen concentrations u and v reach their equilibrium. The number of iterations at each step, the number of time steps needed for the morphogen concentration to reach its inhomogeneous steady state, and the running time for each case are recorded in the following tables. Notice that when the spatial step size in each direction is refined by half (see Table 1 and Table 2) in order to keep the same ratio of %, time step size is changed into a factor of %. Comparing the running time between these two cases it does not show the power of the V-cycle algorithm. Since the fully implicit finite difference discretization is unconditionally stable, we could exploit this advantage by increasing the step size in time to reduce the running time. Thus in Table 3, we double the time step size. Then only one V-cycle is needed, so the running time for V-cycle algorithm is reduced by about half. But the direct iterative method suffers from its limitations because of the change of the ratio %. Thus the number of iterations at each time step must be increased in order to have convergence. Even if the fully implicit finite difference scheme is unconditionally stable, the direct iterative method seems to be unable to take advantage of it. In Table 4, we can see when the spatial grid size is even finer — the running time for direct iterative method is almost 10 times that needed for V-cycle algorithm. 52 Figure 5.9: Numerical simulation of Turing pattern on a square. h = 0.05 number of iterations number of CPU At = 0.01 at each time step time steps time Red-black Gauss-Seidel iteration 25 3000 39.78 seconds V—cycle 1 3000 10.8 seconds Table 5.1: Comparison between single grid and multigrid with h = 0.05 and At = 0.01. h 2 0.025 number of iterations number of CPU At = 0.0025 at each time step time steps time Red-black Gauss-Seidel iteration 22 12000 607.54 seconds V-cycle 1 . 15000 219.18 seconds Table 5.2: Comparison between single grid and multigrid with h = 0.025 and At = 0.0025. h = 0.025 number of iterations number of CPU At = 0.005 at each time step time steps time Red-black Gauss-Seidel iteration 45 5600 569.86 seconds V-cycle 1 7000 103.43 seconds Table 5.3: Comparison between single grid and multigrid with h = 0.025 and At = 0.005. 53 h = 0.0125 number of iterations number of CPU At = 0.0025 at each time step time steps time Red-black Gauss-Seidel iteration 85 11000 8205.23 seconds V-cycle 1 14000 857.52 seconds Table 5.4: Comparison between single grid and multigrid with h = 0.0125 and At = 0.0025. 54 BIBLIOGRAPHY 55 BIBLIOGRAPHY [1] K. E. Atkinson, An Introduction to Numerical Analysis, John Wiley, New York, 1978. [2] R. E. Bank and C. C. Douglas, Sharp estimates for multigrid rates of convergence with general smoothing and acceleration, SIAM J. Numer. Anal., 22(1985), pp. 617-633. [3] J. Bard, A model for generating aspects of zebra and other mammalian coat patterns, J. Theoret. Biol., 93(1981), pp. 363-385. [4] D. Braess and W. Hackbusch, A new convergence proof for the multigrid method including the V cycle, SIAM J. Numer. Anal., 20(1983), pp. 967—975. [5] J. H. Bramble, Multigrid methods, Longman Scientific & Technical, 1993. [6] J. H. Bramble and J. E. Pasciak, The analysis of smoother for multigrid algo- rithm, Math. Comp., 58(1992), pp. 467—488. [7] J. H. Bramble J. E. Pasciak, New convergence estimates for multigrid algorithms, Math. Comput., 49(1987), pp. 311-329. [8] J. H. Bramble, J. E. Pasciak, J. Wang, and J. Xu, Convergence estimates for multigrid algorithms without regularity assumption, Math. Comput., 57(1991), pp. 23-45. [9] J. H. Bramble, J. E. Pasciak, and J. Xu, The analysis of multigrid algorithm for nonsymmetric and indefinite problems, Math. Comput., 51(1988), pp. 389-414. [10] J. H. Bramble, J. E. Pasciak, and J. Xu, Guide to multigrid development, W. Hackbusch and U. Trottenberg, eds., Springer-Verlag, New York, 1982, pp. 220- 312. v [11] N. F. Britton, Reaction-Diffusion Equations and Their Application to Biology, Academic Press, New York, 1986. 56 [12] R. Bulirsch and J. Stoer, Introduction to Numerical Analysis, Springer-Verlag, 1992 e" [13] C. Chiu, F. C. Hoppensteadt, and W. Jager, Analysis and computer simulation of accretion patterns in patterns in bacterial culture, J. of Math. Biology, (1994) 32, pp. 841-855. [14] C. Chin and N. Walkington, An ADI method for Reaction-diflusion systems and Hysteric Reaction-Diffusion Torpor-Systems, SIAM J. Numer. Anal., (1997) 34, No.3, pp. 1185-1206. [15] C. Chiu and N. Walkington, Analysis of hysteric reaction-diffusion systems, Quarterly of Appl. Math., (1998) LVI, No.1, pp. 89-106. [16] C. C. Douglas, Abstract multi-grid with application to elliptic boundary-value problems, in Elliptic Problem Solvers II, G. Birkhoff and A. Schoenstadt, eds., Academic Press, New York, 1984, pp. 453-466. [17] C. C. Douglas, Multi-grid with application to elliptic boundary-value problems, SIAM J. Numer. Anal., 21 (1984), pp. 236-254. [18] C. C. Douglas and J. Douglas, Jr., A unified convergence theory for abstract multigrid or multilevel algorithms, serial and parallel, SIAM J. Numer. Anal., 30 (1993), pp. 136-158. .1" [19] A. Gierer and H. Meinhardt, A theory of biological pattern formation, Kyber- netik, 12(1972), pp. 30-39. [20] A. Greenbaum, Analysis of a multigrid method as an iterative technique for solv— ing linear systems, SIAM J. Numer. Anal., 21 (1984), pp. 473-485. ‘1 [21] D. Hoff A Stability and Convergence of Finite Difference Method for Systems of Nonlinear Reaction-Difi'usion Equations, SIAM J. Numer. Anal., 12 (1978), pp. 1161-1177. [22] G. Horton and S. Vandewalle, A space-time multigrid method for parabolic partial differential equations, SIAM J. Sci. Comput., 16 (1995), pp. 848-864. . [23] Y. Keshet and L. A. Segel, Pattern formation in aspects, Modeling of Patterns in Space and Time, W. Jager, ed., Springer-Verlag, Berlin, pp. 188-197. 57 [24] CC. Kuo and BC. Levy, Two-Color Fourier Analysis of the Multigrid Method with Red-Black Gauss-Seidel Smoothing, Appl. Math. Comput., 29 (1989), pp. 69-87. [25] S. A. Levy, Non-uniform stable solutions to reaction-diffusion equations: appli- cations to ecological pattern formation, Pattern Formation by Dynamic Systems and Pattern Recognition, H. Haken, ed., Springer-Verlag, Berlin, pp. 210-222. [26] S. A. Levin and L. A. Segel, Pattern Generation in Space and Aspect, SIAM Review, 27 (1985), pp. 45-67. [27] N. Li, J. Steiner and S. Tang, Convergence and stability analysis of an explicit fi- nite difference method for 2—dimensional reaction-diffusion equations, J. Austral. Math. Soc. Ser. B, 35 (1993), pp. 223-243. [28] P. Lotstedt and B. Gustafsson, Fourier analysis of multigrid methods for general systems of PDEs, Math. Comp. 60 (1993) pp. 473-493. [29] J. F. Maitre and F. Musy, Multigrid methods: convergence theory in a variational framework, SIAM J. Numer. Anal., 21 (1984), pp. 657-671. [30] J. F. Maitre and F. Musy, Multigrid methods for Symmetric Variational Prob- lems: A General Theory and Convergence Estimates for Usual Smoothers , Appl. Math. Comput, 21 (1987), pp. 21-43. [31] J. Mandel, Multigrid convergence for nonsymmetric, indefinite variational prob- lems and one smoothing step, Appl. Math. Comput., 19 (1986), pp. 201-216. [32] J. Mandel, S. McCormick, and R. Bank, Variational multigrid theory, Multi- grid Methods, (S. McCormick, ed.), SIAM (Society for industrial and Applied Mathematics), Philadelphia, PA, 1987, pp. 131-178. [33] J. Mandel, Algebraic Study of Multigrid methods for symmetric, definite problems , Appl. Math. Comput., 25 (1988), pp. 39-56. [34] J. Mandel, S. McCormick, and J. Ruge, A algebraic theory for multigrid method for variational problem, SIAM J. Numer. Anal., 25 (1988), pp. 91-110. [35] S. F. McCormick, An algebraic interpretation of multigrid methods, SIAM J. Numer. Anal., 19 (1982), pp. 548-560. [36] S. F. McCormick, Multigrid methods for variational problems: further results, SIAM J. Numer. Anal., 21 (1984), pp. 255-263. 58 [37] S. F. McCormick, Multigrid methods for variational problems: general theory for the V-cycle, SIAM J. Numer. Anal., 22 (1985), pp. 634-643. [38] S. McCormick and J. Ruge, Multigrid methods for variational problems, SIAM J. Numer. Anal., 19 (1982), pp. 924-929. \l[39] H. Meinhardt, Models of Biological Pattern Formation, Academic Press, New York, 1982. Methods Partial Differential Eq, 15 (1999), pp. 201-214. 'J’ [40] R. E. Mickens, Nonstandard finite difference schemes for reaction-diffusion equa- tions, Numer Methods Partial Differential Eq, 15 (1999), pp. 201-214. [41] J. D. Murray, Mathematical Biology, Biomathematics Texts, Springer-Verlag, 1989. [42] J. D. Murray, A Per-pattern Formation Mechanism for Animal Coat Markings, J. Theor. Biol., 88 (1981), pp. 161-199. [43] R. A. Nicolaides, 0n multigrid convergence in the indefinite case, Math., Comp. 32 (1978), pp. 1082-1086. [44] J. M. Ortega, Matrix Theory: A Second Course, Plenum Press, New York, 1987. [45] C. V. Pao, Accelerated monotone iterative methods for finite difference equations of reaction-diffusion, Numer. Math, 79 (1998), pp. 621-281. [46] C. V. Pao, Numerical analysis of coupled systems of nonlinear parabolic equa- tions, SIAM J. Numer Anal, 36 (1999), pp. 393-416. [47] S. V. Parter, Estimates for multigrid methods based on red-black Gauss-Seidel smoothings, Numer. Math., 52 (1988), pp. 701-723. [48] J. I. Ramos, Linearization methods for reaction-difiusion equation: I-D problems, Appl. Math. and Comput., 88 (1997), pp. 199-244. [49] J. I. Ramos, Linearization methods for reaction-diffusion equation: multidimen- sional problems, Appl. Math. and Comput., 88 (1997), pp. 225—254. [50] R. D. Reitz, A Study of Numerical Methods for Reaction-Diffusion Equations, SIAM. J. Sci. Stat. Comput., 2 (1981), pp. 95-106. " [51] R. D. Richtnyer and K. W. Morton, Difference methods for initial value problems, Interscience, New York, 1967. 59 [52] J. Smoller, Shock Waves and Reaction-Difiusion Equations, Berlin Heidelberg New York Tokyo, Springer 1983. ‘ [53] J. C. Strikwerd, Finite Difference Schemes and Partial Differential Equations, Wadsworth & Brooks/Cole Advanced Books & Software, Pacific Grove, Calif, 1989. [54] U. Trottenberg and K. Stuben, Multigrid Methods: Fundamental Algorithms, Model Problem Analysis and Applications. Multigrid Methods (Proceedings) Lec- ture Notes in Mathematics. Ed. W. Hackbusch and U. Trottenberg. Spring-Verlag, New York, 1982, pp. 1-176. [55] J. Wang, Convergence analysis of multigrid algorithms for nonselfadjoint and indefinite elliptic problems, SIAM J. Numer. Anal., 30 (1993), pp. 275-285. [56] Y.M. Wang, Petrov-Galerkin methods of nonlinear reaction-diflusion equations, Appl. Math. Comput., 96 (1998), pp. 209-236. [57] P. Wesseling, Theoretical and practical aspects of a multigrid method, SIAM J. Sci. Statist. Comput., 4 (1982), pp. 387-407. [58] D. D. Young, Iterative Solution of Large Linear Systems, Academic, 1971. [59] F. Zhang, Matrix Theory: Basic Results and Techniques, Springer-Verlag, New York, 1999. [60] J. Zhang, Acceleration of Five-point Red-Black Gauss-Seidel in Multigrid for Poisson Equation, Appl. Math. Comput., 80 (1996), pp. 73-93. 60