On Modified Second Paine–de Hoog–Anderssen Boundary Value Problem

Abstract

This article deals with a special case of the Sturm–Liouville boundary value problem (BVP), an eigenvalue problem characterized by the Sturm–Liouville differential operator with unknown spectra and the associated eigenfunctions. By examining the BVP in the Schrödinger form, we are interested in the problem where the corresponding invariant function takes the form of a reciprocal quadratic form. We call this BVP the modified second Paine–de Hoog–Anderssen (PdHA) problem. We estimate the lowest-order eigenvalue without solving the eigenvalue problem but by utilizing the localized landscape and effective potential functions instead. While for particular combinations of parameter values that the spectrum estimates exhibit a poor quality, the outcomes are generally acceptable although they overestimate the numerical computations. Qualitatively, the eigenvalue estimate is strikingly excellent, and the proposal can be adopted to other BVPs.

1. Introduction

The Sturm–Liouville boundary value problem (BVP) is an active research area in mathematics and its applications [1,2]. While the formulation originally comes from applied problems [3,4], the mathematical theory has reached maturity, both in terms of analysis and its numerical solutions [5,6,7,8,9].

One of the intriguing problems in the Sturm–Liouville BVP is finding eigenvalues and their corresponding eigenfunctions. Some attempts to estimate eigenvalues were performed more than a half-century ago [10,11]. Various numerical techniques have been developed, including the Rayleigh–Ritz method [12], matrix-variational method [13], shooting method [14], Pruess’ coefficient approximation method [15], spectral parameter power series method [16] and modified integral series methods among others [17]. A body of literature is continuously being updated by the modification and improvement of the existing techniques, as well as with applications of the Sturm–Liouville problem in numerous physical situations.

The history of the Sturm–Liouville problem goes back to the 19th century when the mathematicians Jacques Charles François Sturm (1803–1855) from Geneva, then part of France, and Joseph Liouville (1809–1882) from France, investigated particular problems of second-order linear differential equations under appropriate boundary conditions and the properties of the corresponding solutions. Their results were published in a sequence of papers in 1836–1837.

The problems arise frequently not only as governing equations of motion based on Newtonian mechanics but also when dealing with the method of separation of variables in tackling linear partial differential equations (PDEs). There are also applications in quantum mechanics where the time-independent Schrödinger equation for a particle with fixed angular momentum quantum numbers moves in a spherically symmetric potential at the energy level associated with spectra [18,19].

The classical Sturm–Liouville problem is composed of a linear second-order ordinary differential equation (ODE) and a set of boundary conditions at the endpoints of the interval where the problem is defined. It can be written in the following self-adjoint (or canonical) form

$$\begin{array}{cc}\hfill -\frac{d}{dz}\left(u\left(z\right)\frac{dy}{dz}\right)+v\left(z\right)y& =\lambda w\left(z\right)y,z_0<z<z_1,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill a_{0}y\left(z_0\right)-a_{1}u\left(z_0\right)\frac{dy}{dz}\left(z_0\right)=0,& b_{0}y\left(z_1\right)+b_{1}u\left(z_0\right)\frac{dy}{dz}\left(z_1\right)=0.\hfill \end{array}$$

(2)

An eigenvalue (or spectrum) is a value of $\lambda$ such that the ODE (1) possesses a nontrivial solution $y\left(z\right)$ subject to the prescribed boundary conditions. This solution is called the corresponding eigenfunction associated with each $\lambda$ and is unique up to scalar multiplications. Note that while some authors impose the classical derivatives $y^{\prime }$ at the boundaries, we opt for quasi-derivatives $\left(u{y}^{\prime }\right)$ instead, where prime denotes the differentiation with respect to the variable z. Furthermore, neither $a_0$ and $a_1$ nor $b_0$ and $b_1$ are both zero.

The regular Sturm–Liouville problem operates under the assumption that the coefficient functions are well-behaved. In this case, $u\left(z\right)$ and $w\left(z\right)$ are strictly positive, and $u\left(z\right)$, $u^{\prime }\left(z\right)$, $v\left(z\right)$, and $w\left(z\right)$ are continuous functions on the bounded closed interval $[z_0,z_1]$. As a consequence, there exists an infinite sequence of eigenvalues ${\lambda }_0<{\lambda }_1<{\lambda }_2<\cdots$ and eigenfunctions $y_0\left(z\right)$, $y_1\left(z\right)$, $y_2\left(z\right)$, ⋯ such that $y_n\left(z\right)$ has only n zeros on the open interval $(z_0,z_1)$. In addition, distinct eigenfunctions are orthogonal with respect to the weight function $w\left(z\right)$, i.e.,

$${\int }_{z_0}^{z_1}w\left(z\right)y_i\left(z\right)y_j\left(z\right)dz=0,\mathrm{whenever}i\ne j.$$

