I am trying to solve the Eigenvalue problem $$ x^2 y''+ x y' + x^2 y = \lambda^2 y\,,\quad x\in(0,1)\,,\quad y(0)=0\,,\quad y'(1)=y(1) $$ The differential equation is the Bessel equation. The solution is given by the Bessel function of the 1st kind $y(x)=J_\lambda(x)$. Bessel functions of the 2nd kind is omitted because of their singularity at $x=0$ . The Eigenvalue $\lambda$ has to be determined using the boundary condition at $x=1$. Using the identity $$J'_\lambda(x)=\frac{1}{2}\left(J_{\lambda-1}\left(x\right)-J_{\lambda+1}\left(x\right)\right)$$ the boundary condition at $x=1$ can be rewritten as $$J_{\lambda-1}\left(1\right)-2J_{\lambda}\left(1\right)-J_{\lambda+1}\left(1\right)=0\,.$$ Now I would like to determine the $\lambda$ with the equation above using the fzero-routine of MATLAB/ GNU Octave. Code below.

%% initial guess
%% function handle
f=@(l) besselj(l-1,1)-2*besselj(l,1)+besselj(l+1,1)
%% set tolerance
%% determine zero

If you plot the function f, click here, you observe that most of the zeros are in $-\infty<x<0$.

To get a sorted sequence $\lambda_{k+1}<\lambda_k$, I would like to search from the "last" zero only in negative $x$-direction to the "next" one. Is there any way to implement this in the fzero-routine?

Thank in advance!

  • $\begingroup$ It's not clear for me what your question is. What do you mean with "from the "last" zero only in negative x-direction to the "next" one"? Normally, one just need a finite set of eigenvalues, then you can make a progressive search for this problem. $\endgroup$ – nicoguaro Jan 13 '15 at 19:51
  • $\begingroup$ Let me clarify the problem. I would like the first 10 eigenvalues. Let's say I found $\lambda_1$. Then $\lambda_2$ should have the property that there is no zero of the function f between $\lambda_1$ and $\lambda_2$. $\endgroup$ – sebastian_g Jan 13 '15 at 20:53
  • $\begingroup$ How can I program a progressive search effectively? How to give a step size? If I detect a change if the sign I could use bisection method, right? $\endgroup$ – sebastian_g Jan 13 '15 at 21:00
  • $\begingroup$ You can specify an interval in place of an initial guess for fzero to search in, e.g. xzero=fzero(f,[-1e-9,xzero_old],opts). $\endgroup$ – Christian Clason Jan 13 '15 at 22:16
  • $\begingroup$ Is there a sign error in your code in front of besselj(l+1,1)? $\endgroup$ – Kirill Jan 13 '15 at 22:35

As you can see in this plot of $\log|f(\lambda)|$, $$ f(\lambda) = J_{\lambda-1}(1) - 2J_\lambda(1) -J_{\lambda+1}(1), $$

enter image description here

the roots $\lambda_k$ are really regular, and are approximately equal to $-k$ (starting from $k\geq0$, $\lambda_0=1.23219$ is an exception).

So the way to get the $k$-th eigenvalue ($k\geq1$), is to bracket the root in $[-k-\frac12,-k+\frac12]$ and use fzero solver in that interval (it accepts a bracketing interval as the initial guess and stays within that interval).

Note that the magnitude of the function grows very quickly (something like $(2\lambda/e)^\lambda\sin\pi\lambda$), so once the function becomes quite large (for large $k$), it may be better to use $\lambda_k=-k$ as an approximation.

| cite | improve this answer | |
  • $\begingroup$ The code in my answer to a related question performs almost exactly this task. $\endgroup$ – Doug Lipinski Jan 14 '15 at 2:41
  • $\begingroup$ @DougLipinski Yes, but your approach requires you to guess a "comb" of locations to separate roots, and you might miss some roots easily. In this exact question it's possible to find all roots guaranteed. $\endgroup$ – Kirill Jan 14 '15 at 3:24
  • 1
    $\begingroup$ Indeed, just choose the intervals in my code to be the ones you suggested. I just wanted to point out some existing code that solves the OP's problem with the choice of intervals you suggested. You've shown the right "comb" for this problem. +1 from me. $\endgroup$ – Doug Lipinski Jan 14 '15 at 3:34
  • $\begingroup$ The bracket you proposed works quite well. Thanks a lot. The log-plot is quite helpful. $\endgroup$ – sebastian_g Jan 14 '15 at 10:15

Here is the code I used now.

%% function handle
f=@(l) besselj(l-1,1)-2*besselj(l,1)-besselj(l+1,1);
%% set tolerance
%% compute first positive eigenvalue
%% compute negative eigenvalues
for k=1:(kmax-1)
  if not(fl==1)
    error('iteration    did     not     convert.')

By the way, I would like to use this to verify my finite difference method.

| cite | improve this answer | |

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.