home

PDF (letter size)

Contents

1 Definitions

Regular Sturm-Liouville ODE is an eigenvalue boundary value ODE. This means the ODE has an eigenvalue in it $\lambda$, where solutions exists only for specific values of eigenvalues. The ODE is$$\frac{d}{dx}\left(p\left(x\right)\frac{dy}{dx}\right)+q\left(x\right)y\left(x\right)+\lambda \sigma \left(x\right)y\left(x\right)=0a<x<b$$

Some books use $-q\left(x\right)y\left(x\right)$ in the above instead of $+q\left(x\right)y\left(x\right)$. The above also can be written as

$$\begin{array}{}\text{(1)}& p{y}^{\prime \prime }+p^{\prime }{y}^{\prime }+(q+\lambda \sigma )y& =0\text{(2)}& {\left(p{y}^{\prime }\right)}^{\prime }+qy& =-\lambda \sigma y\end{array}$$

With the restrictions that $p\left(x\right),q\left(x\right),\sigma \left(x\right)$ are real functions that are continuous everywhere over $a\le x\le b$ and also we need $p\left(x\right)>0,\sigma \left(x\right)>0$. The $\sigma \left(x\right)$ is called the weight function. But this is not all. The boundary conditions must be linear homogeneous, of this form$$\begin{array}{rl}{\beta }_{1}y\left(a\right)+{\beta }_2{y}^{\prime }\left(a\right)& =0\\ {\beta }_{3}y\left(b\right)+{\beta }_4{y}^{\prime }\left(b\right)& =0\end{array}$$

Where ${\beta }_i$ are just real constants. Some of them can be zero but not all. For example, ${\beta }_{1}y\left(a\right)=0,{\beta }_4{y}^{\prime }\left(b\right)=0$ is OK.

Boundary conditions do not have to mixed, but they can be in general. But they must be homogeneous.

Notice that periodic boundary conditions are not allowed. Well, they are allowed, but then the problem is no longer called Sturm-Liouville. The above is just the definition of the equation and its boundary conditions. Below is list of the important properties of this ODE. Each one of these properties have a proof.

The standard way to write S.L. is in the form (2) and not (1) above. This is because the S.L. operator $L$ is defined $$L\left[y\right]=-\lambda \sigma y$$ Where $L\left[y\right]\equiv {\left(p{y}^{\prime }\right)}^{\prime }+qy$. So we break $(q+\lambda \sigma )y$ into two parts, putting the part with the eigenvalue to the right.

$◼$ All eigenvalues ${\lambda }_n$ are real. No complex ${\lambda }_n$.

$◼$ Each eigenvalue ${\lambda }_n$, will have one and only one real eigenfunction ${\varphi }_n\left(x\right)$ associated with it. (in higher dimensions, the eigenvalue problem can have more than one eigenfunction associated with one eigenvalue. For example, heat PDE in rectangle). But we can use Gram-schmidt to make these eigenfunctions orthogonal if we need to). But for exam, just worry about 1D for now.

$◼$ There is smallest $\lambda$, called ${\lambda }_1$ and there are infinite number of eigenvalues. Hence eigenvalues are ordered ${\lambda }_1<{\lambda }_2<{\lambda }_3\cdots$

$◼$ Each eigenfunction ${\varphi }_n\left(x\right)$ will have $n-1$ zeros in $a<x<b$. Note that this does not include the end points. This means, ${\varphi }_3\left(x\right)$ will have two zeros inside the domain. i.e. there are two $x$ locations where ${\varphi }_3\left(x\right)=0$ in $a<x<b$.

$◼$ Eigenfunctions ${\varphi }_n\left(x\right)$ make a complete set of basis. This means any piecewise continuous function $f\left(x\right)$ can be represented as $f\left(x\right)\sim \sum _{n=1}^{\mathrm{\infty }}{a}_n{\varphi }_n\left(x\right)$. This is called the generalized Fourier series.

$◼$ This follows from (5). Each eigenfunction is orthogonal to another eigenfunction, but with the weight in there. This means ${\int }_a^b{\varphi }_n\left(x\right){\varphi }_m\left(x\right)\sigma \left(x\right)dx=0$ if $n\ne m$.

$◼$ Rayleigh quotient relates an eigenvalue to its eigenfunction. Starting with the SL ODE ${\left(p{\varphi }^{\prime }\right)}^{\prime }+q\varphi =-\lambda \sigma \varphi$, then multiplying by$\varphi$ and integrating, we obtain ${\int }_a^b\varphi {\left(p{\varphi }^{\prime }\right)}^{\prime }+q{\varphi }^{2}dx=$ $-{\int }_a^b\lambda \sigma {\varphi }^{2}dx$ and solving for $\lambda$, gives $$\lambda =\frac{{\int }_a^b\varphi {\left(p{\varphi }^{\prime }\right)}^{\prime }+q{\varphi }^{2}dx}{{\int }_a^b{\varphi }^2\sigma \left(x\right)dx}$$ Carrying integration by parts on first part of the integral in numerator, it becomes $$\lambda =\frac{{(-p\varphi {\varphi }^{\prime })}_b^a+{\int }_a^{b}p{\left({\varphi }^{\prime }\right)}^2-q{\varphi }^{2}dx}{{\int }_a^b{\varphi }^2\sigma \left(x\right)dx}$$ But this becomes much simpler when we plug-in the boundary conditions that we must use, making the above $$\lambda =\frac{{\int }_a^{b}p{\left({\varphi }^{\prime }\right)}^2-q{\varphi }^{2}dx}{{\int }_a^b{\varphi }^2\sigma \left(x\right)dx}$$ And if $q=0,p=1$ and $\sigma \left(x\right)=1$, then it becomes $\lambda =\frac{{\int }_a^b{\left({\varphi }^{\prime }\right)}^{2}dx}{{\int }_a^b{\varphi }^{2}dx}.$ Rayleigh quotient is useful to show that $\lambda$ can be positive without solving for $\varphi$ and also used to estimate the value of the minimum $\lambda$ by replacing $\varphi$ by a trial function and actually solving for $\lambda$ using $\lambda =\frac{{\int }_a^b\varphi {\left(p{\varphi }^{\prime }\right)}^{\prime }+q{\varphi }^{2}dx}{{\int }_a^b{\varphi }^2\sigma \left(x\right)dx}$ to obtain a numerical estimate of ${\lambda }_1$.