The singular Sturm–Liouville problem comprises whether the coefficient functions have singularities at the boundary points or whether the interval is unbounded. While the regular Sturm–Liouville problem was introduced by and the focus of Sturm’s and Liouville’s works, the singular Sturm–Liouville problem was initiated by Hermann Weyl from Göttingen, Germany, who investigated some ODEs with singularities and proposed the topic of the essential spectrum in 1910 [20]. Indeed, the progressive development of quantum theory in the 1920s and 1930s and a breakthrough in the general spectral theorem for unbounded self-adjoint operators in the Hilbert space provided an impetus for further examination into the spectral theory of Sturm–Liouville self-adjoint differential operators [6,21].

The body of literature on Sturm–Liouville theory is overwhelmingly plentiful. Thanks to the advancement in computational software and hardware, the numerical aspect of the Sturm–Liouville problem is particularly one of the most active research areas of mathematical physics. The following brief literature review is by no means exhaustive. While a comparison between the Sinc–Galerkin technique and differential transform method reveals that the latter is more efficient than the former [22], a comparative study between the Sinc–Galerkin technique and variational iteration method indicates that the former is better than the latter when dealing with singular problems [23].

During the past decade, both regular and singular fractional derivatives Sturm–Liouville problems and operators have gained a lot of attention [24,25,26]. In particular, Klimet and Agrawal showed for the first time that the eigenvalues of two classes of the fractional Sturm–Liouville operators are real-valued, and the eigenfunctions associated with two distinct eigenvalues are orthogonal [25]. Furthermore, the existence of a countable set of orthogonal eigenfunctions for the regular fractional Sturm–Liouville problem was proven by the methods of fractional variational analysis [27]. Recently, an efficient numerical method for estimating eigenvalues and eigenfunctions of the Caputo-type fractional Sturm–Liouville problem based on the Lagrange polynomial interpolation was proposed [28].

When it comes to applications, the Sturm–Liouville problem features an abundance of them in applied mathematics and physics. The BVP itself originally sprang from the heat conduction problem in a nonhomogeneous thin bar, but equally classic problems from PDE include the vibration of plucked strings, thin membranes, and sound waves. The Sturm–Liouville BVP emerges naturally as an immediate consequence of implementing the method of separation of variables. For these vibration problems, the eigenvalues determine the frequency of oscillation, whereas the associated eigenfunctions correspond to the shape of the vibrating waves at any point in time [29].

The most straightforward applications of the Sturm–Liouville problem bring about the various Fourier series, and more sophisticated applications lead to the generalizations of Fourier series involving Bessel functions, Hermite polynomials, and other special functions. One example is an electrostatic field generated by a spherical capacitor, where Laplace’s equation in spherical coordinates serves as the governing equation and the Legendre polynomials appear in the solution [3]. Another example comes from a quantum-mechanical system, where the wave function of a particle is governed by a one-dimensional time-independent Schrödinger equation. In this case, the eigenvalues represent the energy levels of the atomic system and the eigenfunctions are the wave functions of observable physical quantities [30,31].

The focal point of this article is the so-called modified second Paine–de Hoog–Anderssen (PdHA) BVP [32]. What we refer to as the “first” and “second” PdHA problems deal with the corresponding invariant function of the ODE in the canonical form in the Sturm–Liouville—in this case, the PdHA problem. While the first PdHA problem employs an exponential invariant function, i.e., $q\left(\widehat{z}\right)=e^{\widehat{z}}$, the second “classical” PdHA problem adopts a special case of the reciprocal quadratic binomial function, i.e., $q\left(\tilde{z}\right)={(\widehat{z}+0.1)}^{-2}$. We do not examine the former in this article. Our focus is to investigate the latter by modifying the invariant function to involve several parameters, i.e., $q\left(\tilde{z}\right)=c{(a\widehat{z}+b)}^{-2}$, where $a,b$ and $c>0$. Hence, the name of the modified second PdHA problem. Even though the classical problem appeared more than four decades ago, investigating the modified version still proves to be illuminating, as we will encounter in this article.

Since finding eigenvalues and eigenfunctions is a central aspect of the Sturm–Liouville problem, our motivation is to provide an alternative for estimating the former without actually solving the original modified second PdHA BVP. Our main contribution is a quantitative comparison between the eigenvalues obtained from an estimate with numerical solution results. By combining different parameters, we observe how the eigenvalues behave as we vary the parameters. Although the particular choice of the problem is rather narrow and specific, we anticipate that our contribution may illuminate other problems that can be solved using a similar technique. We also invite comments and debates from expertise for potential approaches to improve the estimated accuracy.

This article is organized as follows. Section 2 deals with the modified second PdHA problem. Using a transformation well-known in the Sturm–Liouville theory, we verify that a particular Sturm–Liouville problem in the canonical form can be transformed into the modified second PdHA BVP, where the corresponding invariant function admits the form of reciprocal quadratic binomial function. Section 3 covers how to estimate the spectrum without solving an eigenvalue problem but by incorporating the landscape function and effective potential instead. Section 4 compares the eigenvalue estimate with numerical simulations for both Dirichlet and Neumann boundary conditions. Finally, Section 5 concludes our discussion.

2. Modified Second Paine–de Hoog–Anderssen Problem

We start with the following definition of the generalized second PdHA problem.

