National University of Kaohsiung Repository System:Item 310360000Q/10828

全文

(1)國立高雄大學應用數學系碩士班碩士論文. 非線性偶合薛丁格方程數值解的行為 Numerical Solution Profiles of the Coupled Nonlinear Schrödinger Equations. 研究生：劉建祥撰指導教授：王偉仲. 中華民國九十七年七月.

(2) Numerical Solution Profiles of the Coupled Nonlinear Schrödinger Equations. by Chien-Hsiang Liu Advisor Weichung Wang. Department of Applied Mathematics, National University of Kaohsiung Kaohsiung, Taiwan 811, R.O.C. July 2008.

(3) Contents 中文摘要. ii. 英文摘要. iii. 1. Introduction. 1. 2. Spike Solution of 2-Coupled Schrödinger System. 2. 3. The Hyperplane-Constrained Continuation Method. 19. 4. The course of the discrete Laplacian 4.1 The Rectangle with Dirichlet Boundary Conditions of Laplacian 4.2 The Disc with Dirichlet Boundary Conditions of Laplacian 4.3 The Disc with Neumann Boundary Conditions of Laplacian 4.4 The Disc with Robin Boundary Conditions of Laplacian 4.5 The Ring with Dirichlet Boundary Conditions of Laplacian 4.6 The Ring with Neumann Boundary Conditions of Laplacian 4.7 The Ring with Robin Boundary Conditions of Laplacian. 36 36 39 42 44 46 49 51. Appendix. 53. References. 89. i.

(4) 非線性偶合薛丁格方程數值解的行為（標楷體 18 號字.置中.粗體）指導教授：王偉仲教授國立台灣大學數學系（標楷體 12 號字，置中，如有兩個以上則如上面格式順序書列）學生：劉建祥國立高雄大學應用數學系碩士班（標楷體 12 號字，置中）摘要（以下皆為標楷體 12 號字，左右對齊）本篇論文利用一些演算法，藉著數值計算，然後模擬非線性偶合薛丁格方程數值解的行為。我們用的演算法分別為固定點迭代法 ( Fixed Point Method ) 以及加平面限制延拓法 ( Hyperplane-Constrained Continuation Method )。在介紹這兩個演算法的精神與過程後，會展示出一些數值實驗的簡介以及結果。最後再提供一些不同定義域跟不同邊界條件下，二階微分的 Laplacian 用有限差分法，離散化成一個矩陣的過程。. 關鍵字：非線性偶合薛丁格方程、固定點迭代法、加平面限制延拓法、有限差分法 (關鍵字為粗體，注意縮排). (本論文摘要上下左右邊緣留白一吋，即 2.54 公分). ii.

(5) Numerical Solution Profiles of the Coupled Nonlinear Schrödinger Equations （Times New Roman 體 18 號字，置中，粗體） Advisor: Dr. Weichung Wang Department of Mathematics National Taiwan University. Student: Chien-Hsiang Liu Department of Applied Mathematics National University of Kaohsiung （Times New Roman 體 12 號字，置中） Abstract In this thesis, we use some algorithms to simulate the solutions type of coupled nonlinear Schrödinger equations. The main algorithms used by us are the Fixed Point Method and the Hyperplane-Constrained Continuation Method. We will summarize the two methods and present the results of numerical experiments designed to showcase some of the more interesting behavior. Finally we arrange the courses of the discrete Laplacian using the Finite Different Method for specific domains and boundary conditions.. Keywords: Coupled nonlinear Schrödinger equations, Fixed Point Method, Hyperplane-Constrained Continuation Method, Finite Different Method (Keywords 為粗體，注意縮排). (本論文摘要上下左右邊緣留白一吋，即 2.54 公分). iii.

(6) 1. Introduction In this thesis, we consider that the numerical solution profiles of following two couple nonlinear Schr¨odinger equation:  2 3 2   ε ∆u − λ1 u + µ1 u + βuv = 0 in Ω, (1) ε2 ∆v − λ2 v + µ2 v 3 + βu2 v = 0 in Ω,   u=v=0 on ∂Ω, where the domain Ω is a subset of RN for N = 1, 2, 3 with parameters λ1 , λ2 , µ1 , µ2 > 0 and β < 0. The variables (ε, λ, µ, and β) will affect the behavior of the ground state solutions. From previous numerical experiments, we find that ε is associated with the breadth of the ground state solutions, µ1 and µ2 have been found to be associated with the height of the ground state solutions, while λ1 and λ2 seem to be associated with the breadth and height of the ground state solutions. Finally, β will affect the attractive or repulsive power between the ground state solutions. When β > 0, the action is attractive. Otherwise, β < 0 represents that the action is repulsive and visa versa. In Section 2, we will introduce the Fixed Point Method used here to compute the solutions to the Equations (1).We then perform some specifically designed experiments. Using Matlab, we will plot various profiles ( highlighting various variable dependencies ) associated with a given numerical solution. We also calculate the energy of the solutions obtained throughout the course of the simulation. We will conclude Section 2 with a discussion of these results. Scale the variables u¯(x) = u(εx) and v¯(x) = v(εx). After reduction of Equation (1), we obtain the following:   u − λ1 u¯ + µ1 u¯3 + β u¯v¯2 = 0 in Ωε ,  ∆¯ (2) ∆¯ v − λ2 v¯ + µ2 v¯3 + β u¯2 v¯ = 0 in Ωε ,   u¯ = v¯ = 0 on ∂Ωε , where Ωε = {x ∈ RN | εx ∈ Ω}. Our model changes to Equation (2), and we will introduce the Hyperplane-Constrained Continuation Method in step-by-step detail in Section 3. We will then use this method to approximate the form of the solutions and plot the corresponding shapes and bifurcations. ∆u can be used as a finite difference in some special domain ( corresponding to a specific metric ) whether using the Fixed Point Method presented in Section 2 or the Hyperplane-Constrained Continuation Method presented in Section 3. The most useful domains for the work here are rectangular and disc shaped. We also work with a hole in the domain i.e. a hole is home to a rectangle or a disc. In Section 4 we describe the inference process that maps the given domain to a matrix for each case.. 1.

(7) 2. Spike Solution of 2-Coupled Schr¨ odinger System Previous work, [8], has shown that the shape of the two ground state solutions will become spike-like if the variable ε2 goes to zero. In this vein, we have designed and carried out a series of numerical experiments. A significant part of my work has, therefore, involved comparing and arranging the results from such numerical experiments. Here we will briefly introduce our method. The domain that we use here will be a disc centered at the origin. Application of Finite Difference Methodology leads to a matrix that has seven diagonal lines. The λ term and the β term are placed into main diagonal line and the µ term is moved to the right hand side of Equation (1). We must now solve a linear system: Ax = b. We Use the Fixed Point Method starting at an initial point on the disc that is defined by the user. If we are lucky enough to obtain a convergent solution from such a choice, we output the solution to the user. Because there are many solutions almost everywhere in the disc, we can boldly use the iterative to converge the numerical solutions from the initial models. One great advantage of using this method is that it is easier to compare cases that differ in the value of one particular variable. As we are also calculating the energy of the solutions during the process, we may hope to determine relations between the variables and the energy. Prior work has given us an explicit knowledge of the energy functional of the system. According to [8] we have: Z Z Z λ1 µ1 ε2 2 2 |∇u| + u − u4 EΩ (u, v) := 2 Ω 2 Ω 4 Ω Z Z Z ε2 λ2 µ2 2 2 + |∇v| + v − v4 (3) 2 Ω 2 Ω 4 Ω Z β − u2 v 2 , 2 Ω for u, v ∈ H01 (Ω). We can then compute the integrals by numerical integration. Suppose that we cute the disc into N radial sections and into M azimuthal sections. Then ∆r = 2N2+1 and ∆θ = 2π . Using the Finite Difference Method on M the disc, we find that ∆u is given as £u as (4), where £ is a matrix whose size is R(N ×M )×(N ×M ) . Further details for tuning the Laplacian into a matrix operator are presented in Section 4.  2 3 2   ε £u − λ1 u + µ1 u + βuv = 0 ε2 £v − λ2 v + µ2 v 3 + βu2 v = 0   u=v=0. 2. in Ω. in Ω. on ∂Ω.. (4).

(8) The discretization of a disc for N = 3 and M = 8. 1 0.8. 19. 0.6. 18 20. 11. 0.4. 10 12. 0.2. 3. y value. 4 0. 17. 8 6. −0.2. 9. 1. 5. 13. 21. 2. 7 16. −0.4. 14 15. −0.6. 24. 22 23. −0.8 −1 −1. Figure 2.1.. −0.5. 0 x value. 0.5. 1. It is the discretization of the Laplacian over a disc. The radial direction is cut into N parts and the angular direction is cut into M parts where N = 3 and M = 8.. Because the angular direction is cut into slices of equal width, we can get away with only describing a sector of angular width ∆θ. We can label variables according to each radial section of width ∆r as follows: d21,j , d22,j , . . . , d2N,j , ∀j = 1, 2, . . . , M . In fact, d21,1 = d21,2 = · · · = d21,M , d22,1 = d22,2 = · · · = d22,M , . . . , d2N,1 = d2N,2 = · · · = d2N,M . We define the variables as squared in anticipation of the simplified version of the equations. The image of d21,1, d22,1, and d23,1. The blue *s are grid points. 0.25 0.2 0.15. y value. 0.1 0.05 d21,1. 0. d22,1. d23,1. −0.05 −0.1 −0.15 −0.2 −0.25 0. 0.2. 0.4. 0.6. 0.8. 1. x value. Figure 2.2.. The image of d21,1 , d22,1 , and d23,1 . The red area is d21,1 , the green area is d22,1 , and the pink area is d23,1 . The other d2i,j are treated similarly for all i and j.. 3.

