Approximate nonradial solutions for the Lane-Emden problem in the ball


 In this paper we provide a numerical approximation of bifurcation branches from nodal radial solutions of the Lane Emden Dirichlet problem in the unit ball in ℝ2, as the exponent of the nonlinearity varies. We consider solutions with two or three nodal regions. In the first case our numerical results complement the analytical ones recently obtained in [11]. In the case of solutions with three nodal regions, for which no analytical results are available, our analysis gives numerical evidence of the existence of bifurcation branches. We also compute additional approximations indicating presence of an unexpected branch of solutions with six nodal regions. In all cases the numerical results allow to formulate interesting conjectures.


Introduction
Consider the classical Lane-Emden problem −Δu = |u| p−1 u in B u = 0 on ∂B, (1.1) when B ⊂ R 2 is the open unit ball centered at the origin and p > 1. It is well known that (1.1) admits only one positive solution and infinitely many sign changing solutions as it can be easily shown by exploiting the oddness of the nonlinearity and standard variational methods.
Among the sign changing solutions (also called nodal solutions) there are infinitely many radially symmetric ones which are obtained by applying variational methods in the space H 1 0,rad (B) which is the subspace of H 1 0 (B) given by radial functions. Actually it could be proved that, for any m ∈ N, m > 1 there exists a unique radial solution to (1.1) with m nodal regions and which is positive at the origin ( [13], [14], [15]). By nodal region of a function u we mean a connected component of the set {x ∈ B : u(x) ≠ 0}. In particular there exists only one radial solution with two nodal regions and positive at the origin. We will recall later other properties of this solution which, from now on, will be denoted by u 2,p for every p > 1. Among nonradial solutions there are infinitely many "obvious solutions" which are invariant by the action of the discrete group of symmetries given by the rotations of angle 2π k , k ≥ 2 and which are obtained by the odd reflections of the unique positive solution in the spherical sector of angle 2π 2k = π k .
Note that the radial symmetry and the discrete 2π k -symmetry correspond to all the possible symmetries of the eigenfunctions of the Laplace operator −Δ in the space H 1 0 (B). Until recently there was no evidence that other kind of nonradial solutions of (1.1) could exist. Then, in the paper [8], studying the Morse index m(u 2,p ) of the radial solution u 2,p with two nodal regions it was proved that m(u 2,p ) = 12, for p sufficiently large. This was derived, among other things, by a previous accurate analysis of the asymptotic behavior of u 2,p as p → +∞, done in [12]. The definition of Morse index of a solution will be recalled in Section 2.
On the other side, as p → 1, the solution u 2,p , suitably normalized, converges (as one would expect) to a radial eigenfunction of −Δ in H 1 0 (B) with two nodal regions, whose Morse index is 6 (see [11] for the details). Then the change of Morse index m(u 2,p ), as p varies from 1 to ∞, allows to claim that the solutions u 2,p become degenerate (i.e., the linearized operator at u 2,p admits zero as an eigenvalue) and this, in turn, allows to prove that some branches of nonradial solutions bifurcate from the solutions u 2,p , for some values of the exponent p. All this is proved in the paper [11] and will be summarized in Section 2. Since the nonradial bifurcating solutions are close to the radial ones at the starting of the branches, it can be shown that they have two nodal regions, the nodal line does not touch ∂B but they are not radial. Moreover it can be proved that they have some discrete symmetry ( [11]). These solutions are called quasiradial in [11] and they are different from the symmetric solutions v k obtained by the odd reflections of the unique positive solutions in the spherical sector of angle π k , k ≥ 2, since the nodal line of v k , obviously, intersects the boundary of B. Thus they are "new" sign changing solutions characterised by the fact of having an interior nodal line. Let us mention that nodal solutions of this type were found in more general symmetric domains in [7] by some energy comparison, but not proving the existence of bifurcation branches.
The bifurcation result is obtained in [11] by analytical methods relying on degree arguments and on a result of [6]. This approach gives a rigorous proof of the existence of three bifurcation branches and also of the discrete symmetry of the nodal solutions along the branches.
However, as typical of analytical methods, it fails to give other quantitative and qualitative information such as the one about the values of the exponent p at which bifurcation occurs or about the comparison between the energy levels of the radial solutions and those of the nonradial solutions along the bifurcating branches. Other questions like the shape of the nodal lines, the possible existence of secondary bifurcation branches or the study of bifurcation from other nodal radial solutions have not yet been treated by analytical methods.
In this paper we use numerical evidence to give a description of the bifurcation from nodal radial solutions with 2 or 3 nodal regions.
The numerical procedure is the following. The first step is to establish numerical approximations to the radially symmetric solutions with two and three nodal regions. Then we look for the values p 0 where the linearized operator at the corresponding radial solution has an eigenvalue zero. Once these values are detected we find numerically approximations of the nonradial solutions along the bifurcating branches by using the eigenfunctions of the linearization corresponding to the zero eigenvalues. This allows to start the nonradial branches and then we use a branch following technique to continue the branch for larger values of p. The Newton iteration behind these computations also provided, by "accidental" convergence, an additional unexpected branch which apparently has no path connection with the other branches found.
Once the approximate solution branches have been obtained we investigate some analytical and geometrical properties of the numerical solutions, in particular their energy. In this way we obtain several numerical results which, on one side complement the analytical results of [11] and on the other side help to formulate conjectures about the exact solutions.
In particular we obtain: 1) some values of p at which bifurcation from the approximate radial solution u 2,p with two nodal regions occurs, 2) the energy levels of the approximate radial and nonradial solutions so to compare them, 3) the pictures of the approximate nonradial solutions, their level sets, their nodal lines, 4) some values of p at which bifurcation from the approximate radial solutions with three nodal regions occurs, 5) the nonexistence of secondary bifurcations, at least up to certain values of p, 6) the bifurcation diagram in the cases 1) and 4), 7) an independent new branch of approximations with 6 nodal regions, which near p = 1 have a similar shape as the "obvious solutions" obtained by odd reflection of the unique positive solution on the spherical sector with opening angle π 5 , but for larger p substantially differ from these, preserving however their discrete rotational symmetry. The numerical results 1), 2), 3), and 6) complement and give a nice description of the analytical results obtained in [11], while 4) and 6) allow to conjecture that bifurcation arises also from other nodal radial solutions, in particular from the ones with three nodal regions. Finally 5) indicates that secondary bifurcation should not occur, while 7) leads to the conjecture that some bifurcation from the "obvious solutions" should occur.
Our results will be stated in Section 2 together with a short survey of some previous analytical results which motivated our numerical investigation. In Section 3 we describe the numerical results.

