Sturm–Liouville theory

Template:Short description In mathematics and its applications, a Sturm–Liouville problem is a second-order linear ordinary differential equation of the form <math display="block">

\frac{\mathrm{d}}{\mathrm{d}x} \left[p(x) \frac{\mathrm{d}y}{\mathrm{d}x}\right] + q(x)y = -\lambda w(x) y

</math> for given functions <math>p(x)</math>, <math>q(x)</math> and <math>w(x)</math>, together with some boundary conditions at extreme values of <math>x</math>. The goals of a given Sturm–Liouville problem are:

To find the Template:Mvar for which there exists a non-trivial solution to the problem. Such values Template:Mvar are called the eigenvalues of the problem.
For each eigenvalue Template:Mvar, to find the corresponding solution <math>y = y(x)</math> of the problem. Such functions <math>y</math> are called the eigenfunctions associated to each Template:Mvar.

Sturm–Liouville theory is the general study of Sturm–Liouville problems. In particular, for a "regular" Sturm–Liouville problem, it can be shown that there are an infinite number of eigenvalues each with a unique eigenfunction, and that these eigenfunctions form an orthonormal basis of a certain Hilbert space of functions.

This theory is important in applied mathematics, where Sturm–Liouville problems occur very frequently, particularly when dealing with separable linear partial differential equations. For example, in quantum mechanics, the one-dimensional time-independent Schrödinger equation is a Sturm–Liouville problem.

Sturm–Liouville theory is named after Jacques Charles François Sturm (1803–1855) and Joseph Liouville (1809–1882), who developed the theory.

Main resultsEdit

The main results in Sturm–Liouville theory apply to a Sturm–Liouville problem Template:NumBlk2{\mathrm{d}x} \left[p(x)\frac{\mathrm{d}y}{\mathrm{d}x}\right] + q(x)y = -\lambda\, w(x)y </math>|1}} on a finite interval <math>[a, b]</math> that is "regular". The problem is said to be regular if:

the coefficient functions <math>p, q, w</math> and the derivative <math>p'</math> are all continuous on <math>[a, b]</math>;
<math>p(x) > 0</math> and <math>w(x) > 0</math> for all <math>x \in [a, b]</math>;
the problem has separated boundary conditions of the form

Template:NumBlk2

The function <math>w = w(x)</math>, sometimes denoted <math>r = r(x)</math>, is called the weight or density function.

The goals of a Sturm–Liouville problem are:

to find the eigenvalues: those Template:Mvar for which there exists a non-trivial solution;
for each eigenvalue Template:Mvar, to find the corresponding eigenfunction <math>y = y(x)</math>.

For a regular Sturm–Liouville problem, a function <math>y = y(x)</math> is called a solution if it is continuously differentiable and satisfies the equation (Template:EquationNote) at every <math>x \in (a, b)</math>. In the case of more general <math>p, q, w</math>, the solutions must be understood in a weak sense.

The terms eigenvalue and eigenvector are used because the solutions correspond to the eigenvalues and eigenfunctions of a Hermitian differential operator in an appropriate Hilbert space of functions with inner product defined using the weight function. Sturm–Liouville theory studies the existence and asymptotic behavior of the eigenvalues, the corresponding qualitative theory of the eigenfunctions and their completeness in the function space.

The main result of Sturm–Liouville theory states that, for any regular Sturm–Liouville problem:

The eigenvalues <math>\lambda_1, \lambda_2, \dots</math> are real and can be numbered so that <math>\lambda_1 < \lambda_2 < \cdots < \lambda_n < \cdots \to \infty.</math>
Corresponding to each eigenvalue <math>\lambda_n</math> is a unique (up to constant multiple) eigenfunction <math>y_n = y_n(x)</math> with exactly <math>n - 1</math> zeros in <math>[a, b]</math>, called the Template:Mvarth fundamental solution.
The normalized eigenfunctions <math>y_n</math> form an orthonormal basis under the w-weighted inner product in the Hilbert space <math>L^2\big([a, b], w(x)\,\mathrm{d}x\big)</math>; that is, <math display="block">

