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?

$\endgroup$
5
  • $\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$ Jul 27, 2017 at 18:09
  • $\begingroup$ @Dr.WolfgangHintze Actually your example works fine for me. Mathematica returns an easily identifiable real solution. $\endgroup$
    – a06e
    Jul 27, 2017 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$ Jul 27, 2017 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$
    – a06e
    Jul 27, 2017 at 20:08
  • $\begingroup$ Ok. I still use version 10.1.0.0. $\endgroup$ Jul 28, 2017 at 7:42

1 Answer 1

10
$\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

$\endgroup$
4
  • $\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$
    – a06e
    Jul 22, 2017 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, 2017 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$
    – a06e
    Jul 22, 2017 at 20:37
  • $\begingroup$ Thank you for bringing DEigensystem to my attention. I had not considered it. $\endgroup$
    – a06e
    Jul 26, 2017 at 15:27

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.