Preliminaries and Statement of the Numerical Results
Let us consider the energy functional: which is related to the variational formulation of (1.1), and the so called Nehari manifold defined by: It is well known that there exists only one positive solution of (1.1) which is obtained by minimizing Ep(u) on N and which has Morse index one. We recall the definition of Morse index for a solution to (1.1).

Definition 1.
Let u be a weak solution of (1.1). The Morse index m(u) is the number of the negative eigenvalues, counted by multiplicity, of the linearized operator L u at u: Concerning sign changing solutions it is proved in [3] that the minimum of the functional E p on the so called nodal Nehari set: is attained (v + and v − are the positive and negative part of v). A functionũ p which attains the minimum of Ep on N ± is therefore called a least-energy nodal solution of (1.1) and it has only two nodal regions and Morse index m(ũ p) = 2, ∀p > 1 (see [3] and also [5]).
While the least-energy positive solution is always radially symmetric by the famous result of Gidas, Ni, and Nirenberg ([10]), the same does not hold for the least-energy nodal solution. Indeed it was proved in [1] and then extended in [9] that the Morse index of the radial solution u 2,p with two nodal regions satisfies: and hence cannot be a least-energy nodal solution. Moreover it is proved in [4] and [16] that a least-energy nodal solutionũ p is even in one Cartesian variable and its nodal line touches ∂B (which gives another reason forũ p not being radial). These properties suggest thatũp should be odd with respect to the other variable and hence coincides with the odd reflection of the unique positive solution on the half-disk {(x, y) ∈ B : x > 0}. This is an interesting question which is still open so far. As recalled in Section 1, problem (1.1) admits a unique nodal radial solution, positive at the origin, having exactly m nodal regions, for every m ∈ N, m ≥ 2. Moreover it has infinitely many sign changing solutions which are invariant under rotations of angle 2π k , k ≥ 2, obtained by the odd reflections of the unique positive solution in the spherical sector of amplitude π k . These solutions are "obvious" nodal solutions to (1.1). Recently in [11], exploiting some results of [8], new types of nodal solutions have been found, by studying the bifurcation from the nodal radial solution u 2,p (having 2 nodal regions) as the exponent p varies.
More precisely the authors introduce the space given by the functions which are even in one Cartesian variable and 2π k -periodic in the angular variable, ∀k ∈ N, k ≥ 2.
Their result is the following. So the solutions bifurcating from the radial solutions u 2,p cannot be one of the "obvious solutions" described before, at least for p close to the bifurcation values p k . Several questions arise from this interesting result. We list some of them. Q1) Is it possible to estimate the values of p at which bifurcation occurs?
Obviously this is very difficult by analytical methods because usually bifurcation theory only allows to prove the existence of the bifurcating branches. This question was one of our first motivations for the numerical investigation. Next, one could consider the least-energy nodal solution in the space H 1 0,k (B) and wonder whether this is radial or not. Q2) Do the solutions in H 1 0,k (B) on the bifurcating branches have lower energy than the radial ones, for the same values of the exponent p? As mentioned before, the radial solutions u 2,p are not least-energy in the whole space H 1 0 (B), but could minimize the energy among the solutions in H 1 0,k (B). In [11] it has been shown that this does not happen for p large.
The bifurcation shown in Theorem 1 is global, but analytically is not proved what shape the solutions have along the branches. Also the behavior of the bifurcation branches as p grows is not clear. Q3) What is the shape of the solutions along the bifurcation branches? Is the nodal line always in the interior of B? Do the solutions exhibit peaks, as p increases? Q4) How does the bifurcation diagram look? Are there secondary bifurcation branches? Finally, one would like to know whether nonradial bifurcation occurs also from other radial solutions with more than two nodal regions. Recently in [2] the Morse index of all radial solutions of (1.1) has been computed for p large, showing that, as in the case of two nodal regions, the Morse index should change as p increases from 1 to ∞. This seems to indicate that some bifurcation could occur in the spaces H 1 0,k (B), for some k. The rigorous analytical proof of such bifurcations presents some difficulty and, as far as we know, has not been carried out. Q5) Is there bifurcation from the radial solution with three nodal regions? What symmetry do the solution have along the branches? In this paper we use numerical evidence to indicate answers to the above questions.
Our results must of course be understood in the sense that the numerical solutions we find are only approximations of possible "true" solutions.
The numerical results can be summarized as follows.  Table 2) (Question Q2), iii) the nodal line of the approximated nonradial solutions on the bifurcation branches is in the interior of B and their negative part exhibits some peaks, (see Figures 4,5,6) (Question Q3), iv) the bifurcation diagram in Figure 7 shows that the bifurcating branches do not intersect and there are no secondary bifurcation branches (Question Q4), v) also from v 3,p there are bifurcating branches approximatively at the values of the exponent q 1 = 1.1, q 2 = 2.0, q 3 = 3.1. The solutions on the branches have symmetry of type 2π j+4 , j = 1, . . . , 3, (see Figure 12) (Question Q5), vi) an additional independent solution branch appears to exist, the solutions on it have symmetry of type 2π 5 , (see Figure 13).
Similar results as ii) -iv) hold for the approximate nonradial solutions bifurcating from v 3,p (see Section 3).