Definition 1. (Generalized second Paine–de Hoog–Anderssen (PdHA) problem). The following Sturm–Liouville BVP in the Schrödinger (or Liouville normal) form is defined as a generalized second Paine–de Hoog–Anderssen (PdHA) problem with the third-kind (Robin) boundary conditions:

$$\begin{array}{cc}\hfill -\frac{d^2\widehat{y}}{d{\widehat{z}}^2}+q\left(\widehat{z}\right)\widehat{y}& =\lambda \widehat{y},{\widehat{z}}_0<\widehat{z}<{\widehat{z}}_1,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\alpha }_0\widehat{y}\left({\widehat{z}}_0\right)-{\alpha }_1\dot{\widehat{y}}\left({\widehat{z}}_0\right)& =0,{\beta }_0\widehat{y}\left({\widehat{z}}_1\right)+{\beta }_1\dot{\widehat{y}}\left({\widehat{z}}_1\right)=0,\hfill \end{array}$$

(4)

where dot denotes the derivative with respect to $\widehat{z}$ and the invariant function q is given by

$$q\left(\widehat{z}\right)=\frac{c}{{\left(a\widehat{z}+b\right)}^n},a,b,c>0,n\in \mathbb{N}.$$

(5)

For $n=2$, we call the above a modified second PdHA problem. Figure 1 displays plots of the invariant function $q\left(\widehat{z}\right)$ for different values of parameters. The left panel of Figure 1 sketches the function for fixed values of a and c but for different values of b in a logarithmic scale. The right panel of Figure 1 illustrates the plots of q when a and b are fixed but the parameter c is varied. In both instances, the function is decreasing for an increasing value of $\widehat{z}$.

For $n\ne 2$, the problem seems to be open. Our conjecture is that the spectrum exhibits a similar qualitative behaviour with the one for $n=2$. Although for the BVP in the Liouville normal form can still be solved numerically, recovering the corresponding Sturm–Liouville problem in the canonical form turns out to be nontrivial.

Corollary 1.

In particular, for the modified second PdHA problem ($n=2$), when $a=1=c$, $b=0.1$, by imposing the Dirichlet conditions ${\alpha }_0=1={\beta }_0$, ${\alpha }_1=0={\beta }_1$, and, selecting the endpoints as ${\widehat{z}}_0=0$ and ${\widehat{z}}_1=\pi$, the above BVP reduces to the well-known classical second PdHA problem [32]:

$$-\frac{d^2\widehat{y}}{d{\widehat{z}}^2}+\frac{1}{{\left(\widehat{z}+\frac{1}{10}\right)}^2}\widehat{y}=\lambda \widehat{y},0<\widehat{z}<\pi ,\widehat{y}\left(0\right)=0=\widehat{y}\left(\pi \right).$$

We need the following theorem on a transformation of the Sturm–Liouville problem from the canonical form into the Schrödinger form [8]:

Theorem 1.

The Sturm–Liouville problem in the canonical form with eigenvalue λ and the corresponding eigenfunction $y\left(z\right)$ is given explicitly as follows:

$$\begin{array}{cc}\hfill -\frac{d}{dz}\left[u\left(z\right)\frac{dy}{dz}\right]=\lambda y,& z_0<z<z_1,\hfill \\ \hfill a_{0}y\left(z_0\right)-a_1\frac{dy}{dz}\left(z_0\right)=0,& b_{0}y\left(z_1\right)+b_1\frac{dy}{dz}\left(z_1\right)=0.\hfill \end{array}$$

can be converted into the following BVP in the Schrödinger (or Liouville normal) form

$$\begin{array}{cc}\hfill -\frac{d^2\widehat{y}}{d{\widehat{z}}^2}+q\left(\widehat{z}\right)\widehat{y}=\lambda \widehat{y},& {\widehat{z}}_0<\widehat{z}<{\widehat{z}}_1,\hfill \\ \hfill {\widehat{\alpha }}_0\widehat{y}\left(\widehat{z_0}\right)-{\widehat{\alpha }}_1\frac{d\widehat{y}}{d\widehat{z}}\left(\widehat{z_0}\right)=0,& {\widehat{\beta }}_0\widehat{y}\left(\widehat{z_1}\right)+{\widehat{\beta }}_1\frac{d\widehat{y}}{d\widehat{z}}\left(\widehat{z_1}\right)=0,\hfill \end{array}$$

by performing Liouville’s transformation

$$\widehat{z}=\int \frac{dz}{\sqrt{u}},v=u^{-1/4},\mathit{and}y\left(z\right)=v\left(z\right)\widehat{y}\left(z\right),$$

(6)

where

$$\begin{array}{ccc}\hfill {\alpha }_0& =a_{0}v\left({\widehat{z}}_0\right)-\frac{a_1}{\dot{z}\left({\widehat{z}}_0\right)}\frac{dv}{d\widehat{z}}\left({\widehat{z}}_0\right),\hfill & \hfill {\alpha }_1=a_1\frac{v\left(\widehat{z_0}\right)}{\dot{z}\left(\widehat{z_0}\right)},\\ \hfill {\beta }_0& =b_{0}v\left({\widehat{z}}_1\right)+\frac{b_1}{\dot{z}\left(\widehat{z_1}\right)}\frac{dv}{d\widehat{z}}\left({\widehat{z}}_1\right),\hfill & \hfill {\beta }_1=b_1\frac{v\left(\widehat{z_1}\right)}{\dot{z}\left(\widehat{z_1}\right)},\end{array}$$

