探討反應擴散方程系統之穩態解

全文

(1)國立高雄大學應用數學系碩士論文. 探討反應擴散方程系統之穩態解 An Investigation into the Steady State Solutions of Reaction Diffusion Systems. 研究生：謝偉勃撰指導教授：陳晴玉. 中華民國一百零四年七月.

(2) 致謝首先要感謝指導教授陳晴玉老師，不管大學時期修過的生物數學或是研究所的專題課在知識跟方向上都給了我很大的幫助，至於甚麼漸進與擾動之類的就暫時先不提它吧。這次論文能夠完成除了歸功於陳老師之外，郭岳承老師在延拓法與線性代數方面的指導也給了我很大的幫助。然後曾昱豪老師在口試技巧與論文中一些錯誤的提點也給了我很大的助益。. 能平安度過研究所這兩年還多虧了很多人的幫助；像是願意耐心回答我粗淺問題的施信宏老師、常常一起討論各式問題的吳士逸同學、提供了很多建議跟方向的陳正時學長以及其他系上的老師、同學、學長姐、學弟妹。也要感謝系辦的雅鳳姐、佩君姐在行政事務上的幫助。. 最後除了感謝家人願意讓我念研究所之外，因為需要感謝的人太多了，就感謝豬豬吧！.

(3) 探討反應擴散方程系統之穩態解學生：謝偉勃國立高雄大學應用數學系. 指導教授：陳晴玉教授國立高雄大學應用數學系. 摘要. 在這篇論文中，我們會專注在反應擴散方程式系統的穩態解並將利用穩態解來解釋一些實際觀察到的生物圖形。我們的目標是研究方程式定義域的大小對從系統所產生的圖形之影響。我們會從圖靈對反應擴散方程式系統的假設出發，並進行了線性穩定性分析，推導出擴散驅動不穩定的必要條件。這個結果提供了一個在線性化系統中定義域的大小與允許的波數間的定量關係。為了評估由線性化系統所提供的預測，我們使用了一個以定義域的大小作為參數的弧長延拓法來定位分岔點。該方法是作用在一個以正交多項式作為基底函數的譜方法離散得來的離散化系統。不論是一維或二維的情況得到的數值結果都會與線性化分析中得到的結論比較。. 關鍵字: 圖形生成、反應擴散方程系統、延拓法、譜方法、正交多項式。.

(4) An Investigation into the Steady State Solutions of Reaction Diffusion Systems. Student: Xie, Wei Bo Department of Applied Mathematics National University of Kaohsiung Advisor: Professor Ching-Yu Chen Department of Applied Mathematics National University of Kaohsiung. ABSTRACT. In this work we focus on steady state solutions of reaction diffusion systems with applications to biological patterns. Our goal is to investigate how domain size affects the patterns generated by the systems. We adopt Turing's idea and carry out linear stability analysis to derive the necessary conditions for diffusion driven instability. This provides for the formulation that relates the domain size to permissible wave numbers for the linearised system. To evaluate the prediction offered by the linearised system, we use an arclength continuation method with domain size as a parameter to locate the bifurcation points. The method is implemented using a spectral collocation scheme for which orthogonal polynomials are used as base functions. Results are compared with those from the linearised system for both one and two dimensional cases.. Keywords:. pattern formation, reaction diffusion equations, continuation method, spectral collocation method, orthogonal polynomials..

(5) TABLE OF CONTENTS Chapter 1. Introduction ............................................................................................... 1. 2. The effect of domain size on wave numbers ............................................... 3. 2.1. Linear analysis ................................................................................. 3. 2.2. Analysis on an discretized system ................................................... 6. 3. The choice of base functions ..................................................................... 10 3.1. Orthogonal polynomials that satisfy zero flux boundary conditions ............................................................................................... 10 3.1.1. Standard p-polynomials ....................................................... 11 3.1.2. Chebyshev p-polynomials ..................................................... 12 3.1.3. Legendre p-polynomials........................................................ 14 3.2. Recurrence relation of p-polynomials ............................................ 15 3.2.1. Three-term recurrence relation............................................. 15 3.2.2. Five-term recurrence relation ............................................... 18 3.3. Numerical result about p-polynomials............................................ 21 3.3.1. Rate of convergence.............................................................. 22 3.3.2. Boundary value problem(BVP)............................................ 26 3.4. Standard p-polynomials, Legendre p-polynomials and Chebyshev p-polynomials are orthogonal basis of L2 ([−1, 1])........................... 29 4. Numerical results ...................................................................................... 32.

(6) 4.1. Numerical methods......................................................................... 33 4.2. Only one wavenumber is permissible .............................................. 35 4.3. Multiple permissible wavenumbers ................................................. 44 5. Conclusion................................................................................................. 52 References .............................................................................................................. 53 REFERENCES............................................................................................... 53.

(7) CHAPTER 1 Introduction In 1952, Turing considered a model to describe the mechanisms of pattern formation [1], the model had successfully produced many biological patterns such as those of leopards, zebra, snakes, lady bugs. The model is based on two diffusible chemicals, morphogenes, which interact with each other. Turing’s model has the following form Ut = DU ∇2 U + F (U, V ), Vt = DV ∇2 V + G(U, V ),. (1.1). where U (x, t) is the activator and V (x, t) is the inhibitor, DU and DV are diffusion coefficients, F (U, V ) and G(U, V ) are reaction kinetics. Turing thought that equations (1.1) generates patterns if the homogeneous steady state is stable to small perturbations in the absence of diffusion but unstable to small perturbations when diffusion is present. The above property is called diffusiondriven instability. Subsequent study on the model showed that domain size plays significant role in the observed patterns[7]. For example, the pattern on the coat of the leopard in the above figure, the spotted body with stripe trail, is often attributed to the difference in dimension during embryonic stage. In this thesis, we will therefore focus on the changes of steady state patterns as a result of varying domain size for both one and two dimensional cases using an arclength continuation method. Although our focus is on the domain size, the method 1.

(8) here can be used for any parameter of interests with suitable rescaling.. 2.

(9) CHAPTER 2 The effect of domain size on wave numbers 2.1. Linear analysis. To discuss the the influence of domain size, we consider the following onedimensional reaction diffusion equations with Schnakenberg’s kinetics and zero flux condition imposed on the boundary . Ut = DU Uxx + F (U, V ), Vt = DV Vxx + G(U, V ),. (2.1). and the zero flux boundary conditions give Ux = Vx = 0 on both x = 0 and x = L, where U (x, t) is the activator and V (x, t) is the inhibitor, x ∈ [0, L], DU and DV are diffusion coefficients, with F (U, V ) = k1 − k2 U + k3 U 2 V , G(U, V ) = k4 − k3 U 2 V . We nondimensionalise (2.1) with the following r r k3 k3 ∗ DU t ∗ x , v=V , t = 2 , x = , u=U k2 k2 L L r r k1 k3 k4 k3 L2 k2 DV , a= , b= , γ= , d= DU k2 k2 k2 k2 DU and obtain the dimensionless system ut = uxx + γf (u, v), vt = dvxx + γg(u, v), where f (u, v) = a − u + u2 v, g(u, v) = b − u2 v, 3. (2.2).

(10) with boundary conditions ux = vx = 0 on both x = 0 and x = 1, Let (u0 , v0 ) be the homogeneous steady state. According to Turing’s idea of diffusion-driven instability, we require (u0 , v0 ) to be stable without diffusion and unstable with diffusion. To examine the stability of the steady state, we consider small perturbations near (u0 , v0 ), u(t, x) = u0 + uˆ(t, x), v(t, x) = v0 + vˆ(t, x), where kˆ u(t, x)k 1 and kˆ v (t, x)k 1. Without diffusion the linearised system is : uˆt uˆ = A vˆ vˆ ,. (2.3). t. where A=γ. . fu fv gu gv. . . (u0 ,v0 ). For stability, we require T r(A) < 0 and det(A) > 0. With diffusion, linearisation near (u0 , v0 ) gives uˆt uˆ uˆxx = A + vˆt vˆ dˆ vxx .. (2.4). To satisfy zero flux boundary conditions, uˆ and vˆ are linear combination of the following form uˆ ∝ eλt cos(kx), vˆ ∝ eλt cos(kx), k = iπ, i ∈ N, 4. (2.5).