Basis functions
For approximately solving equation (1.1) we use a Newton iteration. In the single Newton steps a Ritz-Galerkin procedure with the basis functions defined below is applied. An appropriate start of the Newton iteration is obtained via path-following techniques explained in the next subsection.
By elliptic regularity every solution of the problem is in C 2 (B). In view of this fact, we choose basis functions which also have this property. As a first step we transform equation (1.1) into polar coordinates, i.e., We are finally interested only in approximations to solutions of (1.1) which are symmetric with respect to reflection at some line through the origin. In order to avoid redundant rotational copies, we therefore use only the basis functions F ji (and not the G ji ).
The continuity of the basis functions on B requires that f ji (0) = 0 for j = 1, . . . , m, i = 1, . . . , n j . Let us formulate a sufficient condition for the continuity of the first derivatives of the basis functions. For this purpose we compute ∇F ji (in polar coordinates), which gives Since f ji ∈ C 2 [0, 1], this expression is continuous (without any further assumptions) in all points of B except for r = 0. To achieve continuity also at r = 0, using the continuous differentiability of the functions f ji and l'Hospital's rule we get if j = 0, . . . , m, i = 1, . . . , n j . For continuity of ∇F ji at r = 0, the right-hand-sides of (3.1) need to be φindependent. For j = 1 this requirement is fulfilled, for j ≠ 1 we have to assume that f ′ ji (0) = 0. By similar calculations the continuity of the second derivatives of F ji at r = 0 gives the additional constraints f ′′ ji (0) = 0 for j = 1 and for j ≥ 3. Thus together with the Dirichlet boundary conditions on ∂B required in (1.1) we have the restrictions We choose the functions f ji as follows. Let m ∈ N, n ∈ N and denote by b 1 (r), . . . , b n+3 (r) the B-splines of order 4 corresponding to the knots 0, 0, 0, 0, 1 n , 2 n , . . . , n−1 n , 1, 1, 1, 1. These are piecewise cubic polynomials which are globally C 2 -functions. Taking the above requirements into account and using the properties of the B-splines the appropriate choice for the basis functions is