and q is the corresponding invariant function of the ODE in the canonical form, given by

$$q\left(\widehat{z}\right)=v\frac{d^2}{d{\widehat{z}}^2}\left(\frac{1}{v}\right).$$

A detailed outline of the proof of this theorem is available in a recent preprint, where an application related to perturbed potential temperature distribution in an atmospheric boundary layer was also discussed via the WKB method and numerical simulation [33].

Proposition 1.

The following Sturm–Liouville problem in the canonical form with Dirichlet boundary conditions

$$\begin{array}{cc}\hfill -\frac{d}{dz}\left\{a{\left[(1-2k)(z-d)\right]}^{-4k/(1-2k)}{\displaystyle \frac{dy}{dz}}\right\}& =\lambda y,z_0<z<z_1,\hfill \\ \hfill y\left(z_0\right)=0=y\left(z_1\right),\end{array}$$

where

$$\begin{array}{cccc}\hfill k& =\frac{1}{2}\left(-1\pm \frac{1}{a}\sqrt{a^2+4c}\right),\hfill & \hfill z_0& =\frac{b^{1-2k}}{1-2k}+d,\hfill \\ \hfill a& \ne 0,d\in \mathbb{R},\hfill & \hfill z_1& =\frac{{(\pi +b)}^{1-2k}}{1-2k}+d,\hfill \end{array}$$

reduces to the modified second PdHA problem

$$-\frac{d^2\widehat{y}}{d{\widehat{z}}^2}+\frac{c}{{(a\widehat{z}+b)}^2}\widehat{y}=\lambda \widehat{y},0<\widehat{z}<\pi ,\widehat{y}\left(0\right)=0=\widehat{y}\left(\pi \right).$$

(7)

Proof. 

Since

$$u\left(z\right)={\left[a(1-2k)(z-d)\right]}^{-4k/(1-2k)},$$

then

$$\widehat{z}=\int \frac{dz}{\sqrt{u}}=\frac{1}{a}\left\{{\left[a(1-2k)(z-d)\right]}^{1/(1-2k)}-b\right\},$$

for which we can express z in terms of $\widehat{z}$

$$z=\frac{{\left(a\widehat{z}+b\right)}^{1-2k}}{a(1-2k)}+d.$$

We express the function u in terms of the transformed variable $\widehat{z}$:

$$u\left(\widehat{z}\right)={\left(a\widehat{z}+b\right)}^{-4k},$$

and the function v using the transformation $v=u^{-1/4}$:

$$v\left(\widehat{z}\right)={\left(a\widehat{z}+b\right)}^k.$$

It follows that

$$\begin{array}{cc}\hfill \frac{d}{d\widehat{z}}\left(\frac{1}{v}\right)& =-ak{\left(a\widehat{z}+b\right)}^{-k-1},\hfill \\ \hfill \frac{d^2}{d{\widehat{z}}^2}\left(\frac{1}{v}\right)& =a^{2}k(k+1){\left(a\widehat{z}+b\right)}^{-k-2},\hfill \\ \hfill q\left(\widehat{z}\right)=v\frac{d^2}{d{\widehat{z}}^2}\left(\frac{1}{v}\right)& =\frac{a^{2}k(k+1)}{{\left(\widehat{z}+b\right)}^2}=\frac{c}{{\left(a\widehat{z}+b\right)}^2}.\hfill \end{array}$$

The boundary points ${\widehat{z}}_0=0$ and ${\widehat{z}}_1=\pi$ can be attained by substituting $z_0$ and $z_1$ into the expression for $\widehat{z}$, respectively. The proof is completed. □

Corollary 2.

The classical second PdHA problem can be achieved by taking $a=1=c$, $b=0.1$, and $u\left(z\right)={\left[(2\mp \sqrt{5})(z-d)\right]}^{2(3\pm \sqrt{5})}$.

Example 1.

Consider a regular Sturm–Liouville problem with Dirichlet boundary conditions

$$-\frac{d}{dz}\left\{{\left[\left(2+\sqrt{5}\right)(z-d)\right]}^{2(3-\sqrt{5})}\frac{dy}{dz}\right\}=\lambda y,y\left(0\right)=0=y\left(z_1\right).$$

where $d=-{\left(0.1\right)}^{2+\sqrt{5}}/(2+\sqrt{5})\approx 1.37\times {10}^{-5}$ and $z_1={(\pi +0.1)}^{2+\sqrt{5}}/(2+\sqrt{5})+d\approx 34.4068$. The BVP in the Liouville normal form is given by

$$-\frac{d^2\widehat{y}}{d{\widehat{z}}^2}+\frac{1}{{\left(\widehat{z}+\frac{1}{10}\right)}^2}\widehat{y}=\lambda \widehat{y},\widehat{y}\left(0\right)=0=\widehat{y}\left(\pi \right).$$

