Semilinear Elliptic Problem with Neumann Boundary Conditions
Tsung-Min Hwang
Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan
Weichung Wang
Department of Mathematics Education, National Tainan Teachers College, Tainan 700, Taiwan
Received 3 May 2001; accepted 21 August 2001
A semilinear elliptic equation d∆u − u + up = 0 over the unit ball in R2 with positive solution and the homogeneous Neumann boundary condition is considered. This equation models applications like chemotactic aggregation and biological pattern formation. Focusing on solving the discretized version of the equation, this work proposes an efficient algorithm that combines a newly developed discretization scheme on polar coordinates with a fast Fourier solver. An analysis of the induced matrix structures proves the algorithm converges to positive solutions; the analysis also establishes the q-axial symmetry and monotonicity behavior of the solutions. Numerical experiments were conducted to visualize various solution forms that are new to the best of our knowledge. The experiments also illustrated sensitivity behavior of the solutions. c2002 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 18: 261–279, 2002; Published online in Wiley InterScience (www.interscience.wiley.com). DOI 10.1002/num.10006
Keywords: semilinear Neumann problem; discretization on polar coordinates; q-axial symmetry; mono- tonicity; numerical methods; visualization
I. INTRODUCTION
We consider positive solutions of the following semilinear elliptic equation in B1 ⊆ Rn subject to the homogeneous Neumann boundary condition
Correspondence to: Tsung-Min Hwang, Department of Mathematics, National Taiwan Normal University, 88, Sec. 4 Ting-Chou Road, Taipei 116, Taiwan (e-mail: [email protected])
Contract grant sponsors: National Science Council and National Center for Theoretical Sciences in Taiwan.
c
2002 Wiley Periodicals, Inc.
d∆u − u + up= 0 in B1, u > 0 in B1,
∂u
∂*n = 0 on ∂B1,
(1.1)
where d is a positive constant, ∆ is the Laplace operator, 1 < p < n+2n−2 for n ≥ 3, 1 < p < ∞ for n = 2, and*n is the unit outer normal to ∂B1. In particular, we let B1be the unit ball in R2 throughout this article.
The partial differential equation arises in many applications. For example, Keller and Segel [1]
reduced the steady state problem for a chemotactic aggregation model with logarithmic sensitivity to the form of (1.1). Gierer and Meinhardt [2] studied activator-inhibitor systems modeling biological pattern formation. In the systems, the equation (1.1) is essential when the diffusion rate of the inhibitor is sufficiently large. Moreover, Ni and Takagi [3] treated the equation as the shadow system of some reaction-diffusion system in morphogenesis.
These applications also lead to many recent researches regarding the solution properties to equation (1.1). Let an "energy" function Jd : W1,2(B1) → R associated with (1.1) be defined by
Jd(v) = Z
B1
[1
2(d|∇v|2+ v2) − 1
p + 1vp+1+ ]dx,
where v+= max(v, 0). The mountain-pass lemma [4] suggests that the constant cd= inf
h∈Γ max
0≤t≤1Jd(h(t))
is a positive critical value of Jd, where Γ = {h ∈ C([0, 1], W1,2(B1)) | h(0) = 0, h(1) = e}
and e ∈ W1,2(B1), e 6= 0 in B1with Jd(e) = 0. Furthermore, cd is the least positive critical value as shown in [5]. A critical point ud of Jd with critical value cd is called a least-energy solution of (1.1). In [6], Ni and Takagi showed that, provided d is small enough, udhas only one local maximum point over B1, and it must lie on the boundary. Furthermore, let P0be a local maximum point of a least-energy solution u of (1.1) on B1and O be the origin. Chern and Lin [7] showed that either u is a constant or P06= O and u is axially symmetric with respect to −−→OP0. They also proved that the single peak solution u is monotonically decreasing along the each sphere with same radial coordinate. Lin [8] asserted that, if a least-energy solution is not a constant, the maximum point must be on the boundary. Although the analytical results have provided many important insights to the problem, further understanding toward the equation remains challenging and interesting because of its complexity.
On the other hand, numerical methods provide effective tools for exploring the solution be- haviors. Iterative type of algorithms are common numerical techniques used to solve nonlinear elliptic equations. Chen, Ni, and Zhou [9] discussed mountain-pass, scaling iterative, direct it- eration, and monotone iteration algorithms. They solved nonlinear partial differential equations by iterative algorithms and surveyed related references therein. Although the iterative type algo- rithms may solve the model equation in some cases, some kinds of solutions are not easy to catch by straightforward implementations.
Concentrating on solving the discretized version of the problem, we present following results in this article.
1. We first propose an iterative algorithm that results in the discovery of q-axial symmetric solutions. A solution is q-axial symmetric if the solution is invariant under the reflection with respect to axes ←→OPi, where O is the origin and Pi’s are points evenly distributed on the boundary for i = 1, . . . , q.
2. We seamlessly combine a recently developed discretization scheme on polar coordinates [10] and a variant of fast Fourier solver [11, 12] to efficiently solve the induced linear sys- tems in the iterations. We also assert the scheme converges to q-axial symmetric solutions analytically and numerically.
3. By combining the above two techniques (the new discretization scheme and fast Fourier solver), we can not only fully take advantage of the block diagonal matrix structure to accelerate the solving processes, but prove the discretized solutions are positive, q-axial symmetric, and satisfy a monotonicity behavior.
We use properties of M-matrix to prove the discretized solutions are positive. The main tools used in the proof of q-axial symmetry include sine- and cosine-based orthonormal similar transformation, symmetry, and periodicity of sine and cosine, and block diagonal matrix structure. Applying similar procedures, we show that starting from the q apexes on the axes of symmetry, each peak of the q-axial symmetric solutions is monotonically decreasing along the grid point circles. A grid point circle is formed by the grid points that have same radial coordinate.
4. Several q-axial symmetric solution shapes are visualized. The numerical results demon- strate diversity of the solutions that have never been discovered before to our best knowl- edge. From the numerical viewpoint, we further explore sensitivity of the solutions.
Namely, the single peak solution is quite stable, whereas q-axial symmetric solutions are sensitive to the initial approximations. Besides, we numerically show that the peaks of multiple-peak solutions lie on the boundary of the domain.
The article is organized as follows. Section II. introduces a discrete iterative algorithm, as- sociated with the discretization scheme and a fast Fourier solver. The convergence analysis of the discrete iterative algorithm is also demonstrated in this section. Section III. performs ma- trix analysis to explore discretized solution behaviors. We prove that the discretized solutions are positive and, under mild assumptions, the solutions are q-axial symmetric and retain the monotonicity behavior. Section IV. illustrates visualization results and discusses numerical expe- riences. The numerical experiments were conducted to explore the q-axial symmetry, sensitivity, and monotonicity behaviors of the solutions. Section V. finally concludes the article.
II. NUMERICAL ALGORITHMS
We consider numerical schemes used to solve the partial differential equation (1.1) in this section.
We first introduce an iterative algorithm incorporated with finite difference scheme in Figure 1.
The iterative algorithm solves a sequence of discretized partial differential equation:
d∆euk+1− euk+1= −upk, euk+1> 0 in B1,
∂euk+1
∂*n = 0 on ∂B1, (2.1)
where euk+1and uk are the unknown and known discrete values of u, respectively. The iterative algorithm shown in Figure 1 presents the idea for solving the problem. One essential feature of
FIG. 1. Iterative algorithm.
the algorithm is to normalize the iterates euk at each iteration. Without doing so, the right hand side of the first equation in (2.1) either converge to zero or blow up to infinity and the iterative process tends to trivial solutions.
A. Discretization scheme
To discretize the equation over the unit ball domain, it is convenient to rewrite the equation in the polar coordinates system. Such transformation, however, causes the coordinate singularity at the origin even the solution itself is regular. Numerical schemes such as pole conditions are thus generally used [13, 14, 15]. The accuracy of such methods greatly depends on the choice of the conditions. To overcome the difficulty, Lai [10] recently proposed a simple and elegant discretization scheme. We elaborate our discretization scheme based on [10] as follows.
Let ∆r = 2ν+12 be the radial mesh width and ∆θ = 2πω be the azimuthal mesh width, for positive integers ν and ω. The grid locations are then half-integered in radial direction and integered in azimuthal direction. In other words,
ri = (i −1
2)∆r and θj = (j − 1)∆θ,
for i = 1, . . . , ν and j = 1, . . . , ω. The standard central finite difference method discretizes equation (2.1) and forms the linear system to be solved in Step 1 of Figure 1 as
Au = b, (2.2)
where the unknown vector u in (2.2) is represented by
u =
u1
u2
u...ν
and ui=
ui1
ui2
u...iω
, for i = 1, . . . , ν, (2.3)
which stands for ˜uk+1in (2.1) with uij ≈ u(ri, θj). Furthermore, the matrix
A =
δI − β1Ψ (1 + µ1)I
(1 − µ2)I δI − β2Ψ (1 + µ2)I
... ... ...
(1 − µν−1)I δI − βν−1Ψ (1 + µν−1)I (1 − µν)I ηI − βνΨ
, (2.4)
where µi =2i−11 ; βi =(i−1/2)12(∆θ)2, for i = 1, . . . , ν; δ = −2 −(∆r)d 2; η = −1 + µν−(∆r)d 2; and
Ψ =
2 −1 −1
−1 2 ...
... ... −1
−1 −1 2
∈ Rω×ω. (2.5)
The right hand side b in (2.2) is given by
b =
b1 b2 b...ν
with bi=
bi1 bi2 b...iν
= −(∆r)2 d
(uk)pi1 (uk)pi2 (uk...)piω
for i = 1, . . . , ν, (2.6)
where uk ∈ Rωνis obtained in the previous iterative step. Therefore, the Step 1 in Algorithm 1 can be substituted by
Step 1 Solve the linear systems A˜uk+1= b(k).
B. Fast Fourier Solver
A fast Fourier solver is used to solve the linear system (2.2). The main idea of the scheme is to transform the ων × ων matrix A in (2.2) into a ω × ω block diagonal matrix that each block is a tridiagonal ν ×ν matrix. Note that the linear system can also be solved by (block) cyclic reduction methods proposed by [16] and [11]. Although both of the fast Fourier and cyclic reduction methods have similar computational cost, the fast Fourier solver, however, can be easily parallelized and is somewhat easier to be implemented. Moreover, the ω × ω block diagonal matrix generated by the fast Fourier scheme has nice structure and properties. We can then use the structures and properties to prove the existence of positive q-axial symmetric solutions.
It is well known (see e.g., [17]) that the symmetric circulant matrix Ψ in (2.5) can be explicitly diagonalized by an ω × ω orthogonal matrix W , that is,
W>ΨW = diag(λ1, . . . , λω) ≡ Λ, (2.7) where W ≡ [w1, . . . , wω] ≡ [wjk]ωj,k=1with
wjk=
τ cos[(j − 1)(k − 1)∆θ], for 1 ≤ k ≤ ω/2 + 1,
τ sin[(j − 1)(k − ω/2 − 1)∆θ], for ω/2 + 2 ≤ k ≤ ω, (2.8)
in which τ =p
2/ω for 1 < k ≤ ω, k 6= ω/2 + 1, and τ =p
1/ω, for k = 1 and k = ω/2 + 1;
and
λi=
4 sin2[(i − 1)π/ω], for 1 ≤ i ≤ ω/2 + 1,
4 sin2[(i − ω/2 − 1)π/ω], for ω/2 + 2 ≤ i ≤ ω, (2.9) provided that ω is an even number. Note that for the case where ω is an odd number, ω/2 + 1 in (2.8), (2.9) is replaced by (ω + 1)/2, as well as τ =p
2/ω for 1 < k ≤ ω, and τ =p 1/ω, for k = 1.
The equation (2.7) suggests a discrete iterative algorithm with fast Fourier solver in Figure 2.
Let cW = diag(W, . . . , W ) and let
Π =
Π1
Π2
Π...ω
with Πj =
e>j e>ω+j e>2ω+j e>(ν−1)ω+j...
, for j = 1, 2, . . . , ω, (2.10)
be the permutation matrix, where eiis the ith column of the unit matrix. From (2.4), (2.7), and (2.10), we have
ΠcW>AcW Π> = bA = diag( bA1, . . . , bAω), (2.11)
FIG. 2. Discrete iterative algorithm with fast Fourier solver.
where bAiis the tridiagonal matrix of the form
Abi=
δ − β1λi 1 + µ1
1 − µ2 δ − β2λi 1 + µ2
... ... ...
1 − µν−1 δ − βν−1λi 1 + µν−1
1 − µν η − βνλi
, (2.12)
for i = 1, . . . , ω. Therefore, the linear system Au = b in (2.2) becomes
Abu = bb,b (2.13)
where bu = ΠcW>u and bb = ΠcW>b.
C. The Discrete Iterative Algorithm
Figure 2 illustrates the discrete iterative algorithm that combines the numerical schemes discussed in Section A. and B. Theorem 2.1 demonstrates the convergence analysis of the algorithm.
Theorem 2.1. Assume p ≥ 2. Choose any u0(x) > 0 on B1and let αk+1and uk+1be defined as in Figure 2. Then (uk+1, αk+1) has a convergent subsequence (ukj, αkj) → (u∞, α∞) such that
u ≡ α1/(p−1)∞ u∞ is a solution to Ax = −xp.
Proof. By definitions of vector norms, we have k upk k1=
Xν j=1
Xω i=1
(uk)pji=k uk kpp≤k ukkp2= 1. (2.14)
On the other hand,
k uk kp≤k ukkq≤ (νω)(p−q)/(pq)k uk kp, for p ≥ q.
Take q = 2. Then
1 =k ukkp2≤ (νω)(p−2)/2k uk kpp. (2.15) Combining (2.14) and (2.15), we get
(νω)(2−p)/2≤k uk kpp≤ 1.
Representing the solution ˜uk+1by
˜uk+1= −A−1upk, we have
k ˜uk+1kp≥k upkkpmin
y6=0
k A−1ypkp
k ypkp =k upk kp
maxx6=0
k Ax kp
k x kp
−1
≥ (νω)(2−p)/2k A k−1p . On the other hand,
k ˜uk+1kp≤k A−1 kpk upk kp≤k A−1 kpk upk k1≤k A−1 kp.
This implies that
(νω)(2−p)/2k A k−1p ≤k ˜uk+1kp≤k A−1kp.
Thus, the sequence {˜uk} is bounded. There exist subsequences {˜ukj} and {αkj} of {˜uk} and {αk}, respectively, such that (˜ukj, αkj) → (˜u∞, α∞). From the definition of uk follows that ukj → u∞≡ α∞˜u∞. From (2.1), we have that
Au = A(α1/(p−1)∞ u∞) = αp/(p−1)(Aα∞−1u∞)
= −αp/(p−1)∞ up∞= −(α∞1/(p−1)u∞)p= −up.
III. SOLUTION BEHAVIORS
In this section, matrix analysis reveals three main behaviors of the solutions computed by the discrete iterative algorithm with Fast Fourier solver shown in Figure 2. First, Theorem 3.1 shows that the convergent solution is positive provided that the initial approximation is chosen to be positive. Second, Lemma 3.2 and Theorem 3.3 further show that the solutions of the linear system (2.2) preserve the q-axial symmetry if the right hand side possess the q-axial symmetry.
Third, Theorem 3.4 demonstrates the solutions of (2.2) retain the monotonicity behavior.
A. Positive Solutions
To show the algorithm shown in Figure 2 leads to positive convergent solutions, we only need to prove the following theorem.
Theorem 3.1. If the right hand side vector of the linear system (2.2) is negative, then the solution of the system is positive.
Proof. Let b and u be the right hand side vector and the solution of the linear system in (2.2), respectively. Let
A ≡ ΠAΠ˜ >=
T − 2D D D
D T − 2D ...
... ... D
D D T − 2D
,
where
T =
δ 1 + µ1
1 − µ2 δ 1 + µ2
... ... ...
1 − µν−1 δ 1 + µν−1
1 − µν η
(3.1)
and D = diag(β1, . . . , βν). Since 1−µj+1> 0 and 1+µj > 0 for j = 1, . . . , ν −1, there exists D = diag(d˜ 1, . . . , dν), where d1is any positive number and dj+1= djp
(1 − µj+1)/(1 + µj) for j ≥ 1, such that ˜T ≡ ˜D−1T ˜D is a symmetric tridiagonal matrix with off-diagonal entries γj = p
(1 + µj)(1 − µj+1) for j = 1, . . . , ν − 1. Let ˆD = diag( ˜D, . . . , ˜D). Then ¯A ≡
Dˆ−1A ˆ˜D is symmetric. Therefore, the linear system in (2.2) is equivalent to
A ˆ¯D−1Πu = ˆD−1Πb. (3.2)
Because ¯A is symmetric and strictly diagonally dominant with nonnegative off-diagonal entries,
− ¯A is an M-matrix, that is − ¯A−1 ≥ 0. From (3.2) follows that u > 0.
B. The q-axial Symmetry
A solution is called q-axial symmetric if the solution is axially symmetric with respect to q axes that evenly divide the unit ball into q parts. Let f be a grid function on (ri, θj) for i = 1, . . . , ν and j = 1, . . . , ω and denote fi,j= f(ri, θj). Let ω = q(2m + 1), for arbitrary positive integers q and m, be the partition numbers on the azimuthal direction. The function f is called to be q-axial symmetric over the grid points (ri, θj) if f satisfies the following conditions for each 1 ≤ i ≤ ν:
fi,2m(¯j−1)+¯j= fi,1, fi,¯ı+1= fi,ω−¯ı+1,
fi,2m(¯j−1)+¯j+¯ı= fi,2m(¯j−1)+¯j−¯ı= fi,¯ı+1, (3.3) for ¯j = 2, . . . , q and ¯ı = 1, . . . , m.
We illustrate the q-axial symmetry conditions by the discretization shown in Figure 3. The figure shows the grid points, denoted by ∗, ◦, and +, for ν = 4, ω = 21, and q = 3. Let P1, P8, and P15be the points lying on the boundary with 0, 120, and 240 degrees in angular coordinates.
•
The first condition in (3.3) ensures that the corresponding function values on OP1, OP8, and OP15are equal. Namely, fi,1= fi,8= fi,15for i = 1, 2, 3, 4.•
The second condition guarantees that the values at m (which is equal to 3 in the exam- ple) points lying on the upper and lower sides of OP1oppositely are equal to each other correspondingly. That is, fi,2= fi,21, fi,3= fi,20, and fi,4= fi,19for i = 1, 2, 3, 4.FIG. 3. Grid points for ν = 4, ω = 21, and q = 3.
•
The second symmetry property discussed above is guaranteed for OP8and OP15by the third condition. Furthermore, all the corresponding points must be equal. In other words,fi,2= fi,21= fi,9= fi,7= fi,16= fi,14, fi,3= fi,20= fi,10= fi,6= fi,17= fi,13, fi,4= fi,19= fi,11= fi,5= fi,18= fi,12,
for i = 1, 2, 3, 4.
Next we will show the solution of (2.2) satisfy the q-axial symmetric conditions with any positive integer q. In order to prove the main result in Theorem 3.3, we first show the following useful lemma.
Lemma 3.2. Let ω = q(2m+1) for arbitrary given positive integer m and even positive integer q. Suppose v ≡
v1 · · · vω>
∈ Rωand satisfies the q-axial symmetric conditions (3.3). Let W ≡
w1 · · · wω
be defined as in (2.8). Then
v>wk = 0, for k = ω/2 + 2, . . . , ω, (3.4) v>wqk+` = 0, for k = 0, . . . , m, ` = 2, . . . , q,
and wqk+1preserve the symmetric properties:
wqk+1(2m(j − 1) + j) = wqk+1(1),
wqk+1(i + 1) = wqk+1(ω − i + 1),
wqk+1(2m(j − 1) + j + i) = wqk+1(2m(j − 1) + j − i) = wqk+1(i + 1),
for k = 0, . . . , m, j = 2, . . . , q, and i = 1, . . . , m.
Proof. For convenience, the indices 2m(j − 1) + j and qk + ` are denoted by ¯j and ¯k, respectively.
By the definition of wω/2+1+¯k, we have
wω/2+1+¯k(ω − i + 1) = τ sin(2π − i∆θ) = −wω/2+1+¯k(i + 1), (3.5) wω/2+1+¯k(¯j) = τ sin
(` − 1)(j − 1)
q 2π
, (3.6)
wω/2+1+¯k(¯j + i) = τ sin
(` − 1)(j − 1)
q 2π + (¯k − 1)i∆θ
, (3.7)
wω/2+1+¯k(¯j − i) = τ sin
(` − 1)(j − 1)
q 2π − (¯k − 1)i∆θ
, (3.8)
for `, j = 1, . . . , q and i, k = 1, . . . , m. Using results (3.5)-(3.8), the identities ([18, p.193])
Xn j=2
sin
2π(j − 1)`
n
= 0, for n ≥ 2 and all integers `, (3.9) Xn
j=1
cos
2π(j − 1)`
n
=
n, if ` is a multiple of n,
0, otherwise, (3.10)
and the symmetry property of v, we have
v>wω/2+1+¯k
=
Xq j=1
v1wω/2+1+¯k(¯j)
+Xm
i=1
vi+1h
wω/2+1+¯k(i + 1) + wω/2+1+¯k(ω − i + 1)i
+ Xm i=1
Xq j=2
vi+1
hwω/2+1+¯k(¯j + i) + wω/2+1+¯k(¯j − i)i
= τXm
i=1
vi+1 Xq j=2
sin
(` − 1)(j − 1)
q 2π + (¯k − 1)i∆θ
+ sin
(` − 1)(j − 1)
q 2π − (¯k − 1)i∆θ
= τ Xm i=1
2vi+1cos[(¯k − 1)i∆θ]
Xq j=2
sin
(` − 1)(j − 1)
q 2π
= 0.
This implies that
v>wk = 0, for k = ω/2 + 2, . . . , ω.
Similarly, by the definition of w¯k, we have
w¯k(ω − i + 1) = τ cos(2π − i∆θ) = w¯k(i + 1), (3.11) w¯k(¯j) = τ cos
(` − 1)(j − 1)
q 2π
, (3.12)
w¯k(¯j + i) = τ cos
(` − 1)(j − 1)
q 2π + (¯k − 1)i∆θ
, (3.13)
w¯k(¯j − i) = τ cos
(` − 1)(j − 1)
q 2π − (¯k − 1)i∆θ
, (3.14)
for `, j = 1, . . . , q and i, k = 1, . . . , m. Using (3.10), (3.11)-(3.14) and the symmetry of b1, we have
v>w¯k = Xq j=1
v1w¯k(¯j) +Xm
i=1
vi+1[w¯k(i + 1) + w¯k(ω − i + 1)]
+Xm
i=1
Xq j=2
vi+1[w¯k(¯j + i) + w¯k(¯j − i)]
= τv1 Xq j=1
cos
(` − 1)(j − 1)
q 2π
+ 2τXm
i=1
vi+1cos[(¯k − 1)i∆θ]
+τXm
i=1
2vi+1cos[(¯k − 1)i∆θ]
Xq j=2
cos
(` − 1)(j − 1)
q 2π
= τXm
i=1
2vi+1cos[(¯k − 1)i∆θ]
Xq j=1
cos
(` − 1)(j − 1)
q 2π
= 0,
for ` = 2, . . . , q. Substituting ` = 1 into (3.11)-(3.14), we get wqk+1(ω − i + 1) = wqk+1(i + 1),
wqk+1(¯j) = 1 = wqk+1(1),
wqk+1(¯j + i) = wqk+1(¯j − i) = cos(qki∆θ) = wqk+1(i + 1), for i, k = 0, . . . , m and j = 2, . . . , q.
The Lemma 3.2 also holds for the case where q is odd, provided the initial index ω/2 + 2 in (3.4) is replaced by (ω + 1)/2 + 1.
Theorem 3.3. Let ω = q(2m + 1) for arbitrary given positive integer q and m. Let u and b be defined as in (2.3) and (2.6), respectively. Then u satisfies the q-axial symmetric conditions (3.3) provided that b satisfies the conditions.
Proof. Let W ≡ [w1 · · · wω] be defined in (2.8). From (2.11), the solution of (2.2) is given by
u = cW Π>diag
Aˆ−11 · · · ˆA−1ω ΠcW>b
= cW Π>diag
Aˆ−11 · · · ˆA−1ω
˜b1
˜b...ω
,
where
˜bi=
b>1wi
b>ν...wi
, for i = 1, . . . , ω.
Denote ˆA−1i ≡h a(i)kji
∈ Rν×ν. Then
u = cW Π>
Pν
r=1a(1)1rb>rw1
Pν ...
r=1a(1)νrb>rw1
Pν ...
r=1a(ω)1r b>rwω Pν ...
r=1a(ω)νr b>rwω
= cW
Pν
r=1a(1)1rb>rw1
Pν ...
r=1a(ω)1r b>rwω
Pν ...
r=1a(1)νrb>rw1 Pν ...
r=1a(ω)νr b>rwω
≡
u1
u...ν
,
where uj = Xω
k=1
" ν X
r=1
a(k)jr b>rwk
# wk
=
q/2+1X
`=1
Xm k=0
" ν X
r=1
a(qk+`)jr b>rwqk+`
# wqk+`
+ Xq
`=q/2+2 m−1X
k=0
" ν X
r=1
a(qk+`)jr b>rwqk+`
#
wqk+`+ Xω
k=ω/2+2
" ν X
r=1
a(k)jr b>rwk
# wk,
(3.15) for j = 1, . . . , ν provided that ω is an even number. For the case that ω is an odd number, q/2 and ω/2 + 2 in (3.15) are replaced by (q − 1)/2 and (ω + 1)/2 + 1, respectively, and all the rest of proof remains unchanged. By the symmetry of b and Lemma 3.2, ujcan be simplified by
uj= Xm k=0
" ν X
r=1
a(qk+1)jr b>rwqk+1
# wqk+1.
By symmetric properties of wqk+1, the q-axial symmetry conditions (3.3) for u is obtained.
Combining the symmetry result with the convergence result in Theorem 2.1, we conclude that the convergent solutions of the problem are also q-axial symmetric.
We now show that the q-axial symmetric solutions have a monotone behavior in the following sense. Let a grid point circle be the grid points that have same radial coordinate. For a q-axial symmetric solution, the highest solution value of each grid point circle are located at the q apexes.
The solution values on two sides of the apexes are monotonically decreasing along the grid point circles. Note that Chern and Lin [7] proved the monotone behavior of continuous solution for the single peak axial symmetric solution.
Theorem 3.4. Let ω = q(2m + 1) for arbitrary given positive integer q and m. Let u and b be defined as in (2.3) and (2.6), respectively. Suppose b satisfies the q-axial symmetric conditions (3.3) and
bi,j+1− bi,j> 0, for i = 1, . . . , ν, j = 1, . . . , m.
Then u satisfies
ui,j+1− ui,j < 0, for i = 1, . . . , ν, j = 1, . . . , m. (3.16)
Proof. Let ˆu =
ˆu3/2 ˆu5/2 ˆum+1/2...
and ˆuj+1/2=
ˆu1,j+1/2 ˆu2,j+1/2 ˆuν,j+1/2...
≡
u1,j+1− u1,j
u2,j+1− u2,j
uν,j+1...− uν,j
and
ˆb =
ˆb3/2 ˆb5/2 ˆbm+1/2...
and ˆbj+1/2=
ˆb1,j+1/2 ˆb2,j+1/2 ˆbν,j+1/2...
≡
b1,j+1− b1,j b2,j+1− b2,j bν,j+1...− bν,j
> 0.
From Theorem 3.3, u satisfies the q-axial symmetric conditions. Thus
ui,ω = ui,2 and ui,m+2= ui,m+1 for i = 1, . . . , ν. (3.17) On the other hand, for j = 1, . . . , m,
ˆuν+1,j+1/2= uν+1,j+1− uν+1,j = uν,j+1− uν,j= ˆuν,j+1/2. (3.18) Based on the linear system in (2.2) and the results in (3.17) and (3.18), the linear system with unknown vector ˆu and right hand side vector ˆb is
C ˆu = ˆb, where
C =
T − 3D D
D T − 2D D
... ... ...
D T − 2D D
D T − 2D
,
in which D = diag(β1, . . . , βν) and T is defined in (3.1). As in the proof of Theorem 3.1, there exists a diagonal matrix ˆD with positive diagonal entries such that − ˆD−1C ˆD is an M-matrix. It follows that ˆu < 0, that is, the inequality in (3.16) holds.
IV. NUMERICAL RESULTS
We conducted numerical experiments in computing solutions for the problem with different initial approximations. Our experiments focused on q-axial symmetry, sensitivity, and monotonicity behaviors of the solutions. The discrete algorithm shown in Figure 2 was used in our experiments.
To determine whether the iteration is convergent, we checked two criteria. First, the difference between consecutive iterates uk must be less than 10−13 in 2-norm, i.e., k uk+1− uk k2<
10−13. Second, the residual of the linear system (2.2) must be less than 10−13. Throughout the experiments, the constant d = 0.01 and p = 2.
In the first experiment set, we discovered that q-axial symmetric solutions can be obtained if the q-axial symmetry conditions hold for the initial approximations and all intermediate iterates.
FIG. 4. The convergence solutions with 1-, 2-, 3-, and 8-peaks.
The initial approximations were chosen so that different number of peaks lie on the boundary and the approximations conform to the q-axial symmetry conditions (3.3). After solving the linear system (2.2) in Step 1 of Figure 2 in each iteration, we corrected some of the solution components so the modified iterates satisfy the q-axial symmetry conditions exactly in numerical. Various multiple peak solutions were then obtained. Note that such modification can be justified by Theorem 3.3. In contrast, without performing the modifications, the iterations would approach q- axial symmetric solutions and then converged to single peak solutions eventually. Such behavior occurred since the computed solutions were not q-axial symmetric numerically, even they should be in theory. The slight difference between the computed and theoretical solutions led to the single peak solutions. These sensitive behavior were further explored in the second experiment set. Figure 4 illustrates the convergent solutions. Although only the solutions with single, two, three, and eight peaks are plotted in the figure, many other q-axial symmetric solutions were obtained in our tests.
In the second experiment set, we constructed initial approximations by perturbing q-axial symmetric solutions obtained in the first experiment set. We observed that the perturbations result in convergence to single peak solution. The following highlight the findings in two- and three-peak solution cases.
•
Small perturbation on peak heights leading to one peak solution.As shown in part (b) of Figure 4, two (equivalent) largest solution values occur in the grid points that are closest to (1, 0) and (−1, 0). We slightly lowered an apex so that the height ratio of the lower one to the higher one is 0.9999. One or two peaks in the three-peak
FIG. 5. The single peak moves from interior of the domain to the boundary.
solution were slightly lowered similarly. All of the initial approximations converged to the single peak solution in the form of part (a) in Figure 4.
•
Small perturbations of the peak locations leading to one peak solution.Two of the three peaks in the three-peak solution (part (c) of Figure 4) are slightly shifted by 2π60 so that the two peaks become closer to each other. Again, this initial approximate converged to the one peak solution.
The discussion above suggests that the single peak solution has a large "convergence area." As an extreme example, random initial approximations converge to single peak solutions in our tests.
In contrast, the q-axial symmetric solutions can be very sensitive to the initial approximations.
In the third experiment set, we focused on the locations of the peaks. The following results suggested that peaks occur on the boundary of the domain, except a trivial solution that has a peak at the origin.
•
Initial solutions contain peaks on the boundary.All the initial solutions used in the first and second experiment set contained peaks lying on the boundary. All of them converged to the solutions with peaks on the boundary.
FIG. 6. The two peaks move from the interior of the domain to the boundary.
FIG. 7. (a) The solution with "infinity" peaks. (b) The solution with a peak in the origin.
•
Initial solutions contain peaks in the interior of the domain.Part (a) of Figure 5 shows an initial solution containing one peak located at (12, 0). The peak gradually moved to the boundary as iteration proceeding. Part (b) of the figure shows the solution shape after 690 iterations. Finally, the solution converged to the single peak solution shown in Figure 4. Figure 6 visualizes similar convergence process for the two- peak initialization. Part (a) shows the initial solution with peaks located at (−12, 0) and (12, 0). Part (b) is the result after 490 iterations. The solution converged to the two-peak solution shown in Figure 4.
•
Initial solutions contain peaks lying on the boundary and interior of the domain.We constructed initial solution sets with peaks located at (i) (13, 0) and (1, 0), (ii) (−13, 0) and (1, 0), and (iii) (1, 0), (12cos2π3 ,12sin2π3), and (12cos4π3,12sin4π3 ). Whether the heights of these peaks are equal to each other or not, these initial approximations all led to the single-peak solution shown in Figure 4.
•
Figure 7 illustrates two radial symmetric solutions. In part (a), the solution has ‘‘infinity’’peaks on the boundary. In part (b), the solution has a single peak in the origin.
Figure 8 illustrates the monotone behavior of the three-peak q-axial symmetric solution that was discretized with ν = 64. In Figure 8, the solid, dashed, and dotted lines were produced by the solution values on the grid points with radial coordinate equal to127129,12985, and 12943, respectively.
As shown in Figure 8, the highest points of the three lines occurred when the azimuthal coordinate equal to 0, 120, and 240, all located on the axes of symmetry. The solution values, plotted by the three lines, decreased from degree 0 to 60, and from 0 to 300 for the first peak. Similar behaviors were found for the second and third peak located at degrees 120 and 240, respectively.
We observed the monotone behavior on all of the q-axial symmetric solutions computed in our experiments. Furthermore, it is also clear from the figure that no consistent monotone decreasing phenomenon occur radially. For example, when the azimuthal coordinate equals to 0, the solid line is higher than other two lines; however, the dashed line is the highest when the azimuthal coordinate equals to 60.
FIG. 8. Monotone behavior of the three-peak solution. The solution values are plotted against the azimuthal coordinate of the grid points.
V. CONCLUSION
Our main goals in this article are to develop an efficient numerical algorithm for solving the equation (1.1), to analyze and justify the algorithm, to explore the solution characteristics, to visualize the solutions, and, hopefully, to motive future analytical results on the equation. To achieve these goals, we have presented a discrete iterative algorithm combining a discretization scheme and a fast Fourier solver. We have proved the convergence of the algorithm. By exploring the structure of the matrices and introducing the q-axial symmetry, we have demonstrated that the computed solutions are positive and that q-axial symmetric solutions with the monotonicity behaviors exist. Numerical results not only consisted with the theoretical inferences in [7]; the results had further presented diverse solution shapes and characters that are new to the best of our knowledge. The numerical techniques leading to the solutions have been discussed in detail.
Critical challenges and further developments remain. We study the problem from the discrete viewpoint, parallel results in continuous perspective do not seem trivial. Implementation of the algorithm in parallel architectures is natural and may greatly benefit time consuming applications.
Furthermore, we believe that our numerical schemes can be suitably and efficiently applied to other equations with the same or higher order dimension domains.
The authors are grateful to Wen-Wei Lin and Wei-Cheng Wang for many helpful discussions and suggestions. We also thank the anonymous referees for carefully reading the manuscript and their comments.
References
1. E. F. Keller and L. A. Segel, Initiation of slime mold aggregation viewed as an instability, J Theor Biol, 26 (1970), 399–415.
2. A. Gierer and H. Meinhardt, A theory of biological pattern formation, Kybernetik (Berlin), 12 (1972), 30–39.
3. W.-M. Ni and I. Takagi, Locating the peaks of least-energy solutions to a semilinear Neumann problem, Duke Math, 70 (1993), 247–281.
4. A. Ambrosetti and P. H. Rabinowitz, Dual variational methods in critical point theory and applications, J Funct Anal, 14 (1973), 349–381.
5. C.-S. Lin and W.-M. Ni, On the diffusion coefficient of a semilinear Neumann problem, S. Hildebrandt, D. Kinderlehrer, and M. Miranda, editors, Calculus of Variations and Partial Differential Equations, Lecture Notes in Mathematics 1340, Springer-Verlag, Berlin-Heidelberg-New York-Tokyo, 1988, 160– 174.
6. W.-M. Ni and I. Takagi, On the shape of least-energy solutions to a semilinear Neumann problem, Commun Pure Appl Math, 70 (1991), 819–851.
7. J.-L. Chern and C.-S. Lin, The symmetry of least-energy solutions for semilinear elliptic equations, Preprint, 2000.
8. C.-S. Lin, Locating the peaks of solutions via the maximum principle, Comm Pure Appl Math, 54 (2001), 1065–1095.
9. G. Chen, W.-M. Ni, and J. Zhou, Algorithms and visualization for solutions of nonlinear elliptic equation, part i: Dirichlet problem, Int J Bifurcation Chaos, 7 (2000), 1565–1612.
10. M.-C. Lai, A note on finite difference discretizations for Poisson equation on a disk, Numer Meth Part Differ Eq, 17 (2001), 199–203.
11. B. L. Buzbee, G. H. Golub, and C. W. Nielson, On direct methods for solving Poisson’s equations, SIAM J Numer Anal, 7(4), (1970), 627–656.
12. R. W. Hockney, A fast direct solution of Poisson’s equation using Fourier analysis, J Assoc Comput Mach, 12 (1965), 95–113.
13. H. Eisen, W. Heinrichs, and K. Witsch, Spectral collocation methods and polar coordinate singularities.
J Comput Phys, 96(2) (1991), 241–257.
14. W. Huang and D. M. Sloan, Pole condition for singular problems: the pseudospectral approximation, J Comput Phys, 107(2) (1993), 254–261.
15. P. N. Swarztrauber and R. A. Sweet, The direct solution of the discrete Poisson equation on a disk, SIAM J Numer Anal, 10 (1973), 200.
16. Oscar Buneman, A compact non-iterative Poisson solver, Technical Report Rep. 294, Stanford Univer- sity Institute for Plasma Research, Stanford, CA, 1969.
17. P. J. Davis, Circulant matrices, Wiley, New York, 1979.
18. K. E. Atkinson, An introduction to numerical analysis, Wiley, New York, 2nd Ed., 1989.