\langle y_n, y_m\rangle = \int_a^b y_n(x) y_m(x) w(x)\,\mathrm{d}x = \delta_{nm},

</math> where <math>\delta_{nm}</math> is the Kronecker delta.

Reduction to Sturm–Liouville formEdit

The differential equation (Template:EquationNote) is said to be in Sturm–Liouville form or self-adjoint form. All second-order linear homogenous ordinary differential equations can be recast in the form on the left-hand side of (Template:EquationNote) by multiplying both sides of the equation by an appropriate integrating factor (although the same is not true of second-order partial differential equations, or if Template:Mvar is a vector). Some examples are below.

Bessel equationEdit

<math display="block">x^2y + xy' + \left(x^2-\nu^2\right)y = 0</math> which can be written in Sturm–Liouville form (first by dividing through by Template:Mvar, then by collapsing the first two terms on the left into one term) as <math display="block">\left(xy'\right)'+ \left (x-\frac{\nu^2} x \right )y=0.</math>

Legendre equationEdit

<math display="block">\left(1-x^2\right)y-2xy'+\nu(\nu+1)y=0</math> which can be put into Sturm–Liouville form, since Template:Math, so the Legendre equation is equivalent to <math display="block">\left (\left(1-x^2\right)y' \right )'+\nu(\nu+1)y=0</math>

Example using an integrating factorEdit

Divide throughout by Template:Math: <math display="block">y-\frac{1}{x^2}y'+\frac{2}{x^3}y=0</math>

Multiplying throughout by an integrating factor of <math display="block">\mu(x) =\exp\left(\int -\frac{dx}{x^2}\right)=e^{{1}/{x}},</math> gives <math display="block">e^{{1}/{x}}y-\frac{e^{{1}/{x}}}{x^2} y'+ \frac{2 e^{{1}/{x}}}{x^3} y = 0</math> which can be put into Sturm–Liouville form since <math display="block">\frac{d}{dx} e^{{1}/{x}} = -\frac{e^{{1}/{x}}}{x^2} </math> so the differential equation is equivalent to <math display="block">\left (e^{{1}/{x}}y' \right )'+\frac{2 e^{{1}/{x}}}{x^3} y = 0.</math>

Integrating factor for general second-order homogenous equationEdit

Multiplying through by the integrating factor <math display="block">\mu(x) = \frac 1 {P(x)} \exp \left(\int \frac{Q(x)}{P(x)} \, dx\right),</math> and then collecting gives the Sturm–Liouville form: <math display="block">\frac{d}{dx} \left(\mu(x)P(x)y'\right) + \mu(x)R(x)y = 0,</math> or, explicitly: <math display="block">\frac{d}{dx} \left(\exp\left (\int \frac{Q(x)}{P(x)} \,dx\right)y' \right )+\frac{R(x)}{P(x)} \exp \left(\int \frac{Q(x)}{P(x)}\, dx\right) y = 0.</math>

Sturm–Liouville equations as self-adjoint differential operatorsEdit

The mapping defined by: <math display="block">Lu = -\frac{1}{w(x)} \left(\frac{d}{dx}\left[p(x)\,\frac{du}{dx}\right]+q(x)u \right)</math> can be viewed as a linear operator Template:Mvar mapping a function Template:Mvar to another function Template:Mvar, and it can be studied in the context of functional analysis. In fact, equation (Template:EquationNote) can be written as <math display="block">Lu = \lambda u.</math>

This is precisely the eigenvalue problem; that is, one seeks eigenvalues Template:Math and the corresponding eigenvectors Template:Math of the Template:Mvar operator. The proper setting for this problem is the Hilbert space <math>L^2([a,b],w(x)\,dx)</math> with scalar product <math display="block"> \langle f, g\rangle = \int_a^b \overline{f(x)} g(x)w(x)\, dx.</math>

