Special Issue in honor of Alan C. Lazer Electronic Journal of Differential Equations, Special Issue 01 (2021), pp. 301–314. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu or https://ejde.math.unt.edu INFINITELY MANY SOLUTIONS AND ASYMPTOTICS FOR RESONANT OSCILLATORY PROBLEMS PHILIP KORMAN, DIETER S. SCHMIDT Dedicated to the memory of Alan C. Lazer, a great mathematician and a dear friend Abstract. For a class of oscillatory resonant problems, involving Dirichlet problems for semilinear PDE’s on balls and rectangles in Rn, we show the existence of infinitely many solutions, and study the global solution set. The first harmonic of the right hand side is not required to be zero, or small. We also derive asymptotic formulas in terms of the first harmonic of solutions, and illustrate their accuracy by numerical computations. The numerical method is explained in detail. 1. Introduction We study multiplicity of solutions for semilinear equations with linear part at resonance, a direction of research initiated by the classical paper of Landesman and Lazer [13]. Several classes of oscillatory resonant problems on balls and rectangles in Rn are considered. Our focus is on the existence of infinitely many solutions, which we establish by studying global solution curves. Our results are supported by asymptotic analysis and numerical computations, and they extend related research in [1, 2, 15, 16]. Next we describe one of our main results, and the approach used. Let B be the unit ball in R2, x2 + y2 < 1. For a class of oscillatory resonant problems, with h(u) = √ u sin[ln(u3/2 + 1)] and r = √ x2 + y2, ∆u+ λ1u+ h(u) = g(x, y) = µ1ϕ1(r) + e(x, y) for x ∈ B , u = 0 on ∂B, (1.1) we show the existence of infinitely many solutions for any g(x, y) ∈ L2(B)∩Cα(B), α > 0. Here ∆u = uxx(x, y) + uyy(x, y), while (λ1, ϕ1(r)) is the principal eigenpair of the Laplacian on B, with zero boundary conditions, µ1 ∈ R, e(x, y) ∈ ϕ⊥1 in L2(B), and e(x, y) ∈ L2(B) ∩ Cα(B), for some α ∈ (0, 1). The function e(x, y) is not assumed to be radially symmetric. Decompose solutions of (1.1) as u(x, y) = ξ1ϕ1(r) + U(x, y), with ξ1 ∈ R and U(x, y) ∈ ϕ⊥1 in L2(B) (here ϕ⊥1 denotes the orthogonal complement of ϕ1 in L2(B); similar decomposition is used throughout the paper). We prove that the solution set of (1.1) is exhausted by a continuous 2020 Mathematics Subject Classification. 35J25, 35J61, 65N25. Key words and phrases. Global solution curves; asymptotic distribution of solutions. ©2021 This work is licensed under a CC BY 4.0 license. Published February 23, 2022. 301 302 P. KORMAN, D. SCHMIDT EJDE/SI/01 solution curve (u(x, y), µ1)(ξ1) parameterized by ξ1 ∈ R, and a section of this curve µ1 = µ1(ξ1) oscillates toward ±∞ as ξ1 → ∞, see Figure 1 below. We find this result to be rather surprising. A more typical behavior is that µ1(ξ1) → 0 as ξ1 → ∞. The choice of nonlinearity h(u) was dictated by rather restrictive conditions that we needed to impose to obtain a continuous curve of solutions. For more natural nonlinearities we can only assert the existence of a continuum of solutions (as in the case of (1.3) below), and in some cases we have to rely on the numerically produced solutions curves (e.g., for (1.2) below), for which we can still derive accurate asymptotic formulas. In dimension n = 1 a similar result was proved in [11] for h(u) = up sinu, 1/2 < p < 1. We remark that solution curves of the form (u(x, y), µ1)(ξ1) appeared previously in Schaaf and Schmitt [16]. For a model resonant problem ∆u+ λ1u+ u sinu = µ1ϕ1(r) + e(x, y) for x ∈ B , u = 0 on ∂B, (1.2) we provide asymptotical and computational evidence, as well as heuristic justifi- cation of the following conjecture: there exist two numbers 0 < a < A so that the problem (1.2) has infinitely many solutions for µ1 ∈ (−a, a), there are at most finitely many solutions for µ1 outside of (−a, a), and no solutions exist for |µ1| > A. A similar situation occurs for rectangles in two dimensions. The existence of in- finitely many solutions for an interval of µ1’s (bounded or unbounded) is a new, and actually a rare phenomenon (as evidenced by the results of this paper, in- cluding numerical computations). For the problem (1.2) on balls and rectangles in dimensions higher than two, and for “most” other nonlinear terms, infinitely many solutions occur only at µ1 = 0. The restriction to the case µ1 = 0 is common in the literature, see e.g., [1, 3, 4, 15, 16]. The nonlinear term in (1.2) occurred previously in [7], in connection with the oscillatory bifurcation from infinity. Problem (1.2) can be seen as a limiting case of another model problem ∆u+ λ1u+ up sinu = µ1ϕ1(r) + e(x, y) for x ∈ B , u = 0 on ∂B , (1.3) with p ∈ (0, 1), to which the well known results of Costa et al [1] apply at µ1 = 0 (the results of [1] hold for more general domains), see also [15]. It follows from [1] that for µ1 = 0 the problem (1.3) has infinitely many solutions, and moreover u/maxB u → ϕ1(r) for large solutions. We derive a rather precise asymptotic formula for µ1 = µ1(ξ1) in case |ξ1| is large, and this formula tends to be accurate for small |ξ1| as well. The nonlinear term u sinu in (1.2) has linear growth at infinity, while (1.3) requires sublinear growth and is related to the result in [1]. In addition to balls in R2, we obtain asymptotic formulas and perform compu- tations of solution curves for rectangular domains, and for radial solutions on balls in any dimension. 2. Global solution set for a ball in R2 Let J0(z) be the Bessel function of order zero, with J0(0) = 1, and denote by ν1 > 0 its first root, ν1 ≈ 2.405. The principal eigenpair of the Laplacian on the unit ball B ∈ R2 is λ1 = ν2 1 ≈ 5.78, ϕ1(r) = c0J0(ν1r) with r = √ x2 + y2, with c0 chosen so that ∫ B ϕ2 1(r) dx dy = 2πc20 ∫ 1 0 J2 0 (ν1r)r dr = 1 , EJDE-2021/SI/01 INFINITELY MANY SOLUTIONS 303 which is c0 = 1√ 2π ∫ 1 0 J2 0 (ν1r)r dr ≈ 1.09 . (2.1) Observe that ϕ1(0) = c0. We shall also need the second eigenvalue λ2. Recall (see e.g. [14]) that the eigenvalues of the Laplacian on B with zero boundary condition are λn,m = α2 n,m (n = 0, 1, 2, . . .; m = 1, 2, . . .) with the corresponding eigenfunctions Jn(αn,mr)(α cosnθ+β sinnθ), where αn,m is the m-th root of Jn(x), the n-th Bessel function (α and β are arbitrary constants). One calculates λ2 = α2 1,1 ≈ 14.62, with α1,1 ≈ 3.83, and ϕ2 = J1(α1,1r)(α cos θ+ β sin θ). The principal eigenvalue is simple, while all other eigenvalues have multiplicity two, because any two Bessel functions with indices differing by an integer do not have any roots in common, see Watson [17]. Let us now recall the following result from [9] and [11]. It deals with PDE’s on a general domain Ω ⊂ Rn, ∆u+ h(u) = µ1ϕ1(r) + e(x) for x ∈ Ω, u = 0 on ∂Ω . (2.2) Here x ∈ Rn, r = |x| and (λ1, ϕ1(r)) is the principal eigenpair of the Laplacian on Ω, with zero boundary conditions, µ1 ∈ R, e(x) ∈ ϕ⊥1 in L2(Ω), and e(x) ∈ Cα(Ω) ∩ L2(Ω), for some α ∈ (0, 1). Theorem 2.1. For problem (2.2) assume that h(u) ∈ C2(R), f(x) ∈ L2(Ω), and h′(u) < λ2 − λ1, for all u ∈ R , (2.3) |h(u)| < γ|u|+ c, with 0 < γ < λ2 − λ1, c ≥ 0, and u ∈ R . (2.4) Then the solution set of (2.2) consists of a continuous curve (u(x), µ1)(ξ1) param- eterized by ξ1 ∈ R. The following lemma extends a similar result in [11]. Recall that solutions of (2.2) are decomposed as u(x) = ξ1ϕ1(r) + U(x). Lemma 2.2. In the conditions of Theorem 2.1 assume that lim|u|→∞ h(uz)/u = 0, uniformly in z ∈ R. Then as ξ1 → ±∞, the solutions of (2.2) satisfy u(x)/ξ1 → ϕ1(x) in H1(Ω) (and also in L2(Ω)). Moreover, if h(ξ1u) = O (|ξ1|p) as |ξ1| → ∞ uniformly in u ∈ R, then ‖U(x)‖H1(Ω) = O ( |ξ1|p/2 ) as |ξ1| → ∞. Proof. By Theorem 2.1 we have a solution curve (u(x), µ1)(ξ1). From (1.1), ∆U + λ1U + h(ξ1ϕ1 + U) = µ1ϕ1 + e . Letting here U = ξ1V , obtain ∆V + λ1V = −h(ξ1(ϕ1 + V )) ξ1 + µ1 ξ1 ϕ1 + e ξ1 . Multiplying by V and integrating, we conclude that ∫ Ω |∇V |2 dx = O ( |ξ1|p−1 ) , as ξ1 → ±∞, and the lemma follows (since V ⊥ ϕ1 in L2(Ω), obtain − ∫ Ω ∆V V dx =∫ Ω |∇V |2 dx ≥ λ2 ∫ Ω V 2 dx). � Next we present one of our main results. Theorem 2.3. There exist h(u) ∈ C2(R) for which problem (1.1) has infinitely many positive solutions or any g(x, y) ∈ L2(B) ∩ Cα(B). Moreover, all solutions of (1.1) lie on a continuous solution curve (u(x, y), µ1)(ξ1), and µ1(ξ1) oscillates toward ±∞, as ξ1 → +∞. 304 P. KORMAN, D. SCHMIDT EJDE/SI/01 Proof. We exhibit h(u) satisfying the conditions of Theorem 2.1, for which µ1(ξ1) oscillates toward ±∞, as ξ1 → +∞. Take h(u) = √ u sin[ln(u3/2 + 1)] for u ≥ 0 (observe that h(u) ∈ C2[0,∞) and h′(u) < 1 for u ∈ [0,∞)), then extend h(u) arbitrarily to (−∞, 0) so that h(u) ∈ C2(R) and the conditions (2.3) and (2.4) hold. Calculate (for u ≥ 0) h′(u) = sin[ln(u3/2 + 1)] 2 √ u + 3u cos[ln(u3/2 + 1)] 2(u3/2 + 1) , and one of the anti-derivatives of h(u): H(u) = 1 3 (u3/2 + 1 ( sin[ln(u3/2 + 1)]− cos[ln(u3/2 + 1)] ) = √ 2 3 (u3/2 + 1) sin[ln(u3/2 + 1)− π 4 ] . Multiply the equation in (1.1) by ϕ1 and integrate over B, µ1 = ∫ B h (ξ1ϕ1 + U(x, y))ϕ1 dx dy . (2.5) Since ‖U(x, y)‖L2(B) = O(ξ 1 4 1 ) by Lemma 2.2, and ‖h′(ξ1ϕ1 +U)‖L2(B) = O(1/ξ 1/2 1 ) as ξ1 →∞, we have by the mean-value theorem ‖h (ξ1ϕ1 + U)− h (ξ1ϕ1) ‖L2(B) = o(1) . Then (2.5) becomes µ1(ξ1) = 2π ∫ 1 0 h(ξ1ϕ1)ϕ1r dr + o(1) . (2.6) Denoting f(r) = rϕ1(r) ϕ′1(r) , and integrating by parts (using ϕ1(0) = c0)∫ 1 0 h (ξ1ϕ1)ϕ1r dr = 1 ξ1 ∫ 1 0 f(r) d dr H (ξ1ϕ1) dr = − 1 ξ1 ∫ 1 0 f ′(r)H (ξ1ϕ1) dr − 1 ξ1 f(0)H (c0ξ1) . (2.7) The second term on the right oscillates toward ±∞ as ξ1 → ∞ with the am- plitude approaching √ 2 3 c 3/2 0 |f(0)|ξ1/2 1 . We show next that it dominates the first term on the right in (2.7). A computation shows that f ′(r) > 0 for all r ∈ (0, 1), and hence the first term is estimated in absolute value by a quantity approaching√ 2 3 c 3/2 0 ξ 1/2 1 ∫ 1 0 f ′(r)J 3/2 0 (ν1r) dr. Calculate f(0) = ϕ1(0) ϕ′′1 (0) = − 2 ν2 1 , so that |f(0)| = 2 ν2 1 ≈ 0.34. (From ϕ′′1 + 1 rϕ ′ 1 + ν2 1ϕ1 = 0, it follows that 2ϕ′′1(0) = −ν2 1ϕ1(0).) An- other computation shows that ∫ 1 0 f ′(r)J 3/2 0 (ν1r) dr ≈ 0.1, so that the second term dominates in (2.7). Since g(x, y) ∈ Cα(B) the convergence u(x,y) ξ1 → ϕ1(r) is in C2(B) by the elliptic regularity. It follows that u(x, y) > 0 for large ξ1, and hence it satisfies (1.1) with the original h(u) (before the extension). We conclude that for any g(x, y) ∈ L2(B) ∩ Cα(B) the problem (1.1) has infinitely many positive solutions. � In Figure 1 we present an approximation of the solution curve µ1 = µ1(ξ1) for the problem (1.1), with h(u) = √ u sin [ ln(u3/2 + 1) ] , computed using the formula (2.6). We performed computations on very large intervals (along both ξ1 and µ1 EJDE-2021/SI/01 INFINITELY MANY SOLUTIONS 305 -Log[- ] for <-1 Log[ ] for >1 Log[ξ] 10 20 30 40 -20 -10 10 20 Figure 1. Solution curve µ1 = µ1(ξ1) of (1.1), oscillating to ±∞. Values with |µ1| < 1 are not shown. axes), and to make the resulting picture manageable a logarithmic scale is used for both ξ1 and µ1. (Log denotes the natural logarithm in Mathematica. Note that the curve in this presentation is not continuous. It only appears so, since all values with |µ| < 1 have been omitted.) The result is the solution curve oscillating toward ±∞. In Figure 2 we show the same curve in the original (ξ1, µ1) coordinates. It is not apparent from that picture that µ1(ξ1) → ±∞ as ξ1 → ∞, since µ1(ξ1) keeps the same sign on large, and ever increasing intervals. 200 400 600 800 1000 ξ -15 -10 -5 5 10 15 20 μ Figure 2. Solution curve µ1 = µ1(ξ1) of (1.1) for ξ1 ∈ (0, 1000) There are other h(u) which can be handled similarly (including Mathematica being able to calculate the integral H(u) in elementary functions). We mention h(u) = u sin[ln(u2 + 1)] and h(u) = sin [ln(u+ 1)]. Next we turn to more “natural” nonlinearities h(u). Theorem 2.4. If 0 < p < 1 there is a continuum of solutions of (1.3) (u(x), µ1)(ξ1) parameterized by the first harmonic ξ1 ∈ (−∞,∞). Along this continuum, limξ1→±∞ µ1(ξ1) = 0. Moreover, the asymptotic formula (2.10) below holds. 306 P. KORMAN, D. SCHMIDT EJDE/SI/01 Proof. The existence of a continuum of solutions of (1.3) (u(x), µ1)(ξ1) parameter- ized by the first harmonic ξ1 ∈ (−∞,∞) follows by the result of Schaaf and Schmitt [16], which was based on Dancer [2], see also Korman [11]. By the Lemma 2.2 and elliptic regularity it follows that u/ξ1 → ϕ1 in C2(B), as |ξ1| → ∞. We now derive an asymptotic formula for µ1(ξ1). Multiply the equation (1.3) by ϕ1, and integrate over B, then use Lemma 2.2 and elliptic regularity µ1 = ∫ B up sinuϕ1 dx dy = ∫ 2π 0 ∫ 1 0 up sinuϕ1r dr dθ = 2πξp1 ∫ 1 0 ϕp+1 1 sin (ξ1ϕ1) r dr + o(ξp1) . (2.8) Integration by parts gives∫ 1 0 ϕp+1 1 sin (ξ1ϕ1) r dr = ∫ 1 0 ϕp+1 1 r ξ1ϕ′1 d dr [− cos (ξ1ϕ1)] dr = − 1 ξ1 g(r) cos (ξ1ϕ1) ∣∣1 0 + 1 ξ1 ∫ 1 0 g′(r) cos (ξ1ϕ1) dr , (2.9) where g(r) = ϕp+1 1 r/ϕ′1. Observe that g(1) = 0, while g(0) = cp+1 0 lim r→0 r ϕ′1(r) = cp+1 0 ϕ′′1(0) = −2cp0 ν2 1 . (From ϕ′′1 + 1 rϕ ′ 1 + ν2 1ϕ1 = 0, it follows that 2ϕ′′1(0) = −c0ν2 1 .) Hence −g(r) cos (ξ1ϕ1) ∣∣1 0 = −2cp0 cos c0ξ1 ν2 1 . One checks that g′(r) ∈ C[0, 1] is a bounded function. It follows that the oscillating integral ∫ 1 0 g′(r) cos (ξ1ϕ1) dr is o(1). From (2.8), µ1 = −4πξp−1 1 cp0 cos c0ξ1 ν2 1 + o ( ξp−1 1 ) . (2.10) It follows that µ1(ξ1)→ 0 as ξ1 →∞, concluding the proof. � In limiting case p = 1 the formula (2.10) (while not rigorously justified) indicates that µ1(ξ1) is asymptotic to a multiple of cos c0ξ1, suggesting that there is a µ0 > 0 so that for |µ| < µ0 the problem (1.3) has infinitely many solutions. The same conclusion is supported by our numerical computations, including the following example. Example 2.5. We computed the solution curve µ1 = µ1(ξ1) for the following example, with the linear part at resonance, ∆u+ λ1u+ u sinu = µ1ϕ1(r) + xy for (x, y) ∈ B , u = 0 on ∂B . (2.11) Observe that ∫ B xyϕ1(r) dx dy = 0. The solution curve µ1 = µ1(ξ1) (solid line) is presented in Figure 3. Notice a close agreement with the asymptotic formula (2.10) (dashed line). EJDE-2021/SI/01 INFINITELY MANY SOLUTIONS 307 5 10 15 20 25 ξ1 -2 -1 1 2 μ1 Figure 3. Solution curve µ1 = µ1(ξ1) of (2.11), compared with (2.10) Example 2.6. We computed the solution curve µ1 = µ1(ξ1) for an example of (1.3), with p = 1/2, ∆u+ λ1u+ u1/2 sinu = µ1ϕ1(r) + x2y − 3xy4 for x ∈ B , u = 0 on ∂B . (2.12) Observe that ∫ B (x2y−3xy4)ϕ1(r) dx dy = 0. The solution curve µ1 = µ1(ξ1) (solid line) is presented in Figure 4. The solutions u(x, y) are not radially symmetric, although they get arbitrarily close to radially symmetric functions as |ξ1| → ∞, according to our results. Again, computations show a good agreement with the asymptotic formula (2.10) (dashed line). 5 10 15 20 25 ξ1 -1.0 -0.5 0.5 μ1 Figure 4. Solution curve µ1 = µ1(ξ1) of (2.12), compared with (2.10) 3. Asymptotic formula in case of a rectangle Let R = {0 < x < a} × {0 < y < b} be a rectangle in R2. In this section we present computation of the solution curve µ1 = µ1(ξ1) on R for the nonlinearity 308 P. KORMAN, D. SCHMIDT EJDE/SI/01 considered above ∆u+ λ1u+ u sinu = µ1ϕ1(x, y) + e(x, y) for (x, y) ∈ R , u = 0 on ∂R , (3.1) and derive an asymptotic formula for µ1(ξ1). Here the principal eigenfunction ϕ1(x, y) = 2√ ab sin π ax sin π b y satisfies ∫ R ϕ2 1 dx dy = 1, and the corresponding princi- pal eigenvalue of −∆ is λ1 = π2 a2 + π2 b2 . It is assumed that ∫∫ R e(x, y)ϕ1(x, y) dx dy = 0. We decompose the solution of (3.1) as u(x, y) = ξ1ϕ1(x, y) + U(x, y), with∫∫ R U(x, y)ϕ1(x, y) dx dy = 0, and ξ1 = ∫∫ R u(x, y)ϕ1(x, y) dx dy. Unlike the case of a ball, we shall use the stationary phase method to derive the following asymptotic formula for ξ1 large: µ1(ξ1) ∼ 4 √ ab π sin ( 2√ ab ξ1 − π 2 ) , (3.2) after we recall the following known lemmas. Lemma 3.1. Assume that f(x) ∈ C2(x0−a0, x0+a0) for some constants α , a0 > 0. Then as µ→∞∫ x0+a0 x0−a0 f(x)e−iαµ(x−x0)2 dx = e−i π 4 √ π αµ f(x0) +O ( 1 µ ) . This lemma follows from the part (i) of the following more general lemma, which we will also need, see e.g., [8, p. 83] for the proof. Lemma 3.2. (i) Assume that f(x) and g(x) are of class C2[a, b] and g(x) has a unique critical point x0 on [a, b], and moreover x0 ∈ (a, b) and g′′(x0) 6= 0 (so that x0 gives a global max or global min on [a, b]). Then as µ → ∞ the following asymptotic formula holds∫ b a f(x)eiµg(x) dx = ei[µg(x0)±π4 ] √ 2π µ|g′′(x0)| f(x0) +O ( 1 µ ) , where one takes “plus” if g′′(x0) > 0 and “minus” if g′′(x0) < 0. (ii) Assume that the functions f(x) and g(x) > 0 are of class C2[0, 1], and satisfy g′(x) < 0 for allx ∈ (0, 1], and g′(0) = 0, g′′(0) < 0 . Then, as µ→∞,∫ 1 0 f(x)eiµg(x)dx = ei(µg(0)−π4 ) √ π 2µ|g′′(0)| f(0) +O ( 1 µ ) . Turning to the derivation of (3.2), multiply (3.1) by ϕ1 and integrate over R. Then use that u ∼ ξ1ϕ1 for ξ1 large (sinu ∼ sin ξ1ϕ1 by the mean-value theorem) µ1 = ∫∫ R u sinuϕ1 dx dy ∼ ξ1 ∫∫ R ϕ2 1 sin ξ1ϕ1 dx dy = ξ1Im ∫∫ R ϕ2 1e iξ1ϕ1 dx dy = 4ξ1 ab Im ∫∫ R sin2 π a x sin2 π b ye i 2√ ab ξ1 sin π a x sin π b y dx dy . EJDE-2021/SI/01 INFINITELY MANY SOLUTIONS 309 For ξ1 large there are fast oscillations around zero, except near the stationary point (x0, y0) = (a2 , b 2 ). The approximation of the integral near the stationary point provides the dominant contribution to this integral. Using Taylor’s formula near (x0, y0) sin π a x sin π b y ≈ 1− π2 2a2 (x− x0)2 − π2 2b2 (y − y0)2 , because the second mixed partials vanish at (x0, y0). Then µ1 ∼ 4ξ1 ab Im ( e i 2√ ab ξ1 ∫ a 0 sin2 π a xe−iξ1α1(x−x0)2 dx ∫ b 0 sin2 π b ye−iξ1α2(y−y0)2 dy ) , where α1 = π2 a2 √ ab and α2 = π2 b2 √ ab . Using Lemma 3.1, for large ξ1, µ1 ∼ 4ξ1 ab Im ( e i 2√ ab ξ1e−i π 4 √ π α1ξ1 e−i π 4 √ π α2ξ1 ) = 4 √ ab π sin ( 2√ ab ξ1 − π 2 ) . This formula, as well as our numerical calculations, suggests that there exist two constants 0 < a < A so that the problem (3.1) has infinitely many solutions for µ1 ∈ (−a, a), there are at most finitely many solutions for µ1 outside of (−a, a), and no solutions exist for |µ1| > A. Example 3.3. On the rectangle R1 = {0 < x < 1} × {0 < y < 2} we computed the solution curve µ1 = µ1(ξ1) for the problem ∆u+ λ1u+ u sinu = µ1ϕ1(x, y) + (x− 1 2 )(y − 1) for (x, y) ∈ R1 , u = 0 on ∂R1 , (3.3) with λ1 = 5π2/4 and ϕ1(x, y) = √ 2 sinπx sin π 2 y. Observe that ∫∫ R (x − 1 2 )(y − 1)ϕ1(x, y) dx dy = 0. The solution curve µ1 = µ1(ξ1) (solid line) is presented in Figure 5. Once again, we have a close agreement with the asymptotic formula (3.2) (dashed line). As in the case of balls, oscillations of µ1(ξ1) are decaying in the dimensions n > 2, as will follow from the asymptotic formula that we present next. Consider the n- dimensional rectangle R = (0, a1) × (0, a2) × · · · × (0, an), and the problem (3.1) on R. The principal eigenfunction of the Laplacian on R with Dirichlet boundary conditions, and satisfying ∫ R ϕ2 1 dx = 1, is ϕ1 = 2n/2 √ a1a2 · · · an sin π a1 x1 sin π a2 x2 · · · sin π an xn . As in the above derivation, one shows that µ1(ξ1) ∼ 2 n 2 (3−n2 ) (a1a2 · · · an) n/4 πn/2 ξ 1−n2 1 sin ( 2n/2 √ a1a2 · · · an ξ1 − n π 4 ) → 0 , as ξ1 →∞, for n > 2. 4. Numerical computations We now describe the Mathematica program that was used to produce the solution curves µ1 = µ1(ξ1), presented above. To avoid ambiguity of notation, in this section we shall write ξ instead of ξ1 and µ instead of µ1, so that u = ξϕ1+U , and µ = µ(ξ). 310 P. KORMAN, D. SCHMIDT EJDE/SI/01 5 10 15 20 ξ1 -2 -1 1 2 μ1 Figure 5. Solution curve µ1 = µ1(ξ1) of (3.3), compared with (3.2) Our program handles semilinear equations on rather general domains solving ∆u+ f(u) = µϕ1(x, y) + e(x, y) for (x, y) ∈ Ω , u = 0 on ∂Ω (4.1) on bounded domains Ω ∈ R2, including rectangles and ellipses in R2, and radially symmetric solutions on balls in Rn. Here ∆u = uxx(x, y) +uyy(x, y), while (λ1, ϕ1) is the principal eigenpair of the Laplacian on Ω, with zero boundary conditions and∫ Ω ϕ2 1 dx dy = 1, µ ∈ R, e(x, y) ∈ ϕ⊥1 in L2(Ω). Choosing a step size h and an initial value ξ = ξ0, let ξn = ξ0 + nh. We are looking for µ = µn for which the problem (4.1) has a solution u(x) ≡ un(x, y), with∫ Ω u(x, y)ϕ1(x, y) dx dy = ξn . (4.2) We utilize Mathematica’s ability (the NDSolve command) to solve linear Dirichlet problems of the type ∆u+ a(x, y)u = b(x, y) for (x, y) ∈ Ω , u = 0 on ∂Ω on some domains Ω ∈ R2, including ellipses around the origin, and rectangles. Assuming that (µn, un(x, y))(ξn) is already computed, we use Newton’s method to calculate (µn+1, un+1(x, y))(ξn+1). We calculate (µn+1, un+1(x, y)) using a se- quence of iterates (µk, uk(x, y)) beginning with (µ0, u0(x, y)) = (µn, un(x, y)). As- suming that (µk, uk(x, y)) is already computed, we approximate f(u) ≈ f(uk) + f ′(uk)(u− uk), and solve the linear problem ∆w + f ′(uk)w = µ1ϕ1(x, y) + f ′(uk)uk − f(uk) + e(x, y) for (x, y) ∈ Ω , w = 0 on ∂Ω for µ and w by the algorithm that is described next. One can decompose the solution in the form w(x, y) = µw1(x, y) +w2(x, y), where w1 and w2 are solutions of ∆w1 + f ′(uk)w1 = ϕ1(x, y) , for (x, y) ∈ Ω , w1 = 0 on ∂Ω , ∆w2 + f ′(uk)w2 = f ′(uk)uk − f(uk) + e(x, y) , for (x, y) ∈ Ω , w2 = 0 on ∂Ω . (4.3) EJDE-2021/SI/01 INFINITELY MANY SOLUTIONS 311 After calculating w1 and w2, we look for µ such that ∫ Ω w(x, y)ϕ1(x, y) dx dy = ξn+1, and declare that value of µ to be our new iterate µk+1, so that µk+1 = ξn+1 − ∫ Ω w2(x, y)ϕ1(x, y) dx dy∫ Ω w1(x, y)ϕ1(x, y) dx dy . (4.4) The corresponding w(x, y) is our next iterate uk+1(x, y) = µk+1w1(x, y) +w2(x, y). The iterations are stopped once the relative error µk+1−µk µk is small enough. Controlling the accuracy of iterates is the major improvement of the present algorithm, compared with the one we used in [12]. Remark 4.1. There is a better way of choosing the initial iterate u0(x, y) at ξn+1 (corresponding to the “predictor” of the predictor-corrector method). Write u = u(x, y, ξ). Approximate u(x, y, ξn+1) ≈ u(x, y, ξn) + uξh = un(x, y) + uξh. Differentiate the equation (4.1) in ξ, and compare the result with the first formula in (4.3) to get: uξ = µ′(ξn)w1 ≈ µn+1−µn h w1. So that we take u0(x, y) = un(x, y) + (µn+1 − µn)w1, using the last function w1 computed at ξn. Our experiments show considerably faster convergence. 5. Asymptotics and numerics for radial solutions Let B denote the unit ball around the origin in Rn, x ∈ Rn and r = |x|. We present computations of the solution curve µ1 = µ1(ξ1) for the radial solutions u = u(r) of the model problem u′′(r) + n− 1 r u′(r) + λ1u+ sinu = µ1ϕ1(r) + e(r) , for 0 < r < 1 u′(0) = u(1) = 0 , (5.1) and derive an asymptotic formula for µ1(ξ1). By Theorem 2.1 there exists a con- tinuous solution curve (u(r), µ1)(ξ1) that exhausts the solution set of (5.1), and moreover u(r) ξ1 → ϕ1 in C2(B) as ξ1 → ∞. (The same result holds if one replaces e(r) by e(x) ∈ ϕ⊥1 .) Restricting to radial solutions we can perform numerical com- putations for n ≥ 3.) Recall that the principal eigenfunction of the Laplacian on B is ϕ1(r) = c0r −n−2 2 Jn−2 2 (ν1r), where ν1 denotes the first root of the Bessel function Jn−2 2 (r), and c0 is chosen so that∫ B ϕ2 1(r) dx = ωnc 2 0 ∫ 1 0 J2 n−2 2 (ν1r)r dr = 1 , which is c0 = 1√ ωn ∫ 1 0 J2 n−2 2 (ν1r)r dr . Here ωn = nπn/2 Γ(n2 +1) gives the area of the unit ball in Rn, where Γ denotes the gamma function. The corresponding principal eigenvalue is λ1 = ν2 1 . It is assumed that∫∫ B e(r)ϕ1(r) dx = ωn ∫ 1 0 e(r)ϕ1(r)rn−1 dr = 0. We decompose the solution of (3.1) as u(r) = ξ1ϕ1(r) + U(r), with ∫∫ B U(r)ϕ1(r) dx = 0, and ξ1 = ∫∫ B u(r)ϕ1(r) dx. We now derive an asymptotic formula for the function µ1(ξ1), by using that u(x) ∼ ξ1ϕ1 for large ξ1. Multiplying the PDE version of equation (5.1) by ϕ1 and 312 P. KORMAN, D. SCHMIDT EJDE/SI/01 integrating over B µ1 = ωn ∫ 1 0 (sinu(r))ϕ1(r)rn−1 dr ∼ ωn ∫ 1 0 (sin ξ1ϕ1)ϕ1(r)rn−1 dr = ωn ξ1 ∫ 1 0 ϕ1(r)rn−1 ϕ′1(r) d dr (− cos ξ1ϕ1) dr . (5.2) Integrating by parts, and denoting f1(r) ≡ ϕ1(r)rn−1/ϕ′1(r), obtain (ϕ1(1) = 0) µ1 ∼ ωn ξ1 ( f1(0) cos ξ1ϕ1(0) + ∫ 1 0 f ′1(r) cos ξ1ϕ1(r) dr ) . (5.3) The next steps depend on the dimension n. Assume that n = 2. Then ω2 = 2π, ϕ1(r) = c0J0(ν1r), where ν1 is the first root of the Bessel’s function J0(x), c0 is given by (2.1), f1(r) = ϕ1(r)r ϕ′1(r) , and f1(0) = ϕ1(0) ϕ′′1(0) = − 2 ν2 1 . We conclude that for large ξ1 (using that ϕ1(0) = c0, and that the oscillatory integral in (5.3) is o(1)) µ1 ∼ − 4π ξ1ν2 1 cos c0ξ1 , (5.4) which is consistent with (2.10) at p = 0. Assume now that n = 3. Then ω3 = 4π. Since J 1 2 (x) = sin x√ x , we have ν1 = π2, ϕ1(r) = 1√ 2π sinπr r , ϕ1(0) = √ π 2 . Also, f1(r) = ϕ1r 2 ϕ′1 . It follows that f1(0) = 0, so that the first term in (5.3) vanishes. Then by Lemma 3.2 (ii), µ1 ∼ 4π ξ1 Re ∫ 1 0 f ′1e iξ1ϕ1 dr = 4π ξ1 Re ei(ξ1ϕ1(0)−π4 ) √ π 2ξ1|ϕ′′1(0)| f ′1(0) = 4π ξ1 cos ( ξ1 √ π 2 − π 4 )√ π 2ξ1|ϕ′′1(0)| f ′1(0) . As above, calculate 3ϕ′′1(0) = −π2ϕ1(0), so that |ϕ′′1(0)| = π 5 2 3 √ 2 . A short calcula- tion shows that f ′1(0) = ϕ1(0) ϕ′′1(0) = − 3 π2 . We conclude that for large ξ1, µ1 ∼ − 12 √ 3 √ 2 √ 2ξ 3 2 1 π 7 4 cos ( ξ1 √ π 2 − π 4 ) . (5.5) For numerical computations of radial solutions we used the algorithm described in a previous section. Computations can be performed accurately in any dimension, EJDE-2021/SI/01 INFINITELY MANY SOLUTIONS 313 because the linearization of (5.1) involves a boundary value problem for an ODE, readily handled by Mathematica. The formula (4.4) takes the form µk+1 = ξn+1 − ωn ∫ 1 0 w2(r)ϕ1(r)rn−1 dr ωn ∫ 1 0 w1(r)ϕ1(r)rn−1 dr . Example 5.1. For n = 3 we solved the problem (5.1) with e(r) = cosπr r (observe that ∫ B e(r)ϕ1(r) dx = 0), see Figure 6. The solution of (5.1) (solid line) is in a good agreement with the asymptotic formula (5.5) (dashed line). 5 10 15 20 25 30 ξ1 -0.05 0.05 0.10 0.15 μ1 Figure 6. Solution curve µ1 = µ1(ξ1) of (5.1), n = 3, compared with the asymptotic formula (5.5) References [1] D. Costa, H. Jeggle, R. Schaaf K. Schmitt; Oscillatory perturbations of linear problems at resonance, Results in Mathematics, 14 (1988), 275–287. [2] E. N. Dancer; On the use of asymptotics in nonlinear boundary value problems, Ann. Mat. Pura Appl., 131 (1982), (4), 167–185. [3] D. G. de Figueiredo, W.-M. Ni; Perturbations of second order linear elliptic problems by nonlinearities without Landesman-Lazer condition, Nonlinear Anal., 3 (1979) no. 5, 629– 634. [4] R. Iannacci, M. N. Nkashama, J. R. Ward, Jr; Nonlinear second order elliptic partial differ- ential equations at resonance, Trans. Amer. Math. Soc., 311 (1989), no. 2, 711–726. [5] A. Galstian, P. Korman, Y. Li; On the oscillations of the solution curve for a class of semilinear equations, J. Math. Anal. Appl., 321 (2006), no. 2, 576–588. [6] B. Gidas, W.-M. Ni, L. Nirenberg; Symmetry and related properties via the maximum prin- ciple, Commun. Math. Phys., 68 (1979), 209–243. [7] P. Korman; An oscillatory bifurcation from infinity, and from zero, NoDEA Nonlinear Dif- ferential Equations Appl. 15, (2008) no. 3, 335–345. [8] P. Korman; Global Solution Curves for Semilinear Elliptic Equations, World Scientific, Hack- ensack, NJ, 2012. [9] P. Korman; Curves of equiharmonic solutions, and problems at resonance, Discrete Contin. Dyn. Syst., 34 (2014) no. 7, 2847–2860. [10] P. Korman; Nonlinear elliptic equations and systems with linear part at resonance, Electron. J. Differential Equations, Paper No. 67 (2016). [11] P. Korman; Global solution curves in harmonic parameters, and multiplicity of solutions, J. Differential Equations, 296 (2021), 186–212. 314 P. KORMAN, D. SCHMIDT EJDE/SI/01 [12] P. Korman, D.S. Schmidt; Continuation of global solution curves using global parameters, Preprint, arXiv:2001.00616. [13] E. M. Landesman, A. C. Lazer; Nonlinear perturbations of linear elliptic boundary value problems at resonance, J. Math. Mech., 19 (1970), 609–623. [14] Y. Pinchover, J. Rubinstein; An Introduction to Partial Differential Equations. Cambridge University Press, Cambridge, 2005. [15] R. Schaaf, K. Schmitt; A class of nonlinear Sturm-Liouville problems with infinitely many solutions, Trans. Amer. Math. Soc., 306 (1988) no. 2, 853–859. [16] R. Schaaf, K. Schmitt; Asymptotic behavior of positive solution branches of elliptic problems with linear part at resonance, Z. Angew. Math. Phys., 43 (1992) no. 4, 645–676. [17] G. N. Watson; A treatise on the theory of Bessel functions, Cambridge University Press, Cambridge, 1995. Philip Korman Department of Mathematical Sciences, University of Cincinnati, Cincinnati, OH 45221- 0025, USA Email address: kormanp@ucmail.uc.edu Dieter S. Schmidt Department of Computer Science, University of Cincinnati, Cincinnati, OH 45221-0030, USA Email address: schmiddr@ucmail.uc.edu 1. Introduction 2. Global solution set for a ball in R2 3. Asymptotic formula in case of a rectangle 4. Numerical computations 5. Asymptotics and numerics for radial solutions References