Numerical Evaluation of 2D Ground States

A ground state is defined as the positive radial solution of the multidimensional nonlinear problem −Δu(x) + u(x) − f (u(x)) = 0, x ∈ R; lim |x|→∞ |u(x)| = 0, with the function f being either f (u) = a|u|p−1u or f (u) = a|u|pu+ b|u|2pu. The numerical evaluation of ground states is based on the shooting method applied to an equivalent dynamical system. A combination of fourth order Runge-Kutta method and Hermite extrapolation formula is applied to solving the resulting initial value problem. The efficiency of this procedure is demonstrated in the 1D case, where the maximal difference between the exact and numerical solution is ≈ 10−11 for a discretization step 0.00025. As a major application, we evaluate numerically the critical energy constant. This constant is defined as a functional of the ground state and is used in the study of the 2D Boussinesq equations.


Introduction
The main objective of this paper is the numerical evaluation of a ground state which is defined (see, e.g., [1]) as the positive radial solution of with a nonlinearity f of the types f (U) = a|U| p−1 U, where 1 < p < ∞ for n ≤ 2 and 1 < p < n+2 n−2 for n ≥ 3, or f (U) = a|U| p U + b|U| 2p U with 1 < p < ∞ for n ≤ 2.
It is well known that for some functions f the problem (1) possesses an infinite number of distinct non-trivial solutions (see, e.g., [2] for n = 3).We consider here the problems (1) for which there exists precisely one positive solution.In particular, the uniqueness of the positive solution to (1) with f (U) = a|U| p−1 U is established in [3][4][5].Additionally, for problems (1) with double power nonlinearity f (U) = a|U| p U + b|U| 2p U the uniqueness of the positive solution is proved in [5,6] under some additional assumptions on the coefficients a and b.
In practice, the semi-linear Poison equation (1) arises in many areas of applied mathematics, such as nuclear physics, population genetics, nonlinear optics, etc.If one searches solitary waves or standing waves or stationary states in non-linear dispersive equations of the Klein-Gordon, Schrödinger and a e-mail: natali@math.bas.bgBoussinesq type, then one is led to the problem (1).In a different context, the equation ( 1) appears as a stationary state of the nonlinear heat equation.
In the one-dimensional case the solution of (1) can be found exactly, see [7] for f (U) = a|U| p−1 U and [8] for the double power nonlinearity.But in the multidimensional case no explicit solution to (1) is known.
In the present paper we give an algorithm for the numerical evaluation of ground states.Our method is stable with respect to the round-off errors and has a high order of convergence: for 1D problems the computational order of convergence of the numerical solution to the exact one is O(h 3 ), while in the 2D case it is O(h 4 ).Moreover, in the 1D case, the maximal error between the exact and the numerical solution is ≈ 10 −11 for a discretization step 0.00025.As an application, the critical energy constant, important for the study of the 2D Boussinesq equations, is evaluated numerically as a particular functional of the ground state.
The paper is organized as follows: in Section 2 we recall some properties of ground states and give a stable algorithm for the numerical computation of U. Section 3 gives the numerical evaluation of 1D and 2D ground states for the two types of nonlinearities.A comparison with known results is also given there.Finally, an application of the computed ground states to the evaluation of the critical energy constant to the Boussinesq equation is briefly discussed in Section 4.

Numerical method
In [9], [10] it is proved that any positive solution of (1) must be radially symmetric with respect to some point (which is chosen to be the origin of coordinates).Thus, for n = 2, we shift (1) to the one-dimensional ODE problem for u(r) = U(x), r = x 2 1 + x 2 2 , and u is extended to an even function.This u satisfies Eq. ( 2) can be written as the system of ordinary differential equations (ODEs): The positive solution of ( 2) is a monotonically increasing function of r ∈ (−∞, 0] and it decays exponentially at infinity (see, e.g., [9,10]): in particular, there exists a positive constant The difficulties of the evaluation of such types of solutions (called "minimal" solutions in [11] and "principal solutions" in [12, p. 355]) are described for example in [11] and [13, p. 220].These difficulties arise from the stability properties of the numerical procedure -the round-off errors may deteriorate the asymptotic solution at infinity if one computes the solution in the direction of the solution decay.For example, the standard MATLAB 'ode45' program for the 1D problem (4) (see Example 1 below) with the exact initial condition u(0) = 1.5 gives a false solution tending to −∞ (the effect described by Gautschi for Bessel function recurrence).On the other hand, for initial data u(0) = 1.5 + 0.00001 we get a false positive periodic osculating solution with period ≈ 15.This numerical solution has the behaviour of the exact solution only in the interval (0, 7).
A stable procedure, in particular with respect to the round-off error, evaluates the solution in the direction of the solution growth as it is demonstrated in [11].In the literature such kind of algorithms are referred to as "Miller's algorithm" [13, p. 221].