In this space Template:Mvar is defined on sufficiently smooth functions which satisfy the above regular boundary conditions. Moreover, Template:Mvar is a self-adjoint operator: <math display="block"> \langle L f, g \rangle = \langle f, L g \rangle .</math>

This can be seen formally by using integration by parts twice, where the boundary terms vanish by virtue of the boundary conditions. It then follows that the eigenvalues of a Sturm–Liouville operator are real and that eigenfunctions of Template:Mvar corresponding to different eigenvalues are orthogonal. However, this operator is unbounded and hence existence of an orthonormal basis of eigenfunctions is not evident. To overcome this problem, one looks at the resolvent <math display="block">\left (L - z\right)^{-1}, \qquad z \in \Reals,</math> where Template:Mvar is not an eigenvalue. Then, computing the resolvent amounts to solving a nonhomogeneous equation, which can be done using the variation of parameters formula. This shows that the resolvent is an integral operator with a continuous symmetric kernel (the Green's function of the problem). As a consequence of the Arzelà–Ascoli theorem, this integral operator is compact and existence of a sequence of eigenvalues Template:Mvar which converge to 0 and eigenfunctions which form an orthonormal basis follows from the spectral theorem for compact operators. Finally, note that <math display="block">\left(L-z\right)^{-1} u = \alpha u, \qquad L u = \left(z+\alpha^{-1}\right) u,</math> are equivalent, so we may take <math>\lambda = z+\alpha^{-1}</math> with the same eigenfunctions.

If the interval is unbounded, or if the coefficients have singularities at the boundary points, one calls Template:Mvar singular. In this case, the spectrum no longer consists of eigenvalues alone and can contain a continuous component. There is still an associated eigenfunction expansion (similar to Fourier series versus Fourier transform). This is important in quantum mechanics, since the one-dimensional time-independent Schrödinger equation is a special case of a Sturm–Liouville equation.

Application to inhomogeneous second-order boundary value problemsEdit

Consider a general inhomogeneous second-order linear differential equation <math display="block">P(x)y + Q(x)y' +R(x)y = f(x)</math> for given functions <math>P(x), Q(x), R(x),f(x)</math>. As before, this can be reduced to the Sturm–Liouville form <math>Ly = f</math>: writing a general Sturm–Liouville operator as: <math display="block">Lu = \frac{p}{w(x)}u + \frac{p'}{w(x)}u' + \frac{q}{w(x)}u,</math> one solves the system: <math display="block">p = Pw,\quad p' = Qw,\quad q = Rw.</math>

It suffices to solve the first two equations, which amounts to solving Template:Math, or <math display="block"> w' = \frac{Q-P'}{P}w:= \alpha w.</math>

A solution is:

<math display="block">w = \exp\left(\int\alpha \, dx\right), \quad p = P \exp\left(\int\alpha \, dx\right), \quad q = R \exp\left(\int\alpha \, dx\right).</math>

Given this transformation, one is left to solve: <math display="block">Ly = f.</math>

In general, if initial conditions at some point are specified, for example Template:Math and Template:Math, a second order differential equation can be solved using ordinary methods and the Picard–Lindelöf theorem ensures that the differential equation has a unique solution in a neighbourhood of the point where the initial conditions have been specified.

But if in place of specifying initial values at a single point, it is desired to specify values at two different points (so-called boundary values), e.g. Template:Math and Template:Math, the problem turns out to be much more difficult. Notice that by adding a suitable known differentiable function to Template:Mvar, whose values at Template:Mvar and Template:Mvar satisfy the desired boundary conditions, and injecting inside the proposed differential equation, it can be assumed without loss of generality that the boundary conditions are of the form Template:Math and Template:Math.

Here, the Sturm–Liouville theory comes in play: indeed, a large class of functions Template:Mvar can be expanded in terms of a series of orthonormal eigenfunctions Template:Mvar of the associated Liouville operator with corresponding eigenvalues Template:Mvar: <math display="block">f(x) = \sum_i \alpha_i u_i(x), \quad \alpha_i \in {\mathbb R}.</math>

