Special Issue in honor of John W. Neuberger Electronic Journal of Differential Equations, Special Issue 02 (2023), pp. 67–80. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu or https://ejde.math.unt.edu NUMERIC ESTIMATES OF THE PRINCIPAL EIGENVALUE OF THE p-LAPLACIAN USING INTERVAL ARITHMETIC JIŘÍ BENEDIKT, JAN PŮLPÁN Abstract. We present a numerical algorithm for computing rigorous upper and lower estimates of the principal eigenvalue of the p-Laplacian. To control all possible errors including the rounding errors of the computer arithmetic, we use the interval arithmetic. We implement our algorithm in the Julia programming language using IntervalArithmetic.jl package [12]. 1. Introduction We consider the eigenvalue problem −(|u′|p−2u′)′ = λ|u|p−2u in (0, 1), u(0) = u(1) = 0 (1.1) where p > 1 and λ ∈ R is the spectral parameter. Though our algorithm is general enough to treat more general problems, for instance problems with non- constant coefficients or partial differential equations with the p-Laplacian on various domains, we present our algorithm applied to the problem (1.1) for the sake of clarity. Moreover, the dependence of the principal eigenvalue λ1,p of (1.1) on p is expressed explicitly by λ1,p = (p− 1) ( 2π p sin π p )p (1.2) which allows us to track the errors of our numerical estimates easily. Historical remarks on the practical use of p-Laplace equations can be found in [4] and the references therein. Similarly to [2, 3], we estimate the principal eigenvalue λ1,p of (1.1) from above by using the variational characterization λ1,p = min ∫ 1 0 |v′|pdt∫ 1 0 |v|pdt (1.3) where the minimum is taken over all v ∈ W 1,p 0 (0, 1), v 6≡ 0, and from below by using the following consequence of [1, Theorem 2.1]. 2020 Mathematics Subject Classification. 34B09, 34B15, 65L15, 65G30. Key words and phrases. Quasilinear equation; p-Laplacian; principal eigenvalue; interval arithmetic. ©2023 This work is licensed under a CC BY 4.0 license. Published March 27, 2023. 67 68 J. BENEDIKT, J. PŮLPÁN EJDE/SI/02 Theorem 1.1. Assume v ∈ C1(0, 1), v > 0 in (0, 1) and |v′|p−2v′ ∈ C1(0, 1). Then λ1,p ≥ inf t∈(0,1) − ( |v′|p−2v′ )′ vp−1 . (1.4) Obviously, if we choose an arbitrary (but admissible) test function v and compute the Rayleigh quotient in (1.3), we obtain an upper estimate and from (1.4) we obtain a lower estimate. In both cases the optimal choice of the test function v is the principal eigenfunction ϕ1,p, corresponding to the principal eigenvalue. The more similar the test function v to ϕ1,p is, the better estimate we obtain. But from the practical point of view, we have a very limited choice of test functions if we want to compute useful explicit estimates even in the case when the domain is a ball (cf. [2, 3]), and the computation becomes practically impossible on more complex domains. Hence, our approach is to construct the test function from a numerical ap- proximation of the principal eigenfunction. In one dimension we use the shooting method, in higher dimensions there are numerous effective algorithms for various domains – see, e.g., [5, 6, 7, 8, 9, 10, 11, 15]. But, we have to be careful when insert- ing the numerical solution into (1.3) and (1.4). The numerical solution is typically obtained as an interpolation function through some control points. So first we need to make sure that the interpolation function is an admissible function. Second, the integrals in (1.3) cannot be computed symbolically due to the general power p. Consequently, we integrate numerically. Third, all the operations in computer arithmetic (double precision, IEEE 754) necessarily produce rounding errors. We use the interval arithmetic to control both the error of the numerical integration and the rounding errors to get guaranteed bounds for λ1,p. This article is organized as follows. We start with the easier part, which is the upper estimate, in Section 2. Then we deal with the lower estimate in Section 3. Finally, we present numerical results for various values of p in Section 4. 2. Upper estimate Firstly we compute a numerical approximation of the principal eigenvalue λ1,p and the corresponding eigenfunction ϕ1,p by the shooting method. Because of the uniqueness of the solution of the associated initial value problem, it must be ϕ′1,p(0) 6= 0. Since problem (1.1) is homogeneous, we can normalize ϕ1,p assuming ϕ′1,p(0) = 1. Clearly, λ ∈ R is an eigenvalue of (1.1) if and only if the solution (u1, u2) of the initial value problem for a system of two first-order equations u′1(t) = |u2(t)| 1 p−1 · sign(u2(t)) in (0, 1), u1(0) = 0, u′2(t) = −λ|u1(t)|p−1 · sign(u1(t)) in (0, 1), u2(0) = 1, (2.1) satisfies u1(1) = 0. By the bisection method we find λ such that u1(1) = 0 whereas (λ, u1) becomes the principal eigenpair. For this purpose we need to start the bisection with an initial interval which contains λ1,p and no other eigenvalue. In our simple case we can use (1.2). In general we can use analytical estimates, for example [2] or [3] on a ball in RN . The Julia code follows. """ plaplace solve(λinit, p, n; u20=1.0, dom=(0.0, 1.0)) Numericaly solves p-Laplace equation using shooting method. EJDE-2023/SI/02 ESTIMATES OF THE PRINCIPAL EIGENVALUE 69 Arguments: λinit ... initial interval for λ1 p ... p of p-Laplacian n ... number of solution interpolation points u20 ... initial condition for u2 dom ... domain Return: t, tI ... division points where the solution is interpolated, double precision value and its interval representation U1, U1 I ... numerical approximation of u1 at t, tI U2, U2 I ... numerical approximation of u2 at t, tI Λ1 ... numerical approximation of the first eigenvalue λ1 """ function plaplace solve(λinit, p, n; u20=1.0, dom=(0.0, 1.0)) function sl(du, u, P, t) λ, p = P du[1] = abs(u[2])^(1/(p-1)) * sign(u[2]) du[2] = -λ * abs(u[1])^(p-1) * sign(u[1]) end tl, tr = dom u0 = [0.0; u20;] # initial condition a, b = λinit Λ1 = (a + b)/2 ∆t = (tr-tl)/(n-1) e = 1e-12 # stop condition while (b-a) >= e prob = ODEProblem(sl, u0, dom, (Λ1, p)) sol = solve(prob, saveat=∆t, abstol=1e-8, reltol=1e-8) if sol(tr)[1] == 0 break else probA = ODEProblem(sl, u0, dom, (a, p)) solA = solve(probA, saveat=∆t, abstol=1e-8, reltol=1e-8) probS = ODEProblem(sl, u0, dom, (Λ1, p)) solS = solve(probS, saveat=∆t, abstol=1e-8, reltol=1e-8) if solA(tr)[1] * solS(tr)[1] < 0 b = Λ1 else a = Λ1 end Λ1 = (a+b)/2 end 70 J. BENEDIKT, J. PŮLPÁN EJDE/SI/02 end prob = ODEProblem(sl, u0, dom, (Λ1, p)) sol = solve(prob, saveat=∆t, abstol=1e-8, reltol=1e-8) t = LinRange(0, 1, n-1) tI = [@interval(i) for i in t] U1 = [u[1] for u in sol(t).u] U1 I = [@interval(u[1]) for u in sol(t).u] U2 = [u[2] for u in sol(t).u] U2 I = [@interval(u[2]) for u in sol(t).u] return t, tI, U1, U1 I, U2, U2 I, Λ1 end The function plaplace solve returns vectors U1 and U2 of numerical approxi- mations of the values of u1 = ϕ1,p and u2 at n (a parameter) equidistant division points t of the interval [0, 1]. As the test function v for which we compute the Rayleigh quotient in (1.3) we choose the cubic spline interpolation function through the points U1 with zero second derivatives at the endpoints. These boundary con- ditions seem to be natural since for p = 2 we clearly have ϕ′′1,2(0) = ϕ′′1,2(1) = 0. Since we set the first and the last element of U1 to exact zero and the cubic spline interpolation function is of class C2, we have v ∈ W 1,p 0 (0, 1) and so our choice of the test function in (1.3) is admissible. The coefficients of the cubic spline are standardly obtained from the vector of the second derivatives at the division points which is computed as the solution of a system of linear algebraic equations with a tri-diagonal matrix and a right- hand side constructed from the interpolation values and the division points. As it was already mentioned, these computations cannot be carried out in the double precision arithmetic without rounding errors. Hence, instead of t and U1 we use vectors tI and U1 I of the corresponding intervals which are also returned by the Julia function plaplace solve. This makes all the computations to be carried out in the interval arithmetic and we obtain the coefficients of the cubic spline represented as intervals instead of double precision numbers. The properties of the interval arithmetic guarantee that these intervals contain the exact values of the coefficients of the test function v. The following function cubic natural spline returns the interval coefficients csc V. """ cubic natural spline(t, tI, U, UI, Uld2, Urd2; ns=10) IA interpolation of given points by natural cubic spline. Returns spline coefficients ‘csc V‘ as well as interval values ‘V‘ of the spline function. Arguments: t, tI ... division points U, UI ... values to be interpolated Uld2, Urd2 ... left and right boundary values of second derivative ns ... number of division points for each single piece of spline EJDE-2023/SI/02 ESTIMATES OF THE PRINCIPAL EIGENVALUE 71 """ function cubic natural spline(t, tI, U, UI, Uld2, Urd2; ns=10) # A matrix n = length(UI) dv = [4..4 for i in 1:n-2] ev = [1..1 for i in 1:n-3] A = Array(SymTridiagonal(dv, ev)) A−1 = inv(A) # right-hand side h = 1.0/(n-1) rhs = [] for i in 3:length(UI) append!(rhs, 6/h^2 * (UI[i] - 2 * UI[i-1] + UI[i-2])) end rhs[1] = rhs[1]-Uld2 rhs[end] = rhs[end]-Urd2 # second derivatives vector Ud2 = [] append!(Ud2, @interval(Uld2)) append!(Ud2, A−1*rhs) append!(Ud2, @interval(Urd2)) # spline coefficients csc V = Vector[] for i in 1:length(UI)-1 a = b = c = d = 0 a = (Ud2[i+1]-Ud2[i])/(6*h) b = Ud2[i]/2 c = (UI[i+1]-UI[i])/h - h*(2*Ud2[i]+Ud2[i+1])/6 d = UI[i] append!(csc V, [Interval{Float64}[a,b,c,d]]) end V = Interval{Float64}[] for i in 1:length(UI)-1 x dom = t[i]..t[i+1] x int = mince(x dom, ns) f(x) = csc V[i][4] + (x-t[i])*(csc V[i][3] + (x-t[i])*(csc V[i][2] + csc V[i][1]*(x-t[i]))) append!(V, f.(x int)) end return csc V, V end 72 J. BENEDIKT, J. PŮLPÁN EJDE/SI/02 Now we compute the integrals in (1.3). Though the cubic spline is a piecewise polynomial, we integrate numerically because of the general power p. Since the division at n points may be too coarse, we divide each single subinterval by ns sub- division points. On each of the subsubintervals we compute an interval which must contain the range of all possible values of any cubic function with coefficient from the respective interval in csc V. The vector V of these intervals is also returned by the above function cubic natural spline. To compute the numerator in (1.3) we need a vector of intervals which represent the values of v′. The interval coefficients of the piecewise quadratic function v′ are easily obtained from the coefficients of v. A vector of intervals representing the values of v′ on each subsubinterval is returned by the following function der cubic spline. """ der cubic spline(csc, t, tI, ns) Computes the first derivative ‘V‘ of a given interval cubic spline. Arguments: csc ... cubic spline coefficients interval representation t, tI ... division points ns ... number of division points for each single piece of spline """ function der cubic spline(csc, t, tI, ns) V tmp = Interval[] csc Vder = [ [@interval(3) * c[1], @interval(2) * c[2], c[3]] for c in csc ] for i in 1:length(t)-1 x dom = t[i]..t[i+1] x int = mince(x dom, ns) f(x) = csc Vder[i][3] + (x-t[i])*(csc Vder[i][2] + (x-t[i])*csc Vder[i][1]) append!(V tmp, f.(x int)) end V = Interval[] for i in 1:length(Utmp)-1 append!(V, V tmp[i] ∪ V tmp[i+1]) end append!(V, V tmp[end]) return V end Finally, we compute the p-th power of the intervals representing v and v′, multi- ply by the length of the subsubintervals, sum up over all subsubintervals and divide the numerator by the denominator. Everything in the interval arithmetic, of course. We get an interval which is guaranteed to contain the exact value of the Rayleigh quotient for the test function v. Consequently, the right end-point of the interval is the desired rigorous upper estimate of λ1,p. EJDE-2023/SI/02 ESTIMATES OF THE PRINCIPAL EIGENVALUE 73 """ function upper estimate(V1, V1 der, p) Returns guaranteed upper estimate of the first eigenvalue ‘λ1 up‘. Arguments: p ... p of p-Laplacian V1 ... interval values of v1 V1 der ... interval values of v′1 """ function upper estimate(V1, V1 der, p) f(x) = abs(x)^(p) ni = mince(0..1, length(V1)) numerator = 0..0 for i in 1:length(V1 der) numerator = numerator + f(V1 der[i]) * diam(ni[i]) end denominator = 0..0 for i in 1:length(V1) denominator = denominator + f(V1[i]) * diam(ni[i]) end λ1 up = sup(numerator/denominator) return λ1 up end The parameters n and ns are to be experimented with to get a sharp enough estimate. 3. Lower estimate There are two reasons why such a straightforward approach which we used in Section 2 does not work for the lower estimate. First, both the numerator and the denominator in (1.4) go to 0 when t → 0+ or t → 1-. If we represent them again by a vector of intervals, the first and the last interval would have to contain 0 and probably also negative numbers. Consequently, the ratio in (1.4) would become (−∞,+∞) and the lower estimate would be −∞ which is true but useless. Second, if v ∈ C2(0, 1), then for p < 2 the condition |v′|p−2v′ ∈ C1(0, 1) is not satisfied at a point tm where v′(tm) = 0 and v′′(tm) 6= 0. On the other hand, for p > 2 the condition is satisfied but ( |v′|p−2v′ )′ = 0 at t = tm. Hence, the lower estimate would not be even a positive number. Our strategy to avoid the first problem is that we choose the test function v as a numerical approximation of the restriction to (0, 1) of the principal eigenfunction ϕ̃ on (−r, 1 + r) for a small r > 0 instead of on (0, 1). If we insert ϕ̃ into (1.4), we obtain the principal eigenvalue λ̃ on (−r, 1 + r) which is smaller but close to λ1,p if r > 0 is small. Since we normalize by ϕ̃′(−r) = 1, the value ϕ̃(0) is asymptotically 74 J. BENEDIKT, J. PŮLPÁN EJDE/SI/02 equal to r for r → 0+. Consequently, both the numerator and the denominator in (1.4) now go to a positive constant when t → 0+ or t → 1−. A similar approach is applicable even in higher dimensions thanks to the Hopf Maximum Principle (see [13, Prop. 3.2.1 and 3.2.2, p. 801], [14, Theorem 5, p. 200]). To compute the numerical approximation of ϕ̃ we use the function plaplace solve again. We may even decrease the test function by a small positive constant (the test function must not become negative) which does not change the numerator, decreases the denominator, and improves the lower estimate. To deal with the second problem, we still use the cubic spline interpolation function but instead of ϕ̃ we approximate |ϕ̃′|p−1 ·sign(ϕ̃′), i.e., the function u2 from the associated initial value problem. Let us denote the cubic spline interpolation function through the corresponding numerical values of u2 by v2. The natural boundary conditions for the spline are v′2(0) = v′2(1) = −λ̃ϕ̃p−1(0). Instead of exact λ̃ and ϕ̃ we use the corresponding numerical approximations. Therefore, we need the following code to compute the intervals for the coefficient of the spline using the conditions on the end slopes. """ cubic end slope spline(t, tI, U, UI, Uld1, Urd1; ns=10) IA interpolation of given points by end slope cubic spline. Returns spline coefficients ‘csc V‘ as well as interval values ‘V‘ of the spline function. Arguments: t, tI ... division points U, UI ... values to be interpolated Uld1, Urd1 ... left and right boundary values of first derivative ns ... number of division points for each single piece of spline """ function cubic end slope spline(t, tI, U, UI, Uld1, Urd1; ns=10) # A matrix n = length(UI) dv = [4..4 for i in 1:n-2] ev = [1..1 for i in 1:n-3] A = Array(SymTridiagonal(dv, ev)) A[1,1] = 3.5..3.5 A[end,end] = 3.5..3.5 A−1 = inv(A) # right-hand side h = 1.0/(n-1) rhs = [] for i in 3:length(UI) append!(rhs, 6/h^2 * (UI[i] - 2 * UI[i-1] + UI[i-2])) end rhs[1] = rhs[1] - 3/h * ( (UI[2]-UI[1])/h - Uld1) EJDE-2023/SI/02 ESTIMATES OF THE PRINCIPAL EIGENVALUE 75 rhs[end] = rhs[end] - 3/h * (Urd1 - (UI[end]-UI[end-1])/h) # second derivatives vector sol = A−1*rhs Ud2 = [] σ0 = 3/h * ( (UI[2]-UI[1])/h - Uld1) - sol[1]/2 σ1 = 3/h * (Urd1 - (UI[end]-UI[end-1])/h) - sol[end]/2 append!(Ud2, @interval(σ0)) append!(Ud2, sol) append!(Ud2, @interval(σ1)) # spline coefficients csc V = Vector[] for i in 1:length(UI)-1 a=b=c=d=0 a = (Ud2[i+1]-Ud2[i])/(6*h) b = Ud2[i]/2 c = (UI[i+1] - UI[i])/h - h*(2*Ud2[i]+Ud2[i+1])/6 d = UI[i] append!(csc V, [Interval{Float64}[a,b,c,d]]) end V = Interval{Float64}[] for i in 1:length(UI)-1 x dom = t[i]..t[i+1] x int = mince(x dom,ns) f(x) = csc V[i][4] + (x-t[i])*(csc V[i][3] + (x-t[i])*(csc V[i][2] + csc V[i][1]*(x-t[i]))) append!(V, f.(x int)) end return csc V, V end Since the numerator in (1.4) is −v′2, we compute the interval representation of the numerator by the function der cubic spline. We point out that in the denominator we cannot use a cubic spline interpolation function v1 through the numerical values of ϕ̃ (i.e., u1) since we would fail to guarantee |v′1|p−1·sign(v′1) = v2. Instead, we reconstruct v1 from v2 as a primitive function to |v2| 1 p−1 · sign(v2). The additive constant is chosen as the smallest one such that v1 does not become negative. The following Julia function get v1 returns the interval representation of v1. """ get v1(p, V2, t) Rebuilds interval expression of ‘V1‘ by integrating ‘V2‘. Arguments: 76 J. BENEDIKT, J. PŮLPÁN EJDE/SI/02 p ... p of p-Laplacian V2 ... interval values of v2 t ... division points """ function get v1(p, V2, t) f(x) = abs(x)^(1/(p-1))*sign(x) ni = mince(0..1, length(V2)) V1 tmp = Interval[0..0] for i in 1:length(V2) append!(V1 tmp, V1 tmp[end] + f(V2[i]) * diam(ni[i])) end V1 = Interval[] for i in 1:length(V1 tmp)-1 append!(V1, V1 tmp[i] ∪ V1 tmp[i+1]) end V1 = V1 .- inf(minimum(V1)) return V1 end Finally, we find the smallest left end-point of all intervals representing −v′2/v p−1 1 to get the guaranteed lower estimate. """ lower estimate(V2 der, V1, p) Returns guaranteed lower estimate of first eigenvalue ‘λ1 low‘ and interval values of −u′2/u p−1 1 in ‘Flow‘. Arguments: V2 der ... interval values of v′2 V1 ... interval values of v1 """ function lower estimate(V2 der, V1, p) f(x,y) = -x / y^(p-1) λ1 tmp = f.(V2 der, V1) Flow = Interval[] for i in 1:length(λ1 tmp)-1 append!(Flow, λ1 tmp[i] ∪ λ1 tmp[i+1]) end append!(Flow, λ1 tmp[end]) λ1 low = inf(minimum(Flow)) EJDE-2023/SI/02 ESTIMATES OF THE PRINCIPAL EIGENVALUE 77 return λ1 low, Flow end 4. Results We present results for p ∈ {1.5, 1.6, . . . , 3.0}. The parameters which are subject of experimentation are n, ns and r. Our experiments show that r = 10/(n ∗ ns) is a good choice. The first result is for 20 subintervals of [0, 1] and 10 subsubintervals of each subinterval. n = 21 # number of division points ns = 11 # number of subdivision points of each subinterval λ1 exact(P) = (P-1)*(2*(π/P)/(sin(π/P)))^P λ1 lows = Float64[] λ1 lows err = Float64[] λ1 ups = Float64[] λ1 ups err = Float64[] si = mince(0..1, (n-1)*ns) ps = 1.5:0.1:3.0 for p in ps λ1 = λ1 exact(p) λinit = (3.,1.5*λ1) r = 10/(n*ns) dom = (-r, r+1) ### lower estimate t, tI, U1, U1 I, U2, U2 I, Λ1 = plaplace solve(λinit, p, n, dom=dom); ud1 = -Λ1 * U1[end]^(p-1) csc V2, V2 = cubic end slope spline(t, tI, U2, U2 I, ud1, ud1, ns=ns); V1 = get v1(p, V2, t); V2 der = der cubic spline(csc V2, t, tI, ns); λ1 low, Flow = lower estimate(V2 der, V1, p) append!(λ1 lows, λ1 low) append!(λ1 lows err, λ1 low-λ1) ### upper estimate t, tI, U1, U1 I, U2, U2 I, Λ1 = plaplace solve(λinit, p, n); U1[end] = 0 U1 I[end] = 0..0 csc V1, V1 = cubic natural spline(t, tI, U1, U1 I, 0., 0., ns=ns); V1 der = der cubic spline(csc V1, t, tI, ns); λ1 up = upper estimate(V1, V1 der, p); 78 J. BENEDIKT, J. PŮLPÁN EJDE/SI/02 append!(λ1 ups, λ1 up) append!(λ1 ups err, λ1 up-λ1) end The computation time was 18 seconds (in total) on one core of Mac Mini, 3.2GHz 6-core Intel i7, 32GB DDR4 RAM. The results are shown in Figures 1 and 2. The results for n = 201 and n = 101 are shown in Figures 3 and 4, the computation time was 81 seconds on the same machine. Since the graphs in Figure 3 are too close to each other, we give the numerical values in the following table. The lower estimate is rounded down and the upper estimate is rounded up so the estimates in the table are still guaranteed. Figure 1. Estimates for n = 21 and ns = 11. Figure 2. Errors for n = 21 and ns = 11. EJDE-2023/SI/02 ESTIMATES OF THE PRINCIPAL EIGENVALUE 79 Figure 3. Estimates for n = 201 and ns = 101. Figure 4. Errors for n = 201 and ns = 101. p λ1 low λ1 λ1 up 1.5 5.316213 5.318718 5.320712 1.6 6.073440 6.076626 6.078731 1.7 6.898054 6.902030 6.904292 1.8 7.798633 7.803521 7.805978 1.9 8.783792 8.789731 8.792420 2.0 9.862458 9.869604 9.872564 2.1 11.044050 11.052580 11.055849 2.2 12.338611 12.348721 12.352344 2.3 13.756919 13.768830 13.772852 2.4 15.310114 15.324548 15.329019 2.5 17.010835 17.028449 17.033426 2.6 18.872410 18.894140 18.899683 2.7 20.909702 20.936359 20.942537 2.8 23.138381 23.171079 23.177967 2.9 25.575610 25.615619 25.623302 3.0 28.239929 28.288762 28.297335 80 J. BENEDIKT, J. PŮLPÁN EJDE/SI/02 Acknowledgments. J. Benedikt was supported by the Grant Agency of the Czech Republic, Grant No. 22-18261S. References [1] Allegretto, W.; Huang, Y. X.; A Picone’s identity for the p-Laplacian and applications. Nonlinear Anal. 32 (1998), 819–830. [2] Benedikt, J.; Drábek, P.; Estimates of the principal eigenvalue of the p-Laplacian. J. Math. Anal. Appl. 393 (2012), 311–315. [3] Benedikt, J.; Drábek, P.; Asymptotics for the principal eigenvalue of the p-Laplacian on the ball as p approaches 1. Nonlinear Anal. 93 (2013), 23–29. [4] Benedikt, J.; Girg, P.; Kotrla, L.; Takáč, P.; Origin of the p-Laplacian and A. Missbach. Electron. J. Differential Equations (2018), paper no. 16, 17 pp. [5] Biezuner, R. J.; Brown, J.; Ercole, G.; Martins, E. M.; Computing the first eigenpair of the p-Laplacian via inverse iteration of sublinear supersolutions, J. Sci. Comput. 52 (2012) 180–201. [6] Biezuner, R. J.; Ercole, G.; Martins, E. M.; Computing the first eigenvalue of the p-Laplacian via the inverse power method. J. Funct. Anal. 257 (2009), 243–270. [7] Bognár, G.; Szabó, T.; Solving nonlinear eigenvalue problems by using p-version of FEM, Computers and Mathematics with Applications 43 (2003), 57–68. [8] Bueno, H.; Ercole, G.; Zumpano, A.; Positive solutions for the p-Laplacian and bounds for its first eigenvalue, Adv. Nonlinear Stud. 9 (2009) 313–338. [9] Ercole, G.; Esṕırito Santo, J. C. D.; Martins, E. M.; Computing the first eigenpair of the p-Laplacian in annuli. J. Math. Anal. Appl. 422 (2015), 1277–1307. [10] Horák, J.; Numerical investigation of the smallest eigenvalues of the p-Laplace operator on planar domains. Electron. J. Differential Equations 2011, No. 132, 30 pp. [11] Lefton, L.; Wei, D.; Numerical approximation of the first eigenpair of the p-Laplacian using finite elements and the penalty method, Numer. Funct. Anal. Optim. 18 (1997), 389–399. [12] Sanders, D. P.; et al.; JuliaIntervals/IntervalArithmetic.jl: v0.16.7, https://zenodo.org/record/3727070#.YV8ouy RZvI. [13] Tolksdorf, P.; On the Dirichlet problem for quasilinear equations in domains with conical boundary points. Comm. Partial Differential Equations 7 (1983), 773–817. [14] Vázquez, J. L.; A strong maximum principle for some quasilinear elliptic equations. Appl. Math. Optim. 12(3) (1984), 191–202. [15] Yao, X.; Zhou, J.; Numerical methods for computing nonlinear eigenpairs. I. Iso- homogeneous cases. SIAM J. Sci. Comput. 29 (2007), 1355–1374. Jiř́ı Benedikt Department of Mathematics and NTIS, Faculty of Applied Sciences, University of West Bohemia, Univerzitńı 8, CZ–301 00 Plzeň, Czech Republic Email address: benedikt@kma.zcu.cz Jan Půlpán Department of Mathematics and NTIS, Faculty of Applied Sciences, University of West Bohemia, Univerzitńı 8, CZ–301 00 Plzeň, Czech Republic Email address: honza@pulpan.net 1. Introduction 2. Upper estimate 3. Lower estimate 4. Results Acknowledgments References