Le 26/01/2015 10:57, Dang, Christophe a écrit :
Hello,
De Claus Futtrup
Envoyé : samedi 24 janvier 2015 17:53
I've made a small script to play around with Bessel functions of the first kind
[...]
for m=0:19 // ...
powerseries_m = ((-1)^m / (factorial(m) * factorial(m + n))) * (x/2)^(2*m);
powerseries = powerseries + powerseries_m; // sum the powerseries
end
Generally speaking, I think you'd better use the possibilities of vector
calculation as much as possible, to have a faster calculation.
The lines above could look like the following, assuming x is a column vector:
m = 0:19;
[x_mat, m_mat] = ndgrid(x, m); // rectangular matrices for group evaluation
f = 1 ./(factorial(m) .* factorial(m + n));
f_mat = ndgrid(f, x)';
powerseries = sum((-1).^m_mat.*f_mat.*(x_mat/2).^(2*m_mat), "c");
However, this quite naive implementation is probably not accurate when x and m
are high,
for f would have some zeros (if m or m+n > 170) leading to zeros in powerseries,
whereas f*(x/2)^(2m) might not be negligible.
The reason why you should follow Nikolay's advice and use built-in functions,
which usually use strong algorithms.
--
Christophe Dang Ngoc Chan
Mechanical calculation engineer
This e-mail may contain confidential and/or privileged information. If you are
not the intended recipient (or have received this e-mail in error), please
notify the sender immediately and destroy this e-mail. Any unauthorized
copying, disclosure or distribution of the material in this e-mail is strictly
forbidden.
_______________________________________________
users mailing list
[email protected]
http://lists.scilab.org/mailman/listinfo/users
Almost every numerical analysis textbook comprises a chapter on the
evaluation of Bessel functions, as for instance
Abramowitz ans Stegun, Handbook of mathematical functions, chapters
9,10 (elegant use of recurrence relation)
Press et al., Numerical recipes, chapter 6 (rational approximations)
Enjoy!
_______________________________________________
users mailing list
[email protected]
http://lists.scilab.org/mailman/listinfo/users
