# Bessel EVP and fzero

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
xguess=1;
%% function handle
f=@(l) besselj(l-1,1)-2*besselj(l,1)+besselj(l+1,1)
%% set tolerance
opts=optimset('TolX',1e-12);
%% determine zero
xzero=fzero(f,xguess,opts);


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?

• 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. – nicoguaro Jan 13 '15 at 19:51
• 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$. – sebastian_g Jan 13 '15 at 20:53
• 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? – sebastian_g Jan 13 '15 at 21:00
• 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). – Christian Clason Jan 13 '15 at 22:16
• Is there a sign error in your code in front of besselj(l+1,1)? – 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),$$ 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.

• The code in my answer to a related question performs almost exactly this task. – Doug Lipinski Jan 14 '15 at 2:41
• @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. – Kirill Jan 14 '15 at 3:24
• 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. – Doug Lipinski Jan 14 '15 at 3:34
• The bracket you proposed works quite well. Thanks a lot. The log-plot is quite helpful. – sebastian_g Jan 14 '15 at 10:15

Here is the code I used now.

kmax=5;
lambda=zeros(kmax,1);
%% function handle
f=@(l) besselj(l-1,1)-2*besselj(l,1)-besselj(l+1,1);
%% set tolerance
opts=optimset('TolX',1e-12);
%% compute first positive eigenvalue
lambda(1)=fzero(f,[0.5,1.5],opts);
%% compute negative eigenvalues
for k=1:(kmax-1)
[lambda(k+1),~,fl]=fzero(f,[-k-0.5,-k+0.5],opts);
if not(fl==1)
error('iteration    did     not     convert.')
end
end


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