$◼$ There is symmetry relation. In operator form, let $L\equiv \frac{d}{dx}\left(p\left(x\right)\frac{d}{dx}\right)+q\left(x\right)$, then we have ${\int }_a^{b}uL\left[v\right]-vL\left[u\right]dx=0$ where $u,v$ are any two different eigenfunctions.

$◼$ There is also what is called Lagrange identity, which says $uL\left[v\right]-vL\left[u\right]=\frac{d}{dx}\left(p(u{v}^{\prime }-v{u}^{\prime })\right)$. This comes from simply expanding $L[v],L\left[u\right]$ and simplifying things. Just calculus.

$◼$ Green formula, follows from Lagrange identity, which just gives the integral form ${\int }_a^{b}uL\left[v\right]-vL\left[u\right]dx={\left[p(u{v}^{\prime }-v{u}^{\prime })\right]}_a^b$. But form Sturm-Liouville, we know that ${\int }_a^{b}uL\left[v\right]-vL\left[u\right]dx=0$ (from 8). So this really just says that ${\left[p(u{v}^{\prime }-v{u}^{\prime })\right]}_a^b=0$. But we know this already from boundary conditions. So I am not sure why this is useful for Sturm-Liouville now. Since it is just saying the same thing again.

$◼$ If $p\left(x\right)=0$ at the left or right end (boundaries), then the problem is now called singular Sturm-Liouville. This is actually the important case. In regular S.L., $p\left(x\right)$ must be positive everywhere. We only consider $p\left(x\right)=0$ at the ends, not in the middle or any other place. When this happens, then solutions to the ODE generate special functions, such as Bessel, Legendre, Chebyshev. These special functions are solution to the Singular S.L. ODE, not the regular S.L. ODE. In addition, at the end where $p\left(x\right)=0$, the boundary conditions must be bounded. Not the regular boundary conditions above. For example, if $p\left(a\right)=0$ where $x=a$ is say the left end, the the boundary conditions at the left end must be $y\left(a\right)<\mathrm{\infty },y^{\prime }\left(a\right)<\mathrm{\infty }$. Notice that singular can be at left end, right end or at both at same time. At the end where $p\left(x\right)=0$, the boundary conditions must be bounded type.

$◼$ Symmetry relation $${\int }_a^{b}uL\left[v\right]-vL\left[u\right]dx=0$$ where $u,v$ are are two eigenfunction. But remember, this does NOT mean the integrand is identically zero, or $uL\left[v\right]-vL\left[u\right]=0$. Only when $u,v$ happened to have same eigenvalue we can say that $uL\left[v\right]-vL\left[u\right]=0$. But for unique eigenfunctions only the integral form is zero ${\int }_a^{b}uL\left[v\right]-vL\left[u\right]dx=0.$ This is important difference, used in the proof below, that different eigenfunctions have different eigenvalues (for the scalar SL case). For higher dimensions, we can have more than one eigenfunction for same eigenvalue.

Small note added April 18, 2021. The above is the same as showing the operator $L$ is Hermitian. In bra/ket Dirac notation, an operator $L$ is Hermitian when

$$\begin{array}{rl}⟨u|L|v{⟩}^{\ast }& =⟨v|L|u⟩\\ {(\int u^{\ast }\left(x\right)L\left[v\left(x\right)\right]dx)}^{\ast }& =\int v^{\ast }\left(x\right)L\left[u\left(x\right)\right]dx\\ \int u\left(x\right)L\left[v^{\ast }\left(x\right)\right]dx& =\int v^{\ast }\left(x\right)L\left[u\left(x\right)\right]dx\end{array}$$

In the proofs below, I did not do the complex conjugate, since it is assumed all eigenfunctions $u,v$ are real, and I only wrote that the symmetry relation implies

$$\int u\left(x\right)L\left[v\left(x\right)\right]dx=\int v\left(x\right)L\left[u\left(x\right)\right]dx$$

It should be really called the Hermitian relation to be more accurate. But for real eigenfunctions, it is the same.

2 Proofs

2.1 Proof symmetry of operator

Given regular Sturm-Liouville (RSL) ODE$$\begin{array}{rl}\text{(1)}& {\left(p{y}^{\prime }\right)}^{\prime }+qy& =-\lambda \sigma yp{y}^{\prime \prime }+p^{\prime }{y}^{\prime }+qy& =-\lambda \sigma y\end{array}$$

In operator form $$L\equiv \frac{d}{dx}\left(p\left(x\right)\frac{d}{dx}\right)+q\left(x\right)$$ And (1) becomes$$L[y]=-\lambda \sigma y$$ When solving RSL ode, since it is an eigenvalue ODE with associated boundary conditions, we will get infinite number of non-negative eigenvalues.

