On the localization and numerical computation of positive radial solutions for ϕ -Laplace equations in the annulus

. The paper deals with the existence and localization of positive radial solutions for stationary partial differential equations involving a general ϕ -Laplace operator in the annulus. Three sets of boundary conditions are considered: Dirichlet–Neumann, Neumann–Dirichlet and Dirichlet–Dirichlet. The results are based on the homotopy version of Krasnosel’ski˘ı’s fixed point theorem and Harnack type inequalities, first es-tablished for each one of the boundary conditions. As a consequence, the problem of multiple solutions is solved in a natural way. Numerical experiments confirming the theory, one for each of the three sets of boundary conditions, are performed by using the MATLAB object-oriented package Chebfun.


Introduction
In the short but clever paper [22], Hayan Wang solved the problem of existence of positive radial solutions for the semilinear elliptic equation with one of the following sets of boundary conditions, w = 0 on |x| = R 1 and |x| = R 2 , (1.1) w = 0 on |x| = R 1 and ∂w/∂r = 0 on |x| = R 2 , (1.2) ∂w/∂r = 0 on |x| = R 1 and w = 0 on |x| = R 2 , (1.3) where r = |x| and ∂w/∂r denotes differentiation in the radial direction, and 0 < R 1 < R 2 < +∞. The tools were a Krasnosel'skiȋ type fixed point theorem in cones and the property of bilateral boundedness of the corresponding Green functions. The first one is based on the fixed point index theory, while the second, as shown in [16], on Harnack type inequalities. Since then, many authors have considered the problem of radial solutions for equations and systems involving the Laplacian or some of its generalizations, various boundary conditions and domains, by using different topological or variational methods. We refer the interested reader to some of these contributions [1-3, 5, 6, 9] and the references therein.
Most of the works that followed deviated from the spirit of the original ideas. On their line, we mention our recent papers [17], [18] and [19]. It is the scope of the present paper to complement them, as close as possible to paper [22], for the case of equations with a general ϕ-Laplacian. Here in the absence of a Green function we are forced to produce Harnack type inequalities for each set of boundary conditions. More exactly, in this paper, we deal with the existence, localization and multiplicity of positive radial solutions to equations involving ϕ-Laplacian operators: where 0 < R 1 < R 2 < +∞, the functions g : [R 1 , R 2 ] → R + , f : R + → R + are continuous and ψ : (−a, a) → R is such that ϕ (s) := s ψ (s) is an increasing homeomorphism between two intervals (−a, a) and (−b, b) (0 < a, b ≤ +∞).
The following particular cases are of much interest due to their corresponding models arising from physics: Looking for radial solutions of (1.4), that is, functions of the form w(x) = v(r) with r = |x|, (1.4) reduces to the ordinary differential equation

Harnack type inequalities
Originally, Harnack's inequality was introduced in order to give estimates from above and from below for the ratio u(x)/u(y) of two values of a positive harmonic function. It can be generically put under the form where k is a positive constant depending on the subdomain ω. Next it was generalized to nonnegative solutions or supersolutions of a wide class of linear elliptic equations. For the origin of the notion and many references, we refer the reader to [14]. More general, we speak about a Harnack type inequality whenever for a given operator L acting on a space of functions defined on a set Ω and endowed with a norm ∥·∥ , there is a subdomain ω ⊂ Ω and a constant k > 0 such that min x∈ω u (x) ≥ k ∥u∥ for all nonnegative functions u satisfying L (u) ≥ 0 and eventually some additional behavior properties. In [15], Harnack inequalities have been put in connection with the compressionexpansion method of Krasnosel'skiȋ for the localization of positive solutions of nonlinear problems. In case of boundary value problems for ordinary differential equations, when a Green function is known, a Harnack inequality immediately can be derived using the bilateral estimates of the Green function. However, Harnack inequalities can be obtained even for differential operators for which a Green function does not exist. This is the case of the ϕ-Laplace operators. Deduction of such inequalities requires a fine analysis and makes use of priority properties of solutions such as monotony and concavity (see, e.g., [10] and [11]). The analysis is even more difficult in the case of radial solutions. It is the goal of this section to obtain Harnack inequalities for ϕ-Laplace operators subject to each of the three boundary conditions (1.6), (1.7), (1.8).

