7

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?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
a06e
  • 11,327
  • 4
  • 48
  • 108
  • (at)becko There's an even simpler example with the same technical "difficulties": DSolve[{y''[x] + b ^2 y[x] == 0, y[0] == 0, y[1] == 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. – a06e 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] == 0, y[1] == 0}, y[x], x] correctly solves the eigenvalue problem, without having to separate the b.c.s. – a06e Jul 27 '17 at 20:08
  • Ok. I still use version 10.1.0.0. – Dr. Wolfgang Hintze Jul 28 '17 at 07:42

1 Answers1

10

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] == 0, 
   f[1] == 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

Jens
  • 97,245
  • 7
  • 213
  • 499
  • 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). – a06e 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. – a06e Jul 22 '17 at 20:37
  • Thank you for bringing DEigensystem to my attention. I had not considered it. – a06e Jul 26 '17 at 15:27