For each eigenvalue, there is associated with it one eigenfunction (for 1D case). Looking at any two different eigenfunctions, say $u,v$, then the symmetry relation says the following$${\int }_a^{b}uL\left[v\right]dx={\int }_a^{b}vL\left[u\right]dx$$ Now we will show the above is true. This requires integration by parts two times. We start from the LHS expression and at the end we should end up with the integral on the RHS. Let $I={\int }_a^{b}uL\left[v\right]dx$, then $$\begin{array}{rl}\text{(1)}& I& ={\int }_a^{b}uL\left[v\right]dx& ={\int }_a^{b}u(\frac{d}{dx}\left(p\frac{d}{dx}v\right)+qv)dx\end{array}$$

Where we used $L\left[v\right]=(\frac{d}{dx}\left(p\frac{d}{dx}\right)+q)v=(\frac{d}{dx}\left(p\frac{d}{dx}v\right)+qv)$ in the above. Hence$$\begin{array}{}\text{(2)}& I={\int }_a^{b}u\stackrel{dV}{\overbrace{\frac{d}{dx}\left(p\frac{dv}{dx}\right)}}dx+{\int }_a^{b}qvudx\end{array}$$ Now we will do integration by parts on the first integral. Let $I_2={\int }_a^{b}u\frac{d}{dx}\left(p\frac{dv}{dx}\right)dx$. Using $\int UdV=UV-\int VdU$, and if we let $U=u,dV=\frac{d}{dx}\left(p\frac{dv}{dx}\right)$, then $dU=\frac{du}{dx},V=p\frac{dv}{dx}$. Hence $$I_2={\left(up\frac{dv}{dx}\right)}_a^b-{\int }_a^b\stackrel{dV}{\overbrace{\frac{dv}{dx}}}\stackrel{U}{\overbrace{p\frac{du}{dx}}}dx$$ We now apply integration by parts again to the second integral above. But now let $U=p\frac{du}{dx}$ and $dV=\frac{dv}{dx}$, hence $dU=\frac{d}{dx}\left(p\frac{du}{dx}\right)$ and $V=v$, therefore the above becomes$$I_2={\left(up\frac{dv}{dx}\right)}_a^b-[{\left(p\frac{du}{dx}v\right)}_a^b-{\int }_a^{b}v\frac{d}{dx}\left(p\frac{du}{dx}\right)dx]$$ Substituting the above into (2) gives$$\begin{array}{}\text{(3)}& I={\left(up\frac{dv}{dx}\right)}_a^b-[{\left(p\frac{du}{dx}v\right)}_a^b-{\int }_a^{b}v\frac{d}{dx}\left(p\frac{du}{dx}\right)dx]+{\int }_a^{b}qvudx\end{array}$$ Now comes the part where the boundary conditions are important. In RSL, the boundary conditions are such that all the terms ${\left(up\frac{dv}{dx}\right)}_a^b,{\left(p\frac{du}{dx}v\right)}_a^b$ vanish. This is because the boundary conditions are$$\begin{array}{rl}{\beta }_{1}u\left(a\right)+{\beta }_2{u}^{\prime }\left(a\right)& =0\\ {\beta }_{3}u\left(b\right)+{\beta }_4{u}^{\prime }\left(b\right)& =0\end{array}$$

And$$\begin{array}{rl}{\beta }_{1}v\left(a\right)+{\beta }_2{v}^{\prime }\left(a\right)& =0\\ {\beta }_{3}v\left(b\right)+{\beta }_4{v}^{\prime }\left(b\right)& =0\end{array}$$

So now (3) becomes$$\begin{array}{rl}I& =-[-{\int }_a^{b}v\frac{d}{dx}\left(p\frac{du}{dx}\right)dx]+{\int }_a^{b}qvudx\\ & ={\int }_a^{b}v\frac{d}{dx}\left(p\frac{du}{dx}\right)dx+{\int }_a^{b}qvudx\\ & ={\int }_a^{b}v(\frac{d}{dx}\left(p\frac{du}{dx}\right)+qu)dx\\ & ={\int }_a^{b}vL\left[u\right]dx\end{array}$$

Therefore, we showed that ${\int }_a^{b}uL\left[v\right]dx={\int }_a^{b}vL\left[u\right]dx$. The only thing to watch for here, is which term to make $U$ and which to make $dV$ when making integration by parts. To remember, always start with $uL[v]$ and make $dV$ the term with $\frac{dv}{dx}$ in them during integration parts.

2.2 Proof eigenvalues are real

Assume $\lambda$ is complex. The corresponding eigenfunctions are complex also. Let $\varphi$ be corresponding eigenfunction$$\begin{array}{}\text{(1)}& L\left[\varphi \right]=-\lambda \sigma \varphi \end{array}$$ Taking complex conjugate of both sides, and since $\sigma$ is real, we obtain$$\overline{L\left[\varphi \right]}=-\bar{\lambda }\sigma \bar{\varphi }$$ But $\overline{L\left[\varphi \right]}=L\left[\bar{\varphi }\right]$ since all the coefficients of the ODE are real. The above becomes$$\begin{array}{}\text{(2)}& L\left[\bar{\varphi }\right]=-\bar{\lambda }\sigma \bar{\varphi }\end{array}$$ But by symmetry, we know that$${\int }_a^b\varphi L\left[\bar{\varphi }\right]-\bar{\varphi }L\left[\varphi \right]dx=0$$ Substituting  (1),(2) into the above gives$$\begin{array}{rl}{\int }_a^b\varphi (-\bar{\lambda }\sigma \bar{\varphi })-\bar{\varphi }(-\lambda \sigma \varphi )dx& =0\\ {\int }_a^b-\bar{\lambda }\sigma \varphi \bar{\varphi }+\lambda \sigma \bar{\varphi }\varphi dx& =0\\ {\int }_a^b(\lambda -\bar{\lambda })\left(\sigma \bar{\varphi }\varphi \right)dx& =0\\ (\lambda -\bar{\lambda }){\int }_a^b\sigma \bar{\varphi }\varphi dx& =0\end{array}$$

