Numerical continuation and the Gelfand problem

li NOlfl'It - IIOILAND

Numerical

Continuation

and

the

Gelfand

Problem

Jeff S. M c G o u g h

Department of Mathematics 601 BB, University of Nevada-Reno Reno, Nevada 89557-0045

ABSTRACT This paper examines a semilinear elliptic equation with nonlinear forcing known as the Gelfand equation. Information about the solution set in terms of the parameter is determined. We are specifically interested in determining the first two "turning points" of the solution continuum. The primary technique is a numerical method known as pseudo-arclength continuation. This is applied to the problem to determine the "bifurcation" diagram for three dimensional domains. © Elsevier Science Inc., 1998

1.

INTRODUCTION

T h e following elliptic b o u n d a r y value p r o b l e m is k n o w n as t h e Gelfand problem:

A u + Ae ~ = 0,

u = 0,

x • 12 c R ~,

z • ~fl.

(1)

It arises in several contexts, one of which is t h e s t e a d y s t a t e e q u a t i o n for a n o n l i n e a r h e a t c o n d u c t i o n p r o b l e m . This p r o b l e m has also been associated w i t h questions in g e o m e t r y a n d r e l a t i v i t y . W e restrict t h e focus here to positive solutions a n d positive )t. In this p a p e r , we will be p r i m a r i l y i n t e r e s t e d in resolving t h e n u m b e r of positive solutions as a function of t h e p a r a m e t e r A. Specifically, we are i n t e r e s t e d in t h e i n t e r v a l of m u l t i p l i c i t y in

APPLIED MA THEMA TICS AND COMPUTA TION 89:225-239 (1998) © Elsevier Science Inc., 1998 655 Avenue of the Americas, New York, NY 10010

0096-3003/98/$19.00 PII S0096-3003(96)00265-2

226

J. S. MCGOUGH

terms of the parameter. This may be represented as a graph where the norm of the solution is plotted against the parameter; known as a bifurcation diagram. Thus we numerically determine the bifurcation diagram for the Gelfand problem. The Gelfand problem has a long history dating back to the last century [1]. In 1853, Liouville [2] found an explicit solution to the problem in the special case where n = 1. Further progress was made examining radially symmetric solutions on the unit ball, 12 = B 1. On B1, radial solutions satisfy:

d2u

n -

dr 2 +

r

1 du d r + ~t e u

O,

u' ( O) = u(1)

0.

In 1914, Bratu [3] found an explicit solution when n --= 2. In both cases, by using the explicit solution, information about existence, uniqueness and multiplicity as a function of the parameter ~ may be gleaned. A thorough study of equation (1) was first published in 1963 by Gelfand [4]. In 1973, Joseph and Lundgren [5], described the solution dependence on the parameter A, for all n where the domain is the unit ball as illustrated in Figure 1: Understanding the non-radial case was still an open question. In the decade that followed, several results clarified aspects of the problem. Crandall and Rabinowitz [6] were able to prove that for n = 2, the qualitative properties above remained true for II, an arbitrary smooth bounded domain. In 1979, Gidas, Ni and Nirenberg [7] demonstrated that for l l = B 1 all solutions are radially symmetric. Thus Joseph and Lundgren's result captured all the solutions when 12 = B r Another domain of interest is the annular region 12 = B , ~ \ B ~ I , n >1 2. One may look for radial solutions as well. Radial solutions satisfy: n-1 u" +

r

u' + ~e u = 0 ,

fir1)

= u(r2) = 0 -

This case distinguishes itself from the ball in t h a t Gidas, Ni and Nirenberg fails to hold, thus leaving the door open on non-radial solutions. Recently, S. S. Lin showed that for the annular case, the plot qualitatively looks the same as for the n = 2 case for the unit ball domain. The m a x i m u m interval of existence may be easily estimated and gives a coarse bound for the first turning point. Let te be the first eigenvalue for

227