EPJ Web of Conferences
Therefore, in order to compute the solution of (3) in the interval (−∞, 0], we start at −∞, or practically at −A, where A is a sufficiently large number.In this way we come to the following stable numerical procedure. Let be a small number.We introduce a uniform grid r i = −ih, i = 0, 1, . . ., n − 1, n with step h in the interval (−A, 0), such that r n = −nh = −A.We apply the shooting method for the evaluation of the approximations u h and v h to u and v respectively.
A Set a positive constant K.

Numerical results
In this section we present some numerical results concerning the convergence of the method for two types of nonlinearities: f (u) = u 2 and f (u) = 3u 2 − u 3 .The one-dimensional problems are used to validate the numerical procedure -we compare the numerical solution to the exact one.
For the two-dimensional problems (3) we estimate the difference between solutions obtained on nested meshes.
The exact ground state u satisfies the equality The accuracy of the numerical solution u h is also tested by the closeness of the functional I h (u h ) to 0, where I h (u h ) is the discretization to I(u) by the composite Simpson rule.We solve the numerical algorithm for the parameter values A = 40 and = 10 −12 .
Example 1: Single power nonlinearity f (u) = u 2 We use the presented algorithm both for the one-and two-dimensional problems with the nonlinearity f (u) = u 2 .The second and third column in Table 1 demonstrate the results for the 1D casethe maximal error between the exact solution u and the numerical solution u h , as well as the computational rate of convergence.The exact solution to the 1D problem (4) is 3  2 {cosh( x 2 )} −2 .Plots of the 2D ground states for f (u) = u p with p = 2, 3, 4 are given in Fig. 1.
The next two columns (fourth and fifth) demonstrate similar results for the 2D case.For every h the error is evaluated by the Runge method using the calculated solution u h on the mesh with step h.

Example 3
We compare the ground state "best-fit formula" u chr of [14] with our results for f (u) = u 2 .We observe that the maximal difference between the two solutions is ≈ 0.04.We check also the global characteristic I h (u h ) and found that for the 'best fit' formula u chr the functional is I h (u chr ) = 0.0932298941 while for our numerical solution, I h (u h ) = 1.7763568394 • 10 −15 .We conclude that our numerical solution is globally more accurate than the "best-fit formula" of [14].
The global behaviour of the solution (global existence or blow up) of ( 5) with f (w) = a|w| p−1 w in the case of positive initial energy smaller than d depends on the value of the functional I(w 0 ) of the initial data w 0 (see Theorem 4.5 and Theorem 4.8 in [15]).Following results in [16,17], we evaluate the constant d from d = p−1 2(p+1) a R 2 |u(x)| p+1 dx, where u is the ground state.For example, for f (w) = w 2 we use the already computed solution u h to (2) and we found the following value of d: d ≈ 7.750793162654.
The computed magnitude of d allows us to determine a-priori the behavior of the solution with small initial energy -whether the solution will blow up or not -without solving the hyperbolic problem (5).

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 201

Table 2
is constructed similarly to the Table1, in the case of problem (3) with nonlinearity f (u) = 3u 2 − u 3 .The exact solution to the 1D problem (4) is

Table 1 and
Table 2 show that the numerical rate of convergence for the 1D examples is about O(h 3 ), while for the 2D examples the rate of convergence is about O(h 4 ).Moreover, the discrete functional I h (u h ) is equal to zero with the high accuracy 1.776356839400 • 10 −15 .