But $\bar{\varphi }\varphi =\left|\varphi \right|$ which is positive. Also the weight $\sigma \left(x\right)$ is positive by definition. Hence for the above to zero, it must be that $\lambda =\bar{\lambda }$. Which means $\lambda$ is real. QED.

So the main tools to use in this proof: Definition of $L\left[\varphi \right]=-\lambda \sigma \varphi ,$ and symmetry relation and that $\overline{L\left[\varphi \right]}=L\left[\bar{\varphi }\right]$. This might come up in the exam.

2.3 Proof eigenfunctions are unique

Now will show that there is one eigenfunction associated with each eigenvalue (Again, this is for 1D, it is possible to get more than one eigenfunction for same eigenvalue for 2D, as mentioned earlier).  By contradiction, assume that $\lambda$ has two eigenfunctions ${\varphi }_1,{\varphi }_2$ associated with it. Hence$$\begin{array}{rl}L\left[{\varphi }_1\right]& =-\lambda \sigma {\varphi }_1\\ L\left[{\varphi }_2\right]& =-\lambda \sigma {\varphi }_2\end{array}$$

From the first equation, $\lambda =\frac{-L\left[{\varphi }_1\right]}{\sigma {\varphi }_1}$, substituting this into the second equation gives$$\begin{array}{rl}L\left[{\varphi }_2\right]& =\frac{L\left[{\varphi }_1\right]}{{\varphi }_1}{\varphi }_2\\ {\varphi }_{1}L\left[{\varphi }_2\right]-{\varphi }_{2}L\left[{\varphi }_1\right]& =0\end{array}$$

By Lagrange identity, ${\varphi }_{1}L\left[{\varphi }_2\right]-{\varphi }_{2}L\left[{\varphi }_1\right]=\frac{d}{dx}\left(p({\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime })\right)$, hence this means that $$\begin{array}{rl}\frac{d}{dx}\left(p({\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime })\right)& =0\\ p({\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime })& =c_1\end{array}$$

Where $c_1$ is some constant. This is the main difference between the above argument, and between Lagrange identity. This can be confusing. So let me talk more about this. In Lagrange identity, we write $${\varphi }_{1}L\left[{\varphi }_2\right]-{\varphi }_{2}L\left[{\varphi }_1\right]=\frac{d}{dx}\left(p({\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime })\right)$$ And when ${\varphi }_2,{\varphi }_1$ also satisfy the SL boundary condition, only then we say that$$\begin{array}{rl}{\int }_a^b{\varphi }_{1}L\left[{\varphi }_2\right]-{\varphi }_{2}L\left[{\varphi }_1\right]dx& =0\\ {p({\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime })|}_a^b& =0\end{array}$$

But the above is not the same as saying ${\varphi }_{1}L\left[{\varphi }_2\right]-{\varphi }_{2}L\left[{\varphi }_1\right]=0$. This is important to keep in mind. Only the integral form is zero for any two functions with the SL B.C..Now we continue. We showed that $p({\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime })=c$. In SL, this constant is zero due to B.C. Hence$$p({\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime })=0$$ But $p>0$ by definition. Hence ${\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime }=0$ or$$\begin{array}{rl}\frac{d}{dx}\left(\frac{{\varphi }_2}{{\varphi }_1}\right)& =0\\ \frac{{\varphi }_2}{{\varphi }_1}& =c_2\end{array}$$

or ${\varphi }_2=c_2{\varphi }_1$. So the eigenfunctions are linearly dependent. One is just scaled version of the other. But eigenfunction must be linearly independent.

Hence assumption is not valid, and there can not be two linearly independent eigenfunctions for same eigenvalue. Notice also that ${\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime }=0$ is the just the Wronskian. When it is zero, we know the functions are linearly dependent. The important part in the above proof, is that ${\varphi }_1{\varphi }_2^{\prime }-{\varphi }_2{\varphi }_1^{\prime }=0$ only when ${\varphi }_1,{\varphi }_2$ happened to have same eigenvalue.

2.4 Proof eigenfunctions are real

The idea of this proof is to assume the eigenfunction is complex, then show that its real part and its complex part both satisfy the ODE and the boundary conditions. But since they are both use the same eigenvalue, then the real part and the complex part must be linearly dependent. This implies the eigenfunction must be real. (think of Argand diagram)

Assume that $\varphi =U+iV$ is complex eigenfunction with real part $U$ and complex part $V$. Then since$$L\left[\varphi \right]=-\lambda \sigma \varphi$$ The above is just writing the Sturm-Liouville ODE in operator form, where $L$ is the operator as above.  Now we have$$L[U+iV]=-\lambda \sigma (U+iV)$$ By linearity of operator $L$$$L\left[U\right]+iL\left[V\right]=-\lambda \sigma U-i\lambda \sigma V$$ Which implies$$\begin{array}{rl}L\left[U\right]& =-\lambda \sigma U\\ L\left[V\right]& =-\lambda \sigma V\end{array}$$

So we showed that the real and complex part satisfy S.L. Now we need to show they are satisfy S.L. boundary conditions. Since$$\begin{array}{rl}{\beta }_1\varphi \left(a\right)+{\beta }_2{\varphi }^{\prime }\left(a\right)& =0\\ {\beta }_3\varphi \left(b\right)+{\beta }_4{\varphi }^{\prime }\left(b\right)& =0\end{array}$$

Where $x=a,x=b$ are the left and right ends of the domain. Then$$\begin{array}{rl}{\beta }_1(U\left(a\right)+iV\left(a\right))+{\beta }_2(U^{\prime }\left(a\right)+i{V}^{\prime }\left(a\right))& =0\\ {\beta }_3(U\left(b\right)+iV\left(b\right))+{\beta }_4(U^{\prime }\left(b\right)+i{V}^{\prime }\left(b\right))& =0\end{array}$$

