1

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.

DCM
  • 21
  • 2
  • 3
    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?

     – Wolfgang Bangerth Jul 05 '18 at 21:02
  • Did you try the basic Euler method? A more advanced method would be 4th order Runge Kutta – physicsGuy Jul 06 '18 at 07:11
  • 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. – nicoguaro Jul 06 '18 at 15:04

1 Answers1

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.

DCM
  • 21
  • 2
  • 1
    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$? – Bill Greene Jul 07 '18 at 16:59