49 views (last 30 days)

David Goodmanson
on 31 Jan 2020

Edited: David Goodmanson
on 1 Feb 2020

HI Jaroslav

Do roots and zeros mean the same thing? The code below uses Newton's method to find the roots, which has the advantage that you can apply it to a vector of roots and get them all at once. I don't have a tolerance check for the number of times to iterate Newton, but j = 1:8 seems to work (good up to at least the 14th decimal place for all roots in this case).

function x = besselzeroj(n,q,opt)

% first q roots of bessel function Jn(x), integer order.

% if opt = 'd', first q roots of dJn(x)/dx, integer order.

% if opt is not provided, the default is Jn(x).

% includes the root x = 0 for dJ0(x)/dx (zeroth order)

% but not for higher orders.

%

% starting point for for zeros of Jn is borrowed from Cleve Moler,

% but starting points for both Jn and Jn' can be found in

% Abramowitz and Stegun 9.5.12, 9.5.13.

%

% David Goodmanson

%

% function x = besselzeroj(n,q,opt)

k = 1:q;

if nargin==3 & opt=='d'

% find starting points for J'

beta = (k + n/2 - 3/4)*pi;

mu = 4*n^2;

x = beta - (mu+3)./(8*beta) - 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);

% Newton

for j=1:8

xnew = x - besseljdv(n,x)./ ...

(besselj(n,x).*((n^2./x.^2)-1) -besseljdv(n,x)./x);

x = xnew;

end

else

% find starting points for J

beta = (k + n/2 - 1/4)*pi;

mu = 4*n^2;

x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);

% Newton

for j=1:8

xnew = x - besselj(n,x)./besseljdv(n,x);

x = xnew;

end

end

end

function y = besseljdv(n,x)

% derivative of bessel function of integer order,

% plain vanilla version

%

% function y = besseljdv(n,x)

y = -besselj(n+1,x) + n*besselj(n,x)./x;

% get rid of nans, integer case

if n==1

y(x==0) = 1/2;

else

y(x==0) = 0;

end

end

format long

besselzeroj(1,12)'

ans =

3.831705970207512

7.015586669815619

10.173468135062722

13.323691936314225

16.470630050877634

19.615858510468243

22.760084380592772

25.903672087618382

29.046828534916855

32.189679910974405

35.332307550083868

38.474766234771614

sars sabina
on 20 Dec 2020

Could you please explain more about the starting point that you used in your code?

% find starting points for J

beta = (k + n/2 - 1/4)*pi;

mu = 4*n^2;

x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);

David Goodmanson
on 21 Dec 2020

Hi s s,

I can't say much more than what is in the answer, that the expression for x is in Abramowitz and Stegun, Handbood of Mathematical Functions 9.5.12. I used just the first three terms in 9.5.12, which is an asymptotic expansion that gets more accurate the larger the roots are. But it turns out that even for small roots, the three-term expression does well at providing a first guess. In the expression, the first term can be derived fairly easily by evaluating a certain integral by the method of steepest descents. I think I could probably get the second term if I worked at it enough. After that, not happening.

John D'Errico
on 31 Jan 2020

This is easier than you think. :) Of course, that often is true, once you know the answer.

By the way, it can be a dangerous idea in general to use i as a variable name, since i is also sqrt(-1) by default in MATLAB. One day that will be important for you.

So what you have is a start:

x = 10:0.1:20;

n = 1;

J = besselj(n,x);

plot(x,J)

Now, do you know of some tool that can find a root, given a start point? Better yet, what if you have an interval in x that brackets a root? Note that fzero does exactly what you need, if you have an interval that brackets the root.

Is there a way you can locate pairs of points that bracket the root? What does that mean? Essentially, you need to find a pair of points x(m) and x(m+1), such that J(m) and J(m+1) have opposite signs. The sign function can help you here. For example:

m = find(diff(sign(J)) ~= 0)

m =

2 34 65 97

[J(m);J(m+1)]

ans =

0.018395515457572 -0.005177480554671 0.013894680686257 -0.002856572403405

-0.006615743297723 0.016599019864009 -0.005764213735631 0.015100612097755

So it appears there will be 4 roots in the interval of interest, because there are 4 zero crossings. Here is the first root:

fun = @(x) besselj(n,x);

[xvl,fval] = fzero(fun,[x(m(1)),x(m(1) + 1)])

xvl =

10.173468135062722

fval =

2.805130098197981e-16

Do the same call for each of the zero crossings identified in m.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 3 Comments

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/502752-bessel-function-roots-and-zeros-on-interval-12-20#comment_791615

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/502752-bessel-function-roots-and-zeros-on-interval-12-20#comment_791615

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/502752-bessel-function-roots-and-zeros-on-interval-12-20#comment_791617

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/502752-bessel-function-roots-and-zeros-on-interval-12-20#comment_791617

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/502752-bessel-function-roots-and-zeros-on-interval-12-20#comment_792179

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/502752-bessel-function-roots-and-zeros-on-interval-12-20#comment_792179

Sign in to comment.