Numerical Continuation

c FIG. 1.

[[ uPl verses A.

A ¢ + A¢ -- 0, ¢ ~ 12 where ¢ = 0, ¢ ~ a12, then the Gelfand problem has no solution for A ~> re. Again the radial case is amenable to computation, for 12= Bl:n=

1 : / ~ = ~r2/4, n = 2 : / z =

5.8, n - 3 : / z =

9.9.

Schaaf [9] demonstrated that for IIuH > M > 0 the solution continuum is bounded in a "corridor" 0 < A. ~< A < A* < oo when n >~ 3 and for convex domains (generalized to starlike by McGough [10]). In other words for large norm, the "curve" is bounded away from A = 0. The main goal of this paper is to determine the interval in A for which multiple solutions exist for the Gelfand problem. This corresponds with estimating the interval in which the upper portion of the continuum is contained. The theorems of Schaaf and McGough only assert the existence of the interval, they do not give tight bounds (some estimates were made). To get precise bounds on the interval, other methods are necessary.

228

J. S. MCGOUGH

For the remainder of this paper, let n = 3. Getting precise information about the location of the turning points (if they exist) and the shape of the bifurcation diagram is rather difficult in the non-radially symmetric domain case. Numerical techniques must be incorporated at this point. 2.

NUMERICAL FORMULATION

We begin with the investigation of whether or not the results on the ball persist to other domains, say a cube for example. The primary concern in this paper is with coarser geometric and topological properties of the domains. Accurate approximation of the domain (the usual reason for choosing the finite element method) will not be crucial. All that is needed is to distinguish convex, starlike, and annular domains. For this reason, the finite element method is not used and thus avoids the additional problem of tessellating space. Thus, for simplicity, we will use finite differences in the approximation of the problem. The finite dimensional approximation begins with a uniform grid imposed over the domain. Grid points are indexed by a lexicographical ordering and the PDE is discretized using a standard seven-point stencil. Let U = (U1, U2, ... , Un) be the approximate solution values at grid points and then A u is replaced with A h u, 1

Ahu= ~{u(x+

h, y, z) + u(x, y + h, z) + u ( x , y, z + h)

+u(z-

h, y, z) +

y - h, z)

+ u ( x , y, z - h) - 6 (x, y, z)}. Once the identification of the values in Ahu are made with U, a nonlinear system that approximates the Gelfand problem arises: F ( U , A) = A U +

Ah2e v = 0

(2)

where e U represent the vector of entries ( e v', e v~, . . . , eV'~) t and A is a symmetric banded matrix. At A = 0, the solution is u = 0. This provides an easy starting point for the computations. There are several ways to proceed. The most elementary is called parameter continuation and is used to start the computations. By incrementing A and then solving the set of nonlinear equations, part of the (u, A) curve can be obtained. This method fails at turning points in the curve and so will only resolve the "minimal" branch of the curve.

Numerical Continuation

229

This may be overcome by the use of an arclength method, which expresses the bifurcation curve as a function of arclength, s,

F( U( s),A( s) ) : O,

(u, A) -- ( u ( s ) , A ( s ) ) ,

II~(s)ll 2 + A(s) 2 = 1.

The actual implementation is a pseudo-arclength method [11, E. J. Doedel, Numerical Analysis and Control of Bifurcation Problems, Course Notes, December 1988]. First s is discretized into a discrete variable s k. The dependence of U on s~ is indicated by U(k). Pseudo-arclength continuation requires an additional equation to be appended to the system (2):

( U( k) - U( k +(h(k)

-

1 ) ) . g r ( k - 1) h(k-

1)). A(k-

1) - A s k = 0,

U( k) - U( k -

1)

ASk- 1

i(k) -- A(k) -

~( k -

1)

ASk_ 1

This addition completes the system. Appending A to the end of the U array and augmenting an extra column of zeros to A, the system becomes

f l u ( k ) + u~+l(k)e "(k) -~ o

G(U(k)) =

(U(k)- U(k-

1)). gr(k- 1 ) - Ask=0"

(3)

Continuation methods are well described in the literature; for an extensive list of references, see Allgower and Georg [13]. Pseudo-arclength continuation is not "self-starting." It requires a previous point to get the direction vector. Continuation in the variable ~? = maxl UI + a h (where a is a small positive number) is used to get started (the parameter continuation mentioned above). This works away from bifurcation and turning points (which do not occur in a small neighborhood of the origin for the Gelfand problem).

230 3.

J.S. MCGOUGH SOLVING T H E ALGEBRAIC PROBLEM

Newton's method is used to solve the nonlinear equations. In order for each step to find a U so that g(U)= 0, a sequence U i of vectors is constructed (for each step k):

U(k)i+l= U(k)i- [DG(U(k)i)]-IG(U(k)i),

i=0,1,2,...

DG is the Jacobian matrix of G evaluated at U ~. Note that Newton's Method has two stages: (1) DG(U(k)i)w = -G(U(k)~), (the linear system solver, O(n3)); (2) U(k) ~+1 --= U ( k ) i + w (a vector update, O(n)). Small steps in the pseudo-arclength continuation can place the initial guess in the basin of attraction of the Newton method. Thus we are assured of convergence of the arclength method for small steps. Convergence can be very fast; see [14]. The difficulties for the nonlinear equation solver really lie in inverting the resulting Jacobian matrix. The speed of the solver is based on the performance of the linear algebra routines. The linear systems that arise in forming the Jacobian above are large and sparse. Iterative methods are well suited for the large, sparse linear systems arising here. For storage considerations, a modified conjugate gradient method is used to solve the linear systems arising from inverting the Jacobian in Newton's method. Because of the appended equation from the continuation method, the Jacobian is not symmetric. There are a couple of extensions to conjugate gradient CG that are effective for nonsymmetric matrices. GMRES, the generalized minimum residual method [15], is the one used here. GMRES performs well because it is transpose-free and each iteration minimizes the residual. An additional aspect of the solver is that the error of the linear solver is chosen as a function of the Newton iterations. The idea is that starting the Newton iterations, the approximation of the root is coarse, so only a coarse approximation is needed for the linear solver. As the Newton iteration proceeds, the error tolerance for GMRES is decreased. Thus one reduces the wasted time in "oversolving" the linear problem [16]. For any conjugate-gradient algorithm, preconditioning is almost essential. One good choice of preconditioner, M, is the incomplete Cholesky factorization. This factorization preserves the band structure by setting to zero the elements in the factorization that do not correspond to the original bands. This avoids the "fill in" between bands reducing data storage. It does not act as an inverse like the full Cholesky does, but the operation count is smaller. The preconditioning strategy chosen is a generalization of the

231

Numerical Continuation

modified incomplete Cholesky factorization (MIC). MIC only computes the diagonal in the factors. This appears to be sufficient in concentrating the spectrum enough so that G M R E S ' s run time is halved. This preconditioner can more than double the speed of the unpreconditioned code. A comparison of run times between three types of preconditioning can be found in Table 1 and two types in Table 2. The F F T code is a standard Fishpack (Poisson solver) code. It is used for an exact inverse to the discrete difference operator. It is interesting to see that even though it reduces the iteration count, the run time is still greater. The ratio of t i m e ( M I L U ) / t i m e ( F F T ) is .6555 (obtained by averaging over 100 numerical solutions with n = 2 0 3 + 1). With averaging over three numerical solutions but setting n = 253 + 1, the ratio of t i m e ( M I L U ) / time(FFT) is .5583, and for n = 303 + 1, the ratio t i m e ( M I L U ) / t i m e ( F F T ) is .4239. It appears that MILU grows more slowly than FFT. This is to be expected, since MILU is O(n) and F F T is order n~ * ny * n z *(log2(n~) + log2(ny) + 5) where n~, ny, nz are the discretizations in the three directions. All timings were performed on a DECstation 3100 using the unix t i m e command. Finite differences are easy to set up for rectilinear domains. The examples reported here are for domains that are cubes, dumbbells and cubic shells. These were chosen to illustrate regions of theoretical interest like starlike domains, nonstarlike domains, or domains with nontrivial topology.

4.

N U M E R I C A L RESULTS

The following bifurcation diagrams are plots of the sup norm of u against the parameter A (where the sup norm of u is maxkl UE]). The discretization interval h is 1/30. For the solid domains, the cubes, this generates 293 + 1