(11) we call k in equation (2.5) wavenumbers. The existence of an unstable steady state requires at least one λ > 0. Thus, substitute (2.5) into (2.4) we get λ. . uˆ vˆ. . = B(k 2 ). . uˆ vˆ. . ,. (2.6). where 2. B(k ) =. . γfu − k 2 γfv γgu γgv − dk 2. .. Recall that in the absence of diffusion we require tr(A) < 0, which consequently gives tr(B(k 2 )) < 0. Hence, for instability it is necessary det(B(k 2 )) < 0. Notice that det(B(k 2 )) is quadratic in k 2 : det B(k 2 ) = d(k 2 )2 − γ(gv + dfu )k 2 + det(A). Hence, det(B(k 2 )) < 0 requires (gv + dfu ) > 0 and (gv + dfu )2 − 4d · det(A) > 0. Furthermore, equation (2.5) requires the wavenumbers being discrete points satisfying k 2 = (iπ)2 , i ∈ N. Thus, to ensure instability it is necessary to verify that there is permissible k 2 between the roots of det(B(k 2 )), namely p γ {(dfu + gv ) − (dfu + gv )2 − 4d(fu gv − fv gu )}, 2d p γ 2 k2 = {(dfu + gv ) + (dfu + gv )2 − 4d(fu gv − fv gu )}. 2d. k12 =. (2.7). To examine the influence of domain size on the pattern observed, we focus on γ because γ ∝ L2 . Furthermore, equation (2.7) show that by changing γ we can determine which k is permissible. Since permissible wavenumber can effect final patterns, we say domain size effect the patterns. 5.

(12) 2.2. Analysis on an discretized system. In this subsection we we provide an alternative way to verify (2.7). For each k = iπ, i ∈ N, there is an interval of γ, γk− < γ < γk+ , for which the diffusion-driven instability holds, namely γk− = k 2 γk+ = k. 2. (dfu + gv )2 − 4d(fu gv − fv gu ) fu gv − fv gu. !. (dfu + gv )2 − 4d(fu gv − fv gu ) fu gv − fv gu. !. (dfu + gv ) −. p. (dfu + gv ) +. p. , (2.8) .. At γk+ and γk− since eigenvalues of B(k 2 ) in linear analysis change from both being negative to one positive and the other negative, γk+ and γk− are bifurcation points. Thus, we seek bifurcation points on the uniform steady state. To do this, we consider steady state solutions of (2.2) . 0 = uxx + γf (u, v), 0 = dvxx + γg(u, v).. (2.9). Discreting (2.9) on a uniform grid {xi }, i = 1, · · · , n, we get . 4 0 0 d4. . ~u ~v. . +γ. 6. f (~u, ~v ) g(~u, ~v ). = 0,. (2.10).

(13) where   u1 v1 u 2    v2 ~u =  ..  , ~v =  .. . . un vn  −2 2  1 −2  .. . 2 4 = (n − 1)    .    , ui = u(xi ), vi = v(xi ),  1 .. ... .. ... .. .. ... .. 1. ... . −2 1 2 −2.    ,  . .   f (u1 , v1 ) g(u1 , v1 )  f (u2 , v2 )   g(u2 , v2 )  and g(~u, ~v ) =  f (~u, ~v ) =  .. ..    . . f (un , vn ) g(un , vn ).   . . We will use a theorem called bifurcation test to find bifurcation points on uniform steady state of (2.10). The theorem are introduced as follows. Definition 2.1. Let G : Rn+1 → Rn . Consider the problem: G(x, λ) = 0. Given a smooth path of solutions Γ = {(x(λ), λ)|λa < λ < λb }, a point (x0 , λ0 ) ∈ Γ is said to be a bifurcation point for above problem if every ball B ⊂ Rn+1 of radius > 0 contains solutions of G(x, λ) = 0 not on Γ.. 7.

(14) Theorem 2.2 (bifurcation test). Let G : Rn+1 → Rn be C 1 . Let Γ : (x(λ), λ) be a smooth solution path of G(x, λ) = 0 parametrized by λ. If det(Gx (x(λ), λ)) changes sign at λ0 ∈ (λa , λb ), then (x(λ0 ), λ0 ) is a bifurcation point of G(x, λ) = 0. Note that in the above theorem λ0 ∈ (λa , λb ) is an interior point of the interval. We can not apply this theorem if λ0 is a boundary point. In other words, this theorem is not valid for turning points. Now, to find bifurcation points on the uniform steady state of (2.10) we apply the bifurcation test theorem on (2.10) by considering 4 0 ~u f (~u, ~v ) G(~u, ~v , γ) = 0 d4 ~v + γ g(~u, ~v ) ,. (2.11). and Γ{(u0 , u0 , · · · , u0 , v0 , v0 , · · · , v0 , γ)|γ > 0}. According to bifurcation test theorem, we seek γ such that det(Gx (γ)) = 0 on Γ, where Gx (γ) =. . 4 0 0 d4. . +γ. . fu I fv I gu I gv I. . , I is a n × n identity matrix.. Rewrite Gx (γ) as Gx (γ) =. . 1 0 0 d. . ⊗4+γ. . fu fv gu gv. . ⊗ I,. where ⊗ is the tensor product. Notice that ∃~v 6= 0 such that 1 0 fu fv → − 0 d ⊗ 4 + γ gu gv ⊗ I v = 0, 8.

