Numerical solution for eigenvectors and eigenvalues of a Sturm-Liouville problem

Question

I have to deal with the following problem in my research:

$$\left[\frac{1}{D}F_{x}\right]_{x} + \frac{f D_{x}}{c D^{2}}F = 0$$

with boundary conditions

$$F(0) = 0$$ $$F_{x}(L) = 0$$

where $f$ is a known constant and $D(x)$ is a known function of $x$. I should clarify that $D(x)$ is known numerically for all $x$ (I have a vector of values), but it is not known as a closed-form analytical function and is not readily approximated as such. Subscripts denote derivatives.

I note that this is a Sturm-Liouville problem:

$$[p(x)F_{x}(x)]_{x} + q(x)F(x) = -\lambda r(x) F(x)$$

with $p = 1/D$, $q = 0$, and $r = f D_{x}/D^{2}$.

$F$ are the eigenvectors and $\lambda = 1/c$ are the eigenvalues. Is there an easy way to solve for them numerically? Thanks for any help.

Answer 1

D'oh! I didn't realize how easy this was. My practical experience with numerical methods begins and ends with finite differences. I got caught up by the fact that this is an eigenvalue problem. I didn't realize that if I replace the derivatives with centered difference operators this reduces to a simple matrix eigenvalue problem!

Letting $D1$ and $D2$ be centered difference matrices, and with $p$ and $r$ defined as above, I can write:

$$[D2 + diag(p_{x}/p)D1]F = \lambda[diag(-r/p)]F$$

$$AF = \lambda BF$$

The Dirichlet condition at $x=0$ is easily built into the first rows of $A$ and $B$ and the Neumann condition at $x=L$ can be built in as well with a ghost point.

As a test, the eigenvectors and eigenvalues have an analytical form when $D(x) = a \exp(2bx)$, and my finite difference method compares with sufficient accuracy for my needs.

Thanks for the comments. Sorry this was so simple.