Then a solution to the proposed equation is evidently: <math display="block"> y = \sum_i \frac{\alpha_i}{\lambda_i} u_i.</math>

This solution will be valid only over the open interval Template:Math, and may fail at the boundaries.

Example: Fourier seriesEdit

Consider the Sturm–Liouville problem: Template:NumBlk2 for the unknowns are Template:Mvar and Template:Math. For boundary conditions, we take for example: <math display="block"> u(0) = u(\pi) = 0.</math>

Observe that if Template:Mvar is any integer, then the function <math display="block"> u_k(x) = \sin kx</math> is a solution with eigenvalue Template:Math. We know that the solutions of a Sturm–Liouville problem form an orthogonal basis, and we know from Fourier series that this set of sinusoidal functions is an orthogonal basis. Since orthogonal bases are always maximal (by definition) we conclude that the Sturm–Liouville problem in this case has no other eigenvectors.

Given the preceding, let us now solve the inhomogeneous problem <math display="block">L y =x, \qquad x\in(0,\pi)</math> with the same boundary conditions <math>y(0) = y(\pi) = 0</math>. In this case, we must expand Template:Math as a Fourier series. The reader may check, either by integrating Template:Math or by consulting a table of Fourier transforms, that we thus obtain <math display="block">L y = \sum_{k=1}^\infty -2\frac{\left(-1\right)^k} k \sin kx.</math>

This particular Fourier series is troublesome because of its poor convergence properties. It is not clear a priori whether the series converges pointwise. Because of Fourier analysis, since the Fourier coefficients are "square-summable", the Fourier series converges in Template:Math which is all we need for this particular theory to function. We mention for the interested reader that in this case we may rely on a result which says that Fourier series converge at every point of differentiability, and at jump points (the function x, considered as a periodic function, has a jump at Template:Pi) converges to the average of the left and right limits (see convergence of Fourier series).

Therefore, by using formula (Template:EquationNote), we obtain the solution: <math display="block">y=\sum_{k=1}^\infty 2\frac{(-1)^k}{k^3}\sin kx= \tfrac 1 6 (x^3 -\pi^2 x).</math>

In this case, we could have found the answer using antidifferentiation, but this is no longer useful in most cases when the differential equation is in many variables.

Application to partial differential equationsEdit

Normal modesEdit

Certain partial differential equations can be solved with the help of Sturm–Liouville theory. Suppose we are interested in the vibrational modes of a thin membrane, held in a rectangular frame, Template:Math, Template:Math. The equation of motion for the vertical membrane's displacement, Template:Math is given by the wave equation: <math display="block">\frac{\partial^2W}{\partial x^2}+\frac{\partial^2W}{\partial y^2} = \frac 1 {c^2} \frac{\partial^2W}{\partial t^2}.</math>

The method of separation of variables suggests looking first for solutions of the simple form Template:Math. For such a function Template:Mvar the partial differential equation becomes Template:Math. Since the three terms of this equation are functions of Template:Math separately, they must be constants. For example, the first term gives Template:Math for a constant Template:Mvar. The boundary conditions ("held in a rectangular frame") are Template:Math when Template:Math, Template:Math or Template:Math, Template:Math and define the simplest possible Sturm–Liouville eigenvalue problems as in the example, yielding the "normal mode solutions" for Template:Mvar with harmonic time dependence, <math display="block">W_{mn}(x,y,t) = A_{mn} \sin\left(\frac{m\pi x}{L_1}\right) \sin\left(\frac{n\pi y}{L_2}\right)\cos\left(\omega_{mn}t\right)</math> where Template:Mvar and Template:Mvar are non-zero integers, Template:Mvar are arbitrary constants, and <math display="block">\omega^2_{mn} = c^2 \left(\frac{m^2\pi^2}{L_1^2}+\frac{n^2\pi^2}{L_2^2}\right).</math>