(15) is the necessary and sufficient condition of det(Gx (γ)) = 0 and by the propertys of tensor product . fu fv gu gv. . ⊗ I,. is invertible.Hence, eigenvalues of 1 gv −dfv − ⊗ 4, fu gv − fv gu −gu dfu are what we seek. Note that eigenvalues of 4 are 2. (n − 1). . n−j π) , j = 1, · · · , n. − 2 − 2cos( n−1. Hence, by the property of tensor product eigenvalues of (2.12) tend to ! p 2 − 4d(f g − f g ) (df + g ) ± (df + g ) u v u v u v v u k2 , k = iπ, i ∈ N. fu gv − fv gu as n → ∞. This confirms the result of linear analysis.. 9. (2.12).

(16) CHAPTER 3 The choice of base functions In this thesis, we choose spectral collocation method(SCM) as our numerical scheme in the reason of computing cost. The goal of SCM is to approximate a function with a suitable finite sum of base function φ0i s, namely, y ≈ yn (x) =. n X. ai φi (x).. i=1. In this Chapter, to improve our numerical method we will introduce several base functions satisfy the zero flux boundary conditions and some properties of the kind of base functions. 3.1. Orthogonal polynomials that satisfy zero flux boundary conditions In SCM, we often choose Legendre polynomials, Chebyshev polynomials, Trigono-. metric polynomials and other orthogonal functions as trail functions. However, our goal is to apply SCM on Turing’s Mechanisms with zero flux boundary conditions, which means orthogonal functions satisfying zero flux boundary conditions may converge faster than other orthogonal functions. Hence, to verify this idea we introduce the following sequences of orthogonal polynomials that satisfy zero flux boundary conditions. Without lose of generality, we consider those orthogonal polynomials with first derivative being zero at 1 and −1. For the sake of brevity, we will call them piggy-polynomials or simply p-polynomials. 10.

(17) 3.1.1. Standard p-polynomials First, consider pn (x) =. 1, xn+1 −. n+1 n−1 x , n−1. n = 1, n ≥ 2.. (3.1). It can be shown that the above sequence of polynomials is a sequence of polynomials with first derivative at 1 and −1 are both zero because 0. pn (x) =. . 0, n = 1, (n + 1)xn − (n + 1)xn−2 , n ≥ 2.. If we apply Gram-Schmidt orthogonalization on above polynomials with the following inner product Z. 1. hf, gi =. f (x)g(x)dx, −1. we can get a sequence of p-polynomials  i<n P hpn , Pi i   Pi , n is odd,  pn − 2 i is odd kPi k Pn = i<n P hpn , Pi i   Pi , n is even.  pn − 2 i is even kPi k. 11. (3.2).

(18) We call (3.2) standard p-polynomials. The first few standard p-polynomials are listed below.. 1 x − 3x 7 4 x − 2x2 + 15 35 3 x5 − 290 153 x + 51 x 65 4 31 2 43 x6 − 33 x + 33 x − 693 2317 3 329 5 x7 − 357 169 x + 1859 x − 1859 x 1110 4 244 2 73 6 x8 − 364 159 x + 689 x − 689 x + 6201 234 5 3084 3 234 7 x9 − 972 391 x + 115 x − 5083 x + 5083 x 23590 6 8890 4 1155 2 8 x10 − 1485 551 x + 9367 x − 9367 x + 9367 x − 3. 3.1.2. 259 103037. Chebyshev p-polynomials In [3] J. Shen construct a sequence of polynomials with first derivative at 1. and −1 are both zero base on Chebyshev polynomials ψn = (n − 1)2 Tn+1 − (n + 1)2 Tn−1 , n ≥ 1,. (3.3). where Tn represent the Chebyshev polynomials. It can be shown that (3.3) is a sequence of polynomials with first derivative at 1 and −1 are both zero because Chebyshev polynomials have following property 0. (a) Tn (1) = n2 , 0. (b) Tn (−1) = (−1)n n2 .. 12.

(19) Notice that ψn are not orthogonal polynomials. Thus, we apply Gram-Schmidt orthogonalization on ψn and get  n ≤ 3,  ψn , hψ , Ψ i n n−2 Ψn = Ψn−2 , n > 3,  ψn − kΨn−2 k2. (3.4). where (1) hψn , Ψn−2 i = hψn , ψn−2 i, hψn , ψn−2 i2 , (2) hΨn , Ψn i = kψn k2 − kΨn−2 k2  4 2 4 2   (n − 1)2 k Tn+1 2k +(n + 21) k Tn−1 k , if n = m, (n − 3) (n + 1) k Tn−1 k , if n − m = 2, (3) hψn , ψm i = 2 2 2 (m − 3) (m + 1) k T k , if m − n = 2,  m−1  0, other wise. We call (3.4) Chebyshev p-polynomials the first few Chebyshev p-polynomials are listed below.. 1 4x − 12x 8x4 − 16x2 + 5 530 3 16x5 − 1270 41 x + 41 x 576 2 44 4 32x6 − 1104 17 x + 17 x − 17 8904 3 1400 5 64x7 − 14112 101 x + 101 x − 101 x 11040 4 2624 2 97 6 128x8 − 14848 49 x + 49 x − 49 x + 49 146232 5 46470 3 3978 7 256x9 − 24336 37 x + 259 x − 259 x + 259 x 245840 6 32600 4 13600 2 8 512x10 − 252160 177 x + 177 x − 59 x + 177 x − 3. 13. 302 177.

(20) 3.1.3. Legendre p-polynomials We can construct another p-polynomials base on Legendre Polynomials, defin-. ing ωn = n(n − 1)Ln+1 − (n + 1)(n + 2)Ln−1 , n ≥ 1,. (3.5). where Ln represent the Legendre polynomials. Like Chebyshev Polynomials, Legendre Polynomials have the following property n(n + 1) , 2 n(n + 1) 0 (b) Ln (−1) = (−1)n . 2 Hence, (3.5) is a sequence of polynomials with first derivative being zero at 1 and −1. 0. (a) Ln (1) =. Similar to (3.3), equation (3.5) is not an orthogonal polynomials. Thus, by appyling Gram-Schmidt orthogonalization again, we get ( ω , n ≤ 3, n hω , ω i n n−2 Ωn = ωn − Ωn−2 , n > 3, kΩn−2 k2. (3.6). where (1) hωn , Ωn−2 i = hωn , ωn−2 i, hωn , ωn−2 i2 (2) hΩn , Ωn i = kωn k2 − , k Ωn−2 k2  2 2 2 2 2 2   (n − 1) n k Ln+1 k +(n + 1) (n + 2)2 k Ln−1 k , (n − 3)(n − 2)(n + 1)(n + 2) k Ln−1 k , (3) hωn , ωm i = (m − 3)(m − 2)(m + 1)(m + 2) k Lm−1 k2 ,   0,. 14. if n = m, if n − m = 2, if m − n = 2, other wise..

(21) We call (3.6) Legendre p-polynomials.The first few Legendre piggy polynomials are listed below. 1 5 3 15 2x − 2 x 35 4 35 2 49 8 x − 4 x + 24 63 5 1015 3 735 8 x − 68 x + 136 x 231 6 455 4 217 2 43 16 x − 16 x + 16 x − 48 429 7 11781 5 6951 3 987 16 x − 208 x + 208 x − 208 x 195195 6 274725 4 30195 2 4015 6435 8 128 x − 1696 x + 3392 x − 1696 x + 6784 12155 9 173745 7 284427 5 42405 3 13365 128 x − 736 x + 1472 x − 736 x + 2944 x 3610035 8 1686685 6 635635 4 165165 2 3367 46189 10 256 x − 7424 x + 3712 x − 3712 x + 7424 x − 7424. 3.2. Recurrence relation of p-polynomials. The ultimate goal of SCM is to approximate a function with a suitable finite n P sum of base functions. Suppose yn (x) = ai φi (x) is an approximation of a function i=1. y(x). The value of yn at x0 can be obtained by evaluating the value of φi at x0 for i = 1, · · · , n. Thus, the cost of computing all φi (x0 )0 s may be a matter of concern. 3.2.1. Three-term recurrence relation When choosing a sequence of orthogonal polynomials as base functions, three-. term recurrence often help us reduce the cost computing φi (x0 )0 s. Let’s take Legendre. 15.

(22) polynomials for example. The first few Legendre polynomials are. 1 x 1 2 2 (3x − 1) 1 3 2 (5x − 1) 1 4 2 8 (35x − 30x + 3) 1 5 3 8 (63x − 70x + 15x) 1 6 4 2 16 (231x − 315x + 105x − 5) 1 7 5 3 16 (429x − 693x + 315x − 35x) 1 8 6 4 2 128 (6435x − 12012x + 6930x − 1260x + 35) If coefficients are the only thing we know about Legendry polynomials, then the cost of computing Li (x0 ), i = 0, · · · , n is about. n(n+1) 4. flops. Fortunately, Legendre. polynomials satisfy the following recurrence relation Ln =. 2n − 1 n−1 xLn−1 − Ln−2 , L0 = 1, L1 = x. n n. Hence, if we use the above recurrence relation the cost of computing φi (x0 )0 s require only about 2n flops. In fact, it is known that every sequence of ordinary monic orthogonal polynomials satisfies a three-term recurrence relation. Definition 3.1. A sequence of polynomials {φn }∞ n=0 is said to be ordinary if deg(φn ) = n. Theorem 3.2 (Three-term recurrence relation). Let {φn }∞ n=0 be a sequence of ordinary monic orthogonal polynomials, then φn = (x − αn )φn−1 − βn φn−2 , n ≥ 1, φ−1 ≡ 0, φ0 ≡ 1, 16.

(23) where. hφn−1 , φn−1 i , n > 1, hφn−2 , φn−2 i 0, n = 1. Proof. Given a large enough constant k. Since both φk and φk−1 are monic, deg(φk − hxφn−1 , φn−1 i , βn = αn = hφn−1 , φn−1 i. (. xφk−1 ) ≤ k − 1. Since {φn }∞ n=0 is ordinary, φk − xφk−1 is a linear combination of some φn . Combining above results, we have φk − xφk−1 =. k−1 X. a i φi .. (3.7). i=0. Taking inner product with φi to on both sides of (3.7) we get ai = −. hφk−1 , xφi i hxφk−1 , φi i =− . hφi , φi i hφi , φi i. Notice that when i ≤ k − 3, deg(xφi ) ≤ k − 2. Since {φn }∞ n=0 is ordinary, xφi =. k−2 X. bj φj , i ≤ k − 3.. j=1. Thus, when i ≤ k − 3, hφk−1 , xφi i ai = − =− hφi , φi i. hφk−1 ,. k−2 P. bj φ j i. j=1. hφi , φi i. = 0.. Therefore, (3.7) becomes φk − xφk−1 = −. hxφk−1 , φk−1 i hxφk−1 , φk−2 i φk−1 − φk−2 . hφk−1 , φk−1 i hφk−2 , φk−2 i. It remains now to proof hxφk−1 , φk−2 i = hφk−1 , φk−1 i. Notice that the above equation holds for k − 1, in other words, φk−1 − xφk−2 = −. hxφk−2 , φk−2 i hxφk−2 , φk−3 i φk−2 − φk−3 . hφk−2 , φk−2 i hφk−3 , φk−3 i 17.

(24) Hence, by taking inner product with φk−1 on the above equation both sides, we obtain hφk−1 , xφk−2 i = hφk−1 , φk−1 i. Since hxφk−1 , φk−2 i = hφk−1 , xφk−2 i, the proof is complete. 3.2.2. Five-term recurrence relation Notice that all p-polynomials mentioned in Section 3.1 are not ordinary which. means we can not simply apply the well-know three-term recurrence relation mentioned in Subsection 3.2.1. One might expect Theorem 3.2 is an sufficient condition and perhaps the three-term recurrence relation still works somehow. Unfortunately, if a sequence p-polynomials satisfy the three-term recurrence relation, then 0. 0. 0. 0. 0 = φn = φn−1 + xφn−1 − αn φn−1 − βn φn−2 , which implies φn (1) = φn (−1) = 0. Since all p-polynomials mentioned in section 3.1 do not satisfy the above boundary conditions, they do not satisfy the three-term recurrence relation. It is not a good idea to construct another p-polynomials satisfy above boundary conditions because solution of Turing’s Mechanisms may not satisfy above boundary conditions. However, base on the technics used in the proof of three-term recurrence relation, we construct another recurrence relation for p-polynomials. Lemma 3.3. Let E = [−1, 1] and β ={p ∈ P (E)| p is a polynomial with first derivative being zero at 1 and −1.}. If {φn }∞ n=1 is a sequence of polynomials with first 18.

(25) derivative being zero at 1 and −1 and satify deg(φn ) =. n. n + 1, when n ≥ 2, 0, when n = 1,. (3.8). then φn is a basis of β. Proof. We prove β ⊆ span{φn |n ∈ N}. Let p ∈ β. Without loss of generality, consider deg(p) = k. Hence, we can find ak−1 such that deg(p − ak−1 φk−1 ) ≤ k − 1 and ak−2 such that deg(p − ak−1 φk−1 − ak−2 φk−2 ) ≤ k − 2. Continue this process, we get p − ak−1 φk−1 − ak−2 φk−2 · · · − a1 φ1 = ax2 + bx. Since β is a vector space, ax2 +bx ∈ β. Hence, a = b = 0 and p ∈ span{φn |n ∈ N}. Theorem 3.4 (Five-term recurrence relation). Let {φn }∞ n=1 be a sequence of monic p-polynomials which satisfy n n + 1, n ≥ 1, (a) deg(φn ) = 0, n = 1, (b) φn is odd when n is even, (c) φn is even when n is odd. then φn = gφn−3 − αn φn−2 − βn φn−4 − γn φn−6 , n ≥ 5, φ0 ≡ 0, φ−1 ≡ 0,. 19.

(26) where g = x3 − 3x, hgφn−3 , φn−2 i , hφn−2 , φn−2 i hgφn−3 , φn−4 i βn = , hφn−4 , φn−4 i ( hφn−3 , φn−3 i , n ≥ 7, γn = hφn−6 , φn−6 i 0, n = 5 or n = 6.. αn =. Proof. Let E = [−1, 1] and β = {p ∈ P (E)|p is pig}. Given a large enough constant k . Since both g and φk are monic, deg(φk − gφk−3 ) ≤ k. Since first derivative of gφk−3 being zero at 1 and −1, φk − gφk−3 ∈ β. Thus, by lemma 3.3 we have φk − gφk−3 =. k−1 X. ai φ i ,. (3.9). i=1. taking inner product with φi on both sides of (3.9), we get ai = −. hgφk−3 , φi i hφk−3 , gφi i =− . hφi , φi i hφi , φi i. Notice that deg(gφi ) ≤ k − 3 when i ≤ k − 7, thus, by lemma 3.3 again gφi =. k−4 X. bj φj , i ≤ k − 7.. j=1. Thus, when i ≤ k − 7. hφk−3 , gφi i =− ai = − hφi , φi i 20. hφk−3 ,. k−4 P. b j φj i. j=1. hφi , φi i. = 0..

(27) Therefore, (3.9) becomes φk − gφk−3 =. k−1 X. a i φi .. i=k−6. Recall (b) and (c) in Theorem 3.4, observe that when i = k − 1,k − 3,k − 5 since gφi φk−3 is odd, ai = 0 and hence φk − gφk−3 = −. hgφk−3 , φk−4 i hgφk−3 , φk−6 i hgφk−3 , φk−2 i φk−2 − φk−4 − φk−6 . hφk−2 , φk−2 i hφk−4 , φk−4 i hφk−6 , φk−6 i. Now, the rest is to proof hgφk−3 , φk−6 i = hφk−3 , φk−3 i. Notice that above equation hold for k − 3. In other words, φk−3 − gφk−6 = −. hgφk−6 , φk−7 i hgφk−6 , φk−9 i hgφk−6 , φk−5 i φk−5 − φk−7 − φk−9 . hφk−5 , φk−5 i hφk−7 , φk−7 i hφk−9 , φk−9 i. Hence, by taking inner product with φk−3 both side, we obtain that hgφk−3 , φk−6 i = hφk−3 , gφk−6 i = hφk−3 , φk−3 i. The proof is complete. Note that conditions (b) and (c) are satisfied by Legendre polynomials, Chebyshev polynomials, Taylor series of Trigonometric polynomials, and Most of all, by p-polynomials mention in section 3.1. 3.3. Numerical result about p-polynomials. In this Section, we will show several benefits of choosing piggy polynomials as base functions and some property about piggy polynomials numerically.. 21.

(28) 3.3.1. Rate of convergence We will give several examples to show that when applying SCM on functions. with first derivative being zero at 1 and −1, p-polynomias may be the better choices. Example 3.5. Apply SCM on following functions with standard piggy polynomials, Legendre polynomials and Chebyshev polynomials as base functions. a) cosπx.. b) ex. 3 −3x. .. Let y be the function we want to approximate, {φi }ni=1 are the base functions n P ai φi is the approximation of we chose and {xj }nj=1 are collocation points. If yn = i=1. y created by SCM, then {ai }ni=1  φ1 (x1 ) φ2 (x1 )  φ1 (x2 ) φ2 (x2 )  .. ..  . . φ1 (xn ) φ2 (xn ). are given by ··· ··· .. . ···.   φn (x1 )  a1  y(x1 )   a φn (x2 )  2  y(x )  .  =  .2 .. .   .. . . an φn (xn ) y(xn ).   . . In Figure 3.3.1 we plot the target functions and their approximations using standard p-polynomials, Legendre polynomials and Chebyshev polynomials as base functions. As expected, numerical results in Figure 3.3.1 show that standard p-polynomials are better choices of base functions.. 22.

(29) f (x). The numerical result of cosπx.. 0. −1 −1 −0.8−0.6−0.4−0.2 0 x. 0.2 0.4 0.6 0.8. 1. (a). The numerical result of ex. 3 −3x. .. 8. f (x). 6 4 2 −1 −0.8−0.6−0.4−0.2 0 x. 0.2 0.4 0.6 0.8. 1. (b). Figure 3.1: The solid lines are the target functions, the dashed lines are approximation generated by standard p-polynomials, the circle line represent Legendre polynomials and the x-mark lines are Chebyshev polnomials. In this figure n = 4. 23.

(30) Example 3.6. Apply SCM on the following initial value problem(IVP): 2x dy = 2 − x, y(0) = 0, dx x +1 with standard p-polynomials, Legendre polynomials and Chebyshev polynomials as base functions. Suppose {φi }ni=0 are base functions, unlike Example 3.3.1, we only chose n − n P 1 collocation points {xj }n−1 for the reason of uniquness. If y = ai φi is the n j=1 i=1. {ai }ni=1. approximation generated by SCM, then  0 0 0 φ1 (x1 ) φ2 (x1 ) · · · φn (x1 ) .. .. .. ..  . . . .  0 0 0  φ (x ) φ (x ) · · · φ (x ) 1 n−1 2 n−1 n n−1 φ1 (0) φ2 (0) ··· φn (0). are given by     a1 F (x1 ) ..   a2     .  =  . . ..  F (xn−1 ) an 0. In Figure 3.2 we plot the exact solution of above IVP and numerical solutions computed using SCM. Similar to Figure 3.3.1, Figure 3.2 show that for the above IVP standard p-polynomials are better choice of base functions than Legendre polynomials and Chebyshev polynomials numerically.. 24.

(31) The numerical result of Example 3.6.. f (x). 0.15. 0.1. 5 · 10−2. −1. −0.5. 0 x. 0.5. 1. Figure 3.2: The solid lines are the target functions, the dashed lines are approximation generated by standard p-polynomials, the circle line represent Legendre polynomials and the x-mark lines are Chebyshev polnomials. In this figure n = 4 and collocation points are Legendre-Gauss points.. 25.

(32) 3.3.2. Boundary value problem(BVP) By giving the following example, we will show that choosing p-polynomials as. base functions may bring great benefit on boundary value problems. Example 3.7. Use SCM to find eigenvalues of the following BVP. d2 y 0 0 + λy = 0, y (0) = y (π) = 0. 2 dx. (3.10). Before applying SCM on (3.10), since first derivative being zero y may not π be zero at 1 and −1, we change variable first. Let x(t) = (t + 1) and consider 2 u(t) = y(x(t)) then (3.10) can be expressed by u(t) 2 d2 u 0 0 ( )2 2 + λu = 0, u (−1) = u (1) = 0. π dt Write un =. n P. ai φi as the approximation generated by SCM. If we chose some other. i=1. functions as base functions, then the coefficients {ai }ni=1 satisfy 00 λΦ ( π2 )2 Φ ~a = ~a, 0 B where  00. Φ. =. Φ. =. B. =. ~a. =. 00. φ1 (x1 )  φ001 (x2 )  ..  . 00  φ1 (xn−2 ) φ1 (x1 )  φ1 (x2 )  ..  . φ (x ) 1 n−2 0 φ1 (1) 0 φ1 (−1) ( a1.  00 00 φ2 (x1 ) · · · φn (x1 ) 00 00 φ2 (x2 ) · · · φn (x2 )  , .. .. ..  . . . 00 00 φ2 (xn−2 ) · · · φn (xn−2 )  φ2 (x1 ) · · · φn (x1 ) φ2 (x2 ) · · · φn (x2 )  , .. .. ..  . . . φ2 (xn−2 ) · · · φn (xn−2 ) 0 0 φ2 (1) · · · φn (1) , 0 0 φ2 (−1) · · · φn (−1) a2 · · · · · · an )> . 26. (3.11).

(33) Unfortunately, in (3.11) we can’t solve λ directly. However, if we chose piggy polynomials as base functions, the coefficients {ai }ni=1 satisfy 2 00 ( )2 Φ ~a = λΦ~a, π where  00. Φ.  =   . Φ.  =  . ~a. =. 00. φ1 (x1 ) 00 φ1 (x2 ) .. . 00 φ1 (xn ) φ1 (x1 ) φ1 (x2 ) .. . φ1 (xn ) ( a1. 00. φ2 (x1 ) 00 φ2 (x2 ) .. . 00 φ2 (xn ) φ2 (x1 ) φ2 (x2 ) .. . φ2 (xn ) a2 · · ·. ··· ··· .. . ··· ··· ··· .. . ··· ···. (3.12). 00. φn (x1 ) 00 φn (x2 ) .. . 00 φn (xn ) φn (x1 ) φn (x2 ) .. . φn (xn ) an )> ,.   ,    , . which we can solve directly for λ by solving a generalized eigenvalue problem of Φ. 00. and Φ. In Figure 3.3, we plot relative errors in log scale for eigenvalues λ = 1, 4 and 9 using standard p- polynomials as base functions. Note that we did not plot the case 00. λ = 0 for the reason that Φ in (3.12) is not full rank. Numerical result in Figure 3 show that. |λ−λn | |λ|. decays exponentially to zero for λ = 1, 4 and 9. Notice that when. choosing Legendre p-polynomials or Chebyshev p-polynomials as base functions, the numerical result are almost identical.. 27.

(34) n| ) log10 ( |λ−λ |λ|. 0. −5. −10. −15. 6. 7. 8. 9 n. 10. 11. 12. Figure 3.3: The square-line is λ = 1, the circle-line is λ = 4 and the cross-lines is λ = 9. In this figure, we use standard p-polynomials as base functions and collocation points are Legendre-Gauss points.. 28.

(35) 3.4. Standard p-polynomials, Legendre p-polynomials and Chebyshev p-polynomials are orthogonal basis of L2 ([−1, 1]) There are other familiar functions that satisfy zero flux boundary condition,. such as cos(πx), but are not used as base functions in SCM here. The reason for this is that suppose {φn }∞ n=1 are base functions we choose and f is the function we want ∞ P to approximate, then f may not equal to a n φn . n=1. To avoid this kind of problem, in this subsection we will show that the ppolynomials we have introduced above are are orthogonal basis of L2 ([−1, 1]). Lemma 3.8. Let E = [−1, 1] and β ={p ∈ P (E)| p is a polynomial with first derivative being zero at 1 and −1.}, then β is dense in L2 (E). Proof. Observe that β satisfies the following: 0. 0. 0. 0. 0. a) β is an algebra; let f, g ∈ β and α ∈ R, since (f + g) = f + g , (αf ) = αf and 0. 0. 0. (f · g) = f · g + f · g, all f + g, αf and f · g ∈ β. b) The constant function x 7→ 1 lies in β. c) The strictly increasing function f (x) = x3 − 3x, x ∈ E lies in β. Since E is a compact set, by Stone-Weierstrass Theorem β is dense in C(E, R) together with the supnorm. Notice that E has finite measure, thus β is dense in L2 (E). Lemma 3.9. Let H be a Hilbert space and β a total subspace of H, in other words, β is a subspace of H and cl(β) = H. If β has a countable orthonormal Hamel basis {un |n ∈ N}, then {un |n ∈ N} is an orthonormal basis of H.. 29.