Approximate solutions
We are looking for approximate solutions of equation (1.1) in the finite dimensional subspace L spanned by the functions F ji (j = 0, . . . , m, i = 1, . . . , n j ) using Newton's method and transformation into polar coordinates, which requires to solve for v ∈ L, where ω k−1 denotes the previous Newton iterate in L. This amounts to a linear algebraic system. The Newton step is then completed by the update ω k = ω k−1 + v. The choice of the initial approximation ω 0 for Newton's method is crucial. We apply path-following techniques to find such initial approximations first to the desired sign changing radially symmetric solutions and then to the desired non-symmetric ones.
This means that for p = 1 each nontrivial solution v is an eigenfunction and λ is the corresponding eigenvalue of −Δ. This suggests to use appropriate scalar multiples of radially symmetric eigenfunctions of −Δ as initial approximation for Newton's iteration for p close to 1. More precisely, for finding the radially symmetric branch with k nodal regions we can start with the k th radially symmetric eigenfunction of −Δ, simultaneously choosing λ (close to) the associated eigenvalue. Supposing that, by this technique, we found an approximation to a solution on one of the sign changing radially symmetric branches, we can use a path-following method to follow the branch in order to find solutions for larger values of p.
With this method we established numerical approximations to the radially symmetric branches with two and three nodal regions, respectively, see Figures 1 and 2.  2. Non-radial solutions: we look for non-symmetric solution branches bifurcating from a radially symmetric branch. Candidates for such bifurcation points arise at values p 0 where the linearization of the given problem at the corresponding radially symmetric solution has an eigenvalue zero. If we choose a parameter value p 1 close to p 0 , (trying both p 1 < p 0 and p 1 > p 0 ), we can apply Newton's method with an initial approximation of the form ω 0 = αu 1 + βΦ, where u 1 denotes an approximation to the radially symmetric solution at p 1 , Φ denotes an approximation to the eigenfunction corresponding to the zero eigenvalue, and α and β are real coefficients (to be determined appropriately). If our initial guess ω 0 was successful, i.e., if it is close enough to an exact non-symmetric solution, then Newton's method converges to this solution on the non-symmetric branch. Once we have found an approximation to one solution on this new branch, we can find further approximations by the branch following technique.
We applied the above procedure to find first approximations on and then bifurcations from the radially symmetric branches with two and three nodal regions, respectively. In both cases we found three nonsymmetric branches bifurcating from the symmetric ones.
In case of the symmetric branch with two nodal regions the approximate values for the bifurcation points are p 1 ≈ 1.7, p 2 ≈ 3.4, p 3 ≈ 10.2. At the first bifurcation point p 1 the eigenfunction corresponding to the zero eigenvalue has a discrete rotational symmetry with angle 2π 3 (3-symmetry), see Figure 3. Correspondingly the non-symmetric branch that bifurcates at this point also has 3-symmetry, see Figure 4. The situation is analogous at the two remaining bifurcation points we found, i.e., at the points p 2 and p 3 . At these points a non-symmetric branch with 4-and 5-symmetry, respectively, bifurcates, see Figures 5 and 6.  Analogously, in case of the symmetric branch with three nodal regions the approximate values for the bifurcation points are q 1 ≈ 1.1, q 2 ≈ 2.0 and q 3 ≈ 3.1. At q 1 a non-symmetric branch with 5-symmetry bifur-   cates, at q 2 a non-symmetric branch with 6-symmetry and at q 3 a non-symmetric branch with 7-symmetry, see Figures 8,9 and 10. Finally, the approximations on the aforementioned independent branch, which we found by "accidental" convergence of the Newton iteration, have 5-symmetry, see Figures 11 and 17.  By applying analogous techniques we could not find any secondary bifurcations on the new nonsymmetric branches.

Analysis of the approximate solution branches
We investigate some analytical properties of the numerical solutions, in particular the dependence of their central peaks on p and their energy. The conclusions of these investigations help us to formulate corresponding conjectures about the exact solutions.
When we plot the level sets of the numerical solutions of the branches with 2 nodal regions, we can observe that the nodal lines around the origin form closed curves, see Figures 14,15 and 16. Moreover the    function values at these central peaks appear to remain bounded and even to converge decreasingly as p → ∞, see the above figures, Figure 18 and Table 1 with central peak values of the solutions with two nodal regions. Analogous results can be achieved for the case of the branches with 3 nodal regions.
We computed approximately the energy B 1 2 |∇ω| 2 − 1 p + 1 |ω| p+1 dx for some of the numerical solutions ω on each of the branches, with the results contained in Table 2. We can observe that for each fixed p the least energy solutions are not the radially symmetric ones, but instead the 3-symmetric solutions with two nodal regions. The energy of the other branches mostly increases as the symmetry increases. The solutions on the radially symmetric branch with three nodal regions and the corresponding non-symmetric branches, as well as the solutions on the additional independent branch, have all higher energy values than the solutions with two nodal regions.