(9) For i = 1, . . . , N and j = 1, . . . , M , the area of d2i,j is calculated as d2i,j = (π(i∆r)2 − π((i − 1)∆r)2 ). ∆θ 1 = (2i − 1)∆θ∆r2 . 2π 2. (5). In fact, the area is independent of j. Assume that £ = D−1 AD, i.e. A = D£D−1 where     D12 d2i,1         D22 d2i,2 2 N M ×N M 2     ∈ RM 2 . D = ∈R , where Di =  .. ..   . .     2 2 DN di,M Then the matrix A is symmetric and negative by definition. Since Z. Z. −. Z. u∆u = Ω. Z. 2. |∇u| + Ω. u∇u ∂Ω. and. u∇u = 0,. (6). ∂Ω. we can use the discretization of the Laplacian to get the first term of the energy equation (3) i.e. Z Z 2 |∇u| = − u∆u ≈ uT D2 £u = uT D2 D−1 ADu = uT DADu. Ω. Ω. Therefore, we define E : RN × RM → R by ε2 T λ1 µ1 E(u, v) ∼ u DADu + uT D2 u − (d21 u41 + · · · + d2k u4k + · · · + d2N M u4N M ) = 2 2 4 ε2 T λ2 T 2 µ2 2 4 4 + v DADv + v D v − (d1 v1 + · · · + d2k vk4 + · · · + d2N M vN M ) (7) 2 2 4 β 2 2 2 2 − (d1 u1 v1 + · · · + d2k u2k vk2 + · · · + d2N M u2N M vN M ), 2 where the lower index of d2k , uk , and vk ( i.e k ) is a transform of N and M , defined as k = M (i − 1) + j for i = 1, 2, . . . , N and j = 1, 2, . . . , M , so that k = 1, 2, . . . , M N. Further details are also presented in Section 4. We minimize E(u, v) defined in (7) on a Nahari manifold by the Method of Lagrange Multipliers to get (8). See attached file for detail. We can transform the matrix £ into A by the origan Equations (4) the vector u and v are also present and transform synchronously.  2 −1 ° 3 ° 2   ε D ADu − λ1 u + µ1 u + βu ◦ v = 0 in Ω. 3 2 (8) ε2 D−1 ADv − λ2 v + µ2 v ° + βv ◦ u° = 0 in Ω.   u=v=0 on ∂Ω.. 4.

