Serving the Quantitative Finance Community

 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 6th, 2011, 4:12 pm

What is the best/easiest way to compute the roots x of:J_n(x) = 0for n = 0, 1, ...
 
User avatar
Alan
Posts: 3050
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

Computing the zeroes of Bessel functions

May 7th, 2011, 12:36 am

Truly, the easiest way is to buy Mathematica, where it is a built-in.If I _had_ to code it myself, I would attempt something like this.First, look up the asymptotic large x expression for J_n(x) inAbramowitz and Stegun; it is proportional to a simple cos(...); the k-th zero of the cos(...) is trivially found at(*) x_k = [(2 k -1)/2 + n/2 + 1/4] pi, (k=1,2,...)For n small, these can be remarkably accurate, even for low order zeros.For example, for J_0(x), the first 5 approximate zeros are at{2.35619, 5.49779, 8.63938, 11.781, 14.9226} Now, the exact first 5 zeros are(**) {2.40483, 5.52008, 8.65373, 11.7915, 14.9309}For n=0, the asymptotic zero is so close to the exact, that a simple Newton's method root finder,started at any element of (*) will rapidly converge to the corresponding exact value of (**). For larger n, the lowest order roots are not so well-approximated, so this may require somefine-tuning to make sure the Newton's method routine picks up all the correct roots.
Last edited by Alan on May 6th, 2011, 10:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 7th, 2011, 9:18 am

That's the answer I was hoping for, i.e. use results from other sources. There is no analytic solution.The case J_0 is of relevance for me here. The poor man's solution is a combination of lookup tables and asymptotic expansions for large x. Another way would be to use Bessel's integral form for J_n(x) if I had energy to work it out.Since we have an analytical J_n(x) as a series of Gamma functions we could write an offline(let it run as along as needed) program to crawl on the positive x axis, compute zeros and store them in permanent memory for later use.
Last edited by Cuchulainn on May 6th, 2011, 10:00 pm, edited 1 time in total.
 
User avatar
Alan
Posts: 3050
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

Computing the zeroes of Bessel functions

May 7th, 2011, 3:24 pm

I just confirmed that the percentage error in formula I postedfor x_k (and n=0) drops to less than 0.01% by k=12.So, if you just need J_0(x) zeros, perhaps the very simplest solution is to use the (6 digit) exact values: {2.40483, 5.52008, 8.65373, 11.7915, 14.9309, 18.0711, 21.2116, 24.3525, 27.4935, 30.6346, 33.7758}for the first 11 and the x_k formula I posted for k >= 12. (Or, whatever cutoff k satisfies your error tolerance). Then, you are done. p.s. Another easy solution is to just tell me how many zeros you want and to what precision, I ask Mathematica, post thelist and you simply copy it.
Last edited by Alan on May 6th, 2011, 10:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 9th, 2011, 4:43 am

This Bessel theory is cool stuff! It is possible to find the roots using properties of Bessel functions and Newton's method:d/dx J(0,x) = - J(1,x)d/dx J(n,x) = (J(n-1,x) - J(n+1,x))/2I get convergence in [2,5] iterations depending on the desired accuracy. It's a few lines of code (I call Boost Bessel functions).Initial estimates: for any x, the interleavng of roots of J(n+1) with those of J(n) could be used and the fact that the roots differ by pi = 3.1415... asymptotically.How many roots needeed? If we have exponentially small terms we probably only need a few terms. A good stress test of this might be this?
Last edited by Cuchulainn on May 8th, 2011, 10:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 9th, 2011, 5:34 am