(36) Proof. For a Hilbert space H, it is known that if {en |n ∈} is an orthonormal set of H, then {en |n ∈ N} is an orthonormal basis of H and hx, en i = 0 ∀n ∈ N implies x = 0, where x ∈ H are equivalent. On the other hand, if β is a subspace of H, then an orthonormal set of β is an orthonormal set of H. Let {un |n ∈ N} be an orthonormal Hamel basis of β. To prove {un |n ∈ N} is an orthonormal basis of H we show that hx, un i = 0 ∀n ∈ N implies x = 0, x ∈ H. Let x ∈ H, assume hx, un i = 0 ∀n ∈ N. Given > 0, since cl(β) = H, ∃ y ∈ β such that kx − yk < . Thus, hx, xi = = ≤ <. hx − y + y, xi hx − y, xi + hy, xi kx − ykkxk + hy, xi kxk + hy, xi.. Notice that y ∈ β, since {un |n ∈ N} is a Hamel basis, y is a finite linear combination of vector in {un |n ∈ N}. Therefore, hx, xi < kxk. Since is arbitrary, hx, xi = 0 and hence x = 0. Theorem 3.10. E = [−1, 1] and β ={p ∈ P (E)| p is a polynomial with first derivative being zero at 1 and −1.}, then an orthonormal Hemal basis of β is an orthonormal basis of L2 (E). Proof. Since L2 (E) is a Hilbert space, by Lemma 3.8 β is a total subspace of L2 (E). Since β is a subspace of P (E), any orthonormal Hamel basis of β is countable. Thus, by Lemma 3.9 an orthonormal Hamel basis of β is an orthonormal basis of L2 (E). Lemma 3.3 proves that the standard p-polynomials, Legendre p-polynomials 30.

