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.

  • 3
    $\begingroup$ That's not quite, but close to a question of the form "I need to do math. Where do I start?" What have you already tried? What kind of methods have you already found and what kind of question do you have specifically? $\endgroup$ Commented Jul 5, 2018 at 21:02
  • $\begingroup$ Did you try the basic Euler method? A more advanced method would be 4th order Runge Kutta $\endgroup$
    – physicsGuy
    Commented Jul 6, 2018 at 7:11
  • $\begingroup$ There are several questions in this site about eigenvalue problems, see for example 1 and 2. That being said, the easiest method is probably finite differences, although finite element methods are also used. $\endgroup$
    – nicoguaro
    Commented Jul 6, 2018 at 15:04

1 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.

  • 1
    $\begingroup$ So how do you take the derivative of $D$ in the case when all you have is a vector of values at (I assume) specific values of $x$? $\endgroup$ Commented Jul 7, 2018 at 16:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.