# How to solve a Sturm-Liouville problem with Mathematica (or, how to go from the complex to the general *real* solution)?

Let $a>0$ be a constant positive number. I am stuck trying to solve the following regular Sturm-Liouville problem:

$$\frac{\mathrm d}{\mathrm d x}((a+x)f'(x)) = -v f(x),\qquad f(0)=f(1)=0$$

where $v$ is the eigenvalue.

According to Mathematica, the general solution of the ODE is (found by running DSolve[{D[(x + a) f'[x], x] + v f[x] == 0}, f[x], x]):

$$f(x) = c_1 I_0(2\sqrt{-(a+x)v}) + c_2 K_0(2\sqrt{-(a+x)v})$$

However, no combination of $c_1,c_2,v$ can give a non-trivial function $f(x) \not\equiv 0$ while satisfying $f(0)=f(1)=0$. However according to Sturm Liouville theory, eigenvalues and non-trivial eigenfunctions must exist.

So how can I use Mathematica to solve this problem? Or, how can I obtain the general real solution of an ODE with real coefficients?

• (at)becko There's an even simpler example with the same technical "difficulties": DSolve[{y''[x] + b ^2 y[x] == 0, y == 0, y == 0}, y[x], x] – Dr. Wolfgang Hintze Jul 27 '17 at 18:09
• @Dr.WolfgangHintze Actually your example works fine for me. Mathematica returns an easily identifiable real solution. – becko Jul 27 '17 at 18:27
• the solution returned in my example is f = 0. Therefore I adopt the procedure to request only the b.c. at x = 0, find the solution, and implement the second b.c. afterwards. This then leads to the eigenvalue equation. – Dr. Wolfgang Hintze Jul 27 '17 at 19:52
• @Dr.WolfgangHintze In Mathematica v11.1.1, DSolve[{y''[x] + b ^2 y[x] == 0, y == 0, y == 0}, y[x], x] correctly solves the eigenvalue problem, without having to separate the b.c.s. – becko Jul 27 '17 at 20:08
• Ok. I still use version 10.1.0.0. – Dr. Wolfgang Hintze Jul 28 '17 at 7:42

Here is the numerical approach, showing that the eigenvalues and eigenfunctions are real:

{eigN, funcN} =
With[{a = 1},
N@DEigensystem[{D[(x + a) f'[x], x],
DirichletCondition[f[x] == 0, True]}, f[x], {x, 0, 1}, 3]];

Plot[funcN, {x, 0, 1}, PlotLegends -> eigN] The eigenvalues in the plot above are the negative of v.

For the analytical solution starting from your code, you'll have to find the roots of a transcendental equation in v defining the eigenvalue condition:

PiecewiseExpand[
DSolveValue[{D[(x + a) f'[x], x] + v f[x] == 0, f == 0,
f == 0}, f[x], x]]


The output of the above will give tell you the necessary condition for a non-trivial solution: it is

condition =
BesselI[0, 2 Sqrt[-a v]]-(BesselI[0,2 Sqrt[-(1 + a) v]] BesselK[0, 2 Sqrt[-a v]])/
BesselK[0, 2 Sqrt[-(1 + a) v]] == 0


This is the transcendental equation for the eigenvalues v. Solving it will require assigning a numerical value to a and then using FindRoot. At this stage the problem becomes similar to this one: Labeling solutions of an Eigenvalue equation involving Bessel functions.

The results of the root finding agree with the DEigensystem results shown earlier:

Block[{a = 1}, Chop@First@FindRoot[condition, {v, 10.}]]


v->14.3377

• Thank you! Is there any hope of obtaining an analytical expression? Perhaps related to Bessel functions (the eigenvalues should be related to zeros of these functions). – becko Jul 22 '17 at 20:26
• The condition in the output of DSolveValue (which I didn't bother copying) looks too complicated to reduce to straight Bessel zeros. I think it will have to be left in implicit form, especially since a is an arbitrary parameter. – Jens Jul 22 '17 at 20:29
• Yes, the condition can be obtained by evaluating $f(0),f(1)$ (the one that I obtained using DSolve, in terms of Bessel funs) and demanding that they are zero). Is there no connection between these eigenfunctions and the general form of $f(x)$ found by DSolve? I am curious about this, because I tried to make $f(0) = f(1) = 0$ and all I got was the trivial zero function. – becko Jul 22 '17 at 20:37
• Thank you for bringing DEigensystem to my attention. I had not considered it. – becko Jul 26 '17 at 15:27