Another root-finding problem is:f(x) = x d/dx J(0,x) - a J(0,x) = 0If I use Newton then we need a 'good' way to compute d/dx f(x)? (J(0)" = -J(1)' = (J(2) - J(0))/2)
Last edited by Cuchulainn on May 8th, 2011, 10:00 pm, edited 1 time in total.
 
User avatar
Alan
Posts: 3050
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

Computing the zeroes of Bessel functions

May 9th, 2011, 2:21 pm

Given a very good starting point (say my x_k formula for J0), then you canbastardize Newton's method with a crude denominator without much loss of efficiency.
 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 9th, 2011, 6:33 pm

QuoteOriginally posted by: AlanGiven a very good starting point (say my x_k formula for J0), then you canbastardize Newton's method with a crude denominator without much loss of efficiency.NR has a mind of its own regarding which root it converges to. Keep it on a leash.
 
User avatar
AVt
Posts: 90
Joined: December 29th, 2001, 8:23 pm

Computing the zeroes of Bessel functions

May 9th, 2011, 7:13 pm

Trust Alan [ even here :-) ]As BesselJ ~ cos a good guess will give success, as zeroes are like findingzeroes for cos by Newton: you only need the regime.BTW: why you need them and in what generality (your OP was wider)?And do not forget: Bessel is often computed by recursion from asymptotics.
 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 9th, 2011, 7:49 pm

Did I say that I did not trust Alan's tables of numbers?I am using using C++ so I have just J(n) from Boost (for me it is a black box). Rationale is to stress it by computing the roots which are the same as MM. The Boost authors have documented which algos they use e.g. asymptotics etc.QuoteAs BesselJ ~ cos a good guess will give success, as zeroes are like findingzeroes for cos by Newton: you only need the regime.J ~ sin??But these are less accurate than NR? The o(1/xsqrt(x))) remainder term is a bit scary but I have not compared it with NR yet, so maybe it's OK to use it. NR just uses function calls. If I had a choice between asymptics and iteration, I would choose the latter.QuoteBTW: why you need them and in what generality (your OP was wider)?Yes, e.g. exact solutions of equations and stress tests of BesselJ (and j,y, I,K).Almost forgot: curiosity and for fun as well.
Last edited by Cuchulainn on May 8th, 2011, 10:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 10th, 2011, 11:26 am

Last edited by Cuchulainn on May 9th, 2011, 10:00 pm, edited 1 time in total.
 
User avatar
nwhitehe
Posts: 0
Joined: March 3rd, 2006, 6:57 am

Computing the zeroes of Bessel functions

May 10th, 2011, 10:45 pm

You might also look at John Harrison, "Fast and Accurate Bessel Function Computation", ARITH2009. He has asymptotic expansions for J_0, Y_0, J_1, Y_1 that are better than Abramowitz & Stegun's. He uses them to get high precision function values, but they can also be used to get high precision zeros.For example, here is J_0(x) to three terms (from the paper):
 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 11th, 2011, 9:44 am

QuoteOriginally posted by: nwhiteheYou might also look at John Harrison, "Fast and Accurate Bessel Function Computation", ARITH2009. He has asymptotic expansions for J_0, Y_0, J_1, Y_1 that are better than Abramowitz & Stegun's. He uses them to get high precision function values, but they can also be used to get high precision zeros.For example, here is J_0(x) to three terms (from the paper):This will probably be fast as well.Will these results allow us to put bounds on Fourier-Bessel series? (i.e. when to truncate it).Looking at Table 1 of the paper for J(0) and J(1) I compare with the Math Toolkit with Newton Raphson on top of it (with 'naive' start point of being a distance of pi away). I get:
Last edited by Cuchulainn on May 10th, 2011, 10:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Topic Author
Posts: 23029
Joined: July 16th, 2004, 7:38 am

Computing the zeroes of Bessel functions

May 11th, 2011, 12:24 pm

QuoteOriginally posted by: nwhiteheYou might also look at John Harrison, "Fast and Accurate Bessel Function Computation", ARITH2009. He has asymptotic expansions for J_0, Y_0, J_1, Y_1 that are better than Abramowitz & Stegun's. He uses them to get high precision function values, but they can also be used to get high precision zeros.For example, here is J_0(x) to three terms (from the paper):I coded this formula (hopefully no error..):