(10) Multiplying by D in both sides of the equations in (8), we have  2 ° 2 ° 2   ε ADu − λ1 Du + µ1 Du ◦ u + βDu ◦ v = 0 in Ω. 2 2 ε2 ADv − λ2 Dv + µ2 Dv ◦ v ° + βDv ◦ u° = 0 in Ω.   u=v=0 on ∂Ω. Moving the nonlinear term to the right hand side and rearranging we get:  2 ° 2 ° 2   (ε A − λ1 I + βv ) ◦ Du = −µ1 u ◦ Du in Ω. 2 2 (ε2 A − λ2 I + βu° ) ◦ Dv = −µ2 v ° ◦ Dv in Ω.   u=v=0 on ∂Ω.. (9). We thus obtain the above equations where u ◦ v = (u1 v1 , . . . , uN vN )> denotes the r Hadamard product of u and v, u° = u ◦ · · · ◦ u denotes the r-times Hadamard prodN ×M r uct of u. If u is a vector in R , then u° is also in RN ×M for all r. We will summarize the spirit of the Fixed Point Method in this paragraph. First, lets take a vector ue0 and a vector ve0 in Ω, and normalize them as u0 = keu10 k2 u e0 and 1 1 1 v0 = kev0 k2 ve0 . Let ζk = keuk k2 and ηk = kevk k2 for all k. The purpose of the normalization is in order to prevent the solutions converging to trivial solutions ( i.e. converging to zero ). Second, we can find a form of the equations (9) which looks similar to a linear system Ax = b. That is, rewrite the equations (9) into the type (10).  ° 2 2 2  uk+1 = −µ1 u° in Ω.  (ε A − λ1 I + βηk vk ) ◦ De k ◦ Duk ° 2 ° 2 2 (10) (ε A − λ2 I + βζk uk ) ◦ De vk+1 = −µ2 vk ◦ Dvk in Ω.   uk = vk = 0 on ∂Ω. We now put u0 and v0 into (10), and solve them iteratively to get ue1 and ve1 , etc. until k k2 | |kvk+1 k2 −kvk k2 | we reach convergence. Once the solutions satisfy max( |kuk+1kuk2k−ku , )< k2 kvk k2 tolerance or once we have performed 10000 iterations, we output the solutions. Suppose we get a convergent solution u∞ for k going to ∞. Then we will output the solutions 1 1 (11) u = (ζ∞ ) 2 u∞ and v = (η∞ ) 2 v∞ . The reason is as following. From the first item of (10), we get: ° 2 )◦D (ε2 A − λ1 I + βη∞ v∞. 1 2 u∞ = −µ1 u° ∞ ◦ Du∞ . ζ∞. (12). Now plugging (11) into Equation (12), Equation (12) will become −1 ° v2)◦D (ε2 A − λ1 I + βη∞ η∞. −1 −3 1 2 (ζ∞ ) 2 u = −µ1 (ζ∞ ) 2 u° ◦ Du. ζ∞. (13). After doing Similar things in v, we find that Equation (13) is equivalent to Equation (9). We therefore output the solutions giving by (11).. 5.

(11) We now turn to a series of numerical experiments that were designed specifically to investigate some interesting phenomena. A full presentation of the results obtained will be presented after their descriptions. Experiment 1. In this experiment, We will consider an ensemble of systems all of which share certain parameter values but differ by other parameter values. Take N = 60, M = 180, and a convergence tolerance= 10−10 . Fix u on the (N, M ) = (30, 1) and v on the (N, M ) = (30, 90).. ε 2 = 0.01, λ = 30, µ = 1, engergy = 2.3686. ε 2 = 0.01, λ = 30, µ = 10, engergy = 0.23686. 1. 11. 1. 0.8. 10. 0.8. 0.6. 9. 0.6. 0.4. 8. 0.4. 3.5. 3. 7. 0.2 2. 6. 0. y. y. 0.2. 2.5. 0. 5 −0.2. 1.5. −0.2 4. −0.4. −0.4. 1. 3 −0.6. 2. −0.8 −1 −1. 1 −0.8 −0.6 −0.4 −0.2. Figure 2.3.. 0 x. 0.2. 0.4. 0.6. 0.8. −0.6. −1 −1. 1. 0.5. −0.8 −0.8 −0.6 −0.4 −0.2. 0 x. 0.2. 0.4. 0.6. 0.8. 1. Cases of the differing µ values. The upper left picture is the lateral view of the solutions for µ = 1. The upper right picture is the lateral view of the solutions for µ = 10. The lower pictures are their respective vertical views.. We will first discuss cases with differing µ values. We immediately see a difference in the height of the solutions. In Figure 2.3, the µ value in the left side is equal to 1, and that of the right side is equal to 10. The change in the breadth of the curves can hardly be seen in the vertical views. Thus, we may conjecture that the variable µ acts heavily on the height of the solutions. We may also note that as µ decreases, the height increases.. 6.

(12) 2. 2. ε = 0.03, λ = 30, µ = 1, engergy = 9.6907. ε = 0.007, λ = 30, µ = 1, engergy = 1.5207. 1. 1. 0.8. 0.8. 9. 0.6. 8. 0.4. 7. 0.2. 6. 0. 5. −0.2. 4. −0.4. 3. −0.6. 2. −0.8. 1. 12 0.6 10. 0.4. 8. 0. y. y. 0.2. 6. −0.2 −0.4. 4. −0.6 2. −0.8 −1 −1. −0.8 −0.6 −0.4 −0.2. Figure 2.4.. 0 x. 0.2. 0.4. 0.6. 0.8. −1 −1. 1. 10. −0.8 −0.6 −0.4 −0.2. 0 x. 0.2. 0.4. 0.6. 0.8. 1. Cases of the differing ε2 values. The upper left picture is the lateral view of the solutions for ε2 = 0.03. The upper right picture is the lateral view of the solutions for ε2 = 0.007. The lower pictures are their vertical views, respectively.. Next, we discuss cases of differing ε2 value. As ε2 decreases to zero, the solutions approach spike-like shapes. If ε2 is equal to 0.03, then ε should be equal to 0.1732. Similarly, if the ε2 is equal to 0.007, then ε should be equal to 0.0837. In fact, the ε2 value is the key variable controlling the spike-like nature of the solutions. Figure 2.4 shows that change in ε2 value affects the solutions’ somatotype and height, but that the change in height is not as clear as that found by varying µ. Because µ is fixed at -1 fixedly in this case, we can take the left hand case of Figure 2.3 as a reference i.e. the case of ε2 equals 0.01 (ε equals 0.1). We see that if the ε2 value is larger, the solutions are wider and taller.. 7.

(13) 2. 2. ε = 0.01, λ = 10, µ = 1, engergy = 1.0767. ε = 0.01, λ = 50, µ = 1, engergy = 3.5169. 1. 1. 8. 0.8. 12. 0.8. 7. 0.6. 10. 0.6 6. 0.4. 0.4. 0. 4. −0.2. 3. −0.4. 0.2 y. y. 8. 5. 0.2. 0. 6. −0.2 4. −0.4 2. −0.6. −0.6 2 1. −0.8 −1 −1. −0.8 −0.6 −0.4 −0.2. Figure 2.5.. 0 x. 0.2. 0.4. 0.6. 0.8. −0.8 −1 −1. 1. −0.8 −0.6 −0.4 −0.2. 0 x. 0.2. 0.4. 0.6. 0.8. 1. Cases of the diffring λ values. The upper left picture is the lateral view of the solutions for λ = 10. The upper right picture is the lateral view of the solutions for λ = 50. The lower pictures are their vertical views, respectively.. The last variable to consider is λ. The effect of changing the λ value is not similar to the effect of changing ε2 . We can see from Figure 2.5 and from the left hand part of Figure 2.3 that the solutions will be fatter and shorter when the λ value is smaller. We will further compare the effect of varying λ values on the vectors u and v a little later on this work. In the above case, comparisons are made only on spike shaped solutions i.e. situations where the variable ε2 is sufficiently small. The variable β turns out to affect the attraction of u and v, so we will not compare cases of various β values in this experiment.. 8.

(14) Experiment 2. After observing the solution shape here, we will be interested in the energy variation. Thus we will present calculated values of the energy for solutions over a disc domain. Take N = 60, M = 180, ε2 = 0.005, λ1 = λ2 = 30, µ1 = µ2 = 1, β = −108 , and the convergent tolerance= 10−10 . We set the spaces of the solution of u on N = 10, 20, 30, 40, and 50, and M = 1. Consider a solution of v which coils a circuit from M = 3 to M = 177. We then observe the energy variation. N = 60, M = 180, ε2 = 0.005, λ = 30, µ = 1, β = −100000000. N = 60, M = 180, ε2 = 0.005, λ = 30, µ = 1, β = −100000000. 1.4. 0. 1.35. r10. r10. r20. r20 Energy relative differences (log10). r30 r40. 1.3. r50 Energy. 1.25 1.2 1.15 1.1. r30 r40 −5. r50. −10. 1.05 1. 0. 0.5. Figure 2.6.. 1. 1.5. θj. 2. 2.5. 3. 3.5. −15. 0. 1. 2. 3. 4. θj. 5. 6. The energy variation according to the space of u and v. Here we present data for N = 10, 20, 30, 40, 50. The right figure presents of the log10 of the relative error along the y axis.. N = 60, M = 180, ε = 0.005, λ = 10, µ = 1, β = −100000000. N = 60, M = 180, ε = 0.005, λ = 10, µ = 1, β = −100000000. 1. 1.346. 1.346. 0.8. 1.310. 1.310. 0.6 0.4. 1.275. 1.4. 1.275. 1.239. 1.3. 1.239. 0.2 1.204 0. 1.204. 1.2. 1.169 −0.2. 1.169 1.1. 1.133. −0.4. 1.098. −0.6. 1.133 1 −1. 1. −0.5 1.063. −0.8 −1 −1 −0.5 0 0.5 1 The five radius are 0.15702, 0.32231, 0.4876, 0.65289, 0.81818.. Figure 2.7.. 1.027. 0.5 0. 0 0.5. −0.5. 1 −1 The five radius are 0.15702, 0.32231, 0.4876, 0.65289, 0.81818.. Here the color bar is used to show the energy variation from three different aspects. The left picture is a vertical view.. In this experiment, we find that the lower energy solutions u and v are close to the center of the disc, and higher energy costs occur close to the center of the disc. We find that irregular energy distributions occur in the neighborhood of M = 1. In these cases, it may be that the solutions u and v are close to each other in order to exclude the opposite side, and in doing so produce a larger energy configuration. We can also see that the energy variation decreases as the solution of v toward an angle equal to π in the counterclockwise direction. The energy will of course increase. 9. 1.098 1.063 1.027.

(15) again after going through angle equal to π. We may, therefore, set the solution for u on the angular origin and the solution for v on the angle equal to π to see what radius space solutions will produce a minimum energy. The last result of note is that the minimum energy will occur on the N = 30 circle. Experiment 3. We will address the energy variation according to different λ values of a system, and find the value that produces a minimum of energy. We take N = 60, M = 180, ε2 = 0.005, µ1 = µ2 = 1, β = −108 , and the convergent tolerance= 10−10 . We then observe the cases of (λ1 , λ2 ) = (30, 30), (30, 40), and (30, 50). We will not present this in point-by-point detail because doing so would take too much time. We first set an interval in which to observe the energy change. We then subdivide the interval and then determine which space has minimum energy. We further subdivide this space and continue the computation recursively. Figures 2.8, 2.9, and 2.10 are set in an interval by interval in order in order to see the energy variation by area. N = 60, M = 180, ε2 = 0.005, λ1 = 30, λ2 = 30, µ = 1, β = −100000000, the min. energy location for peak u is on (35,0), and peak v is on (35,90).. 1.4. 1.4. 1.35. 1.35. 1.3. 1.3. 1.25. 1.25 Energy. Energy. N = 60, M = 180, ε2 = 0.005, λ1 = 30, λ2 = 30, µ = 1, β = −100000000, the min. energy location for peak u is on (35,0), and peak v is on (35,90).. 1.2. 1.2. 1.15. 1.15. 1.1. 1.1. 1.05. 1.05. 1 −1. −0.9. −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 The location of peak v on the left radius.. Figure 2.8.. −0.2. 1 0.1. −0.1. 1.6 1.55. 1.5. 1.5. 1.45. 1.45. 1.4. 1.4. Energy. Energy. 1.6. 1.35. 1. 1.35. 1.3. 1.3. 1.25. 1.25. 1.2. 1.2. −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 The location of peak v on the left radius.. Figure 2.9.. 0.9. N = 60, M = 180, ε2 = 0.005, λ1 = 30, λ2 = 40, µ = 1, β = −100000000, the min. energy location for peak u is on (35,0), and peak v is on (30,90).. 1.55. −0.9. 0.3 0.4 0.5 0.6 0.7 0.8 The location of peak u on the left radius.. The energy variation of λ1 = 30 and λ2 = 30 on N = 10, 15, 20, 25, 30, 35, 40, 45, 50, and 55.. N = 60, M = 180, ε2 = 0.005, λ1 = 30, λ2 = 40, µ = 1, β = −100000000, the min. energy location for peak u is on (35,0), and peak v is on (30,90).. 1.15 −1. 0.2. −0.2. −0.1. 1.15 0.1. 0.2. 0.3 0.4 0.5 0.6 0.7 0.8 The location of peak u on the right radius.. 0.9. 1. The energy variation of λ1 = 30 and λ2 = 40 on N = 10, 15, 20, 25, 30, 35, 40, 45, 50, and 55. The left picture is a basic on the solution v, and arrange energy for each solution u. The right is correspondingly. 10.

(16) 2. 2. N = 60, M = 180, ε = 0.005, λ = 30, λ = 50, µ = 1, β = −100000000, 1 2 the min. energy location for peak u is on (35,0), and peak v is on (30,90).. 1.75. 1.75. 1.7. 1.7. 1.65. 1.65. 1.6. 1.6. 1.55. 1.55. Energy. Energy. N = 60, M = 180, ε = 0.005, λ = 30, λ = 50, µ = 1, β = −100000000, 1 2 the min. energy location for peak u is on (35,0), and peak v is on (30,90).. 1.5. 1.5. 1.45. 1.45. 1.4. 1.4. 1.35. 1.35. 1.3 −1. −0.9. −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 The location of peak v on the left radius.. Figure 2.10.. −0.2. 1.3 0.1. −0.1. 0.2. 0.3 0.4 0.5 0.6 0.7 0.8 The location of peak u on the left radius.. 0.9. 1. The energy variation of λ1 = 30 and λ2 = 50 on N = 10, 15, 20, 25, 30, 35, 40, 45, 50, and 55. The left picture is a basic on the solution v, and arrange energy for each solution u. The right is corresponds to the left. It views from the solution of u.. After such recursive calculation, We show the minimum energy space in Figure 2.11. In the cases of (λ1 , λ2 ) = (30, 30), (30, 40), and (30, 50), the minimum energy spaces are located at (34, 34), (34, 31), and (34, 28), respectively. The energy will be larger values of λ in the system.. N = 60, M = 180, ε2 = 0.005, µ = 1, β = −100000000. N = 60, M = 180, ε2 = 0.005, µ = 1, β = −100000000. 1.4. 1.4. 1.35. 1.35 1.3. 1.25. 1.25. 1.2. 1.2. 1.15. Energy. Energy. 28 1.3. 31. 1.1. 1.15. 34. 1.1. 1.05 1 0.95. 34. 30−30 30−40 30−50 min sp. 34 −0.65. −0.6 −0.55 −0.5 −0.45 The location of peak v on the left radius.. Figure 2.11.. 30−30 30−40 30−50 min sp −0.4. 1.05 1. −0.35. 0.95 0.4. 34 0.45. 0.5 0.55 0.6 0.65 The location of peak u on the left radius.. 0.7. 0.75. The minimum energy spaces about (λ1 , λ2 ) = (30, 30), (30, 40) and (30, 50), after subdividing the interval. We fix the peak of u on M = 1 and the peak of v on the M = 90. The minimum energy spaces are located on the radio space of the peak of are located on the (the radio space of peak of u, the radio space of the peak v)= (34, 34), (34, 31), and (34, 28) for (λ1 , λ2 ) = (30, 30), (30, 40), and (30, 50).. 11.

(17) Experiment 4. We take N = 60, M = 180 and observe the solution for u and for v variously locate on the M = 1 and M = 90. We list the different variables, corresponding minimum energy, and the location of the minimum energy configuration as follows. Case (λ1 , λ2 ) ε2 µ β Min. Eng. (ur , vr ). 01 (30, 50) 0.005 1 −1014 0.134566081689497D+01 (34, 28). 02 (30, 50) 0.005 1 −1000 0.134566081689497D+01 (34, 28). 03 (30, 50) 0.005 1 −1 0.134566081689497D+01 (34, 28). Case (λ1 , λ2 ) ε2 µ β Min. Eng. (ur , vr ). 04 (30, 50) 0.005 1 −0.001 0.134566081689497D+01 (34, 28). 05 (30, 50) 0.005 10 −1 0.134566081689497D+00 (34, 28). 06 (30, 50) 0.005 0.001 −1 0.134566081689497D+04 (34, 28). Case (λ1 , λ2 ) ε2 µ β Min. Eng. (ur , vr ). 07 (30, 50) 1 1 −1000 0.467974836047890D+03. 08 (50, 30) 0.005 1 −1000 0.134566081689497D+01 (28, 34). 09 (30, 40) 0.005 1 −1000 0.117934943480590D+01 (31, 34). Case (λ1 , λ2 ) ε2 µ β Min. Eng. (ur , vr ). 10 (30, 30) 0.01 1 −1000 0.223767278528725D+01 (45, 45). 11 (30, 40) 0.01 1 −1000 0.254427449864588D+01 (45, 40). 12 (30, 50) 0.01 1 −1000 0.284999932025183D+01 (45, 36). 12.

(18) Case (λ1 , λ2 ) ε2 µ β Min. Eng. (ur , vr ). 13 (30, 30) 0.003 1 −1000 0.601612204905536D+00 (28, 28). 14 (30, 40) 0.003 1 −1000 0.707591934264562D+00 (28, 25). 15 (30, 50) 0.003 1 −1000 0.821203750097474D+00 (28, 23). Case (λ1 , λ2 ) ε2 µ β Min. Eng. (ur , vr ). 16 (30, 30) 0.0005 1 −1000 0.143088685313069D+00 (14, 14). 17 (30, 40) 0.0005 1 −1000 0.177526999280233D+00 (14, 12). 18 (30, 50) 0.0005 1 −1000 0.216148141406401D+00 (14, 11). From the several cases above, We can generalize some results as follows. First, the value of β is expected to have a smaller effect on the energy. We can see it from the cases of 01, 02, 03, and 04. The minimum energies are the same as that for β = −1014 and β = −10−3 . If the value of µ is smaller (i.e. µ approaches to 0 ), the energy will become larger. In contrast, we can find it from cases 03, 05, and 06. Additionally, there is an episode in these cases. The digits of energy are the same in the three cases, and the difference among them resides only in the decimal point location. As the values of λ1 and λ2 change together, the minimum energy space of u and v is also change. We can see this clearly in cases 02 and 08. If the value of ε2 is grows away from zero, then the solution’s shape would no longer be a spike. Take ε2 = 1 as an example ( i.e. case 07 ). The solution approaches a bell shape as seen in Figure 2.12. There is also a vexing situation in this case: no matter where we start u and v initially, they always converge to the same space, that seen in Figure 2.12. The computing process has also proven so elusive, that we could not observe the energy variation in this case. Each energy is also closed to 0.467974836048D+03. Since we know that the variables β and µ are not very influential on the energy variation, we have designed a group of case experiments which will help to resolve the minimum energy space. We start by choosing ε2 values such that the solution shape will be a spike. Cases 10, 11, and 12 are generalized as a group, while cases 13, 14,. 13.

(19) 15, and the cases 16, 17, 18 represent two more groups. We can learn something by comparing the group. The first thing is that if we find solutions for different systems that have the same λ value ( i.e. the solutions u in the cases 10, 11, and 12 all have λ1 = 30 ), then the various solutions will share the same space of minimal energy. My advisor and Professor Wu consider that this may be coincidental and thus, the big question remains why the minimal spaces do not change. The second observation is that when the value of ε2 is smaller, the space of the solutions of u will be supported near the center of the disc. We can see this in cases 02, 12, 15, and 18. The third thing to observe is about the space of the solutions of v. If the λ2 value for the solution is large, then the solution for v will also move toward the center of the disc. We can see this in the three groups mentioned above. ε 2 = 1, λ = 30, µ = 1, engergy = 467.9748 1 14. 0.8 0.6. 12. 0.4 10 0.2 8. −0.2. 6. y. 0. −0.4 4 −0.6 2. −0.8 −1 −1. Figure 2.12.. −0.8 −0.6 −0.4 −0.2. 0 x. 0.2. 0.4. 0.6. 0.8. 1. While the variable ε2 is equal to 1, the solutions shape as above. This data from case 07.. Experiment 5. In this paragraph, we present the results of experiments completed during my master’s work. Here, I show only pictures of these results and give simple descriptions of the experiments in the captions. Some of these results will be used to check or test the above four experiments. N = 30, M = 90, ε = 0.005, λ = 10, µ = 1, u on the (15, 0), v on the (15, 22). N = 30, M = 90, ε = 0.005, λ = 10, µ = 1, u on the (15, 0), v on the (15, 11) 0.3351 0.3351. 0.3351. 0.3351. 0.3351. 0.3351 Energy. Energy. 0.3351. 0.3351. 0.3351. 0.3351 0.3351 0.3351 0.335. 0.3351 0.3351. 0. 2. 4. 6 log(−β). 8. 10. 12. 14. 0.335. 0. 2. 4. 6 log(−β). 8. 10. 12.

(20) N = 30, M = 90, ε = 0.005, λ = 10, µ = 1, u on the (15, 0), v on the (15, 45). 0.3351. 0.3351. 0.3351. 0.3351. 0.3351. 0.3351 Energy. Energy. N = 30, M = 90, ε = 0.005, λ = 10, µ = 1, u on the (15, 0), v on the (15, 33). 0.3351. 0.3351. 0.3351. 0.3351. 0.335. 0.335. 0.335. 0. 2. 4. Figure 2.13.. 6 log(−β). 8. 10. 12. 0.335. 0. 2. 4. 6 log(−β). 8. 10. 12. A plot of the energy variation associated with variable β values.. We can see that the β is less relation about the energy change from Figure 2.13. We take N = 30, M = 90, ε2 = 0.005, λ = 10, µ = 1, and tolerance= 10−10 . We fix the solution for u at (N, M ) = (15, 0), and set the solution for v angles of 11, 22, 33, and 45 degrees. We then observe that the energy variation is the same as that found in Experiment 4. N = 30, M = 90, ε2 = 0.005, µ = 1, β = −100000000, u on the (15, 0). N = 30, M = 90, ε2 = 0.005, µ = 1, β = −100000000, u on the (15, 0). 1. 1.5 (λ1,λ2)=(10,10). (λ1,λ2)=(30,10). 1.4. (λ1,λ2)=(10,20). 0.9. (λ1,λ2)=(10,30). (λ1,λ2)=(30,20) (λ1,λ2)=(30,30). 1.3. 0.8. Energy. Energy. 1.2 0.7. 0.6. 1.1 1. 0.5. 0.9. 0.4. −0.9. 0.8. −0.8. −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 The location of peak v on the left radius.. Figure 2.14.. −0.1. 0. 0.7 −0.9. −0.8. −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 The location of peak v on the left radius.. −0.1. 0. With a fixed solution for u at (N, M ) = (15, 0), and set the solution for v at M = 45 and various of N , We present the corresponding variation of energy.. Figure 2.14 tells us that the energy would increase as the λ2 value increases. Similarly as the value of λ1 increases, the energy also increases. Take N = 30, M = 90, ε2 = 0.005, µ = 1, and tolerance= 10−10 . Fix the solution u at (N, M ) = (15, 0), and set the solution for v at M = 45 with various N values in order to observe the energy variation. Conjugate different λ values and choose six cases to compare the energy variation for (λ1 , λ2 ) = (10, 10), (10, 20), (10, 30), (30, 10), (30, 20), and (30, 30). What we see is in accord with the statistics in Experiment 4 and Figure 2.4. Moreover, the energy variation seen in the left picture is like that seen in the right picture. 15.

(21) 2. 2. N = 30, M = 90, ε = 0.005, λ = 10, λ = 10, µ = 1, β = −100000000, 1 2 the min. energy location for peak u is on (16,0), and peak v is on (16,45).. 0.46. 0.46. 0.44. 0.44. 0.42. 0.42. 0.4. 0.4. Energy. Energy. N = 30, M = 90, ε = 0.005, λ = 10, λ = 10, µ = 1, β = −100000000, 1 2 the min. energy location for peak u is on (16,0), and peak v is on (16,45).. 0.38. 0.38. 0.36. 0.36. 0.34. 0.34. 0.32 −1. −0.8. −0.6 −0.4 −0.2 The location of peak v on the left radius.. Figure 2.15.. 0.32. 0. 0.2. 0.4 0.6 0.8 The location of peak u on the left radius.. 1. The energy variation of λ1 = 10 and λ2 = 10 for u at M = 0, and v at M = 45. Their routes are from N = 2 to N = 29 dividedly. The left picture is a basic on the peak v and arrange energy for each peak u. The right picture is corresponds with the left. It is viewed from the peak u. Thus, there are 28 × 28 = 784 points in the figure.. N = 30, M = 90, ε2 = 0.005, λ1 = 10, λ2 = 20, µ = 1, β = −100000000, the min. energy location for peak u is on (16,0), and peak v is on (12,45).. N = 30, M = 90, ε2 = 0.005, λ1 = 10, λ2 = 20, µ = 1, β = −100000000, the min. energy location for peak u is on (16,0), and peak v is on (12,45).. 0.7. 0.7. 0.65. 0.65 Energy. 0.75. Energy. 0.75. 0.6. 0.6. 0.55. 0.55. 0.5 −1. −0.8. −0.6 −0.4 −0.2 The location of peak v on the left radius.. Figure 2.16.. 0.5. 0. 1.05. 1.05. 1. 1. 0.95. 0.95 Energy. 1.1. 0.9. 0.85. 0.8. 0.8. 0.75. 0.75. −0.6 −0.4 −0.2 The location of peak v on the left radius.. Figure 2.17.. 0.4 0.6 0.8 The location of peak u on the left radius.. 1. 0.9. 0.85. −0.8. 0.2. N = 30, M = 90, ε2 = 0.005, λ1 = 10, λ2 = 30, µ = 1, β = −100000000, the min. energy location for peak u is on (16,0), and peak v is on (10,45).. 1.1. 0.7 −1. 0. The energy variation of λ1 = 10 and λ2 = 20. Other ways of doing things are the same as Figure 2.15.. N = 30, M = 90, ε2 = 0.005, λ1 = 10, λ2 = 30, µ = 1, β = −100000000, the min. energy location for peak u is on (16,0), and peak v is on (10,45).. Energy. 0. 0. 0.7. 0. 0.2. 0.4 0.6 0.8 The location of peak u on the left radius.. 1. The energy variation of λ1 = 10 and λ2 = 30. Other ways of doing things are the same as Figure 2.15. 16.

(22) 2. 2. N = 30, M = 90, ε = 0.005, λ = 20, λ = 10, µ = 1, β = −100000000, 1 2 the min. energy location for peak u is on (12,0), and peak v is on (16,45).. N = 30, M = 90, ε = 0.005, λ = 20, λ = 10, µ = 1, β = −100000000, 1 2 the min. energy location for peak u is on (12,0), and peak v is on (16,45).. 0.7. 0.7. 0.65. 0.65 Energy. 0.75. Energy. 0.75. 0.6. 0.6. 0.55. 0.55. 0.5 −1. −0.8. −0.6 −0.4 −0.2 The location of peak v on the left radius.. Figure 2.18.. 0.5. 0. 0.95. 0.95. 0.9. 0.9. 0.85. 0.85. Energy. Energy. 1. 0.8. 0.75. 0.7. 0.7. −0.6 −0.4 −0.2 The location of peak on the left radius.. 0. 0.65. 1.4. 1.4. 1.3. 1.3. 1.2. 1.2. Energy. Energy. 1.5. 1.1. 0.4 0.6 0.8 The location of peak u on the left radius.. 1. 1.1. 1. 1. 0.9. 0.9. −0.6 −0.4 −0.2 The location of peak v on the left radius.. Figure 2.20.. 0.2. N = 30, M = 90, ε2 = 0.005, λ1 = 20, λ2 = 30, µ = 1, β = −100000000, the min. energy location for peak u is on (12,0), and peak v is on (10,45).. 1.5. −0.8. 0. The energy variation of λ1 = 20 and λ2 = 20. Other ways of doing things are the same as Figure 2.18.. N = 30, M = 90, ε2 = 0.005, λ1 = 20, λ2 = 30, µ = 1, β = −100000000, the min. energy location for peak u is on (12,0), and peak v is on (10,45).. 0.8 −1. 1. 0.8. 0.75. Figure 2.19.. 0.4 0.6 0.8 The location of peak u on the left radius.. N = 30, M = 90, ε2 = 0.005, λ1 = 20, λ2 = 20, µ = 1, β = −100000000, the min. energy location for peak u is on (12,0), and peak v is on (12,45).. 1. −0.8. 0.2. The energy variation of λ1 = 20 and λ2 = 10 for u at M = 0, and v at M = 45. Their routes are from N = 2 to N = 29 dividedly. The left picture is a basic on the peak v and arrange energy for each peak u. The right picture is corresponds with the left. It is viewed from the peak u. Thus, there are 28 × 28 = 784 points in the figure.. N = 30, M = 90, ε2 = 0.005, λ1 = 20, λ2 = 20, µ = 1, β = −100000000, the min. energy location for peak u is on (12,0), and peak v is on (12,45).. 0.65 −1. 0. 0. 0.8. 0. 0.2. 0.4 0.6 0.8 The location of peak u on the left radius.. 1. The energy variation of λ1 = 20 and λ2 = 30. Other ways of doing things are the same as Figure 2.18.. 17.

(23) We have also performed tests with other values of λ1 and λ2 and arrange the numerical results as follows. (λ1 , λ2 ) (10, 06) (10, 08) (10, 10) (10, 12) (10, 14) (10, 20) (10, 30) (10, 40) (10, 50) (10, 60). Real location of the Pu 5.081967213114754e-001 5.081967213114754e-001 5.081967213114754e-001 5.081967213114754e-001 5.081967213114754e-001 5.081967213114754e-001 5.081967213114754e-001 5.081967213114754e-001 5.081967213114754e-001 5.081967213114754e-001. Real location of the Pv -6.065573770491803e-001 -5.409836065573771e-001 -5.081967213114754e-001 -4.754098360655738e-001 -4.426229508196722e-001 -3.770491803278689e-001 -3.114754098360656e-001 -2.786885245901640e-001 -2.786885245901640e-001 -2.459016393442623e-001. (N, M ) for Pu , Pv (16, 0), (19, 45) (16, 0), (17, 45) (16, 0), (16, 45) (16, 0), (15, 45) (16, 0), (14, 45) (16, 0), (12, 45) (16, 0), (10, 45) (16, 0), (09, 45) (16, 0), (09, 45) (16, 0), (08, 45). (20, 10) (20, 20) (20, 30). 3.770491803278689e-001 3.770491803278689e-001 3.770491803278689e-001. -5.081967213114754e-001 -3.770491803278689e-001 -3.114754098360656e-001. (12, 0), (16, 45) (12, 0), (12, 45) (12, 0), (10, 45). (30, 10) (30, 20) (30, 30). 3.114754098360656e-001 3.114754098360656e-001 3.114754098360656e-001. -5.081967213114754e-001 -3.770491803278689e-001 -3.114754098360656e-001. (10, 0), (16, 45) (10, 0), (12, 45) (10, 0), (10, 45). Table 2.1.. For N = 30 and M = 90, we use various values of λ and calculate the location of their minimum energy support. The values in second column are computed by ∆r × (NPu − 0.5) = (2/(1 + 2M )) × (NPu − 0.5) where NPu is the location of the peak of u. We can find the location of the peak values of u and v in the forth column. The third column values are obtained just as in the second column, but they need a minus sign and use NPv as a substitute for NPu in the computation.. We can obtain some important facts from Table 2.1. For a fixed N and M , the location of the minimum energy configuration is certainly related to λ. If the λ2 value is larger than λ1 , then the minimum energy location moves toward the center of the disc. If the values of λ1 and λ2 change, then the location of minimum energy would change as well. The results are the same as those in Experiment 4, but the division of the domain is not the same.. 18.

(24) 3. The Hyperplane-Constrained Continuation Method In the introduction, we noted that the ε2 term will be absorbed by a change of variables. Despite this, the two equations are essentially the same. In this section, we will offer a summarized description of the Hyperplane-Constrained Continuation Method, and we will show some figures exhibiting specially shaped solutions in some bifurcations. Suppose that we solve for the function g(x) in a difficult way, but that there exists a function f (x) which is similar to g(x), but that can be obtained more easily. Now, we define a new function H(f, g, τ ) and add a variable τ which connects f (x) and g(x) as: H(f, g, τ ) = (1 − τ )f (x) + τ g(x). When τ = 0, H(f, g, τ ) equal to f (x) i.e. H(f, g, τ ) = f (x). If we then increase τ from 0 to 1 somehow, then H(f, g, τ ) will become g(x) at τ = 1 i.e. H(f, g, τ ) = g(x) at τ = 1. In the end, we will obtain the solution of g(x). Now we may apply the above ideas to solve our equation. We rename the variables u and v to u1 and u2 . The difficult part of this particular solution is the β term, i.e. the coupling term, because this term contains both u1 and u2 . In the presence of this term, the solution of u1 is affected by the value of u2 . Similarly, u2 is controlled by u1 , too. We therefore assign a new variable τ in the β term as follows.  ° 3 ° 2   £u1 − λ1 u1 + µ1 u1 + β21 τ u1 ◦ u2 = 0 3 ° 2 H(u(s), τ (s)) = £u2 − λ2 u2 + µ2 u° 2 + β12 τ u1 ◦ u2 = 0   u1 = u2 = 0. in Ω. in Ω. on ∂Ω.. We assume the solutions move along a curve C which is our supposition. Let a virtual variable s be an arc-length along the curve, so u = (u1 , u2 ) and τ are related to s. Since ∆ will be transfer a matrix £, the variable β can be expressed in matrix form. u1 ◦ u2 denotes the Hadamard product of u1 and u2 as described in Section 2, and r u° 1 = u1 ◦ · · · ◦ u1 denotes the r-time Hadamard product of u1 . First, we compute the initial solution, called u0 (which holds at τ = 0), as → − (u0 , τ0 )T = (u10 , u20 , τ0 )T . Second, we introduce a target vector t0 which comes from the u and τ respective different by s and then suitably normalized. Third, we − → consider a hyperplane P1 whose normal vector is t0 and contains a point (¯ u1 , τ¯1 )T = − → (u0 , τ0 )T + δ t 0 where δ is a small step length. Thus, the hyperplane P1 will be intersected at a point by the solution curve C. This point is called (u1 , τ1 )T . Forth, we find (u1 , τ1 )T using Newton’s Method, iteratively converging to (u1 , τ1 )T from (¯ u1 , τ¯1 )T . Once we find (u1 , τ1 )T , we continue using the same steps to find (u2 , τ2 )T and so on. The main idea of the Continuation Method is contained in the first through forth steps, and, in practice, there are some problems and computing difficulties that arise 19.

(25) in going further. We will save an explanation of such details for later on in this paper. − → We may ask how we get t0 in the first place. In fact, the variable s does not exist, but we can find it in another way. We note that H(u(s), τ (s)) is different than s to zero at every point, since H(u(s), τ (s)) = 0 for all (u(s), τ (s)) on the curve C. According to the Chain Rule, we can obtain " # du du dτ ds Hu + Hτ = 0 ⇐⇒ [Hu , Hτ ] dτ = 0. ds ds ds Let x0 = (u0 , τ0 ) and the vector ( du (u0 , τ0 ), dτ (u0 , τ0 ))T = (us (x0 ), τs (x0 ))T be called ds ds t¯0 . Since we can compute the value of the matrix [Hu (x0 ), Hτ (x0 )], the vector t¯0 can be gotten by solving a linear system like Ax = 0. This presents a new problem that must be solved. The problem is that, since the vector t¯0 maybe solved a zero solution using a computer, we must consider solving the case of a nonzero component in the t¯0 . Without a loss of generality, let the last element be 1 in the t¯0 . After we solve for − → a nonzero t¯0 , the vector t0 may quickly be obtained. " # − → us (x0 ) t¯0 − → = − . t0 = → kt¯ k τ (x ) 0. s. 0. Eventually, we will give explicit forms for Hu and Hτ is, in fact, derived from two equations, so we suppose that H1 = £u1 − λ1 u1 + µ1 u31 + β21 τ u1 u22 = 0 and H2 = £u2 − λ2 u2 + µ2 u32 + β12 τ u21 u2 = 0. Rewrite these in matrix form, we have Hu and Hτ as follows. " Hu =. ∂H1 ∂u1 ∂H2 ∂u1. ∂H1 ∂u2 ∂H2 ∂u2. #. # 2 ° 2 2β τ u ◦ u £ − λ1 I + 3µ1 u° + β τ u 21 21 2 1 2 1 , = 2 ° 2 2β12 τ u1 ◦ u2 £ − λ2 I + 3µ2 u° 2 + β12 τ u1 " # 2 β21 u° ◦ u 1 2 Hτ = . 2 β12 u° 1 ◦ u2 − → − → Once we have t0 theb, given a δ, we can easily calculate (¯ u1 , τ¯1 )T = (u0 , τ0 )T +δ t0 . − → We will find the hyperplane P1 via (¯ u1 , τ¯1 )T and t0 . # " # " # " # " # " # " − → − → us (x0 ) u us (x0 ) u¯1 − → u − → u¯1 i i ⇐⇒ h − , i=h − , h t0 , i = h t0 , → → τ¯1 τs (x0 ) τ τs (x0 ) τ¯1 τ ". where u and τ are the coordinates of the hyperplane P1 . Given that we want to get (u1 , τ1 ), we can put u1 and τ1 into P1 and H(u(s), τ (s)) to obtain the following:    H(u1 (s), τ1 (s)) = 0 # " # " £− ¤ £ ¤ u (s) u ¯ 1 1 → → → → T −  =0 − − us (x0 )T , − τs (x0 )  P1 (u1 (s), τ1 (s)) = us (x0 ) , τs (x0 ) τ¯1 τ1 (s) 20.

(26) We now solve for (u1 , τ1 ) using Newton’s Method. In R1 , there is a differentiable function f : R → R for which there exists an x such that f (x) = 0. We want to find x. First, we choose an x0 and compute both f (x0 ) and f 0 (x0 ). We use these values to find x1 as: x1 = x0 −. f (x0 ) f 0 (x0 ). i.e.. f (xk ) for k = 0, 1, 2, · · · . f 0 (xk ). xk+1 = xk −. A similar mechanics holds for performing Newton’s Method in Rn for n ≥ 2. Newton’s Method for x to x 0. 1. 12 f(x) 10. y value. 8. f(x0). 6. 4. 2. 0. −2 −1. x1 0. 1. x0 2. 3. 4. x value. Figure 3.1.. The image of Newton’s Method in R1 for x0 to x1 .. In Rn for n ≥ 2, we suppose that there is a differentiable function f : Rn → Rn , and an x such that f (x) = 0. We will find x ∈ Rn . Then: xk+1 = xk − J −1 (xk ) · f (xk ) for k = 0, 1, 2, · · · , where J(xk ) is the Jacobian for f (xk ). In fact, we do not solve J −1 (xk ), because it is computationally expensive. Let y(xk ) = J −1 (xk ) · f (xk ). Then J(xk ) · y(xk ) = f (xk ). Therefore, we solve J(xk )·y(xk ) = f (xk ), then iteratively compute xk+1 = xk −y(xk ). Once f (xk+1 ) < tolerance which is given by user, we stop the method. In our case, we already have a variable x, so we change this to a variable a instead of using x. Our equation will thus become ( H(u1 (s), τ1 (s)) = 0 f (a1 ) = f (u1 (s), τ1 (s)) = . P1 (u1 (s), τ1 (s)) = 0 Thus, the next step is finding J(x0 ) = J(u0 (s), τ0 (s)). " # Hu (a0 ) Hτ (a0 ) J(a0 ) = J(u0 (s), τ0 (s)) = . (us (a0 ))T τs (a0 ) 21.

(27) We can find the lower part of the Jacobian J(a0 ) to be equal to t¯0 T . We match up the previous states of Hu and Hτ , inserting x0 to obtain the Jacobian J(a0 ). The process proceeds iteratively to finding x1 etc., but there are inevitable problems associated with this process.. Figure 3.2.. The image of Hyperplane-Constrained Continuation Method on the curve C. The pink curve is represented as the process of Newton’s Method.. From the arguments of Section 2, we should be able to detect that the ground state solutions almost everywhere in this system. If the computing accuracy has a little error, then the solutions will occur at a small distance along the original curve C. However, these solutions may not be directed along the curve C. Therefore, we will add another two hyperplane to force the direction of the solutions to lie along the original curve C. Let the two hyperplanes be P2 and P3 and their normal vectors be n2 and n3 , respectively. Considering iterations that take x0 to x1 , we can also rewrite the normal vectors as n20 and n30 . We have that n20 is perpendicular to n30 , n20 is perpendicular − → − → to t0 , and n30 is perpendicular to t0 . What are n20 and n30 ? Indeed, n20 and n30 are also the vectors of H which differ by an amount s from zero. In an earlier paragraph we solved a homogeneous linear system as Ax = 0. In that instance however, we → − supposed the last component of t0 be 1. We now loosen this restriction, and suppose that the above states, n20 and n30 will satisfy with as follows:  # " 3 events  − →T n20   # ( " =0 t0   0 Hu (x0 ) · n20 = 0 n 30 # " =0 (1) (3) [nT20 , 0] (2)  0 Hu (x0 ) · n30 = 0 − →T n30   =0   t0 0 We can make a table to show that the normal vectors and the points pass through. 22.

(28) each of the three hyperplanes when x0 goes to x1 . P1 P2 P3. normal direction pass through − → t0 (u¯1 , τ¯1 ) and (u1 , τ1 ) [nT20 , 0] (u¯1 , τ¯1 ) and (u1 , τ1 ) T [n30 , 0] (u¯1 , τ¯1 ) and (u1 , τ1 ). The equations that are used to develop this solution are listed as follows: H : H(u1"(s), τ1 (s)) #=0 " # h− i h− u1 (s) → →i u¯1 =0 P1 : t0 − t0 τ¯1 τ1 (s) " # " # £ T ¤ u1 (s) £ T ¤ u¯1 P2 : n20 , 0 − n20 , 0 =0 τ1 (s) τ¯1 " # " # £ T ¤ u1 (s) £ T ¤ u¯1 =0 P3 : n30 , 0 − n30 , 0 τ1 (s) τ¯1 Then we find (u1 , τ1 )T by the process of the Newton’s Method. That is solving ak+1 = ak − J −1 (ak ) · f (ak ) and J(ak ) · y(ak ) = f (ak ), where k is the number of steps in Newton’s Method. Since we add two hyperplanes P2 and P3 the Jacobian matrix will be changed.     Hu (ak ) Hτ (ak )) H(ak )  u (a )T  P (a )  τs (ak )   s k   1 k    · y(ak ) =  , T  n20 (ak )   P2 (ak )  0 n30 (ak )T. 0. P3 (ak ). where k = 0, 1, 2, · · · , the initial a0 is (u¯1 , τ¯1 ), and iterations stop when the right hand side less than tolerance. If we assume above system has the form Ay = b, then the dimension of A is (2n + 3) × (2n + 1), the dimension of y is (2n + 1) × 1, and the dimension of b is (2n + 3) × 1. That is, an over determined system. We can therefore add extra elements to the equations so that they are not over determined. We expand the matrix A and the vector y as follows:      Hu (ak ) Hτ (ak ) n20 (ak ) n30 (ak ) y(ak )   u (a )T τ (a )  0 0    s k   s k = ·    bk     n20 (ak )T  0 0 0 dk T n30 (ak ) 0 0 0. H(ak ) P1 (ak ) P2 (ak ) P3 (ak ).    . . Then the dimension of A is (2n + 3) × (2n + 3), and the dimension of [y(ak ), bk , dk ]T is (2n + 3) × 1. In fact, bk and dk should equal to zero for all k. The two systems characterized above are in common use and one can find expressions for (u1 , τ1 )T . We expect that when we have (u1 , τ1 )T , we can perform the same steps to find (u2 , τ2 ), (u3 , τ3 ) ... etc. iteratively. In fact, there are some hidden problems in this 23.

(29) process. It is the main problem occurs when we find a singular vector which is equal to zero in solving the linear system Ax = b. If there exists y such that Ay = 0, then there exists an α ∈ R such that A(x + αy) = b. Therefore, bifurcations are expected to occur along the curve C. Assume we have a matrix B ∈ Rm×n . Using the Singular Value Decomposition (SVD), we have B = U ΣV T , where U ∈ Rm×m , Σ ∈ Rm×n , and V ∈ Rn×n . The matrices U and V satisfy U T U = I, V T V = I. If an element of Σ is written as σij , then σkk are the singular values of B for all k = 1, 2, · · · , min(n, m), and σij = 0 as i 6= j, where i = 1, 2, · · · , n, and j = 1, 2, · · · , m. We consider the system as follows: BT B = = = T ⇒ (B B)V = = T ⇒ (B B)vi =. (U ΣV T )T (U ΣV T ) V ΣT U T U ΣV T V (ΣT Σ)V T V (ΣT Σ)V T V V (ΣT Σ) σi2 vi. where vi ’s are eigenvectors of B T B, σi2 are the corresponding eigenvalues, and i = 1, 2, · · · , n. We call (σi2 , vi ) ”eigenpairs” for all i. We can use the above method to find the eigenvalues of B T B. We just take their positive square roots, which are the singular values of B. What kinds of situations lead to such bifurcations? Suppose that there are three points ua , ub and uc on the curve C and suppose they are conjoint points of our method. Without a loss of generality, we let the first point be ua , the second point be ub , and the last point be uc . Suppose that the σ value is the smallest singular value of the point u on the curve C, and the numbers 1, 2, and 3 are the placings of σa , σb , and σc . 1 denotes the largest value, 3 the smallest value of them, and 2 an intermediate value. In our cases, the matrix A whether overdetermined or expanded, can be used to find the smallest singular values. There are also other cases which are possible as follows. σa σb σc 1 2 3 1 3 2 ∗ 2 1 3 2 3 1 ∗ 3 1 2 3 2 1 There may be bifurcations in the above two * shapes i.e. (1, 3 ,2) and (2, 3, 1) corresponding to the possibility that ∃ σi = 0. Connecting the three points, we note a V shape, as in Figure 3.3. We know that the singular value may be located at the red star or the green star. We see in the picture on the left of Figure 3.4 that the 24.

(30) distribution of points is linear. Suppose that either ud or ue is the singular value, where the point ud is between the points ua and ub , and the point ue is between the points ub and uc , if they exist. Denote a ratio for calculating the point ud or ue as follows: ud =. σb σa ua + ub , σa + σb σa + σb. ue =. σc σb ub + uc . σb + σc σb + σc. This is seen more clearly in the right hand picture of Figure 3.4.. V shape of (1,3,2). σa. 0.8. σc σ value. σb 0.2. 0.4 σb 0.2. 0. 0. Figure 3.3.. |. |. |. 1. 2 x value. 3. 0 −0.2. 4. 0. |. |. |. 1. 2 x value. 3. 4. The V shape of σa , σb , and σc as (1, 3, 2) and (2, 3, 1). The coordinate values are inexpressive.. V shape of (1,3,2). We give a ratio to find ud or ue.. V shape of (1,3,2), if it is linear. 1. 1 σa. 0.8. σc σb. 0.2 |. 0. |. σc. 0.6 σ value. 0.4. σa. 0.8. 0.6 σ value. σa. 0.6. 0.4. −0.2. σc. 0.8. 0.6 σ value. V shape of (2,3,1). |. 0.4. σb. 0.2 |. 0. −0.2. −0.2. −0.4. −0.4. |. |. −σb 0. Figure 3.4.. 1. 2 x value. 3. 4. 0. 1. 2 x value. 3. 4. We want to find the point which has a zero ( singular ) value point. The left picture shows that the zero singular value maybe happened with a linear distribution of points is linear. The right picture shows the ratio described in the text.. When we find ud or ue , we can calculate the associated smallest singular value, called σd or σe . We can decide on our the next step according to this value. I divide up the area and explain various cases in Figure 3.5. First, we divide the space in two with x = σb . According to σa and σb , we can 25.

(31) further divide the domain into three areas. Lets observe the singular values. The interval [0, σb ] is called the F area, the interval [σb , σa ] is called the E area, and the D area contains a singular value which is larger than σa . The L area, M area, and N area are divided up similarly. When σd falls in the D area, we call it d1 . When σd falls instead into the E area, we call it as d2 . Finally, when σd is located in the F area, we call it as d3 . e1 , e2 , and e3 are named similarly. When we connect the points of σa and d1 , and also the points d1 and σb , ( i.e. σa d1 and d1 σb ), we find that with the resulting shape is not a V shape. We find the same result with σa d2 and d2 σb , but the segments σa d3 and d3 σb do form a V shape. Therefore, we will not be able to judge whether the or not a bifurcation exists as the σd falls into the D or E areas. The case of σe falling in the L or M areas is similar. We will thus continue to iterate while σd or σe lie in the F or N areas. There is also another special case in which we must continue iterations. When σd falls into the E area, we should continue iterating with ue . Similarly if σe falls into the M area, then the segments d2 σb and σb e2 would form a V shape. This is another case continued iterations are necessary. Some distributions of σd and σe. σa. D. 0.8. d1. e1. σ value. 0.6. L σc. d2. E 0.4. M e2 0.2. |. 0. −0.2. σb. d3. F. 0. Figure 3.5.. 0.5. 1. e3. |. 1.5. 2 x value. N |. 2.5. 3. 3.5. 4. Some distributions of σd and σe . All written words are defined as in the above paragraph.. In numerical computations, the user should always specify a tolerance, here called ε1 . When three successive points have smallest singular values less then ε1 , we only need further check whether V shapes are formed. If the solution evolves to a V-shape, 26.

(32) the user will need to give another tolerance called ε2 . If the solution further evolves to having a smallest singular value less then ε2 , that solution is defined as a bifurcation point. After identifying a bifurcation point, we need to find its associated direction. If ui is a bifurcation point for some i, then we would take the penultimate smallest singular value of matrix A for the corresponding ui . If we call this value σps , then one direction of the bifurcation point for ui is the singular vector for the corresponding σps . Another direction is its anti-parallel singular vector. There are therefore three directions in which a bifurcation point may advance. In fact, the σps for A is equal to the penultimate smallest eigenvalue of AT A, and the singular vector for A is the associated eigenvector for AT A. We will now name these bifurcations according to specific numerical experiments. We identify them in order along τ starting at zero so the line of zero to the first bifurcation point is named 1. Suppose the first bifurcation point is called b1 , then there are three directions available to b1 from the above discussion. The three directions will be named 1 1 , 1 2, and 1 3. Then we can see that 1 1 is the original direction of 1. If we let the direction of solutions walk along 1 2, and then we find a new bifurcations, then we name them as 1 21, 1 22, and 1 23, and so forth. The origin direction may be called the main shaft, i.e. the path following i 1, for all i = 1, 2, 3.... We note that the solutions’ shapes are symmetric in i 2 and i 3 for all i, and the resulting shapes can be very interesting. − → Consider a set of information, for example, the domain’s range, the length of t , the stopping tolerances of Newton’s method, the V-shape...etc. We denote variables as λ1 = λ2 = 1, µ1 = µ2 = 1, β12 = β21 = −1, the domain is a square whose coordinates are (−5, −5), (−5, 5), (5, 5) and (5, −5), and the edge is cut at 51 ( we need an odd number. ). Then we will obtain several results. Since there are dozens of related figures we will present them completely in the appendix. We present crude results from the vertical views and only typical shapes and bifurcations as follows.. 27.

(33) 28.

(34) 29.

(35) 30.

(36) 31.

(37) 32.

(38) 33.

(39) We first consider the energy variation. Figure 3.6 shows the relation between τ value and the energy. We find that the energy of the bifurcation is less than that on the main perch. Because the bifurcations happen so frequently, we focus only on the beginning curve. 25. 20 Energy. 1 11 12 111 112 121 122 1121 1122. 15. 10. 0. Figure 3.6.. 0.2. 0.4. τ value. 0.6. 0.8. 1. The relation between τ and energy. They are on the beginning curve.. In order to see the details clearly, we also magnify two parts as follows. One of them goes 12 to 121 and to 122, another goes 112 to 1121 and to 1122. Whether or not τ increases or decreases from the perch to the bifurcation, the energy of the bifurcation is less than the energy on the perch. 22.6 19 22.55 22.5. 1 11 12 111 112 121 122 1121 1122. 18. 17.5. Energy. Energy. 18.5. 22.45 1 11 12 111 112 121 122 1121 1122. 22.4 22.35 22.3. 17. 22.25 0.32. 0.34. Figure 3.7.. 0.36. 0.38 0.4 τ value. 0.42. 0.44. 0.46. 0.59. 0.595. 0.6 τ value. 0.605. 0.61. The relation of τ and energy. We magnify two parts from Figure 3.6 in order to see clearly.. 34.

(40) We will describe this behavior in the following paragraphs. We can see that the solutions will be bell-shaped, so I can portray them using the numbers of bells. From the main picture of he 1 1 bifurcation, we see that the solutions for u1 and u2 will separate to be one-and-two bells or one-and-four bells. It is interesting that the process which separates eventually into four bells happens as a bell becomes a ring in one-in-four cases. It is also only owning one simple bell for all bifurcations at present. The bifurcations of 2 1 and 3 1 are both two-and-two cases. In 2 1, one of the solutions is separated fore and aft, and another is separated right and left. Because we don’t set the rotated plane, the solutions will eddy itself. In 3 1, their bells will likely align with four points. The behaviors of 5 1 and 9 1 are likely. 5 1 is three-and-three case, and 9 1 is four-and-five case. In rudimentary processes they will not separate noticeably, so they behave, respectively, like two-Y shapes and two-+ shapes. Although their processes are likely, their final behaviors are not same. In 9 1, one solution will separate to five bells and another only to become four bells. The behaviors of 6 1 and 7 1 are also likely. They are both four-and-four case and they will become two-line shapes or spheral-cross shapes in their final states. The behavior of 11 1 is like them, but it is a six-and-six case. Some parts of 12 1 are like 11 1, and the two-line shape is transformative as an X in 12 1. 22 1 also has two-line shape and spheral-cross shape, but I am unsure of the number of of bells present. Most types in finding shapes are cross shapes. Some contain many bells in the center of the cross as in 22 1. They also typically become ”skew” in their bifurcations. For example see: 11 1, 16 1, 19 1, 22 1, 24 1, and 29 1.. 35.

(41) 4. The Course of the Discrete Laplacian In this section, we will introduce the process for using Finite Different Method of the discrete Laplacian in different domain. There are rectangle and disc in our interesting domain. We also edit some different boundary conditions in the disc. In the end, we also arrange some domains about holey disc and rectangle. 4.1. The Rectangle with Dirichlet Boundary Conditions of Laplacian   ∆u = f (x, y)       u(x, c) = m(x) , y = c (14) u(b, y) = n(y) , x = b    u(x, d) = p(x) , y = d     u(a, y) = q(y) , x = a. Let the domain of rectangle be a ≤ x ≤ b and c ≤ y ≤ d where a, b, c, d ∈ constant. Suppose ν, ω are the partition numbers of the two sides. Then we can get that the step size h = (b − a)/ν, and k = (c − d)/ω. Therefore we define grid points as following xi = a + i × h,. for i = 0, 1, 2, · · · , ν,. yj = c + j × k,. for j = 0, 1, 2, · · · , ω.. Now, we have ∆u = uxx + uyy . The second differential item which is fixed respectively x and y can be obtained from Taylor series. They are u(xi+1 , yj ) − 2u(xi , yj ) + u(xi−1 , yj ) h2 u(xi , yj+1 ) − 2u(xi , yj ) + u(xi , yj−1 ) = k2. uxx = uyy. We rewrite u(xi , yj ) = ui,j , f (xi , yj ) = fi,j . We put above equations into (14), and use the finite difference method, we can get ∆u = uxx + uyy u(xi+1 , yj ) − 2u(xi , yj ) + u(xi−1 , yj ) u(xi , yj+1 ) − 2u(xi , yj ) + u(xi , yj−1 ) + = h2 k2 = fi,j . 2. Let ρ = hk2 , k = i + (j − 1)(ν − 1) with fk = fi,j , for i = 1, · · · , ν − 1, and j = 1, · · · , ω − 1. We suppose ν = 9 and ω = 4, then we can arrange some equations as following. We can see the discretization index from Figure 4.1.. 36.

(42) The discretization of a rectangle for ν = 9 and ω = 4. 4.5 p(x) 4 3.5 3 17. 18. 19. 20. 21. 22. 23. 24. 9. 10. 11. 12. 13. 14. 15. 16. 1. 2. 3. 4. 5. 6. 7. 8. y value. 2.5 2. q(x). n(x). 1.5 1 0.5 0 m(x) −0.5. Figure 4.1.. 0. 2. 4 x value. 6. 8. 10. It is the discretization of this laplacian for a rectangle. The partitions of two sides are ν = 9 and ω = 4.. (i, j) = (1, 1) : (i, j) = (2, 1) : (i, j) = (3, 1) : (i, j) = (4, 1) : (i, j) = (5, 1) : (i, j) = (6, 1) : (i, j) = (7, 1) : (i, j) = (8, 1) :. 1 1 (u2 − 2u1 + u0,1 ) + 2 (u9 − 2u1 + u1,0 ) 2 h k 1 1 (u3 − 2u2 + u1 ) + 2 (u10 − 2u2 + u2,0 ) h2 k 1 1 (u4 − 2u3 + u2 ) + 2 (u11 − 2u3 + u3,0 ) 2 h k 1 1 (u5 − 2u4 + u3 ) + 2 (u12 − 2u4 + u4,0 ) 2 h k 1 1 (u6 − 2u5 + u4 ) + 2 (u13 − 2u5 + u5,0 ) 2 h k 1 1 (u7 − 2u6 + u5 ) + 2 (u14 − 2u6 + u6,0 ) 2 h k 1 1 (u8 − 2u7 + u6 ) + 2 (u15 − 2u7 + u7,0 ) 2 h k 1 1 (u9,1 − 2u8 + u7 ) + 2 (u16 − 2u8 + u8,0 ) 2 h k. = f1 = f2 = f3 = f4 = f5 = f6 = f7 = f8. ================================================ (i, j) = (1, 2) : (i, j) = (2, 2) : (i, j) = (3, 2) : (i, j) = (4, 2) : (i, j) = (5, 2) : (i, j) = (6, 2) : (i, j) = (7, 2) : (i, j) = (8, 2) :. 1 1 (u10 − 2u9 + u0,2 ) + 2 (u17 − 2u9 + u1 ) 2 h k 1 1 (u11 − 2u10 + u9 ) + 2 (u18 − 2u10 + u2 ) h2 k 1 1 (u12 − 2u11 + u10 ) + 2 (u19 − 2u11 + u3 ) 2 h k 1 1 (u13 − 2u12 + u11 ) + 2 (u20 − 2u12 + u4 ) 2 h k 1 1 (u14 − 2u13 + u12 ) + 2 (u21 − 2u13 + u5 ) 2 h k 1 1 (u15 − 2u14 + u13 ) + 2 (u22 − 2u14 + u6 ) 2 h k 1 1 (u16 − 2u15 + u14 ) + 2 (u23 − 2u15 + u7 ) 2 h k 1 1 (u9,2 − 2u16 + u15 ) + 2 (u24 − 2u16 + u8 ) 2 h k 37. = f9 = f10 = f11 = f12 = f13 = f14 = f15 = f16.

(43) ================================================ 1 1 (i, j) = (1, 3) : (u − 2u + u ) + (u1,4 − 2u17 + u9 ) = f17 18 17 0,3 h2 k2 1 1 (i, j) = (2, 3) : (u − 2u + u ) + (u2,4 − 2u18 + u10 ) = f18 19 18 17 h2 k2 1 1 (i, j) = (3, 3) : (u − 2u + u ) + (u3,4 − 2u19 + u11 ) = f19 20 19 18 h2 k2 1 1 (i, j) = (4, 3) : (u − 2u + u ) + (u4,4 − 2u20 + u12 ) = f20 21 20 19 h2 k2 1 1 (i, j) = (5, 3) : (u − 2u + u ) + (u5,4 − 2u21 + u13 ) = f21 22 21 20 h2 k2 1 1 (i, j) = (6, 3) : (u − 2u + u ) + (u6,4 − 2u22 + u14 ) = f22 23 22 21 h2 k2 1 1 (i, j) = (7, 3) : (u24 − 2u23 + u22 ) + 2 (u7,4 − 2u23 + u15 ) = f23 2 h k 1 1 (i, j) = (8, 3) : (u9,3 − 2u24 + u23 ) + 2 (u8,4 − 2u24 + u16 ) = f24 2 h k ================================================ Put boundary conditions m(xi ) = mi,0 , n(yj ) = nν,j , p(xi ) = pi,ω , and q(yj ) = q0,j . For example u1,0 and u0,1 : u1,0 = m1,0 and u0,1 = q0,1 , then we get 1 1 1 1 1 1 −2( 2 + 2 )u1 + 2 u2 + 2 u9 = f1 − 2 q0,1 − 2 m1,0 . h k h k h k Use the same step for other condition items about u2 , · · · , u8 , u9 , u16 , and u17 , · · · , u24 . We will solve the linear system question for Au = b. Furthermore, the matrix   1 −2 k12 I − h12 Ψ I 2 k   1 1 1 A =  I −2 I − h12 Ψ I  k2 k2 k2 1 I −2 k12 I − h12 Ψ k2   −2ρI − Ψ ρI 1   24×24 = 2 , ρI −2ρI − Ψ ρI ∈R h ρI −2ρI − Ψ where ρ =. h2 , k2. and. . 2. −1.   −1 2 Ψ=  .. .  The right hand side b is  b1   b =  b2  b3 . .    where b1 =    .     ∈ R8×8 .  .. . −1  −1 2 .... f1 − h12 q0,1 − k12 m1,0 f2 − k12 m2,0 .. . f7 − k12 m7,0 f8 − h12 n9,1 − k12 m8,0 38. . .        , b2 =       . f9 − h12 q0,2 f10 .. . f15 f16 − h12 n9,2.     ,   .

(44)     and b3 =    . f17 − h12 q0,3 − k12 p1,4 f18 − k12 p2,4 .. . f23 − k12 p7,4 f24 − h12 n9,3 − k12 p8,4.     .   . After we arrange a example for ν = 9 and ω = 4, we can write down the general case for any ν and ω. The matrix   −2ρI − Ψ ρI   ρI −2ρI − Ψ ρI     1  . . .  ∈ R(ν−1)×(ω−1) , .. .. .. A= 2  h     ρI −2ρI − Ψ ρI ρI −2ρI − Ψ where ρ =. h2 , k2. and . 2. −1.   −1 2 Ψ=  .. . .     ∈ R(ν−1)×(ν−1) .  .. . −1  −1 2 .... The right hand side b is     b1 f1,1 − h12 q0,1 − k12 m1,0     f2,1 − k12 m2,0   b2     .   .  .. ..  where b1 =  b=         1   bω−2   fν−2,1 − k2 mν−2,0 fν−1,1 − h12 nν,1 − k12 mν−1,0 bω−1    f1,ω−1 − h12 q0,ω−1 − k12 p1,ω f1,j − h12 q0,j    f2,j f2,ω−1 − k12 p2,ω       .. ..  for j = 2, 3, · · · , ω−2, and bω−1 =  bj =  . .          fν−2,ω−1 − k12 pν−2,ω fν−2,j fν−1,ω−1 − h12 nν,ω−1 − k12 pν−1,ω fν−1,j − h12 nν,j 4.2. The Disc with Dirichlet Boundary Conditions of Laplacian ( ∆u = f (x, y) u(x, y) = c(x, y) , (x, y) ∈ Ω. (15). 2 be the radial mesh width and ∆θ = 2π be the azimuthal mesh Let ∆r = 2ν+1 ω width, where ν and ω is positive integer. Then we define grid points as following. 1 ri = (i − )∆r, 2 θj = (j − 1)∆θ,. for i = 1, 2, · · · , ν, for j = 1, 2, · · · , ω. 39.     .   .