MATLAB Answers

Bessel function - roots and zeros on interval <12,20>

42 views (last 30 days)
Jaroslav Zadnik
Jaroslav Zadnik on 30 Jan 2020
Edited: David Goodmanson on 1 Feb 2020
Hi, I'm beginner here. I have Bessel function of the first kind J1. I need find roots and zeros on interval <10,20>. Could you help me plase? :/

  3 Comments

Walter Roberson
Walter Roberson on 30 Jan 2020
I suggest that you plot the function and use what you see for estimates of starting points for fzero
Jaroslav Zadnik
Jaroslav Zadnik on 30 Jan 2020
I have only:
x=10:0.1:20;
i=1;
J=besselj(i,x)
plot(x,J)
John D'Errico
John D'Errico on 31 Jan 2020
You are inconsistent, since in the title you ask about the interval [12,20], yet in your question, you say [10,20].

Sign in to comment.

Answers (2)

David Goodmanson
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

  0 Comments

Sign in to comment.


John D'Errico
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.

  0 Comments

Sign in to comment.

Sign in to answer this question.