(37) and Chebyshev p-polynomials are orthogonal Hemal basis of β in Theorem 3.10, consequently all three p-polynomials are orthogonal basis of L2 (E).. 31.

(38) CHAPTER 4 Numerical results In this section, we will verify several assumptions based on the linear analysis in Chapter 2 numerically by considering one-dimensional and two-dimensional dimensionless case of equation (1.1). The one-dimensional case is same as equation (2.2) in Chapter 2. For the two-dimensional case, consider the following system ut = 52 u + γf (u, v) , (x, y) ∈ [0, 2] × [0, π], (4.1) vt = d 52 v + γg(u, v) and zero flux boundary conditions gives ux = vx = uy = vy = 0 on the boundary. Similar to linear analysis in Chapter 2, solutions of linearise system of Equation (4.1) are linear combination of the following form uˆ ∝ eλt cos(kx x)cos(ky y), vˆ ∝ eλt cos(kx x)cos(ky y),. (4.2). iπ , i ∈ N and ky = j, j ∈ N. Furthermore, for each pair of wavenumbers 2 (kx , ky ), there is an interval γ− < γ < γ+ such that the above pair of wavenumbers is where kx =. permissible, where p. (dfu + gv )2 − 4d(fu gv − fv gu ) ), p fu gv − fv gu (dfu + gv ) + (dfu + gv )2 − 4d(fu gv − fv gu ) γ+ = (kx2 + ky2 )( ) fu gv − fv gu. γ− = (kx2 + ky2 )(. (dfu + gv ) −. (4.3). Notice that we choose the above rectangle domain rather than square domain so that the bifurcation points remain simple singular points; see [4] for details. 32.

(39) 4.1. Numerical methods. In this subsection we will show how to use p-polynomials discussed in Section 3 as base functions for two-dimensional reaction diffusion equations. First, we change π y + 1), and dropping the hats for variables in equations (4.1) . Let x = xˆ + 1, y = (ˆ 2 convenience. Equation (4.1) becomes    ut = uxx + ( 2 )2 uyy + γf (u, v) π , (x, y) ∈ [−1, 1] × [−1, 1], (4.4)  v = d(v + ( 2 )2 v + γg(u, v)  t xx yy π and the zero flux boundary conditions give ux = vx = uy = vy = 0 on the boundary. To Apply SCM on equation (4.4), we assume un =. n X. αij (t)φi (x)φj (y), vn =. i=1,j=1. n X. βij (t)φi (x)φj (y),. (4.5). i=1,j=1. and get  . α ~0 β~ 0. .  (Φy ⊗ =. Φ00x.  2 2 00 + ( ) Φy ⊗ Φx 0 ~  α π  ~ 2 β 0 d (Φy ⊗ Φ00x + ( )2 Φ00y ⊗ Φx π (4.6). +γ. ~ ~a − Φy ⊗ Φx α ~ + (Φy ⊗ Φx α ~ ) ◦ (Φy ⊗ Φx α ~ ) ◦ (Φy ⊗ Φx β) ~b − (Φy ⊗ Φx α ~ ~ ) ◦ (Φy ⊗ Φx α ~ ) ◦ (Φy ⊗ Φx β). ,. where ⊗ is the tensor product and ◦ denote the element-by-element multiplication; a c ac for example ◦ = b d bd , {(xi , yj )|1 ≤ i ≤ n, 1 ≤ j ≤ n.} are the collocation points, 33.

(40) αnn )> , β~ = ( β11 β21 · · · βnn )> , > dα21 (t) dαnn (t) dβ11 (t) dβ21 (t) 0 ~ ,β = ··· dt dt dt dt   φ1 (y1 ) φ2 (y1 ) · · · φ2 (x1 ) · · · φn (x1 )   φ (y ) φ2 (y2 ) · · · φ2 (x2 ) · · · φn (x2 ) , Φy =  1 . 2 .. .. .. .. ..   .. . . . . . φ1 (yn ) φ2 (yn ) · · · φ2 (xn ) · · · φn (xn )   00 00 00 00 φ2 (x1 ) · · · φn (x1 ) φ1 (y1 ) φ2 (y1 ) · · · 00 00 00  φ1 (y2 ) φ002 (y2 ) · · · φ2 (x2 ) · · · φn (x2 )  00  , Φy =  .. .. .. .. .. ..   . . . . . . 00 00 00 00 φ1 (yn ) φ2 (yn ) · · · φ2 (xn ) · · · φn (xn ) ~a = ( a a · · · a )> , ~b = ( b b · · · b )> . α ~ = ( α11 α21 dα11 (t) 0 α ~ = dt  φ1 (x1 )  φ1 (x2 ) Φx =  ..  . φ1 (xn )  00 φ1 (x1 )  φ001 (x2 ) 00 Φx =  ..  . 00 φ1 (xn ). ···. dβnn (t) ··· dt  φn (y1 ) φn (y2 )  , ..  . φn (yn )  00 φn (y1 ) 00 φn (y2 )  , ..  . 00 φn (yn ). For the one-dimensional case, we assume un =. n X. αi (t)φi (x), vn =. i=1. n X. βi (t)φi (x),. (4.7). i=1. Since our main concern is on how γ effects the final patterns, we will focus on the state solutions of equation (4.1), namely 0 = 52 u + γf (u, v) , (x, y) ∈ [0, L1 ] × [0, L2 ], 0 = d 52 u + γg(u, v). (4.8). and the zero flux boundary conditions give ux = vx = uy = vy = 0 on the boundary. We rescale x and y as before and equation (4.8) becomes  2 2   0 = ( )2 uxx + ( )2 uyy + γf (u, v) L1 L2 , (x, y) ∈ [−1, 1] × [−1, 1], 2 2 2   0 = d ( ) vxx + ( )2 vyy + γg(u, v) L1 L2 34. (4.9). > ,.

(41) with zero flux boundary conditions. To apply SCM on equation (4.9) we assume un =. n X. αij φi (x)φj (y), vn =. i=1,j=1. n X. βij φi (x)φj (y),. (4.10). i=1,j=1. and get  2 2 00 0 + ( ) Φy ⊗ Φx ~  Φy ⊗  α π 0=  ~ 2 β 0 d Φy ⊗ Φ00x + ( )2 Φ00y ⊗ Φx π ~ ~a − Φy ⊗ Φx α ~ + (Φy ⊗ Φx α ~ ) ◦ (Φy ⊗ Φx α ~ ) ◦ (Φy ⊗ Φx β) , +γ ~b − (Φy ⊗ Φx α ~ ~ ) ◦ (Φy ⊗ Φx α ~ ) ◦ (Φy ⊗ Φx β) . Φ00x. (4.11). the notations are similar to those in equation (4.6) and are omitted here. For the one-dimensional case, we assume un =. n X. αi φi (x), vn =. i=1. 4.2. n X. βi φi (x),. (4.12). i=1. Only one wavenumber is permissible. In this subsection we will show that when there is only one permissible wavenumber, we can use it to predict the waveform and the occurrence of heterogeneous patterns generated by reaction diffusion equations. For the one-dimensional case, consider a = 0.1, b = 1.5, d = 19.2. The above parameters are chosen so that the intervals (γk+ , γk− ) in equation (2.8) do not overlap. In Figure 4.1 the notation λk represents the eigenvalue of B(k 2 ) from equation (2.6). Figure 4.1 show that (γk+ , γk− ) corresponding to π, 2π, 3π, 4π are non-overlapping intervals. Hence, the above wavenumbers are the only permissible wavenumber when γ in corresponding (γk+ , γk− ). 35.

(42) 0.6. 0.5. 0.4. k. 0.3. 0.2. 0.1. 0 0. 10. 20. 30. 40. 50. 60. Figure 4.1: In this Figure we plot eigenvalues of B(k 2 ) from equation (2.6) for onedimensional case; only positive λk is shown here. From left to right, k = π, 2π, 3π, 4π.. 36.

(43) In the upper half of Figure 4.2 u0n represents un (0) =. n P. αi φi (0) in equation. i=1. (4.12), in other words, we plot the approximate steady state solutions of equation (2.2) at x = 0 as γ changes, the solutions are generated by an arclength continuation method. The bottom of Figure 4.2 is similar to Figure 4.1. Figure 4.2 show that wavenumber can predict the occurrence heterogeneous patterns precisely. In Figure 4.3 we plot the approximate steady state solutions of equation (2.2) un as γ changes. Figure 4.3 show that wavenumbers can predict the waveform of patterns.. 37.

(44) 0 1.8 n. u. 1.6 1.4 1.2 0. 10. 20. 30. 40. 50. 30. 40. 50. (a). 0.6. k 0.4 0.2 0 0. 10. 20. (b). Figure 4.2: In (a) we plot un (0) as given by equation (4.12) for one-dimensional case as γ changes. The solid lines represent stable steady states. The dashed lines represent unstable steady states. In (b) we plot eigenvalues of B(k 2 ) from equation (2.6).. 38.

(45) 2. 2. 1.8. 1.8. u. u. n 1.6. n 1.6. 1.4. 1.4. 3.5. 14 1. 3 0. 2.5 2. −0.5. 1. 12. 0.5. 0.5. x. −1. 0. 10 8. (a). −0.5. x. −1. (b). 2. 2. 1.8. 1.8. u. u. n 1.6. n 1.6. 1.4. 1.4. 30. 60 1 0.5. 25. 0 20. −0.5 −1. 1. 50. 0.5 0. 40. x. 30. (c). −0.5 −1. x. (d). Figure 4.3: In this figure we present one dimensional spatial patterns given by equation (2.2) for one-dimensional case as γ changes. Waveform of patterns basically looks like cos(kx).. 39.

(46) For the two-dimensional case we choose a = 0.2, b = 1.2, d = 20.55 for the same reasom in one-dimensional case. The λk in Figure 4.4 is similar to the λk in Figure 4.1. Figure 4.4 and table 4.1 show that (γk+ , γk− ) in equation (4.3) corresponding to π π , 0 , (0, 1), , 1 , (π, 0), (π, 1), (π, 2) and (2π, 0) are non-overlapping (kx , ky ) = 2 2 intervals. In two-dimensional case we will focus on the longer intervals corresponding to (kx , ky ) = (π, 2) and (kx , ky ) = (2π, 0) for the ease of observation. 0.1 0.09 0.08 0.07 0.06. k. 0.05 0.04 0.03 0.02 0.01 0 0. 20. 40. 60. 80. 100. Figure 4.4: In this Figure we plot eigenvalues of B(k 2 ) from linear analysis for twodimensional case; only positive λk is shown here. We will focus on (γk+ , γk− ) = (42.1477, 47.8531) and (48.6217, 55.2034). In the upper half of Figure 4.5 u0,0 n represents un (0, 0) in equation (4.10). For the bottom of Figure 4.5 only λk corresponding to (kx , ky ) = (π, 2) and (2π, 0) are shown here. The numerical result in the upper half of Figure 4.5 is generated 40.

(47) pair ofwavenumbers interval of γ pair of wavenumbers π ,0 . 3.0389 ∼ 3.4502. (0, 1). 2 π ,1 . 10.5369 ∼ 11.9633. (π, 0). 2 3π (π, 0). 19.6535 ∼ 22.3139. ,0 . 2 π (0, 2). 29.9923 ∼ 34.0522. ,2 . 2 3π 34.8478 ∼ 39.5650. (π, 2). ,1 . 2 48.6217 ∼ 55.2034. (2π, 1). (2π, 0). 3π ,2 . 57.3420 ∼ 65.1041. (0, 3). 2 π 5π ,3 . 70.5215 ∼ 80.0677. ,0 . 2 2. interval of γ 7.4981 ∼ 8.5131. 12.1554 ∼ 13.8009. 27.3497 ∼ 31.0519. 33.0311 ∼ 37.5021. 42.1477 ∼ 47.8531. 56.1197 ∼ 63.7165. 67.4827 ∼ 76.6175. 75.9711 ∼ 86.2553.. Table 4.1: Wavenumbers and corresponding (γk+ , γk− ). by an arclength continuation method. Figure 4.5 show that pair of wave numbers can predict the occurrence heterogeneous patterns in two-dimensional case. Notice that in the two-dimensional case our patterns are 3D graphs, thus, we only present several graphs of fix γ here. Figure 4.6 are patterns generated form equation (4.1) choosing γ = 45 and γ = 51, in left permissible pair of wavenumbers is (π, 2) in right permissible pair of wavenumbers is (2π, 0). As expect, waveform of Figure 4.11(a) is the same with cos(πx)cos(2y) and waveform of Figure 4.11(b) is the same with cos(2πx).. 41.

(48) 1.6. 0,0 un. 1.4 1.2 40. 45. 50. 55. 50. 55. (a). 0.1. k. 0.05 0 40. 45. (b). Figure 4.5: In (a) we plot un (0, 0) as given by equation (4.10) for two-dimensional case as γ changes. The solid lines represent stable steady states. The dashed lines represent unstable steady states. We focus on (kx , ky ) = (π, 2) and (2π, 0). In (b) we plot eigenvalues of B(k 2 ) from linear analysis.. 42.

(49) =45. 1.6 1.5 1.4 1.3. 4 2 1.5. 2 y. 1 0. 0.5 0. x. (a) =51. 1.5 1.45 1.4 1.35 1.3 4 2 1.5. 2 y. 1 0. 0.5 0. x. (b). Figure 4.6: Spatial patterns given by equation (4.1) for two-dimensional case at γ = 45 which corresponding to (kx , ky ) = (π, 2) and γ = 51 which corresponding to (kx , ky ) = (2π, 0). 43.

(50) 4.3. Multiple permissible wavenumbers. When multiple wavenumbers are permissible, we expect the waveform to be determined the corresponding wavenumber of the largest λ in linear analysis. However, the following numerical results show that when we choose γ such that multiple wavenumbers are permissible, different waveforms are predicted. For the one-dimensional case, Consider a = 0.1, b = 1.2, d = 14.5. We chose this parameters to make sure that only two wavenumbers are permissible. In Figure 4.7 the notation λk represents the eigenvalue of B(k 2 ) from equation (2.6). Figure 4.7 show that (γk+ , γk− ) corresponding to 3π and 4π are overlap with each other. Thus, we focus on the above two intervals. In the upper half of Figure 4.8 n P u0n represents un (0) = αi φi (0) from equation (4.12), Figure 4.8 show that largest i=1. λ may not determine the final pattern. Moreover, Figure 4.9 is a counterexample of our assumption.. 44.

(51) 2.5. 2. k. 1.5. 1. 0.5. 0 0. 10. 20. 30. 40. 50. 60. Figure 4.7: In this Figure we plot eigenvalues of B(k 2 ) from equation (2.6) for onedimensional case; only positive λk is shown here. We will focus on (γk+ , γk− ) = (18.2422, 38.0969) and (32.4305, 67.7279).. 45.

(52) 2. 0 u 1.5 n 1 0.5. 20. 25. 30. 35. 40. 30. 35. 40. (a). k. 2 1 0. 20. 25. (b). Figure 4.8: In (a) we plot un (0) as given by equation (4.12) for one-dimensional case as γ changes. The solid lines represent stable steady states. The dashed lines represent unstable steady states. The circle represent bifurcation points. In (b) we plot eigenvalues of B(k 2 ) from equation (2.6) for one-dimensional case.. 46.

(53) =35.5 1.8 1.7 1.6 1.5. un. 1.4 1.3 1.2 1.1 1 0.9 −1. −0.5. 0. x. 0.5. 1. 0.5. 1. (a) =35.5 1.7 1.6 1.5. un. 1.4 1.3 1.2 1.1 1 −1. −0.5. 0. x (b). Figure 4.9: Two different spatial patterns given by equation (2.2) for one-dimensional case at γ = 35.5.. 47.

(54) For the two-dimensional case, we choose a = 0.2, b = 1.5, d = 27.1 to ensure only two pair of wavenumbers are permissible. In Figure 4.10 the notation λk represents the eigenvalue of B(k 2 ) from linear analysis. By Figure 4.10 we chose π , 1 and (kx , ky ) = (0, 2). In the top of Figure 4.11 the notation un0,0 (kx , ky ) = 2 represents un (0, 0) in equation (4.10). Like Figure 4.8 in one-dimensional case, Figure 4.11 show that largest λ may not determine the final pattern, too. As expect, Figure π 4.12 show that both , 1 and (0, 2) are predicted. Figure 4.13 is an example show 2 that on the same curve, patterns have the same waveform.. 0.35 0.3 0.25. k 0.2 0.15 0.1 0.05 0 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. Figure 4.10: In this Figure we plot eigenvalues of B(k 2 ) from linear analysis for two-dimensional case; only positive λk is shown here. We will focus on (γk+ , γk− ) = (9.3896, 12.0069) and (10.8319, 13.8512).. 48.

(55) 2. 0,0 u 1.8 n 1.6 1.4 9. 10. 11. 12. 13. 14. 12. 13. 14. (a) 0.3. k 0.2 0.1 0 9. 10. 11. (b). Figure 4.11: In (a) we plot un (0, 0) as given by equation (4.10) for two-dimensional case as γ changes. The solid lines represent stable steady state. The solid lines represent stable steady states. The dashed lines represent unstable steady states. The circle represent bifurcation points. In (b) we plot eigenvalues of B(k 2 ) from linear analysis for two-dimensional case.. 49.

(56) =11.1. 1.9 1.8 1.7 1.6 1.5 1.4 4 2 1.5. 2 y. 1 0. 0.5 0. x. (a) =11.1. 1.9 1.8 1.7 1.6 1.5 4 2 1.5. 2 y. 1 0. 0.5 0. x. (b). Figure 4.12:. Two different spatial patterns given by equation (4.1) for two-. dimensional case at γ = 11.1.. 50.

(57) =10. 2. 1.8. 1.6. 1.4 4 2 1.5. 2 y. 1 0. 0.5 0. x. (a) =13. 1.9 1.8 1.7 1.6 1.5 4 2 1.5. 2 y. 1 0. 0.5 0. x. (b). Figure 4.13: Spatial patterns given by equation (4.1) for two-dimensional case at γ = 10 and γ = 13.. 51.

(58) CHAPTER 5 Conclusion In this thesis we derived a formulation that describes the relation between domain size and permissible wavenumbers. The analysis on the discretized system in Chapter 2 and numerical results in Chapter 4 confirm the formulation about the occurrence of heterogeneous patterns. For the final spatial pattern when there is only one permissible wavenumber prediction from the formulation is reliable. In the future, we may apply similar analysis to more complex domains. Numerical results in Chapter 4 are generated using an arclength continuation scheme choosing γ as the varying parameter. Technically, we can choose any other parameters in Turing’s reaction diffusion equation as the parameter in the arclength continuation scheme. Moreover, bifurcation test in Chapter 2 can also be applied on uniform steady state with different parameter from the reaction diffusion equations. Thus, as well as considering different domain, we can also discuss the influence of any other parameters in the equations to the final patterns. In addition to discussing the final patterns generated from Turing’s equations, we also improved our numerical method by considering p-polynomials in Chapter 3 as base functions and proposed a five-term recurrence relation to reduce the cost of computation. However, p-polynomials generated from a five-term recurrence relation in Chapter 3 may have error problems when n is large. In the future we may provide an alternative way to compute the coefficients in the five-term recurrence relation.. 52.

(59) REFERENCES [1] A. M. Turing, The Chemical Basis of Morphogenesis. Philos. Trans. R. Soc. London, Ser. B237, 37(1952). [2] J. Shen, Efficient Spectral-Galerkin method I. Direct solvers of second and forth order equations using Legendre Polynomials, SIAM J. Sci. Compute. 15(1995) 74-87. [3] J. Shen, Efficient Spectral-Galerkin method II. Direct solvers of second and forth order equations using Chebyshev Polynomials, SIAM J. Sci. Compute. 15(1995) 74-87. [4] H.B. Keller, Lectures on numerical methods in bifurcation problems.(Springer, Berlin, 1987) [5] J. D. Murray, Mathematical Biology (Springer, Berlin, 1990). [6] Kronecker Products and Matrix Calculus: with Application, Alexander Graham. [7] Edmund J. Crampin, Eamonn A. Gaffney and Philip K. Maini, Reaction and Diffusion on growing Domains: Scenarios for Robust Pattern Formation(1999). [8] J.C. Eilbeck, The pseudo-spectral method and path following in reaction-diffusion bifurcation studies, SIAM J. Stat. Comput. 7(1986) 599-610. [9] D. Gottlieb, M.Y. Hussaini, S.A Orszag, Theory and application of spectral methods, in: R.G. Voigt, D. Gottlieb, M.Y. Hussaini (Eds.), Spectral Methods fo Partial Differential Equations, SIAM Publications, Philadelphia, 1984, pp. 1-54. 53.

(60)