Case of the boundary conditions (1.7)
, then v is nonnegative, increasing and concave. In addition, for any c ∈ (R 1 , Proof. Let h := L (v) . Integrating from r to R 2 and taking into account that v ′ ( A new integration, this time from R 1 to r gives the expression of the solution, namely and of the associated solution operator Since h ≥ 0, these formulas show that v is nonnegative and increasing. Also v ′ is decreasing, i.e., v is concave. Finally, the concavity implies that the graph of v is over the line joining the point (R 1 , 0) and (R 2 , v (R 2 )) , whose equation is Note that under the assumptions of Theorem 2.1, in (2.1), one has v (c) = min r∈[c,

Case of the boundary conditions (1.8)
, then v is nonnegative, decreasing and concave. In addition, for any c ∈ (R 1 , Proof. If we let h = L (v) , then by integration we obtain Since h is nonnegative, these formulas immediately imply that v is nonnegative and decreasing.
To show that v is concave we need to prove that v ′ is decreasing, equivalently, that the function η (r) = r 1−N r R 1 τ N−1 h (τ) dτ is increasing. Indeed, using the monotonicity of h, one has Finally, since the graph of the concave function v is over the line joining the points (R 1 , v (R 1 )) and (R 2 , 0) , if c is any point in (R 1 , R 2 ) , we have (2.2).

Proof.
Let v be a nonnegative solution. Since h is not identically zero, v is nonzero and since it vanishes at R 1 and R 2 , any maximum point R is interior and so v ′ (R) = 0. Integrating from R to r then gives . First assume that R ≤ R m . The concavity of v implies that the graph of v restricted to [R, R 2 ] is over the line joining the points (R, v (R)) and (R 2 , 0) which at its turn is over the line joining the points (R 1 , v (R)) , (R 2 , 0) , of equation In addition the graph of v on [R 1 , R m ] is over the line joining the points (R 1 , 0) , (R 2 , v (R)) . Then As a result The proof of the case R ≥ R m is similar. Next, integration in (2.4) gives the representation formulas (for the solution operator) To prove the existence of a solution, in virtue of (2.5) and (2.4) it is enough to prove the existence of a number R ∈ (R 1 , R 2 ) such that This immediately follows since the continuous function takes values of opposite sign at the ends R 1 and R 2 .
To prove the uniqueness of the solution, assume that v 1 and v 2 are two nonnegative solutions and let R ′ , R ′′ be two of their maximum points, respectively. Using the representation formula (2.4) it is easy to see that (v 2 − v 1 ) ′ preserves its sign on the whole interval (R 1 , R 2 ) , positive or negative depending on the ordering between R ′ and R ′′ . Thus v 2 − v 1 is monotone and being zero at the ends of the interval it must be identically zero. Hence v 1 = v 2 .

Existence and localization
As mentioned above, the key ingredient together with Harnack inequalities to obtain positive solutions in this paper will be the fixed point index in cones. In particular, we recall the well-known homotopy version of Krasnosel'skiȋ fixed point theorem in cones.

Theorem 3.1 (Krasnosel'skiȋ)
. Let X be a Banach space, K a cone of X and Ω 1 and Ω 2 two relatively open and bounded subsets of K with 0 ∈ Ω 1 ⊂ Ω 1 ⊂ Ω 2 . Let N : K → K be a completely continuous operator satisfying one of the following two conditions: (i) λu ̸ = N(u) for all u ∈ ∂ K Ω 1 and all λ ≥ 1; and there exists h ∈ K \ {0} such that u ̸ = N(u) + λ h for all u ∈ ∂ K Ω 2 and all λ ≥ 0.
Then N has a fixed point u ∈ K such that u ∈ Ω 2 \ Ω 1 .
In the sequel, consider the Banach space of continuous functions C[R 1 , R 2 ] endowed with the usual maximum norm ∥v∥ = max r∈[R 1 ,R 2 ] |v (r)| and denote by P the cone of nonnegative

Case of the boundary conditions (1.7)
By a solution of (1.5)-(1.7) we mean a function and satisfies (1.5). We will look for nonnegative nontrivial solutions on [R 1 , It is clear that v is a nonnegative solution of (1.5)-(1.7) if and only if v is a fixed point of the operator for all s ∈ R + , then T 1 is defined on the whole cone P and T 1 (P) ⊂ P. Moreover, T 1 is completely continuous as follows from the Arzelà-Ascoli theorem.
Here, for a fixed c ∈ (R 1 , R 2 ), we will look for fixed points of the operator T 1 in a subcone of P, namely, . By the Harnack inequality given by Theorem 2.1, it is easy to check that the operator T 1 maps the cone K 1 into itself. Now, for any numbers α, β > 0, consider the open (in K 1 ) sets and We are in the position to apply Theorem 3.1 in order to obtain existence and localization results for problem (1.5)-(1.7). In this way, we localize a solution in the set We will use the following notations: Also, for any α, β > 0, we denote In addition assume that there exist α, β > 0 such that has a positive solution v such that α < ∥v∥ < β/k 1 .
Proof. We shall apply Theorem 3.1. First, let us see that On the other hand, let us prove that and thus Theorem 3.1 implies that the operator T 1 has at least a fixed point located in Note that condition (3.3) trivially holds if b = +∞. Obviously, if ϕ is a classical or a bounded homeomorphism, i.e., if a = +∞, then conditions (3.4) and (3.5) can be rewritten as with suitable positive constants C 1 , C 2 , C 3 and C 4 as come from (3.4)-(3.5).
Hence, if we are only interested on the existence and not on the localization of the solutions, we can establish sufficient conditions for the existence of the numbers α and β satisfying the inequalities above. They are given by asymptotic conditions on the ratio f /ϕ at 0 and at infinity. Theorem 3.3. Assume that the following conditions are satisfied: a = +∞, and Then problem (1.5)-(1.7) has at least one positive solution.
Similarly, an existence result can be obtained if f is sublinear at 0 and superlinear at infinity with respect to ϕ. Theorem 3.4. Assume that ϕ is a classical homeomorphism such that and f satisfies Then problem (1.5)-(1.7) has at least one positive solution.

Remark 3.5.
Note that if ϕ is bounded, then condition f ∞ = +∞ is not possible, since lim x→+∞ ϕ(x) = b and f must be bounded.
Note that if ϕ is singular (i.e., a < +∞, b = +∞), then condition (3.4) is trivially satisfied for α large enough and so the existence of a positive solution for problem (1.5)-(1.7) is ensured provided that there exists a positive number β satisfying (3.5). This holds if f is superlinear at 0 with respect to ϕ, i.e., f 0 = +∞. Thus we have Theorem 3.6. Assume that ϕ is a singular homeomorphism such that < +∞ for all τ > 0 (3.10) and f satisfies f 0 = +∞.
Obviously, the localization of solutions given by Theorem 3.2 allows us to derive multiplicity results provided that there exist several couples of positive numbers (α, β) satisfying assumptions (3.4)-(3.5). Some conclusions in this line are collected in the following (1) Let (α i ) 1≤i≤k , (β i ) 1≤i≤k (k ∈ N) be sets of positive numbers with α i < β i ≤ k 1 α i+1 for each i. If the assumptions of Theorem 3.2 hold for each couple (α i , β i ), then problem (1.5)-(1.7) has k different solutions v i such that α i < ∥v i ∥ < β i /k 1 .
(2) Let (α i ) 1≤i≤k , (β i ) 1≤i≤k (k ∈ N) be sets of positive numbers with α i < β i < k 1 α i+1 for each i. If the assumptions of Theorem 3.2 hold for each couple (α (3) Let (α i ) i∈N , (β i ) i∈N be two sequences of positive numbers with α i < β i ≤ k 1 α i+1 for each i. If the assumptions of Theorem 3.2 hold for each couple (α i , β i ), then problem (1.5)-(1.7) has infinitely many different solutions v i such that α i < ∥v i ∥ < β i /k 1 .
(1) For each i, since α i < β i , Theorem 3.2 ensures that problem (1.5)-(1.7) has a positive solution v i such that α i < ∥v i ∥ < β i /k 1 . Now, it suffices to remark that β i /k 1 ≤ α i+1 implies that ∥v i ∥ < ∥v i+1 ∥, so there exist at least k different such solutions.
(2) For each i, since α i < β i , we can derive from the proof of Theorem 3.2 a better localization result: the solution v i belongs to the set W β i \ V α i , that is, On the other hand, for each j ∈ {1, . . . , k − 1}, since β j < k 1 α j+1 , Theorem 3.2 also implies that problem (1.5)-(1.7) has a positive solution w j located in the set Since α j+1 < β j+1 , one has that w j < w j+1 . Finally, the estimations w j (r) and w j < α j+1 < ∥v m ∥ , for n ∈ {1, . . . , j} and m ∈ {j + 1, . . . , k}, show that w j is also distinct from any v i and so problem (1.5)-(1.7) has at least 2k − 1 different solutions.
The proof of case (3) is analogous and thus we omit it.

Case of the boundary conditions (1.8)
By a solution of (1.5)-(1.8) we mean a function v ∈ C 1 [R 1 , and satisfies (1.5). We will look for nonnegative nontrivial solutions on [R 1 , R 2 ]. It is clear that v is a nonnegative solution of (1.5)-(1.8) if and only if v is a fixed point of the operator T 2 : P → P defined as which is a completely continuous operator.
Let us assume that the functions f and g satisfy the following monotonicity assumptions: For a fixed c ∈ (R 1 , R 2 ), we consider the following subcone of P: Note that the operator T 2 maps the cone K 2 into itself. Indeed, take v ∈ K 2 and let us show that w := T 2 (v) belongs to K 2 . Since f and g are nonnegative, then w is nonnegative and decreasing. Moreover, the monotonicity assumptions on f and g given by (H f ) and (H g ) together with the fact that v is decreasing imply that the function r → r N−1 g(r) f (v(r)) is increasing. Thus, r → L(w)(r) = r N−1 g(r) f (v(r)) For any numbers α, β > 0, define the sets V α and W β as in (3.1) and (3.2), with K 2 instead of K 1 . Then the following existence and localization result for problem (1.5)-(1.8) can be proved as an application of Theorem 3.1, which guarantees the existence of a fixed point of We will use the following notation: A := It is obvious that the following result can be proved in a similar way to Theorem 3.2, so we omit the proof here.  (3.12) (1 0 ) If α < β, then problem (1.5)-(1.8) has a positive solution v such that α < ∥v∥ < β/k 2 .
Remark 3.9. If we take into account that f is decreasing, then conditions (3.11) and (3.12) can be rewritten as (3.14) Note that condition (3.13) is always satisfied for α sufficiently large since the left-hand side in the inequality is independent of α. Furthermore, from the fact that f is continuous with f (0) > 0, it follows that condition (3.14) holds for any β close enough to 0. Remark 3.11. Observe that multiplicity results cannot be derived from Theorem 3.8. Indeed, since A ≥ B and f is decreasing, one has

In view of Theorem 3.8 and Remark 3.9, it is clear that problem (1.5)-(1.8) is always solvable under assumptions (H f ) and (H g ). Thus we have
Therefore, any α satisfying (3.13) must be bigger than any β for which (3.14) holds.

Remark 3.12.
Observe that the results contained in Section 3.2 remain valid for R 1 = 0, i.e., in the ball. Note that problem (1.5)-(1.8) with R 1 = 0 and R 2 = 1, that is, in the unit ball, was considered in [19], but the results are not comparable since there f was assumed to be nondecreasing.

Case of the boundary conditions (1.6)
By a solution of (1.5)-(1.6) we mean a function v ∈ C 1 [R 1 , and satisfies (1.5). We will look for nonnegative nontrivial solutions on [R 1 , R 2 ].
To construct the fixed point operator, we need the following technical result, similar to Lemma 1 in [4]. Denote Moreover, the function Q ϕ : Proof. The existence of R with the desired property follows from the proof of Theorem 2.3. Note that for any h ∈ D b , one has For uniqueness, assume that there exist γ i ∈ R (i = 1, 2) such that Now, by the mean value theorem for integration, there exists s 0 ∈ [R 1 , R 2 ] such that This clearly implies that γ 1 = γ 2 . Finally, for the continuity of Q ϕ , let {h n } n∈N ⊂ D b such that h n → h 0 ∈ D b in C[R 1 , R 2 ]. We may assume that Q ϕ (h n ) → γ 0 . Passing to limit we find that and so γ 0 = Q ϕ (h 0 ), as wished.
In addition, the solution operator is monotone as shows the next lemma. The proof follows similar ideas to those in [11].
Proof. Assume to the contrary that v 1 ̸ ≥ v 2 . Then there exists an interval [t 1 , Hence, by the mean value theorem, there exists R ∈ (t 1 , If f satisfies condition (3.3), then for each v ∈ P, the function h := g f (v) ∈ D b and since h ≥ 0, one has S (h) ≥ S (0) = 0. Hence the operator is well-defined. In addition, thanks to the continuity of Q ϕ and the Arzelà-Ascoli theorem, it is completely continuous.
Notice that v is a nonnegative solution of (1.5)-(1.6) if and only if v is a fixed point of the operator T 3 . Here, for a fixed c ∈ (0, (R 2 − R 1 )/2), we shall look for fixed points of the operator T 3 in a subcone of P, namely, By the Harnack inequality given by Theorem 2.3, it follows that the operator T 3 maps the cone K 3 into itself. Now, for any numbers α, β > 0, consider the relatively open sets V α := {v ∈ K 3 : ∥v∥ < α} and W β := v ∈ K 3 : min r∈I c v (r) < β .
Proof. We shall apply Theorem 3.1. First, let us show that ).
Since f (v(s)) ≤ M α for every s ∈ [R 1 , R 2 ] and S is monotone, we have as wished.
On the other hand, let us prove that v ̸ = T 3 (v) + λ h for all v ∈ ∂ K 3 W β and all λ ≥ 0 with h ≡ 1. Notice that for v ∈ K 3 with min r∈I c v(r) = β, we have that β ≤ v(r) ≤ β/k 3 for all r ∈ I c , and thus m β ≤ f (v(r)) for all r ∈ I c . Hence, f (v(r)) ≥ m β χ I c (r) for all r ∈ [R 1 , R 2 ] (where χ I c denotes the characteristic function of I c ). Then Lemma 3.14 implies that Note that there is R ∈ (R 1 , R 2 ) such that Analogously, if R ≤ R m , then we may prove that Therefore, v ̸ = T 3 (v) + λ for all v ∈ ∂ K 3 W β and all λ ≥ 0. The conclusion follows from Theorem 3.1.

Remark 3.16.
If ϕ is odd then the two conditions (3.16) and (3.17) on β reduce to the unique inequality (3.18) We emphasize that if ϕ is a classical or bounded odd homeomorphism, then conditions (3.15) and (3.18) can be rewritten as for certain positive constants C 1 , C 2 , C 3 and C 4 . Therefore, existence results for sublinear and superlinear nonlinearities can be proven exactly as in Section 3.1.

Theorem 3.17.
Assume that ϕ is odd and that one of the following conditions holds: (i) f 0 = +∞, f ∞ = 0 and ϕ is a classical or bounded homeomorphism satisfying (3.6).

Numerical examples
From numerical point of view we will consider three distinct boundary value problems. In order to solve them we make use of the new and powerful MATLAB package Chebfun which is a product of the numerical analysis group at Oxford University led by Professor Trefethen (see for instance [20] and [21] to quote but a few). The philosophy behind this package is non-standard in numerical analysis and can be summed up in the words of its initiator as "Feel symbolic but run at the speed of numerics".
In short, the method implemented by Chebfun is a Chebyshev type collocation one. Chebfun tries to solve a BVP by using successively to approximate the solution Chebyshev polynomials on grids of size 17, 33, 65 . . . until the spectral convergence is reached. The relative accuracy of each computation carried out by a Chebfun algorithm is usually about 16 digits, and in principle the user need have no knowledge of the underlying algorithms. However, when solving a nonlinear BVP, Chebfun provides useful information on the convergence of the Newtonian method used to solve nonlinear algebraic systems obtained by discretization. In addition, the behavior of the solution coefficients can be visualized (the way in which they decrease to the machine accuracy). We will display these two outputs for each of the three issues considered. In fact, we must emphasize that we have used Chebfun with excellent results in our previous works [7] and [17].
Moreover, in order to observe the behavior of the Chebfun system in solving genuinely nonlinear boundary value problems, we have reported in our recent work [8] the solutions of eight non-linear problems, some of them even singular. In the vast majority of cases, the asymptotic rate of convergence of the Chebyshev collocation implemented by Chebfun is exponential (geometric). Only in the case of the singular problem was this reduced to an algebraic one.

First example: a Dirichlet-Neumann problem
Consider the Dirichlet-Neumann problem for an equation involving a singular homeomor- where g (r) = r + 1 2r 2 + 1 , f (v) = v 2 + 1.  The residual Chebfun satisfies the operator is of order 10 −10 and the boundary conditions are satisfied exactly. From the left panel of Fig. 4.2 it is very clear that Newton method converges with an order of at most 2. From the right panel of the same figure one can observe that a Chebyshev polynomial of order 16, with highly and smoothly decreasing coefficients is the solution of this problem and the asymptotic rate of convergence is exponential.

Second example: a Neumann-Dirichlet problem
We now solve numerically the following problem  The residual Chebfun satisfies the operator is of order 10 −11 and the boundary conditions are satisfied exactly. From the left panel of Fig. 4.4 it is very clear that Newton method converges with an order of at least 2. From the right panel of the same figure one can observe that a Chebyshev polynomial of order 17, with highly decreasing coefficients is the solution of the problem and the asymptotic rate of convergence is again exponential.

Third example: a Dirichlet problem
The last example is giving by the Dirichlet problem where g (r) = 1, The residual Chebfun satisfies the operator is of order 10 −10 and the boundary conditions are  satisfied exactly. From the left panel of Fig. 4.6 it is very clear that Newton method converges with an order of at least 2 and from the right panel of the same figure one can observe that a Chebyshev polynomial of order 24, with highly decreasing coefficients is the solution of the problem and the asymptotic rate of convergence continues to be exponential. We must make an important remark at the end of these three examples. Exponential convergence occurs for solutions represented by Chebyshev polynomials of relatively small order (of the order of a few tens). Moreover and more important, the convergence is so fast that no rounding off plateau appears (see the right panels of the Figures 4.2, 4.4 and 4.6).
From this point of view, the problems in the present work, compared to those in [8], appear to be only slightly nonlinear.
In the papers [12] and [13], the authors solve numerically similar problems. They exclusively use shooting type methods, i.e., they transform a nonlinear boundary value problem into a Cauchy problem and then solve it by finite difference schemes.
Variants of the shooting method have produced remarkable results over time, but we consider that the Chebyshev collocation implemented by Chebfun, through the information it provides, is very reliable. Unfortunately, a direct comparison of our results with the numerical results from the last two cited works is almost impossible.