The functions Template:Mvar form a basis for the Hilbert space of (generalized) solutions of the wave equation; that is, an arbitrary solution Template:Mvar can be decomposed into a sum of these modes, which vibrate at their individual frequencies Template:Mvar. This representation may require a convergent infinite sum.

Second-order linear equationEdit

Consider a linear second-order differential equation in one spatial dimension and first-order in time of the form: <math display="block"> f(x) \frac{\partial^2 u}{\partial x^2} + g(x) \frac{\partial u}{\partial x} + h(x) u= \frac{\partial u}{\partial t} + k(t) u,</math> <math display="block"> u(a,t)=u(b,t)=0, \qquad u(x,0)=s(x). </math>

Separating variables, we assume that <math display="block">u(x,t) = X(x) T(t). </math> Then our above partial differential equation may be written as: <math display="block">\frac{\hat{L} X(x)}{X(x)} = \frac{\hat{M} T(t)}{T(t)}</math> where <math display="block"> \hat{L}=f(x) \frac{d^2}{dx^2}+g(x) \frac{d}{dx}+h(x), \qquad \hat{M} = \frac{d}{dt} + k(t).</math>

Since, by definition, Template:Mvar and Template:Math are independent of time Template:Mvar and Template:Mvar and Template:Math are independent of position Template:Mvar, then both sides of the above equation must be equal to a constant: <math display="block"> \hat{L} X(x) =\lambda X(x),\qquad X(a)=X(b)=0,\qquad \hat{M} T(t) =\lambda T(t).</math>

The first of these equations must be solved as a Sturm–Liouville problem in terms of the eigenfunctions Template:Math and eigenvalues Template:Math. The second of these equations can be analytically solved once the eigenvalues are known.

<math display="block"> \frac{d}{dt} T_n (t)= \bigl(\lambda_n -k(t)\bigr) T_n (t) </math> <math display="block"> T_n (t) = a_n \exp \left(\lambda_n t -\int_0^t k(\tau) \, d\tau\right) </math> <math display="block"> u(x,t) =\sum_n a_n X_n (x) \exp \left(\lambda_n t -\int_0^t k(\tau) \, d\tau\right) </math> <math display="block"> a_n =\frac{\bigl\langle X_n (x), s(x)\bigr\rangle}{\bigl\langle X_n(x),X_n (x)\bigr\rangle}</math>

where <math display="block"> \bigl\langle y(x),z(x)\bigr\rangle = \int_a^b y(x) z(x) w(x) \, dx, </math> <math display="block"> w(x)= \frac{\exp \left(\int \frac{g(x)}{f(x)} \, dx\right)}{f(x)}. </math>

Representation of solutions and numerical calculationEdit

The Sturm–Liouville differential equation (Template:EquationNote) with boundary conditions may be solved analytically, which can be exact or provide an approximation, by the Rayleigh–Ritz method, or by the matrix-variational method of Gerck et al.<ref>Ed Gerck, A. B. d'Oliveira, H. F. de Carvalho. "Heavy baryons as bound states of three quarks." Lettere al Nuovo Cimento 38(1):27–32, Sep 1983.</ref><ref>Augusto B. d'Oliveira, Ed Gerck, Jason A. C. Gallas. "Solution of the Schrödinger equation for bound states in closed form." Physical Review A, 26:1(1), June 1982.</ref><ref>Robert F. O'Connell, Jason A. C. Gallas, Ed Gerck. "Scaling Laws for Rydberg Atoms in Magnetic Fields." Physical Review Letters 50(5):324–327, January 1983.</ref>

Numerically, a variety of methods are also available. In difficult cases, one may need to carry out the intermediate calculations to several hundred decimal places of accuracy in order to obtain the eigenvalues correctly to a few decimal places.

Shooting methods<ref>Template:Cite book</ref><ref>Template:Cite journal</ref>
Finite difference method
Spectral parameter power series method<ref name="KP">Template:Cite journal</ref>

