Bessel EVP and fzero

Question

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?

Thank in advance!

Answer 1

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.

Answer 2

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.