0

im writing a programm for a fitting routine and am currently optimizing the code for faster calculations. The weakes point is a part, where i have to calculate a big amount of bessel functions, which takes around 0.7 s. In my case q has 177 entries, th 100 and R 400.

  Js = zeros(numel(th),numel(q)); tR=sin(th')*R;
  for k = 1:numel(q)
  Js(:,k) = sum(tn.*besselj(0,q(k)*tR),2);
  end

I also tried to make a 3D-Matrix, but it takes slightly longer to calculate.

[Q,T,RR]=meshgrid(q,sin(th),R);
Js1 = besselj(0,Q.*T.*RR);

So, i'm wondering, is there a way to calculate these besselfunctions faster? thanks in advance, kuy

kuy
  • 13
  • 6
  • I don't think there are. You're limited by the use of the built-in `besselj`. – rayryeng Aug 03 '16 at 14:03
  • Have you tried using `bsxfun`? – flawr Aug 03 '16 at 14:06
  • What are the sizes of the inputs? How many dimensions do each of the inputs have? Which ones are vectors and again are they row or column vectors and which ones are higher dim matrices? Telling us the count of entries doesn't give us that info. – Divakar Aug 03 '16 at 14:50
  • @Divakar the matrix, i need to calculate besselfunctions from is of the size 177x100x400 (as an order of magnitude - it will vary a bit), so the vectors are all onedimensional. – kuy Aug 04 '16 at 08:02
  • @flawr, i'm not sure how to use bsxfun in this case, do you mean the same approach as rahnema1 (in the answer) – kuy Aug 04 '16 at 08:02

1 Answers1

0

Wolfram shows a special case of Bessel function where the first argument of the function is 0. So we can pre-compute some variables like sgn_cum and km2 to simplify the computation:

n = 10;
k = 0 : n;
sgn_cum = ((-1 * 0.25) .^ k ./ cumprod([1 1:n]).^2).';
km2 = 2*k;

Given a column vector z we can obtain Bessel function as:

bsxfun(@power,z,km2 ) * sgn_cum

Example:

z= (1:5).';
result = bsxfun(@power,z,km2 ) * sgn_cum;

We can decrease n to speed up the computation but with the cost of lower precision.

rahnema1
  • 15,264
  • 3
  • 15
  • 27
  • hi, thanks for your answer, i tried your function, but was not very successful. for n = 10 the built in besselj is around twice as fast. also the values agree with besselj just until around z = 8, then the function diverges.. the link you sent is for wolfram, i'm using matlab, do you think they use the same algorithm? – kuy Aug 04 '16 at 07:58
  • @kuy , I edited the code to run faster, please check it. the link from wolfram explains theoretical points, I do not know what is the implementation of both softwares, I implemented it directly from formula. if n increased may precision improved – rahnema1 Aug 04 '16 at 10:34
  • yes, i know what you mean (actually it doesnt really matter, what matlab uses..), anyway its some kind of power series expansion and every term adds another turn of the function, but eventually it diverges. since i need to calculate values up to 2e3, i think this will be not very efficient (and its still slower than besselj even for small values..) anyway thanks for your effort! :) – kuy Aug 04 '16 at 11:14