Hence$$\begin{array}{rl}{\beta }_{1}U\left(a\right)+{\beta }_2{U}^{\prime }\left(a\right)& =0\\ {\beta }_{1}V\left(a\right)+{\beta }_2{V}^{\prime }\left(a\right)& =0\\ {\beta }_{3}U\left(b\right)+{\beta }_4{U}^{\prime }\left(b\right)& =0\\ {\beta }_{3}V\left(b\right)+{\beta }_4{V}^{\prime }\left(b\right)& =0\end{array}$$

So the above means both $U$ and $V$ satisfy the boundary conditions of S.L. But since both $U,V$ have the same eigenvalue, then they must be linearly dependent, since we know with S.L. each eigenfunction (or one linearly dependent to it) have only one eigenvalue. This means$$V=cU$$ Where $c$ is some constant. In other words, $$\begin{array}{rl}\varphi & =U+iV\\ & =U+icU\\ & =U(1+ic)\\ & =c_{0}U\end{array}$$

Where $c_0$ is new constant. (OK it happens to be complex constant, but it is OK to do so, we always do this trick in other places, if it will make me feel better, I could take the magnitude of the constant). So all what the above says, is that we assumed $\varphi$ to be complex, and found that it is real. So it can’t be complex.

2.5 Proof eigenfunctions are orthogonal with weight

Given two different eigenfunctions ${\varphi }_1,{\varphi }_2$. Hence$$\begin{array}{}\text{(1)}& L\left[{\varphi }_1\right]& =-{\lambda }_1\sigma {\varphi }_1\text{(2)}& L\left[{\varphi }_2\right]& =-{\lambda }_2\sigma {\varphi }_2\end{array}$$

From symmetry integral relation, since these eigenfunctions also satisfy S.L. boundary conditions, we can write$${\int }_a^b{\varphi }_{1}L\left[{\varphi }_2\right]-{\varphi }_{2}L\left[{\varphi }_1\right]dx=0$$ Replacing (1,2) into the above$$\begin{array}{rl}{\int }_a^b{\varphi }_1(-{\lambda }_2\sigma {\varphi }_2)-{\varphi }_2(-{\lambda }_1\sigma {\varphi }_1)dx& =0\\ {\int }_a^b-{\lambda }_2\sigma {\varphi }_1{\varphi }_2+{\lambda }_1\sigma {\varphi }_2{\varphi }_{1}dx& =0\\ {\int }_a^b({\lambda }_1-{\lambda }_2)\left(\sigma {\varphi }_1{\varphi }_2\right)dx& =0\\ ({\lambda }_1-{\lambda }_2){\int }_a^b\sigma {\varphi }_1{\varphi }_{2}dx& =0\end{array}$$

But $({\lambda }_1-{\lambda }_2)\ne 0$ since there are different eigenvalues for different eigenfunctions. Hence $${\int }_a^b\sigma {\varphi }_1\left(x\right){\varphi }_2\left(x\right)dx=0$$ Which means ${\varphi }_{1,}{\varphi }_2$ are orthogonal to each others with weight $\sigma \left(x\right)$

3 Special functions generated from solving singular Sturm-Liouville

When S.L. is singular, meaning $p=0$ at one or both ends, we end with important class of ODE’s, whose solutions are special functions (not $\sin$ or $\cos$) as the case with the regular S.L. Recall that S.L. is$$\begin{array}{rl}\text{(1)}& \frac{d}{dx}\left(p\left(x\right)\frac{dy}{dx}\right)+q\left(x\right)y\left(x\right)+\lambda \sigma \left(x\right)y\left(x\right)& =0p{y}^{\prime \prime }+p^{\prime }{y}^{\prime }+qy& =-\lambda \sigma y\end{array}$$

With the regular S.L., we say that $p\left(x\right)>0$ over the whole domain, including end points. But with singular, this is not the case.  Here are three important S.L. ODE’s that are singular.

Bessel equation: $$\begin{array}{rl}{x}^2{y}^{\prime \prime }+xy+\lambda x^{2}y& =00<x<1\\ x{y}^{\prime \prime }+y+\lambda xy& =0\end{array}$$

Or in standard form$${\left(x{y}^{\prime }\right)}^{\prime }=-\lambda xy$$ Comparing to SL form, then $p\left(x\right)=x,q\left(x\right)=0,\sigma \left(x\right)=x$. At $x=0$, then $p\left(x\right)=0$, which is what makes it singular. (also the weight happens to be zero), but we only care about $p$ being zero or not, at one of the ends. So to check if S.L. is singular or not, just need to check if $p\left(x\right)$ is zero or not at one of the ends.

As mentioned before, when $p=0$ at one of the ends, we can’t use the standard B.C. for the regular S.L. instead, at the end where $p=0$, we must use what is called bounded boundary conditions, which is in this case $y\left(0\right)<\mathrm{\infty },y^{\prime }\left(0\right)<\mathrm{\infty }$. The solution to this ODE will be in terms of Bessel functions. Notice that this happened to be singular, due to the domain starting at $x=0$. If the domain happened to be from $x=0.5$ to say $x=1$, then it is no longer singular S.L. but regular one.

Legendre equation $${\left((1-x^2)y^{\prime }\right)}^{\prime }=-\lambda y0<x<1$$ Or it can be written as$$(1-x^2)y^{\prime \prime }-2x{y}^{\prime }=-\lambda y$$ We see that $p=1-x^2$. And now it happens to be that $p=0$ at $x=1$. So it is singular at the other end compared to Bessel. In this case $q=0,\sigma =1$. Again, the boundary conditions at $x=1$ must now be bounded. i.e. $y\left(1\right)<\mathrm{\infty },y^{\prime }\left(1\right)<\mathrm{\infty }$. On the other end, where $p$ is not zero, we still use the standard boundary conditions  ${\beta }_{1}y\left(0\right)+{\beta }_2{y}^{\prime }\left(0\right)=0$.