TABLE 1 PDE SOLVER TIMES, 15 3 -t- 1 EQUATIONS Precond. None FFT MILU

Time 1260.9u (21:09) 983.4u (16:31) 629.7u (10:36)

Nonlinear iterations

Krylov iterations

Function calls

Precond. solves

3

26

30

0

2

5

8

7

3

11

15

14

232

J. S. MCGOUGH TABLE 2 PDE SOLVER TIMES, 203 -4- 1 EQUATIONS

Precond. FFT MILU

Time 3593.2u (1:00:12) 2355.4u (30:29)

Nonlinear iterations

Krylov iterations

Function calls

Precond. solves

2

5

8

7

4

15

20

19

equations. Note that all of the problems with supercritical growth have a sharp decay (in A) above the second turning point. This is a numerical effect due to the discretization interval. Experimentation with the discretization interval showed a strong relation between the height of the decay ( y ) and h ( y = c in l / h ) when f ( u ) = e u. Unfortunately, the current computing power available did not permit pushing this cutoff above a third turning point (assuming it exists). For any finite dimensional approximation, the bifurcation plot will be asymptotic to the vertical axis. Thus this cutoff is to be expected. This is easily seen by dividing equation (2) by the norm of U and defining v = U/II Ull, ellUII v Av = -A-IIuII