Figure 2 shows the exponential behaviour in z of the coefficient function $u\left(z\right)$ (left panel) and the corresponding invariant function $q\left(\widehat{z}\right)$ in a logarithmic scale (right panel), which is considered in the classical second PdHA problem.

3. Estimating the Smallest Eigenvalue

Arnold et al. proposed a novel technique to approximate the smallest eigenvalue of a Sturm–Liouville problem by means of a localized landscape and effective potential functions [34,35]. Let $\widehat{w}$ be a landscape function, then it satisfies the ODE

$$-\frac{d^2\widehat{w}}{d{\widehat{z}}^2}+q\left(\widehat{z}\right)\widehat{w}=1,$$

(8)

with appropriate boundary conditions. An effective potential $\widehat{W}$ is defined as $\widehat{W}=1/\widehat{w}$, i.e., the reciprocal of the landscape function. Then, the lowest eigenvalue of the Sturm–Liouville problem is given by

$${\lambda }_0\approx \frac{5}{4}{\widehat{W}}_{\min },$$

(9)

where ${\widehat{W}}_{\min }$ is the minimum value of the effective potential on the corresponding interval of the problem.

Without losing any generality, we can rewrite the invariant function q in the following form:

$$q\left(\widehat{z}\right)=\frac{\widehat{c}}{{\left(\widehat{z}+\widehat{b}\right)}^2},\mathrm{where}\widehat{c}=\frac{c}{a^2},\mathrm{and}\widehat{b}=\frac{b}{a}.$$

Hence, the corresponding landscape function for the modified second PdHA problem satisfies the following BVP:

$$-\frac{d^2\widehat{w}}{d{\widehat{z}}^2}+\frac{\widehat{c}}{{\left(\widehat{z}+\widehat{b}\right)}^2}\widehat{w}=1,\widehat{w}\left(0\right)=0=\widehat{w}\left(\pi \right).$$

In what follows, we will pursue generalized conditions at the endpoints, i.e., the third-kind (Robin) boundary conditions, formulated as follows:

$${\alpha }_0\widehat{w}\left(0\right)+{\alpha }_1\frac{d\widehat{w}}{d\widehat{z}}\left(0\right)=w_0,{\beta }_0\widehat{w}\left(\pi \right)+{\beta }_1\frac{d\widehat{w}}{d\widehat{z}}\left(\pi \right)=w_1.$$

Observe that the first and second kinds or Dirichlet and Neumann boundary conditions can be acquired by taking ${\alpha }_1=0={\beta }_1$ and ${\alpha }_0=0={\beta }_0$, respectively.

Let $\widehat{w}\left(\widehat{z}\right)={(\widehat{z}+\widehat{b})}^{\varphi }$, $\varphi \in \mathbb{R}$, be an Ansatz for the homogeneous ODE, and then, upon substitution, we obtain ${\varphi }^2-\varphi -\widehat{c}=0$, which solves

$$\varphi ={\varphi }_{1,2}=\frac{1}{2}\left(1\pm \sqrt{1+4\widehat{c}}\right),$$

where the subscripts 1 and 2 correspond to the positive and negative signs, respectively. Observe that, for a special case of $\widehat{c}=1$, the positive root ${\varphi }_1=\phi$, the well-known golden ratio, and the conjugate root ${\varphi }_2=1-\phi =-1/\phi$, also known as the negative of the silver ratio. Thus, a complementary solution to the homogeneous ODE is given by

$${\widehat{w}}_c\left(\widehat{z}\right)=C_1{(\widehat{z}+\widehat{b})}^{{\varphi }_1}+C_2{(\widehat{z}+\widehat{b})}^{{\varphi }_2},C_1,C_2\in \mathbb{R}.$$

We seek the particular solution using the variation of parameters technique. Writing the Ansatz as

$${\widehat{w}}_p\left(\widehat{z}\right)=r_1\left(\widehat{z}\right){(\widehat{z}+\widehat{b})}^{{\varphi }_1}+r_2\left(\widehat{z}\right){(\widehat{z}+\widehat{b})}^{{\varphi }_2},$$

we obtain

$$\begin{array}{cc}\hfill r_1\left(\widehat{z}\right)& =\int \frac{{(\widehat{z}+\widehat{b})}^{1-{\varphi }_1}}{({\varphi }_2-{\varphi }_1)}d\widehat{z}=\frac{{(\widehat{z}+\widehat{b})}^{2-{\varphi }_1}}{({\varphi }_2-{\varphi }_1)(2-{\varphi }_1)},\hfill \\ \hfill r_2\left(\widehat{z}\right)& =\int \frac{-{(\widehat{z}+\widehat{b})}^{1-{\varphi }_2}}{({\varphi }_2-{\varphi }_1)}d\widehat{z}=\frac{-{(\widehat{z}+\widehat{b})}^{2-{\varphi }_2}}{({\varphi }_2-{\varphi }_1)(2-{\varphi }_2)}.\hfill \end{array}$$