Chebyshev equation$${\left(\sqrt{1-x^2}{y}^{\prime }\right)}^{\prime }=-\lambda \sqrt{1-x^2}y-1<x<1$$ Or$$\sqrt{1-x^2}{y}^{\prime \prime }-\frac{1}{2}\frac{2x}{\sqrt{1-x^2}}{y}^{\prime }=-\lambda \sqrt{1-x^2}y$$ So we see that $p=\sqrt{1-x^2},q=0,\sigma =\sqrt{1-x^2}$. So where is $p=0$ here? At $x=-1$ then $p=\sqrt{1-1}=0$ and at $x=-1$ then $p=0$ also. So this is singular at both ends. So we need to use bounded boundary conditions at both ends now to solve this.$$\begin{array}{rl}y(-1)& <\mathrm{\infty }\\ y^{\prime }(-1)& <\mathrm{\infty }\\ y(+1)& <\mathrm{\infty }\\ y^{\prime }(+1)& <\mathrm{\infty }\end{array}$$

The solution to this singular S.L. is given in terms of special function Chebyshev.

4 General notes on S-L

$◼$ Solution to $x{y}^{\prime \prime }+y^{\prime }=-\lambda xy$. Set in S-L form, we get ${\left(p{y}^{\prime }\right)}^{\prime }+qy=-\lambda \sigma y$, then $p=x,q=0,\sigma =x$. For $0<x<b$ with $y\left(b\right)=0,y\left(0\right)<\mathrm{\infty }$. This is singular S-L since $p=x$ is zero at one end. Using Frobenius series, the solution comes out to be $y\left(x\right)=\sum _{n=1}^{\mathrm{\infty }}{c}_n{J}_0\left(\sqrt{{\lambda }_n}x\right)$, where $J_0$ is Bessel function of first kind. $\sqrt{{\lambda }_n}$ are roots of $0=J_0\left(\sqrt{\lambda }b\right)$. (check if this should be $0=J_0\left(\sqrt{\lambda }(b-a)\right)$. $c_n$ is constant, which still needs to be found.

$◼$ We can find approximate solution to S-L ODE for large eigenvalue using WKB so we do not have to solve the ode using series method. Given ${\left(p{y}^{\prime }\right)}^{\prime }+qy=-\lambda \sigma y$ which is standard form S-L. Using physical optics correction, we obtain the solution as$$y\left(x\right)\sim {\left(\sigma \left(x\right)p\left(x\right)\right)}^{\frac{-1}{4}}(c_1\cos \left(\sqrt{\lambda }{\int }_0^x\sqrt{\frac{\sigma \left(t\right)}{p\left(t\right)}}dt\right)+c_1\sin \left(\sqrt{\lambda }{\int }_0^x\sqrt{\frac{\sigma \left(t\right)}{p\left(t\right)}}dt\right))$$ Where $c_1,c_2$ are found from boundary conditions now. The above is valid for large $\lambda$, and is found by first letting ${\epsilon }^2=\frac{1}{\lambda }$ and then assuming $y\left(x\right)=\exp \left(\frac{1}{\epsilon }\sum _{n=0}^{\mathrm{\infty }}{\epsilon }^n{S}_n\left(x\right)\right)$. And working though the WKB method. Remember, WKB only works for linear homogeneous ODE and is used to estimate solution for large $\lambda$ (or small $\epsilon$).

$◼$ in 1D, regular S-L is ${\left(p{y}^{\prime }\right)}^{\prime }+qy=-\lambda \sigma y$ with $0<x<L$, with B.C. given by ${\beta }_{1}y\left(0\right)+{\beta }_2{y}^{\prime }\left(0\right)=0$ and ${\beta }_{1}y\left(L\right)+{\beta }_2{y}^{\prime }\left(L\right)=0,$ and in higher dimensions, this problem becomes $\mathrm{\nabla }\cdot \left(p\mathrm{\nabla }\varphi \right)+q\varphi =-\lambda \sigma \varphi$ with boundary conditions now written as ${\beta }_1\varphi +{\beta }_2(\mathrm{\nabla }\varphi \cdot \hat{n})=0$ on boundary $\mathrm{\Gamma }$. In higher dimensions, $\varphi \equiv \varphi \left(\overline{x}\right)$.

$◼$ Lagrange’s identity in 1D is $$uL\left(v\right)-vL\left(u\right)=\frac{d}{dx}\left(p(u\frac{dv}{dx}-v\frac{du}{dx})\right)$$ where $L$ above is the S-L operator, as in $L\equiv \frac{d}{dx}\left(p\frac{d}{dx}\right)+q$ so that we write $L\left(u\right)=-\lambda \sigma y$.  

$◼$ When it says eigenfunctions are normalized, it should mean that  ${\int }_a^b{\varphi }_n^2\left(x\right)\sigma \left(x\right)dx=1$ where $\sigma \left(x\right)$ is the weight function (comes from SL equation) and $a<x<b$ is the domain. This is for 1D SL. Remember for example, that ${\int }_0^L{\cos }^2\left(\sqrt{{\lambda }_n}x\right)dx=\frac{L}{2}$. where $\sqrt{{\lambda }_n}=\frac{n\pi }{L}$. Here, it is not normalized since the result is not one. The weight here is just one and here ${\varphi }_n\left(x\right)=\cos \left(\frac{n\pi }{L}x\right)$

$◼$ Heat PDE in cylinder is $\frac{\partial T(r,z,\theta ,t)}{\partial t}=k(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{{\partial }^{2}T}{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\theta }^2}).$ For steady state, it becomes, say cylinder has length $L$ $$\begin{array}{rl}0& ={\mathrm{\nabla }}^{2}T(r,z,\theta )=\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{{\partial }^{2}T}{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\theta }^2}\\ \left|T(0,z,\theta )\right|& <\mathrm{\infty }\\ \left|T^{\prime }(0,z,\theta )\right|& <\mathrm{\infty }\\ T(r,0)& =f\left(r\right)\\ R(r,L)& =g\left(r\right)\end{array}$$