"

Since IIvii = 1, II Avll is bounded. Thus ell uII

Ai-~

~< c,

so A ---) 0 as II u I I - - , ~ .

Figures 2 and 3 are the two supercritical growth examples. Comparing the supercritical growth examples to Figure 9, the third turning point is numerical artifact and the second branch turn may have significant discretization error for coarse discretizations. However, the loss in accuracy by jumping from one to three dimensions is not as great as would be expected. Figures 4, 5, 6 are generated from punctured domain problems. Here, is the unit cube where various smaller cubes are deleted. It is not so clear that the computations here are flawed in the limit as Itull --* oo. Lin [8] has shown that for the radial case, there are at least two positive radial solutions for each Z ~ (0, A*).

Numerical Continuation

233

11 10 9 8

7 6 5 4

3 2

J

1 0

•

,

i

I

I

I

i

I

I

1

2

3

4

5

6

7

8

9

10

FIG. 2. ~ = (0, 1)3.

More complicated geometries t h a n the cube or p u n c t u r e d cube above were also studied. T h e code has been a d a p t e d to solve over u n i o n s of cubes. Three was the m a x i m u m n u m b e r of cubes used for the investigation. T h e d o m a i n used in Figure 7 is a rectilinear version of a dumbbell. Specifically, it is the u n i o n of two cubes a n d a parallelepiped. A g a i n the upper cutoff occurs in a similar fashion. T h e q u a l i t a t i v e aspect of all the curves is very similar.

3.5 3 2.5 2 1.5 1

0.5 I

0 o

I

I

I

I

I

2

4

6

8

10

FIG. 3. ](u) = e~2 with II = (0, 1)3.

12

234

J. S. MCGOUGH 12 10 8 6 4

2 i

0

5

0

i

10

i

15

20

25

FIG. 4. f~ = (0, 1)3 \ [1/3, 2/3] 3.

