7
$\begingroup$

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?

|improve this question
$\endgroup$
  • $\begingroup$ (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] $\endgroup$ – Dr. Wolfgang Hintze Jul 27 '17 at 18:09
  • $\begingroup$ @Dr.WolfgangHintze Actually your example works fine for me. Mathematica returns an easily identifiable real solution. $\endgroup$ – becko Jul 27 '17 at 18:27
  • $\begingroup$ 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. $\endgroup$ – Dr. Wolfgang Hintze Jul 27 '17 at 19:52
  • $\begingroup$ @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. $\endgroup$ – becko Jul 27 '17 at 20:08
  • $\begingroup$ Ok. I still use version 10.1.0.0. $\endgroup$ – Dr. Wolfgang Hintze Jul 28 '17 at 7:42
9
$\begingroup$

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

|improve this answer
$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – becko Jul 22 '17 at 20:26
  • $\begingroup$ 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. $\endgroup$ – Jens Jul 22 '17 at 20:29
  • $\begingroup$ 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. $\endgroup$ – becko Jul 22 '17 at 20:37
  • $\begingroup$ Thank you for bringing DEigensystem to my attention. I had not considered it. $\endgroup$ – becko Jul 26 '17 at 15:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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