And for axis symmetric (no $\theta$ dependency), it becomes $$0={\mathrm{\nabla }}^{2}T(r,z)=\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{{\partial }^{2}T}{\partial z^2}$$ The solution is found by separation of variables. For the above is in terms of Bessel function order zero$$T(r,z)=\sum _{n=1}^{\mathrm{\infty }}{J}_0\left(\sqrt{{\lambda }_0}r\right)(A_n{e}^{\sqrt{{\lambda }_n}z}+B_n{e}^{\sqrt{{\lambda }_n}z})$$ which is found using series method $T=\sum _{n=0}^{\mathrm{\infty }}{a}_n{r}^{n+\alpha }$. The eigenvalues ${\lambda }_n$ are zeros of $0=J_0\left(\sqrt{\lambda }{r}_0\right)$ where $r_0$ is disk radius.

$◼$ Heat PDE in 2D polar is $\frac{\partial T(r,\theta ,t)}{\partial t}=k(\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\theta }^2})$ and for steady state it becomes $0=\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\theta }^2}$. After separating, the $r$ ODE becomes Euler ODE. Use guess $R\left(r\right)=r^p$ for solution. The solution will be$$T(r,z)=A_0+\sum _{n=1}^{\mathrm{\infty }}{A}_n\cos \left(n\theta \right)r^n+\sum _{n=1}^{\mathrm{\infty }}{B}_n\sin \left(n\theta \right)r^n$$ $◼$ To change $p\left(x\right)y^{\prime \prime }+q\left(x\right)y^{\prime }+r\left(x\right)y\left(x\right)=0$ to S.L. form, multiply the ODE by $\mu =\frac{1}{p\left(x\right)}\exp (\int \frac{q\left(x\right)}{p\left(x\right)}dx)$

$◼$ In using method of eigenfunction expansion to solve nonhomogeneous B.C. PDE, the eigenfunction used ${\varphi }_n\left(x\right)$ are the ones from solving the PDE with homogeneous B.C. For example, if given $\frac{\partial u}{\partial t}=k\frac{{\partial }^{2}u}{\partial x^2}$ with nonhomogeneous B.C. , say $u(0,t)=A,u(L,t)=B$, and we want to use $u(x,t)\sim \sum _{n=1}^{\mathrm{\infty }}{b}_n\left(t\right){\varphi }_n\left(x\right)$ on the nonhomogeneous B.C. PDE, then those ${\varphi }_n\left(x\right)$ are from the corresponding homogeneous B.C. PDE. They should be something like $\sin \left(\frac{n\pi }{L}x\right)$ and so on.

See the table at the top. That is why we write $\sim$ above and not $=$ and remember not to do term-by-term differentiation in $x$ when this is the case. We can do term-by-term differentiation only of the PDE itself also has homogeneous B.C. , even if it had a source term there also. The point is, ${\varphi }_n\left(x\right)$ always come from solving the homogeneous B.C. PDE. (This normally means solving a Sturm-Liouville ODE).

$◼$ I found new relation for eigenvalue $\lambda$, but it is probably not useful as the one the book has. Here how to derive it
Given S-L $$\begin{array}{}\text{(1)}& \frac{d}{dx}\left(p\frac{dy}{dx}\right)+qy=-\lambda \sigma y\end{array}$$ Let ${\varphi }_n\left(x\right)$ be the eigenfunction associated with eigenvalue ${\lambda }_n$. Since eigenfunctions satisfy the ODE itself, then we can write, for any arbitrary eigenfunction (subscript removed for clarity in what follows)$$\begin{array}{rl}\frac{d}{dx}\left(p{\varphi }^{\prime }\right)+q\varphi & =-\lambda \sigma \varphi \\ p^{\prime }{\varphi }^{\prime }+p{\varphi }^{\prime \prime }+q\varphi & =-\lambda \sigma \varphi \end{array}$$

Integrating both sides$$\begin{array}{}\text{(2)}& {\int }_a^b{p}^{\prime }{\varphi }^{\prime }dx+{\int }_a^{b}p{\varphi }^{\prime \prime }dx+{\int }_a^{b}q\varphi dx=-\lambda {\int }_a^b\sigma \varphi dx\end{array}$$ Looking at ${\int }_a^{b}p{\varphi }^{\prime \prime }dx$. Integrating by part. Let $u=p,dv={\varphi }^{\prime \prime }\to du=p^{\prime },v={\varphi }^{\prime }$, hence$$\begin{array}{rl}{\int }_a^{b}p{\varphi }^{\prime \prime }dx& ={\left[uv\right]}_a^b-{\int }_a^{b}vdu\\ \text{(3)}& & ={\left[p{\varphi }^{\prime }\right]}_a^b-{\int }_a^b{p}^{\prime }{\varphi }^{\prime }dx\end{array}$$