One notices from the diagrams t h a t the qualitative behavior is not a function of geometry, in fact for the principal branch drawn, it appears independent of the topology of the domain. The critical factors are the nonlinearity and the dimension. Of the several dozen diagrams produced for the Gelfand nonlinearity in three dimensions, no qualitative differences were seen from the graphs presented here.

12 10 8 6 4

2 0

2

4

t

i

6

8

t

10

FIG. 5. ,0. = (0, 1) 3 \ [13/30, 16/30] 3.

i

12

14

Numerical Continuation

235

12 10 8 6 4

2 0

l

l

I

I

15

20

25

30

•

0

5

10

35

FIG. 6. n = (0, 1)3 \ [7/30, 22/3013.

5.

ERROR

T h e m e s h w i d t h was v a r i e d from h = 1 / 1 0 , to h = 1 / 4 0 . T h e s h a p e of t h e g r a p h a t t h e first t u r n i n g p o i n t does n o t change m u c h b e t w e e n h = 1 / 2 0 a n d h = 1 / 4 0 . If only t h e first t u r n i n g p o i n t is desired, t h e n s e t t i n g h = 1 / 2 0 will be a d e q u a t e . Resolving t h e second t u r n i n g p o i n t is a m u c h more difficult m a t t e r . N u m e r i c a l e x p e r i m e n t s (see F i g u r e 9) on t h e unit ball

12 10 8 6 4

2 0

I

0

1

2

|

i

I

I

I

I

I

3

4

5

6

7

8

9

FIG. 7. Dumbbell domain.

I0

236

J.S.

MCGOUGH

11 10 9 8

7 6 5 4 3 2 1

0 0

2

F~G. 8.

4

6

8

10

12

14

Subcritical growth, 12 = (0, 1)3, solid cube.

have demonstrated that h < 1 / 5 0 is necessary for accurate resolution of secondary turning points. This is impractical, so h = 1 / 3 0 was chosen as a typical mesh size for plot generation. This will resolve the first excursion with good accuracy, the second turning point reasonably well and can be produced in a single overnight run on a DEC 3100. The resulting system size depends on the domain. For example, the solid cube with h = 1 / 3 0 resulted in 293 + 1 ( = 24,390) equations. The " + 1" arises from the parameter.

14 12 10 8

6 4 2 0 0

I

,

I

I

I

I

0.5

1

1.5

2

2.5

3

FIG. 9.

Gelfand equation on the unit ball.

3.5

237

Numerical Continuation

11 10 9 8

7 6 5 4 3 2 1 0 0

0.5

,

f

I

I

I

1

1.5

2

2.5

3

3.5

FIG. 10. Gelfandequation on the annulus.

The subcritical growth curve, Figure 8, is exactly the plot expected and agrees with the theoretical results [17]. To gain some insight into the error of discretization, Figure 9 is an overlay of 10 bifurcation plots. The test problem is A u + Ae ~ = 0 with zero Dirichlet data prescribed on the boundary of the unit ball. This problem reduces to n-1

u" + - -

r

du dr

+ Ae ~ = O,

u'(O) = u(1) = O.

This is solved using the previously mentioned code. Since the exact behavior is known for this case, it provides a test of the solution methods. This specific example also allows for very fine grids, which is not possible in the three-dimensional case and therefore gives upper bounds for how much accuracy can be expected from the methods and code available. For Figure 9, the discretization numbers for this plot are n = 20, 40, 6 0 , . . . , 200. The radial symmetric simplification works for problems defined on annulae as well. Taking 12 = B 2 \ B 1 produces a one-dimensional problem on the interval [1, 2] with zero Dirichlet data. Figure 10 Was generated with 200 discretization points. All other discretizations give the exact same diagram making 20 discretization points appear to be sufficient. Picking larger intervals, i.e., choosing 1~ = B r \ B 1 for large r, did not resolve anything new. If any value of r exhibits multiple solutions, the code could not reliably determine this. Very large intervals require increasingly large numbers of discretization points.

238 6.

J . S . MCGOUGH CONCLUSION

We set out to numerically determine an interval for nonuniqueness. The pseudo-arclength method has worked well. In fact, the methods described above have been effective in numerical continuation for a variety of semilinear elliptic PDEs. One natural idea to improve performance is to a t t e m p t a predictor-corrector approach to the path following. Using the previous two points, one extends the line for a A s distance. This works as an improved guess. In practice, there is no change in iterations or times. Efforts to improve the speed and efficiency of the continuation code will probably lie in the linear solvers. The stopping criterion for the linear solver also has some impact. This criterion m a y be modified in variety of ways, however, if the stepsize is small, not a great deal of linear solves will occur. There are still some questions that remain. The asymptotic behavior of the curve as JJuJJ --* ~ is not known. This is a more difficult problem due to the breakdown of the finite difference method. It is well understood how to compute bifurcation points, but doing this for large systems presents numerical problems. Iterative (Krylov) methods for the linear solver in the continuation algorithm have proved successful. Some integration of a Krylov based eigenvalue routine m a y also be successful in the bifurcation detection aspect. REFERENCES 1 J. Bebernes and D. Eberly, Mathematical Problems from Combustion Theory, volume 83 of Applied Matl~ Sci., Springer-Verlag, New York, 1989. 2 J. Liouville, Sur l'equation aux d~riv~es partielles (a 2in ) 0 / ( a ua v) ± 2 ~ a2 = 0, J. de Math, 18:71-72, (1853). 3 G. Bratu, Sur les equations integrales non lineares, Bull Soc. Math. de Franca 42:113-142, (1914). 4 I. M. Gelfand, Some problems in the theory of quasilinear equations, Amer. Math, Soc. TransL, 29:295-381, (1963). 5 D. D. Joseph and T. S. Lundgren, Quasilinear Dirichlet prolems driven by positive sources, Arct~ Rat. Mech, App£, 49:241-269, (1973). 6 M.G. Crandall and P. H. Rabinowitz, Some continuation and variation methods for positive solutions of nonlinear elliptic eigenvalue problems, Arck Rat. Mech. AnaL, 58(3):208-218, (1975). 7 B. Gidas, W. Ni, and L. Nirenberg, Symmetry and related properties via the maximum principle, Comm. Matb~ Phys., 68:209-243, (1979). 8 S. S. Lin, Positive radial solutions and non-radiai bifurcation for semilinear elliptic equations in annular domains, J. Dill, Eq., 86:367-391, (1990). 9 R. Schaaf, Uniqueness for semilinear elliptic problems--supercritical growth and domain geometry, preprint, 1992.

Numerical Continuation

239

10 J.S. McGough, On solution continua of supercritical quasilinear elliptic problems, Dif. Int. Eqns, 7(5/6):1453-1471, (1994). 11 R. Glowinski, H. B. Keller, and L. Reinhart, Continuation-conjugate gradient methods for the least square solution of nonlinear boundary value problems, SIAM J. Sci. Stat. CompuL, 6(4):793-832, (October 1985). 12 Deleted in proof. 13 E.L. Allgower and K. Georg, Numerical Continuation Methoda An Introduction, Springer-Verlag, New York, 1990. 14 J. Dennis and R. Schnabel, NumericalMethodsfor Unconstrained Optimization and Nonlinear Equation5 Prentice-Hall, Englewood Cliffs, N.J., 1983. 15 Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems,SIAMJ. Sci. Stat. CompuL, 7:856-869, (1986). 16 P. Brown and Y. Saad, Hybrid Krylov methods for nonlinear systems of equations, Technical Report ucrl-97645, LLNL, November 1987. 17 D. (~. de Figueiredo, P. L. Lions, and R. D. Nussbaum, A priori estimates and existence of positive solutions of semilinear elliptic equations, J. Math. pures et appl., 61:41-63, (1982).