Numerical verification for asymmetric solutions of the Hénon equation on the unit square

by   Taisei Asai, et al.

The Hénon equation, a generalized form of the Emden equation, admits symmetry-breaking bifurcation for a certain ratio of the transverse velocity to the radial velocity. Therefore, it has asymmetric solutions on a symmetric domain even though the Emden equation has no asymmetric unidirectional solution on such a domain. We numerically prove the existence of asymmetric solutions of the Hénon equation for several parameters representing the ratio of transverse to radial velocity. As a result, we find a set of solutions with three peaks. The bifurcation curves of such solutions are shown for a square domain.



page 9

page 10


Existence and nonexistence results of radial solutions to singular BVPs arising in epitaxial growth theory

The existence and nonexistence of stationary radial solutions to the ell...

Stable blow-up dynamics in the L^2-critical and L^2-supercritical generalized Hartree equation

We study stable blow-up dynamics in the generalized Hartree equation wit...

On the structure of the solutions to the matrix equation G^*JG=J

We study the mathematical structure of the solution set (and its tangent...

Ocean surface radial velocity imaging in the AT-INSAR Velocity Bunching Model. A functional approach

This work is concerned with the estimation of radial velocities of sea s...

A square root velocity framework for curves of bounded variation

The square root velocity transform is a powerful tool for the efficient ...

Binary perceptron: efficient algorithms can find solutions in a rare well-connected cluster

It was recently shown that almost all solutions in the symmetric binary ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The Hénon equation was proposed as a model for mass distribution in spherically symmetric star clusters, which is important in studying the stability of rotating starts [1]. One important aspect of the model is the Dirichlet boundary value problem


where is a bounded domain, is the location of the star, and stands for the stellar density. Particularly, is set at the center of the domain. The parameter ( if and if ) is the polytropic index, determined according to the central density of each stellar type. The parameter

is the ratio of the transverse velocity to the radial velocity. These velocities can be derived by decomposing the space velocity vector into the radial and transverse components.

When , the Hénon equation coincides with the Emden equation . In this case, the transverse velocity vanishes and the orbit becomes purely radial. Gidas, Ni, and Nirenberg proved that the Emden equation has no asymmetric unidirectional solution in a convex domain [2]. However, Breuer, Plum, and McKenna reported some asymmetric solutions obtained with an approximate computation based on the Galerkin method [3], which were called “spurious approximate solutions” caused by discretization errors. This example shows the need to verify approximate computations. By contrast, a theoretical analysis [4] for large (when the orbit tends to be purely circular) found that the Hénon equation admits symmetry-breaking bifurcation, thereby having several asymmetric solutions even on a symmetric domain.

The importance of the Hénon equation has led to active mathematical study on it over the last decade. For example, Amadori [5] analyzed the bifurcation structure of (1) with respect to parameter . Amadori applied an analytical method to the Hénon equation that had worked for the Emden equation. Additionally, several numerical studies have been conducted on the Hénon equation [6, 7, 8, 9]. In particular, we are motivated by the work of Yang, Li, and Zhu [6], who developed an effective computational method to find multiple asymmetric solutions of (1) on the unit square using algorithms based on the bifurcation method. They generated the bifurcation curve of (1) with and numerically predicted bifurcation points around and using approximate computations.

The purpose of our study is to prove the existence of asymmetric solutions of (1) on the same domain, , using the Newton–Kantorovich theorem (see Theorem 2). We prove their existence through the following steps:

  1. We construct approximate solutions using the Galerkin method with polynomial approximations.

  2. Using the Newton–Kantorovich theorem (Theorem 2), we prove the existence of solutions of (1) with nearby approximations while sharply evaluating the error bound between and in terms of the -norm .

Through the steps above, we successfully prove the existence of several solutions for , including those with three peaks, which were not revealed in [6] (see Figure 1). These solutions are proved after the second bifurcation point.

The remainder of this paper is organized as follows. Some notation is introduced in Section 2. Sections 3 and 4 describe numerical verification based on the Newton–Kantorovich theorem together with evaluations of several required constants. Section 5 shows the results numerically proving the existence of several asymmetric solutions of (1). Subsequently, we discuss the bifurcation structure of the problem for .

2 Preliminaries

We begin by introducing some notation. For two Banach spaces and , the set of bounded linear operators from to is denoted by . The norm of is defined by


Let be the function space of -th power Lebesgue integrable functions over a domain with the -norm When , is the Hilbert space with the inner product . Let be the function space of Lebesgue measurable functions over , with the norm . We denote the first-order Sobolev space in as and define

as the solution space for the target equation (1). We endow with the inner product and norm


where is a nonnegative number chosen as


for a numerically computed approximation ; is explicitly constructed in Section 5. Because the norm monotonically increases with respect to , the norm is dominated by the norm for all . Therefore, the error bound in terms of the norm in (4) can be used as that in terms of the norm . The topological dual space of is denoted by with the usual supremum norm defined in (2).

The bound for the embedding is denoted by . More precisely, is a positive number satisfying


Note that , holds for satisfying

. Explicitly estimating the embedding constant

is important for our numerical verification. We use [10, Corollary A.2] to obtain an explicit value of .

Theorem 1 ([10, Corollary A.2]).

Let be a bounded domain, the measure of which is denoted by . Let if , if . We set . Then, (6) holds for

Here, is defined by

where is the gamma function.

When , to which Theorem 1 is inapplicable, the following evaluation is used:


is the first eigenvalue of the Laplacian in the weak sense. For example, when

, we have .

3 Numerical verification method

This section discusses the numerical verification method used in this paper. We first define the operator as

Furthermore, we define the nonlinear operator as and characterize it as

where . The Fréchet derivatives of and at are denoted by and , respectively, and given by

Then, we consider the following problem:


which is the weak form of the problem (1). To conduct the numerical verification for this problem, we apply the Newton–Kantorovich theorem, which enables us to prove the existence of a true solution near a numerically computed “good” approximate solution (see, for example, [11]). Hereafter, and respectively denote the open and closed balls with center approximate solution and radius in terms of norm .

Theorem 2 (Newton–Kantorovich’s theorem).

Let be some approximate solution of . Suppose that there exists some satisfying


Moreover, suppose that there exists some satisfying


where is an open ball depending on the above value for small . If

then there exists a solution of in with

Furthermore, the solution is unique in .

4 Evaluation for and

To apply Theorem 7 to the numerical verification for problem (1), we need to explicitly evaluate and . The left side of (8) is evaluated as

Therefore, we set

Moreover, the left side of (9) is estimated as

Hence, the desired value of is obtained via

where is the Lipschitz constant satisfying


We are left to evaluate the inverse operator norm , the residual norm , and the Lipschitz constant for problem (7).

4.1 Residual norm

Under the condition , we evaluate as follows:


where is the embedding constant satisfying (6) with .

4.2 Inverse operator norm

In this subsection, we evaluate the inverse operator norm . To this end, we use the following theorem.

Theorem 3 ([12]).

Let be the canonical isometric isomorphism; that is, is given by



is positive, then the inverse of exists, and we have


where denotes the point spectrum of .

The eigenvalue problem in is equivalent to

Recall that denotes the inner product defined in (3) that depends on . Because is already known to be in , it suffices to look for eigenvalues . By setting , we further transform this eigenvalue problem into


Because is chosen so that becomes positive (see (5)), (14) is a regular eigenvalue problem, the spectrum of which consists of a sequence of eigenvalues converging to . To compute on the basis of Theorem 3, we need to enclose the eigenvalue of (14) that minimizes the corresponding absolute value of . We consider the approximate eigenvalue problem


where is a finite-dimensional subspace of such as the space spanned by the finite element basis and Fourier basis. For our problem, will be explicitly chosen in Section 5. Note that (15) is a matrix problem with eigenvalues that can be enclosed with verified numerical computation techniques (see, for example, [13, 14, 15]).

We then estimate the error between the -th eigenvalue of (14) and the -th eigenvalue of (15). We consider the weak formulation of the Poisson equation,


given . This equation has a unique solution for each [16]. Let be the orthogonal projection defined by

The following theorem enables us to estimate the error between and .

Theorem 4 ([17, 18]).

Let . Suppose that there exists such that


for any and the corresponding solution of (16). Then,

where the -norm is defined by .

The right inequality is known as the Rayleigh–Ritz bound, which is derived from the min-max principle:

where , and the minimum is taken over all -dimensional subspaces of . The left inequality was proved in [17, 18]. Assuming the -regularity of solutions to (16) (which follows, for example, when is a convex polygonal domain [16, Section 3.3]), Theorem 4 ensures the left inequality. A more general statement that does not require the -regularity is proved in [18, Theorem 2.1].

When the solution of (16) has -regularity, (17) can be replaced with


The constant satisfying (18) is obtained via (see [19, Remark A.4]), where we denote with . For example, when , an explicit value of is obtained for spanned by the Legendre polynomial basis using [20, Theorem 2.3]. This will be used for our computation in Section 5.

Theorem 5 ([20]).

When , the inequality

holds for

4.3 Lipschitz Constant

Hereafter, we denote . The Lipschitz constant satisfying (10), which is required for obtaining , is estimated as follows:


The numerator of (4.3) is evaluated as

Therefore, we have

Choosing from , for small , we can express them as

Hence, we have

For fixed , this is reduced to


Furthermore, when we set with the center , this is further reduced to


where .

Remark 1.

When , the constant decreases as increases. This is a “good” trend for the verification criterion. However, at the same time, a larger raises the solution altitude, leading to larger absolute error bounds (see the numerical results in Section 5).

5 Results

In this section, we present numerical verification proving the existence of asymmetric solutions of (1) with . All computations were implemented on a computer with 2.20 GHz Intel Xeon E7-4830 CPUs 4, 2 TB RAM, and CentOS 7 using MATLAB 2019b with GCC Version 6.3.0. In the following, the existence of all solutions was proved via the Newton–Kantorovich theorem, and all rounding errors were verified with toolboxes kv Library [21] Version 0.4.49 and Intlab Version 11 [14]. Therefore, the accuracy of all results was guaranteed mathematically. We constructed approximate solutions of (1) for from a Legendre polynomial basis [20]. Specifically, we defined a finite-dimensional subspace

as the tensor product

, where each is defined as

For a fixed integer , we constructed as

Tables 1 and 2 show the approximate solutions together with their verification results. Here, was set to the floating point number after zero to satisfy (5). In the tables, , , , , and denote the constants required by Theorem 2. Moreover, and denote an upper bound for absolute error and relative error , respectively. The values in row “Peak” represent upper bounds for the maximum values of the corresponding approximations. We see that error bounds are affected by the number of peaks — fewer peaks tend to lead to larger error bounds. Moreover, as increases, the peaks approach the corners of the domain and become higher. Therefore, a larger makes verification based on Theorem 2 more difficult. However, we succeeded in proving the existence of solutions in all cases in which , including three-peak solutions not found in [6].

Figure 1 displays the solution curves of (1) for ( is always a multiple of 0.05). If the vertical axis scaling is changed, the curves coincide with those in [6, Figure 2] except for that corresponding to the three-peak solutions after the second bifurcation point around . Note that the verified points where lie on the solution curves. In this sense, the reliability of the result is higher than that from just approximate calculations. According to Figure 1, the bifurcation points are expected to be around and . The single-solution curve bifurcates to three at the first bifurcation point around . Then, one of them further bifurcates to three at the second point around .

0 2
40 40 60 60
40 40 40 40
1.17370e-7 3.96407e-7 1.19312e-8 4.22257e-7
1.70326 2.26200 15.19763 36.47472
6.78398e-1 1.64252 1.43209 1.21150
1.99910e-7 8.96672e-7 1.81325e-7 1.54017e-5
1.15549 3.71537 21.76424 44.18887
4.63296e-8 2.55597e-7 1.44557e-7 2.48634e-5
3.76958e-9 3.98528e-9 2.45351e-9 4.63166e-7
Peak 6.62326 24.36528 29.03437 29.20268
Table 1: Verification results for .

: the number of basis functions for constructing approximate solution
: the number of basis functions for calculating
: upper bound for the residual norm estimated by (11)
: upper bound for the inverse operator norm estimated by Theorem 3
: upper bound for Lipschitz Constant satisfying (10)
: upper bound for required in Theorem 2
: upper bound for required in Theorem 2
: upper bound for absolute error
: upper bound for relative error
Peak: upper bound for the maximum values of the corresponding approximation

70 70 70 70 70
80 80 80 80 80
1.88534e-11 7.91070e-6 4.76970e-7 8.47044e-6 3.47384e-8
6.82420 24.18779 78.96665 21.26750 47.44875
2.31308 1.46531 1.55126 1.18832 1.97091
1.28659e-10 1.91343e-4 3.76648e-5 1.80145e-4 1.64830e-6
15.78486 35.44250 1.22498e+2 25.27251 93.51720
4.95952e-11 1.73351e-4 8.76586e-5 1.53306e-4 2.32064e-6
2.35369e-13 9.86681e-7 5.12219e-7 1.20925e-6 1.16657e-8
Peak 62.30489 68.15045 66.28947 69.69524 64.16408
Table 2: Verification results for .

: the number of basis functions for constructing approximate solution
: the number of basis functions for calculating
: upper bound for the residual norm estimated by (11)
: upper bound for the inverse operator norm estimated by Theorem 3
: upper bound for Lipschitz Constant satisfying (10)
: upper bound for required in Theorem 2
: upper bound for required in Theorem 2
: upper bound for absolute error
: upper bound for relative error
Peak: upper bound for the maximum values of the corresponding approximation

Figure 1: Bifurcation curves for (1) on the unit square .

6 Conclusion

We have numerically proved the existence of asymmetric solutions of the Hénon equation (1) for several parameters of using the Newton–Kantorovich theorem (Theorem 2). This ensures the existence of several solutions of (1), including solutions with three peaks not found in [6]. The bifurcation curve of (1) is illustrated for in Figure 1. Future work should verify the existence of solutions for arbitrary real values of , and describe the bifurcation structure for (1) in a strict mathematical sense.

7 Acknowledgments

We express our sincere thanks to Dr. Kouta Sekine (Toyo University, Japan) for his helpful advice. This work was supported by CREST, JST Grant Number JPMJCR14D4. The second author was supported by JSPS KAKENHI Grant Number JP19K14601.


  • [1] M. Hénon, Numerical experiments on the stability of spherical stellar systems, Astronomy and astrophysics 24 (1973) 229–238.
  • [2] B. Gidas, W.-M. Ni, L. Nirenberg, Symmetry and related properties via the maximum principle, Communications in Mathematical Physics 68 (3) (1979) 209–243.
  • [3] B. Breuer, M. Plum, P. McKenna, Inclusions and existence proofs for solutions of a nonlinear boundary value problem by spectral numerical methods, in: Topics in Numerical Analysis, Springer 15, (2001) 61–77.
  • [4] D. Smets, M. Willem, J. Su, Non-radial ground states for the Hénon equation, Communications in Contemporary Mathematics 4 (03) (2002) 467–480.
  • [5] A. L. Amadori, F. Gladiali, et al., Bifurcation and symmetry breaking for the Hénon equation, Advances in Differential Equations 19 (7/8) (2014) 755–782.
  • [6] Z. Yang, Z. Li, H. Zhu, Bifurcation method for solving multiple positive solutions to Henon equation, Science in China Series A: Mathematics 51 (12) (2008) 2330–2342.
  • [7] Z.-x. Li, Z.-h. Yang, H.-l. Zhu, Bifurcation method for computing the multiple positive solutions to p-Henon equation, Applied Mathematics and Computation 220 (2013) 593–601.
  • [8] Z.-X. Li, Z.-H. Yang, H.-L. Zhu, A bifurcation method for solving multiple positive solutions to the boundary value problem of the Henon equation on a unit disk, Computers & Mathematics with Applications 62 (10) (2011) 3775–3784.
  • [9] Z.-x. Li, H.-l. Zhu, Z.-h. Yang, Bifurcation method for solving multiple positive solutions to Henon equation on the unit cube, Communications in Nonlinear Science and Numerical Simulation 16 (9) (2011) 3673–3683.
  • [10] K. Tanaka, K. Sekine, M. Mizuguchi, S. Oishi, Sharp numerical inclusion of the best constant for embedding