Substituting (3) into (2) gives$$\begin{array}{rl}{\int }_a^b{p}^{\prime }{\varphi }^{\prime }dx+{\left[p{\varphi }^{\prime }\right]}_0^L-{\int }_a^b{p}^{\prime }{\varphi }^{\prime }dx+{\int }_a^{b}q\varphi dx& =-\lambda {\int }_a^b\sigma \varphi dx\\ {\left[p{\varphi }^{\prime }\right]}_0^L+{\int }_a^{b}q\varphi dx& =-\lambda {\int }_a^b\sigma \varphi dx\end{array}$$

Hence$$\lambda =-\frac{{\left[p{\varphi }^{\prime }\right]}_a^b+{\int }_a^{b}q\varphi dx}{{\int }_a^b\sigma \varphi dx}$$ compare to the one in the book $\lambda =\frac{{(-p\varphi {\varphi }^{\prime })}_b^a+{\int }_a^{b}p{\left({\varphi }^{\prime }\right)}^2-q{\varphi }^{2}dx}{{\int }_a^b{\varphi }^2\sigma \left(x\right)dx}$.

$◼$ This is about using eigenfunction expansion to approximate $f\left(x\right)$. Given the eigenfunctions ${\mathrm{\Phi }}_n\left(x\right)$ found by solving SL $L\left[y\right]=\lambda ry$ problem, then these eigenfunctions are called complete in the mean sense, if $\sum _{n=1}^{\mathrm{\infty }}{b}_n{\mathrm{\Phi }}_n\left(x\right)$ converges to the mean of $f\left(x\right)$ at each point in $a\le x\le b$.

This requires only that $f\left(x\right)$ be square integrable over the domain. We can also say that the eigenfunctions are complete, in pointwise sense, if $\sum _{n=1}^{\mathrm{\infty }}{b}_n{\mathrm{\Phi }}_n\left(x\right)$ converges to $f\left(x\right)$ at each point in $a\le x\le b$. This however requires that $f\left(x\right)$ be continuous in $a\le x\le b$ and $f^{\prime }\left(x\right)$ be piecewise continuous.

So clearly the first case is less strict than the second. This means if $f\left(x\right)$ is only square integrable (and not necessarily continuos), then convergence to the mean of $f\left(x\right)$ at each point can be obtained. To obtain the more strict point wise convergence, more restrictions are needed on $f\left(x\right)$ as mentioned above. Clearly point wise convergence implies convergence in the mean. But not the other way around.

$◼$ For Sturm-Liouville, write it as $p{y}^{\prime \prime }+p^{\prime }{y}^{\prime }+(q+\lambda \sigma )y=0$. Now it is easy to compare to. For example, $x^{\prime \prime }+\lambda x=0$, then we see that $p=1,q=0,\sigma =1$. And for $x{y}^{\prime \prime }+x^{\prime }+\lambda xy=0$ then $p=1,q=0,\sigma =x$. Do not divide by $x$ first to re-write $x{y}^{\prime \prime }+x^{\prime }+\lambda xy=0$ as $y^{\prime \prime }+\frac{1}{x}{x}^{\prime }+\lambda y=0$ and then compare, as this makes it hard to compare to standard form and will make error.  Books write Sturm-Liouville as $\frac{d}{dx}\left(p{y}^{\prime }\right)+qy+\lambda \sigma y=0$ but I find this confusing to compare to this way. Better to expand it to $p{y}^{\prime \prime }+p^{\prime }{y}^{\prime }+(q+\lambda \sigma )y=0$ first so now it is in a more familiar form and can read it out more directly.

5 Methods of solutions to some Sturm Liouville problems

Problem 1 If the problem gives S-L equations, and asks to find estimate on the smallest eigenvalue, then use Rayleigh quotient for $\lambda$. And write ${\lambda }_{min}={\lambda }_1\le \frac{{(-p\varphi {\varphi }^{\prime })}_b^a+{\int }_a^{b}p{\left({\varphi }^{\prime }\right)}^2-q{\varphi }^{2}dx}{{\int }_a^b{\varphi }^2\sigma \left(x\right)dx}$. Now we do not need to solve the SL to find solution $\varphi$, this is the whole point. Everything in RHS is given, except for, of course the solution $\varphi$.

Here comes the main idea: Come up with any trial ${\varphi }_{trial}$ and use it in place of $\varphi$. This trial function just needs to satisfy the boundary conditions, which is also given. Then all what we need to do is just evaluate the integral. Pick the simplest ${\varphi }_{trial}$ function of $x$ which satisfies boundary conditions. All other terms $p,q,\sigma$ we can read from the given problem.

At the end, we should get a numerical value for the integral. This is the upper limit of the lowest eigenvalue ${\lambda }_1$

Problem 2 We are given SL problem with boundary conditions, and asked to show that $\lambda \ge 0$ without solving the ODE to find $\varphi$:

Use Rayleigh quotient and argue that the denominator can’t be zero (else eigenfunction is zero, which is not possible) and it also can’t be negative. Argue about the numerator not being negative. Do not solve the ODE! this is the whole point of using Rayleigh quotient.

Problem We are given an SL problem with boundary conditions and asked to estimate large $\lambda$ and corresponding eigenfunction.

This is different from being asked to estimate the smallest eigenvalue, where we use Rayleigh quotient and trial function. In this one, we instead use WKB. For 1D, just use what is called the physical optics correction method, given by $\varphi \left(x\right)\approx {\left(\sigma p\right)}^{\frac{-1}{4}}(c_1\cos \left(\sqrt{\lambda }{\int }_0^x\sqrt{\frac{\sigma \left(t\right)}{p\left(t\right)}}dt\right)+c_2\sin \left(\sqrt{\lambda }{\int }_0^x\sqrt{\frac{\sigma \left(t\right)}{p\left(t\right)}}dt\right))$. Where $\sigma ,p$ in this are all known functions (from the problem itself). Notice that $q$ is not there.