Substituting these to the Ansatz for the particular solution yields

$${\widehat{w}}_p\left(\widehat{z}\right)=-\frac{{(\widehat{z}+\widehat{b})}^2}{(2-{\varphi }_1)(2-{\varphi }_2)}=\frac{{(\widehat{z}+\widehat{b})}^2}{\widehat{c}-2},\widehat{c}\ne 2.$$

Combining together both ${\widehat{w}}_c$ and ${\widehat{w}}_p$, we obtain a general solution for the landscape function:

$$\widehat{w}\left(\widehat{z}\right)=C_1{(\widehat{z}+\widehat{b})}^{{\varphi }_1}+C_2{(\widehat{z}+\widehat{b})}^{{\varphi }_2}+\frac{{(\widehat{z}+\widehat{b})}^2}{\widehat{c}-2}.$$

By imposing the Robin boundary conditions, the constant coefficients $C_1$ and $C_2$ satisfy the following matrix equation:

$$\left[\begin{array}{cc}\hfill b^{{\varphi }_1}\left({\alpha }_0+\frac{{\alpha }_1{\varphi }_1}{b}\right)\hfill & b^{{\varphi }_2}\left({\alpha }_0+\frac{{\alpha }_1{\varphi }_2}{b}\right)\\ \hfill {(\pi +b)}^{{\varphi }_1}\left({\beta }_0+\frac{{\beta }_1{\varphi }_1}{\pi +b}\right)\hfill & {(\pi +b)}^{{\varphi }_2}\left({\beta }_0+\frac{{\beta }_1{\varphi }_2}{\pi +b}\right)\end{array}\right]\left[\begin{array}{c}\hfill C_1\hfill \\ \hfill C_2\hfill \end{array}\right]=\left[\begin{array}{c}\hfill w_0\hfill \\ \hfill w_1\hfill \end{array}\right]+\frac{1}{2-\widehat{c}}\left[\begin{array}{c}\hfill b({\alpha }_{0}b+2{\alpha }_1)\hfill \\ \hfill (\pi +b)\left[{\beta }_0(\pi +b)+2{\beta }_1\right]\hfill \end{array}\right],$$

which solves as $C_1=C_1^R={\widehat{C}}_1/\widehat{C}$ and $C_2=C_2^R={\widehat{C}}_2/\widehat{C}$, where

$$\begin{array}{cc}\hfill {\widehat{C}}_1& =\left[w_0(2-\widehat{c})+b({\alpha }_{0}b+2{\alpha }_1)\right]{(\pi +b)}^{{\varphi }_2}\left({\beta }_0+\frac{{\beta }_1{\varphi }_2}{\pi +b}\right)-b^{{\varphi }_2}\left({\alpha }_0+\frac{{\alpha }_1{\varphi }_2}{b}\right)\left\{w_1(2-\widehat{c})+(\pi +b)\left[{\beta }_0(\pi +b)+2{\beta }_1\right]\right\},\hfill \\ \hfill {\widehat{C}}_2& =b^{{\varphi }_1}\left({\alpha }_0+\frac{{\alpha }_1{\varphi }_1}{b}\right)\left\{w_1(2-\widehat{c})+(\pi +b)\left[{\beta }_0(\pi +b)+2{\beta }_1\right]\right\}-{(\pi +b)}^{{\varphi }_1}\left({\beta }_0+\frac{{\beta }_1{\varphi }_1}{\pi +b}\right)\left[w_0(2-\widehat{c})+b({\alpha }_{0}b+2{\alpha }_1)\right],\hfill \\ \hfill \widehat{C}& =(2-\widehat{c})\left\{b^{{\varphi }_1}{(\pi +b)}^{{\varphi }_2}\left({\alpha }_0+\frac{{\alpha }_1{\varphi }_1}{b}\right)\left({\beta }_0+\frac{{\beta }_1{\varphi }_2}{\pi +b}\right)-b^{{\varphi }_2}{(\pi +b)}^{{\varphi }_1}\left({\alpha }_0+\frac{{\alpha }_1{\varphi }_2}{b}\right)\left({\beta }_0+\frac{{\beta }_1{\varphi }_1}{\pi +b}\right)\right\}.\hfill \end{array}$$

The corresponding coefficients $C_1$ and $C_2$ for Dirichlet boundary conditions

$$\widehat{w}\left(0\right)={\widehat{w}}_0\mathrm{and}\widehat{w}\left(\pi \right)={\widehat{w}}_1,\mathrm{where}{\widehat{w}}_0=\frac{w_0}{{\alpha }_0}\mathrm{and}{\widehat{w}}_1=\frac{w_1}{{\beta }_0},$$

are given as follows:

$$\begin{array}{cc}\hfill C_1=C_1^D& =\frac{\left[{\widehat{w}}_0(2-\widehat{c})+b^2\right]{(\pi +b)}^{{\varphi }_2}-b^{{\varphi }_2}\left[{\widehat{w}}_1(2-\widehat{c})+{(\pi +b)}^2\right]}{(2-\widehat{c})\left[b^{{\varphi }_1}{(\pi +b)}^{{\varphi }_2}-b^{{\varphi }_2}{(\pi +b)}^{{\varphi }_1}\right]},\hfill \\ \hfill C_2=C_2^D& =\frac{b^{{\varphi }_1}\left[{\widehat{w}}_1(2-\widehat{c})+{(\pi +b)}^2\right]-{(\pi +b)}^{{\varphi }_1}\left[{\widehat{w}}_0(2-\widehat{c})+b^2\right]}{(2-\widehat{c})\left[b^{{\varphi }_1}{(\pi +b)}^{{\varphi }_2}-b^{{\varphi }_2}{(\pi +b)}^{{\varphi }_1}\right]}.\hfill \end{array}$$

Similarly, we can also derive the coefficients for Neumann boundary conditions

$$\frac{d\widehat{w}}{d\widehat{z}}\left(0\right)={\tilde{w}}_0\mathrm{and}\frac{d\widehat{w}}{d\widehat{z}}\left(\pi \right)={\tilde{w}}_1,\mathrm{where}{\tilde{w}}_0=\frac{w_0}{{\alpha }_1}\mathrm{and}{\tilde{w}}_1=\frac{w_1}{{\beta }_1}.$$

The corresponding coefficients are given as follows:

$$\begin{array}{cc}\hfill C_1=C_1^N& =\frac{\left[{\tilde{w}}_0(2-\widehat{c})+2b\right]{(\pi +b)}^{{\varphi }_2-1}-b^{{\varphi }_2-1}\left[{\tilde{w}}_1(2-\widehat{c})+2(\pi +b)\right]}{{\varphi }_1(2-\widehat{c})\left[b^{{\varphi }_1-1}{(\pi +b)}^{{\varphi }_2-1}-b^{{\varphi }_2-1}{(\pi +b)}^{{\varphi }_1-1}\right]},\hfill \\ \hfill C_2=C_2^N& =\frac{b^{{\varphi }_1-1}\left[{\tilde{w}}_1(2-\widehat{c})+2(\pi +b)\right]-{(\pi +b)}^{{\varphi }_1-1}\left[{\tilde{w}}_0(2-\widehat{c})+2b\right]}{{\varphi }_2(2-\widehat{c})\left[b^{{\varphi }_1-1}{(\pi +b)}^{{\varphi }_2-1}-b^{{\varphi }_2-1}{(\pi +b)}^{{\varphi }_1-1}\right]}.\hfill \end{array}$$

The minimum value of effective potential, or the maximum value of the landscape function, occurs at $\widehat{z}={\widehat{z}}_c$ such that $d\widehat{w}/d\widehat{z}\left({\widehat{z}}_c\right)=0$, i.e.,

$$C_1{\varphi }_1{({\widehat{z}}_c+\widehat{b})}^{{\varphi }_1-2}+C_2{\varphi }_2{({\widehat{z}}_c+\widehat{b})}^{{\varphi }_2-2}=\frac{2}{(2-\widehat{c})}.$$

We can solve this equation numerically using any standard technique of root-finding algorithms, e.g., the Newton–Raphson method or similar derivative-based techniques.

Figure 3 displays the plots of the landscape function $\widehat{w}\left(\widehat{z}\right)$ and effective potential $\widehat{W}\left(\widehat{z}\right)$ for several values of $\widehat{b}$ but for a chosen fixed value of $\widehat{c}$. The landscape function satisfies the BVP for the modified PdHA problem with homogeneous Dirichlet conditions at both endpoints. We observe that as $\widehat{b}$ increases, its slope at $\widehat{z}=0$ also increases, which causes an increase in its maximum value. Consequently, the minimum values for the potential function decrease as $\widehat{b}$ increases, and an immediate consequence is also a decrease in the eigenvalue estimate. See the left panel of Figure 4.

Figure 5 displays the plots of the landscape function and the corresponding effective potential where the former satisfies the BVP for the modified second PdHA problem with nonhomogeneous Neumann boundary conditions. The slopes are chosen to be positive and negative unity at the left- and right-endpoints, respectively. Both panels depict the curves for a fixed value of $\widehat{c}$ and several values of $\widehat{b}$. Observe that although the slopes of the landscape function at both endpoints remain identical, the initial condition for each landscape function increases as the value of $\widehat{b}$ increases. Consequently, the maximum value also increases accordingly. Conversely, the effective potential exhibits a decreasing minimum value for an increasing value of $\widehat{b}$, which translates to a decreasing value of eigenvalue as $\widehat{b}$ is becoming large. See Figure 4 and Figure 6.

4. Numerical Comparison

In this section, we compare the lowest-order eigenvalues obtained through the landscape function with the ones from numerical simulation. We utilize MATLAB’s finite difference code bvp4c for the latter, accompanied by call functions of the system of the ODE, boundary conditions, and initial guess to solve the corresponding BVP with an unknown parameter.

4.1. Eigenvalue Comparison

In this subsection, we provide comparisons of the eigenvalues of the modified second PdHA problem between the ones obtained from the landscape function and numerical simulations. We also examine further where the prescribed boundary conditions are either Dirichlet or Neumann type.