Shooting methodsEdit

Shooting methods proceed by guessing a value of Template:Mvar, solving an initial value problem defined by the boundary conditions at one endpoint, say, Template:Mvar, of the interval Template:Closed-closed, comparing the value this solution takes at the other endpoint Template:Mvar with the other desired boundary condition, and finally increasing or decreasing Template:Mvar as necessary to correct the original value. This strategy is not applicable for locating complex eigenvalues.Template:Clarify

Spectral parameter power series methodEdit

The spectral parameter power series (SPPS) method makes use of a generalization of the following fact about homogeneous second-order linear ordinary differential equations: if Template:Mvar is a solution of equation (Template:EquationNote) that does not vanish at any point of Template:Closed-closed, then the function <math display="block"> y(x) \int_a^x \frac{dt}{p(t)y(t)^2} </math> is a solution of the same equation and is linearly independent from Template:Mvar. Further, all solutions are linear combinations of these two solutions. In the SPPS algorithm, one must begin with an arbitrary value Template:Math (often Template:Math; it does not need to be an eigenvalue) and any solution Template:Math of (Template:EquationNote) with Template:Math which does not vanish on Template:Closed-closed. (Discussion below of ways to find appropriate Template:Math and Template:Math.) Two sequences of functions Template:Math, Template:Math on Template:Closed-closed, referred to as iterated integrals, are defined recursively as follows. First when Template:Math, they are taken to be identically equal to 1 on Template:Closed-closed. To obtain the next functions they are multiplied alternately by Template:Math and Template:Math and integrated, specifically, for Template:Math: Template:NumBlk2 Template:NumBlk2

The resulting iterated integrals are now applied as coefficients in the following two power series in λ: <math display="block"> u_0 = y_0 \sum_{k=0}^\infty \left (\lambda-\lambda_0^* \right )^k \tilde X^{(2k)},</math> <math display="block"> u_1 = y_0 \sum_{k=0}^\infty \left (\lambda-\lambda_0^* \right )^k X^{(2k+1)}.</math> Then for any Template:Mvar (real or complex), Template:Math and Template:Math are linearly independent solutions of the corresponding equation (Template:EquationNote). (The functions Template:Math and Template:Math take part in this construction through their influence on the choice of Template:Math.)

Next one chooses coefficients Template:Math and Template:Math so that the combination Template:Math satisfies the first boundary condition (Template:EquationNote). This is simple to do since Template:Math and Template:Math, for Template:Math. The values of Template:Math and Template:Math provide the values of Template:Math and Template:Math and the derivatives Template:Math and Template:Math, so the second boundary condition (Template:EquationNote) becomes an equation in a power series in Template:Mvar. For numerical work one may truncate this series to a finite number of terms, producing a calculable polynomial in Template:Mvar whose roots are approximations of the sought-after eigenvalues.

When Template:Math, this reduces to the original construction described above for a solution linearly independent to a given one. The representations (Template:EquationNote) and (Template:EquationNote) also have theoretical applications in Sturm–Liouville theory.<ref name="KP"/>

Construction of a nonvanishing solutionEdit

The SPPS method can, itself, be used to find a starting solution Template:Math. Consider the equation Template:Math; i.e., Template:Mvar, Template:Mvar, and Template:Mvar are replaced in (Template:EquationNote) by 0, Template:Math, and Template:Mvar respectively. Then the constant function 1 is a nonvanishing solution corresponding to the eigenvalue Template:Math. While there is no guarantee that Template:Math or Template:Math will not vanish, the complex function Template:Math will never vanish because two linearly-independent solutions of a regular Sturm–Liouville equation cannot vanish simultaneously as a consequence of the Sturm separation theorem. This trick gives a solution Template:Math of (Template:EquationNote) for the value Template:Math. In practice if (Template:EquationNote) has real coefficients, the solutions based on Template:Math will have very small imaginary parts which must be discarded.

ReferencesEdit

Template:Reflist