Figure 4 displays the plots of the lowest-order eigenvalue ${\lambda }_0$ as a function of the parameter $\widehat{b}$. The left panel of Figure 4 depicts a comparison between the eigenvalues obtained from the landscape function and the results from numerical simulations for several values of parameter $\widehat{c}$ and particular Dirichlet boundary conditions. Although the approximation overestimates the numerical calculations, the eigenvalue generally exhibits remarkable quantitative agreement for a wide range of parameter values.

The right panel of Figure 4 presents a comparison of the eigenvalues obtained by an estimate from the landscape function and numerical simulations for particular Neumann boundary conditions with $\widehat{c}=1$. We observe that although the qualitative behavior for both outcomes is similar, i.e., generally decreasing for increasing value of $\widehat{b}$, its quantitative behavior is not. For $\widehat{b}\le 2$, the approximation is far from accurate, whilst for $b>2$, the eigenvalue estimate demonstrates a better approximation with its numerical counterpart despite the former still overestimating the latter.

Figure 6 shows a similar comparison with the right panel of Figure 4 but for a small value of $\widehat{c}$ (left panel) and for the value of $\widehat{c}$ where the landscape function is nearly singular (right panel). A distinct qualitative behavior emerges from this finding. For small $\widehat{c}$, the eigenvalue estimate has a better quantitative agreement with the numerical results. As $\widehat{c}$ increases, it turns out that the accuracy fades away for smaller values of $\widehat{b}$, as we also observed in the right panel of Figure 4. For this particular case, we also noticed that for $\widehat{b}>2$, the estimate improves significantly whereas it was completely inaccurate for $\widehat{b}\le 2$.

4.2. Numerical Eigenfunction

In this subsection, we present the corresponding eigenfunctions of the modified second PdHA BVP for Dirichlet and Neumann conditions. We only focus on the first two lowest-order eigenfunctions since higher-order eigenfunctions can be calculated accordingly by simply modifying the values of the initial guess for the eigenvalue and/or eigenfunction.

Figure 7 displays the eigenfunctions for the modified second PdHA problem with particular Dirichlet boundary conditions. The lowest-order and the second-order eigenfunctions are presented in the left and right panels of Figure 7, respectively. A fixed value of $\widehat{c}$ is prescribed and different curves correspond to various values of $\widehat{b}$. We observe that for increasing values of $\widehat{b}$, both maximum and minimum values of the eigenfunctions decrease. A similar trend also occurs for the eigenvalues, as we discussed in the previous subsection. From these two lowest-order eigenfunctions, even though an increase in the value of $\widehat{b}$ is rather large, e.g., for $\widehat{b}\ge 1$, remarkably, the eigenfunction profiles do not change drastically. The opposite situation occurs when the value of $\widehat{b}$ is relatively small, i.e., $\widehat{b}\le 0.5$, the decrease in the maximum value and quantitative profile of the eigenfunctions is quite dramatic.

Figure 8 illustrates similar eigenfunctions for the modified second PdHA BVP but with specified Neumann conditions. The left and right panels correspond to the first- and second-order eigenfunction profiles for different values of $\widehat{b}$ and a fixed parameter $\widehat{c}$. Compared to the previous case, we recognize that the eigenfunctions exhibit completely different quantitative and qualitative behaviours when we modify the boundary conditions from Dirichlet to Neumann type. While the slopes of the eigenfunction remain identical at both endpoints, the initial values do not. Rather, they rise in tandem with an increase of $\widehat{b}$. The boundary values at the right endpoints, however, do not necessarily follow this trend. An increase in the maximum values of the eigenfunctions is consistent with the decrease in eigenvalues as we confirmed earlier. Different from the previous case, a relatively large change in the parameter value, e.g., from $\widehat{b}=0.5$ to $6.0$, appears to correspond with a substantial change in the eigenfunction profiles quantitatively.

5. Conclusions

We considered a special case of a regular Sturm–Liouville BVP in this article. By transforming the problem into the Liouville normal form, we focus on the second PdHA problem. Although the classical PdHA problem was introduced four decades ago, a modified version of the second PdHA problem appears to be absent from the literature. The exposition above suggests that the modified second PdHA problem might provide additional insights when it comes to the behavior of the eigenvalue and corresponding eigenfunctions with respect to changes in the parameter values.

By incorporating the landscape function and effective potential, we can estimate the lowest-order eigenvalue without solving the corresponding eigenvalue problem. We find that the approach is pretty remarkable even though numerical techniques could deliver the desired output within seconds. In particular, for special cases of the corresponding invariant function that corresponds to the modified second PdHA problem, both landscape and effective potential functions can be solved and expressed analytically.

By prescribing Dirichlet- and Neumann-type boundary conditions, we demonstrated the comparison of eigenvalues between the estimate and numerical simulations. Generally, the eigenvalues obtained from the landscape function overestimate the numerical results quantitatively but they exhibit relatively excellent qualitative agreement. The only exception occurred when we imposed Neumann boundary conditions and selected small values of parameter $\widehat{b}$ but for larger values of $\widehat{c}$.